Actual source code: eige.c

petsc-master 2020-08-13
Report Typos and Errors

  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: }