Actual source code: eige.c

      2: #include <petsc/private/kspimpl.h>
      3: #include <petscdm.h>
      4: #include <petscblaslapack.h>

      6: typedef struct {
      7:   KSP ksp;
      8:   Vec work;
      9: } Mat_KSP;

     11: static PetscErrorCode MatCreateVecs_KSP(Mat A,Vec *X,Vec *Y)
     12: {
     13:   Mat_KSP        *ctx;
     14:   Mat            M;

     18:   MatShellGetContext(A,&ctx);
     19:   KSPGetOperators(ctx->ksp,&M,NULL);
     20:   MatCreateVecs(M,X,Y);
     21:   return(0);
     22: }

     24: static PetscErrorCode MatMult_KSP(Mat A,Vec X,Vec Y)
     25: {
     26:   Mat_KSP        *ctx;

     30:   MatShellGetContext(A,&ctx);
     31:   KSP_PCApplyBAorAB(ctx->ksp,X,Y,ctx->work);
     32:   return(0);
     33: }

     35: /*@
     36:     KSPComputeOperator - Computes the explicit preconditioned operator, including diagonal scaling and null
     37:     space removal if applicable.

     39:     Collective on ksp

     41:     Input Parameter:
     42: +   ksp - the Krylov subspace context
     43: -   mattype - the matrix type to be used

     45:     Output Parameter:
     46: .   mat - the explict preconditioned operator

     48:     Notes:
     49:     This computation is done by applying the operators to columns of the
     50:     identity matrix.

     52:     Currently, this routine uses a dense matrix format for the output operator if mattype == NULL.
     53:     This routine is costly in general, and is recommended for use only with relatively small systems.

     55:     Level: advanced

     57: .seealso: KSPComputeEigenvaluesExplicitly(), PCComputeOperator(), KSPSetDiagonalScale(), KSPSetNullSpace(), MatType
     58: @*/
     59: PetscErrorCode  KSPComputeOperator(KSP ksp, MatType mattype, Mat *mat)
     60: {
     62:   PetscInt       N,M,m,n;
     63:   Mat_KSP        ctx;
     64:   Mat            A,Aksp;

     69:   KSPGetOperators(ksp,&A,NULL);
     70:   MatGetLocalSize(A,&m,&n);
     71:   MatGetSize(A,&M,&N);
     72:   MatCreateShell(PetscObjectComm((PetscObject)ksp),m,n,M,N,&ctx,&Aksp);
     73:   MatShellSetOperation(Aksp,MATOP_MULT,(void (*)(void))MatMult_KSP);
     74:   MatShellSetOperation(Aksp,MATOP_CREATE_VECS,(void (*)(void))MatCreateVecs_KSP);
     75:   ctx.ksp = ksp;
     76:   MatCreateVecs(A,&ctx.work,NULL);
     77:   MatComputeOperator(Aksp,mattype,mat);
     78:   VecDestroy(&ctx.work);
     79:   MatDestroy(&Aksp);
     80:   return(0);
     81: }

     83: /*@
     84:    KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
     85:    preconditioned operator using LAPACK.

     87:    Collective on ksp

     89:    Input Parameter:
     90: +  ksp - iterative context obtained from KSPCreate()
     91: -  n - size of arrays r and c

     93:    Output Parameters:
     94: +  r - real part of computed eigenvalues, provided by user with a dimension at least of n
     95: -  c - complex part of computed eigenvalues, provided by user with a dimension at least of n

     97:    Notes:
     98:    This approach is very slow but will generally provide accurate eigenvalue
     99:    estimates.  This routine explicitly forms a dense matrix representing
    100:    the preconditioned operator, and thus will run only for relatively small
    101:    problems, say n < 500.

    103:    Many users may just want to use the monitoring routine
    104:    KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
    105:    to print the singular values at each iteration of the linear solve.

    107:    The preconditoner operator, rhs vector, solution vectors should be
    108:    set before this routine is called. i.e use KSPSetOperators(),KSPSolve() or
    109:    KSPSetOperators()

    111:    Level: advanced

    113: .seealso: KSPComputeEigenvalues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSPSetOperators(), KSPSolve()
    114: @*/
    115: PetscErrorCode  KSPComputeEigenvaluesExplicitly(KSP ksp,PetscInt nmax,PetscReal r[],PetscReal c[])
    116: {
    117:   Mat               BA;
    118:   PetscErrorCode    ierr;
    119:   PetscMPIInt       size,rank;
    120:   MPI_Comm          comm;
    121:   PetscScalar       *array;
    122:   Mat               A;
    123:   PetscInt          m,row,nz,i,n,dummy;
    124:   const PetscInt    *cols;
    125:   const PetscScalar *vals;

    128:   PetscObjectGetComm((PetscObject)ksp,&comm);
    129:   KSPComputeOperator(ksp,MATDENSE,&BA);
    130:   MPI_Comm_size(comm,&size);
    131:   MPI_Comm_rank(comm,&rank);

    133:   MatGetSize(BA,&n,&n);
    134:   if (size > 1) { /* assemble matrix on first processor */
    135:     MatCreate(PetscObjectComm((PetscObject)ksp),&A);
    136:     if (!rank) {
    137:       MatSetSizes(A,n,n,n,n);
    138:     } else {
    139:       MatSetSizes(A,0,0,n,n);
    140:     }
    141:     MatSetType(A,MATMPIDENSE);
    142:     MatMPIDenseSetPreallocation(A,NULL);
    143:     PetscLogObjectParent((PetscObject)BA,(PetscObject)A);

    145:     MatGetOwnershipRange(BA,&row,&dummy);
    146:     MatGetLocalSize(BA,&m,&dummy);
    147:     for (i=0; i<m; i++) {
    148:       MatGetRow(BA,row,&nz,&cols,&vals);
    149:       MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
    150:       MatRestoreRow(BA,row,&nz,&cols,&vals);
    151:       row++;
    152:     }

    154:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    155:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    156:     MatDenseGetArray(A,&array);
    157:   } else {
    158:     MatDenseGetArray(BA,&array);
    159:   }

    161: #if defined(PETSC_HAVE_ESSL)
    162:   /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */
    163:   if (!rank) {
    164:     PetscScalar  sdummy,*cwork;
    165:     PetscReal    *work,*realpart;
    166:     PetscBLASInt clen,idummy,lwork,bn,zero = 0;
    167:     PetscInt     *perm;

    169: #if !defined(PETSC_USE_COMPLEX)
    170:     clen = n;
    171: #else
    172:     clen = 2*n;
    173: #endif
    174:     PetscMalloc1(clen,&cwork);
    175:     idummy = -1;                /* unused */
    176:     PetscBLASIntCast(n,&bn);
    177:     lwork  = 5*n;
    178:     PetscMalloc1(lwork,&work);
    179:     PetscMalloc1(n,&realpart);
    180:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
    181:     PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_(&zero,array,&bn,cwork,&sdummy,&idummy,&idummy,&bn,work,&lwork));
    182:     PetscFPTrapPop();
    183:     PetscFree(work);

    185:     /* For now we stick with the convention of storing the real and imaginary
    186:        components of evalues separately.  But is this what we really want? */
    187:     PetscMalloc1(n,&perm);

    189: #if !defined(PETSC_USE_COMPLEX)
    190:     for (i=0; i<n; i++) {
    191:       realpart[i] = cwork[2*i];
    192:       perm[i]     = i;
    193:     }
    194:     PetscSortRealWithPermutation(n,realpart,perm);
    195:     for (i=0; i<n; i++) {
    196:       r[i] = cwork[2*perm[i]];
    197:       c[i] = cwork[2*perm[i]+1];
    198:     }
    199: #else
    200:     for (i=0; i<n; i++) {
    201:       realpart[i] = PetscRealPart(cwork[i]);
    202:       perm[i]     = i;
    203:     }
    204:     PetscSortRealWithPermutation(n,realpart,perm);
    205:     for (i=0; i<n; i++) {
    206:       r[i] = PetscRealPart(cwork[perm[i]]);
    207:       c[i] = PetscImaginaryPart(cwork[perm[i]]);
    208:     }
    209: #endif
    210:     PetscFree(perm);
    211:     PetscFree(realpart);
    212:     PetscFree(cwork);
    213:   }
    214: #elif !defined(PETSC_USE_COMPLEX)
    215:   if (!rank) {
    216:     PetscScalar  *work;
    217:     PetscReal    *realpart,*imagpart;
    218:     PetscBLASInt idummy,lwork;
    219:     PetscInt     *perm;

    221:     idummy   = n;
    222:     lwork    = 5*n;
    223:     PetscMalloc2(n,&realpart,n,&imagpart);
    224:     PetscMalloc1(5*n,&work);
    225:     {
    226:       PetscBLASInt lierr;
    227:       PetscScalar  sdummy;
    228:       PetscBLASInt bn;

    230:       PetscBLASIntCast(n,&bn);
    231:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
    232:       PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,array,&bn,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr));
    233:       if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
    234:       PetscFPTrapPop();
    235:     }
    236:     PetscFree(work);
    237:     PetscMalloc1(n,&perm);

    239:     for (i=0; i<n; i++)  perm[i] = i;
    240:     PetscSortRealWithPermutation(n,realpart,perm);
    241:     for (i=0; i<n; i++) {
    242:       r[i] = realpart[perm[i]];
    243:       c[i] = imagpart[perm[i]];
    244:     }
    245:     PetscFree(perm);
    246:     PetscFree2(realpart,imagpart);
    247:   }
    248: #else
    249:   if (!rank) {
    250:     PetscScalar  *work,*eigs;
    251:     PetscReal    *rwork;
    252:     PetscBLASInt idummy,lwork;
    253:     PetscInt     *perm;

    255:     idummy = n;
    256:     lwork  = 5*n;
    257:     PetscMalloc1(5*n,&work);
    258:     PetscMalloc1(2*n,&rwork);
    259:     PetscMalloc1(n,&eigs);
    260:     {
    261:       PetscBLASInt lierr;
    262:       PetscScalar  sdummy;
    263:       PetscBLASInt nb;
    264:       PetscBLASIntCast(n,&nb);
    265:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
    266:       PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&nb,array,&nb,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&lierr));
    267:       if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
    268:       PetscFPTrapPop();
    269:     }
    270:     PetscFree(work);
    271:     PetscFree(rwork);
    272:     PetscMalloc1(n,&perm);
    273:     for (i=0; i<n; i++) perm[i] = i;
    274:     for (i=0; i<n; i++) r[i]    = PetscRealPart(eigs[i]);
    275:     PetscSortRealWithPermutation(n,r,perm);
    276:     for (i=0; i<n; i++) {
    277:       r[i] = PetscRealPart(eigs[perm[i]]);
    278:       c[i] = PetscImaginaryPart(eigs[perm[i]]);
    279:     }
    280:     PetscFree(perm);
    281:     PetscFree(eigs);
    282:   }
    283: #endif
    284:   if (size > 1) {
    285:     MatDenseRestoreArray(A,&array);
    286:     MatDestroy(&A);
    287:   } else {
    288:     MatDenseRestoreArray(BA,&array);
    289:   }
    290:   MatDestroy(&BA);
    291:   return(0);
    292: }

    294: static PetscErrorCode PolyEval(PetscInt nroots,const PetscReal *r,const PetscReal *c,PetscReal x,PetscReal y,PetscReal *px,PetscReal *py)
    295: {
    296:   PetscInt  i;
    297:   PetscReal rprod = 1,iprod = 0;

    300:   for (i=0; i<nroots; i++) {
    301:     PetscReal rnew = rprod*(x - r[i]) - iprod*(y - c[i]);
    302:     PetscReal inew = rprod*(y - c[i]) + iprod*(x - r[i]);
    303:     rprod = rnew;
    304:     iprod = inew;
    305:   }
    306:   *px = rprod;
    307:   *py = iprod;
    308:   return(0);
    309: }

    311: #include <petscdraw.h>
    312: /* collective on ksp */
    313: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp,PetscInt neig,const PetscReal *r,const PetscReal *c)
    314: {
    316:   PetscReal      xmin,xmax,ymin,ymax,*xloc,*yloc,*value,px0,py0,rscale,iscale;
    317:   PetscInt       M,N,i,j;
    318:   PetscMPIInt    rank;
    319:   PetscViewer    viewer;
    320:   PetscDraw      draw;
    321:   PetscDrawAxis  drawaxis;

    324:   MPI_Comm_rank(PetscObjectComm((PetscObject)ksp),&rank);
    325:   if (rank) return(0);
    326:   M    = 80;
    327:   N    = 80;
    328:   xmin = r[0]; xmax = r[0];
    329:   ymin = c[0]; ymax = c[0];
    330:   for (i=1; i<neig; i++) {
    331:     xmin = PetscMin(xmin,r[i]);
    332:     xmax = PetscMax(xmax,r[i]);
    333:     ymin = PetscMin(ymin,c[i]);
    334:     ymax = PetscMax(ymax,c[i]);
    335:   }
    336:   PetscMalloc3(M,&xloc,N,&yloc,M*N,&value);
    337:   for (i=0; i<M; i++) xloc[i] = xmin - 0.1*(xmax-xmin) + 1.2*(xmax-xmin)*i/(M-1);
    338:   for (i=0; i<N; i++) yloc[i] = ymin - 0.1*(ymax-ymin) + 1.2*(ymax-ymin)*i/(N-1);
    339:   PolyEval(neig,r,c,0,0,&px0,&py0);
    340:   rscale = px0/(PetscSqr(px0)+PetscSqr(py0));
    341:   iscale = -py0/(PetscSqr(px0)+PetscSqr(py0));
    342:   for (j=0; j<N; j++) {
    343:     for (i=0; i<M; i++) {
    344:       PetscReal px,py,tx,ty,tmod;
    345:       PolyEval(neig,r,c,xloc[i],yloc[j],&px,&py);
    346:       tx   = px*rscale - py*iscale;
    347:       ty   = py*rscale + px*iscale;
    348:       tmod = PetscSqr(tx) + PetscSqr(ty); /* modulus of the complex polynomial */
    349:       if (tmod > 1) tmod = 1.0;
    350:       if (tmod > 0.5 && tmod < 1) tmod = 0.5;
    351:       if (tmod > 0.2 && tmod < 0.5) tmod = 0.2;
    352:       if (tmod > 0.05 && tmod < 0.2) tmod = 0.05;
    353:       if (tmod < 1e-3) tmod = 1e-3;
    354:       value[i+j*M] = PetscLogReal(tmod) / PetscLogReal(10.0);
    355:     }
    356:   }
    357:   PetscViewerDrawOpen(PETSC_COMM_SELF,NULL,"Iteratively Computed Eigen-contours",PETSC_DECIDE,PETSC_DECIDE,450,450,&viewer);
    358:   PetscViewerDrawGetDraw(viewer,0,&draw);
    359:   PetscDrawTensorContour(draw,M,N,NULL,NULL,value);
    360:   if (0) {
    361:     PetscDrawAxisCreate(draw,&drawaxis);
    362:     PetscDrawAxisSetLimits(drawaxis,xmin,xmax,ymin,ymax);
    363:     PetscDrawAxisSetLabels(drawaxis,"Eigen-counters","real","imag");
    364:     PetscDrawAxisDraw(drawaxis);
    365:     PetscDrawAxisDestroy(&drawaxis);
    366:   }
    367:   PetscViewerDestroy(&viewer);
    368:   PetscFree3(xloc,yloc,value);
    369:   return(0);
    370: }