Actual source code: eige.c

petsc-master 2019-07-18
Report Typos and Errors

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

  6: static PetscErrorCode MatMult_KSP(Mat A,Vec X,Vec Y)
  7: {
  8:   KSP            ksp;
  9:   DM             dm;
 10:   Vec            work;

 14:   MatShellGetContext(A,&ksp);
 15:   KSPGetDM(ksp,&dm);
 16:   DMGetGlobalVector(dm,&work);
 17:   KSP_PCApplyBAorAB(ksp,X,Y,work);
 18:   DMRestoreGlobalVector(dm,&work);
 19:   return(0);
 20: }

 22: /*@
 23:     KSPComputeOperator - Computes the explicit preconditioned operator, including diagonal scaling and null
 24:     space removal if applicable.

 26:     Collective on ksp

 28:     Input Parameter:
 29: +   ksp - the Krylov subspace context
 30: -   mattype - the matrix type to be used

 32:     Output Parameter:
 33: .   mat - the explict preconditioned operator

 35:     Notes:
 36:     This computation is done by applying the operators to columns of the
 37:     identity matrix.

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

 42:     Level: advanced

 44: .seealso: KSPComputeEigenvaluesExplicitly(), PCComputeOperator(), KSPSetDiagonalScale(), KSPSetNullSpace(), MatType
 45: @*/
 46: PetscErrorCode  KSPComputeOperator(KSP ksp, MatType mattype, Mat *mat)
 47: {
 49:   PetscInt       N,M,m,n;
 50:   Mat            A,Aksp;

 55:   KSPGetOperators(ksp,&A,NULL);
 56:   MatGetLocalSize(A,&m,&n);
 57:   MatGetSize(A,&M,&N);
 58:   MatCreateShell(PetscObjectComm((PetscObject)ksp),m,n,M,N,ksp,&Aksp);
 59:   MatShellSetOperation(Aksp,MATOP_MULT,(void (*)(void))MatMult_KSP);
 60:   MatComputeOperator(Aksp,mattype,mat);
 61:   MatDestroy(&Aksp);
 62:   return(0);
 63: }

 65: /*@
 66:    KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
 67:    preconditioned operator using LAPACK.

 69:    Collective on ksp

 71:    Input Parameter:
 72: +  ksp - iterative context obtained from KSPCreate()
 73: -  n - size of arrays r and c

 75:    Output Parameters:
 76: +  r - real part of computed eigenvalues, provided by user with a dimension at least of n
 77: -  c - complex part of computed eigenvalues, provided by user with a dimension at least of n

 79:    Notes:
 80:    This approach is very slow but will generally provide accurate eigenvalue
 81:    estimates.  This routine explicitly forms a dense matrix representing
 82:    the preconditioned operator, and thus will run only for relatively small
 83:    problems, say n < 500.

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

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

 93:    Level: advanced

 95: .seealso: KSPComputeEigenvalues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSPSetOperators(), KSPSolve()
 96: @*/
 97: PetscErrorCode  KSPComputeEigenvaluesExplicitly(KSP ksp,PetscInt nmax,PetscReal r[],PetscReal c[])
 98: {
 99:   Mat               BA;
100:   PetscErrorCode    ierr;
101:   PetscMPIInt       size,rank;
102:   MPI_Comm          comm;
103:   PetscScalar       *array;
104:   Mat               A;
105:   PetscInt          m,row,nz,i,n,dummy;
106:   const PetscInt    *cols;
107:   const PetscScalar *vals;

110:   PetscObjectGetComm((PetscObject)ksp,&comm);
111:   KSPComputeOperator(ksp,MATDENSE,&BA);
112:   MPI_Comm_size(comm,&size);
113:   MPI_Comm_rank(comm,&rank);

115:   MatGetSize(BA,&n,&n);
116:   if (size > 1) { /* assemble matrix on first processor */
117:     MatCreate(PetscObjectComm((PetscObject)ksp),&A);
118:     if (!rank) {
119:       MatSetSizes(A,n,n,n,n);
120:     } else {
121:       MatSetSizes(A,0,0,n,n);
122:     }
123:     MatSetType(A,MATMPIDENSE);
124:     MatMPIDenseSetPreallocation(A,NULL);
125:     PetscLogObjectParent((PetscObject)BA,(PetscObject)A);

127:     MatGetOwnershipRange(BA,&row,&dummy);
128:     MatGetLocalSize(BA,&m,&dummy);
129:     for (i=0; i<m; i++) {
130:       MatGetRow(BA,row,&nz,&cols,&vals);
131:       MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
132:       MatRestoreRow(BA,row,&nz,&cols,&vals);
133:       row++;
134:     }

136:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
137:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
138:     MatDenseGetArray(A,&array);
139:   } else {
140:     MatDenseGetArray(BA,&array);
141:   }

143: #if defined(PETSC_HAVE_ESSL)
144:   /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */
145:   if (!rank) {
146:     PetscScalar  sdummy,*cwork;
147:     PetscReal    *work,*realpart;
148:     PetscBLASInt clen,idummy,lwork,bn,zero = 0;
149:     PetscInt     *perm;

151: #if !defined(PETSC_USE_COMPLEX)
152:     clen = n;
153: #else
154:     clen = 2*n;
155: #endif
156:     PetscMalloc1(clen,&cwork);
157:     idummy = -1;                /* unused */
158:     PetscBLASIntCast(n,&bn);
159:     lwork  = 5*n;
160:     PetscMalloc1(lwork,&work);
161:     PetscMalloc1(n,&realpart);
162:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
163:     PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_(&zero,array,&bn,cwork,&sdummy,&idummy,&idummy,&bn,work,&lwork));
164:     PetscFPTrapPop();
165:     PetscFree(work);

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

171: #if !defined(PETSC_USE_COMPLEX)
172:     for (i=0; i<n; i++) {
173:       realpart[i] = cwork[2*i];
174:       perm[i]     = i;
175:     }
176:     PetscSortRealWithPermutation(n,realpart,perm);
177:     for (i=0; i<n; i++) {
178:       r[i] = cwork[2*perm[i]];
179:       c[i] = cwork[2*perm[i]+1];
180:     }
181: #else
182:     for (i=0; i<n; i++) {
183:       realpart[i] = PetscRealPart(cwork[i]);
184:       perm[i]     = i;
185:     }
186:     PetscSortRealWithPermutation(n,realpart,perm);
187:     for (i=0; i<n; i++) {
188:       r[i] = PetscRealPart(cwork[perm[i]]);
189:       c[i] = PetscImaginaryPart(cwork[perm[i]]);
190:     }
191: #endif
192:     PetscFree(perm);
193:     PetscFree(realpart);
194:     PetscFree(cwork);
195:   }
196: #elif !defined(PETSC_USE_COMPLEX)
197:   if (!rank) {
198:     PetscScalar  *work;
199:     PetscReal    *realpart,*imagpart;
200:     PetscBLASInt idummy,lwork;
201:     PetscInt     *perm;

203:     idummy   = n;
204:     lwork    = 5*n;
205:     PetscMalloc2(n,&realpart,n,&imagpart);
206:     PetscMalloc1(5*n,&work);
207: #if defined(PETSC_MISSING_LAPACK_GEEV)
208:     SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
209: #else
210:     {
211:       PetscBLASInt lierr;
212:       PetscScalar  sdummy;
213:       PetscBLASInt bn;

215:       PetscBLASIntCast(n,&bn);
216:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
217:       PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,array,&bn,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr));
218:       if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
219:       PetscFPTrapPop();
220:     }
221: #endif
222:     PetscFree(work);
223:     PetscMalloc1(n,&perm);

225:     for (i=0; i<n; i++)  perm[i] = i;
226:     PetscSortRealWithPermutation(n,realpart,perm);
227:     for (i=0; i<n; i++) {
228:       r[i] = realpart[perm[i]];
229:       c[i] = imagpart[perm[i]];
230:     }
231:     PetscFree(perm);
232:     PetscFree2(realpart,imagpart);
233:   }
234: #else
235:   if (!rank) {
236:     PetscScalar  *work,*eigs;
237:     PetscReal    *rwork;
238:     PetscBLASInt idummy,lwork;
239:     PetscInt     *perm;

241:     idummy = n;
242:     lwork  = 5*n;
243:     PetscMalloc1(5*n,&work);
244:     PetscMalloc1(2*n,&rwork);
245:     PetscMalloc1(n,&eigs);
246: #if defined(PETSC_MISSING_LAPACK_GEEV)
247:     SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
248: #else
249:     {
250:       PetscBLASInt lierr;
251:       PetscScalar  sdummy;
252:       PetscBLASInt nb;
253:       PetscBLASIntCast(n,&nb);
254:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
255:       PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&nb,array,&nb,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&lierr));
256:       if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
257:       PetscFPTrapPop();
258:     }
259: #endif
260:     PetscFree(work);
261:     PetscFree(rwork);
262:     PetscMalloc1(n,&perm);
263:     for (i=0; i<n; i++) perm[i] = i;
264:     for (i=0; i<n; i++) r[i]    = PetscRealPart(eigs[i]);
265:     PetscSortRealWithPermutation(n,r,perm);
266:     for (i=0; i<n; i++) {
267:       r[i] = PetscRealPart(eigs[perm[i]]);
268:       c[i] = PetscImaginaryPart(eigs[perm[i]]);
269:     }
270:     PetscFree(perm);
271:     PetscFree(eigs);
272:   }
273: #endif
274:   if (size > 1) {
275:     MatDenseRestoreArray(A,&array);
276:     MatDestroy(&A);
277:   } else {
278:     MatDenseRestoreArray(BA,&array);
279:   }
280:   MatDestroy(&BA);
281:   return(0);
282: }

284: static PetscErrorCode PolyEval(PetscInt nroots,const PetscReal *r,const PetscReal *c,PetscReal x,PetscReal y,PetscReal *px,PetscReal *py)
285: {
286:   PetscInt  i;
287:   PetscReal rprod = 1,iprod = 0;

290:   for (i=0; i<nroots; i++) {
291:     PetscReal rnew = rprod*(x - r[i]) - iprod*(y - c[i]);
292:     PetscReal inew = rprod*(y - c[i]) + iprod*(x - r[i]);
293:     rprod = rnew;
294:     iprod = inew;
295:   }
296:   *px = rprod;
297:   *py = iprod;
298:   return(0);
299: }

301:  #include <petscdraw.h>
302: /* collective on ksp */
303: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp,PetscInt neig,const PetscReal *r,const PetscReal *c)
304: {
306:   PetscReal      xmin,xmax,ymin,ymax,*xloc,*yloc,*value,px0,py0,rscale,iscale;
307:   PetscInt       M,N,i,j;
308:   PetscMPIInt    rank;
309:   PetscViewer    viewer;
310:   PetscDraw      draw;
311:   PetscDrawAxis  drawaxis;

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