Actual source code: eige.c

petsc-3.8.2 2017-11-09
Report Typos and Errors

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

  5: /*@
  6:     KSPComputeExplicitOperator - Computes the explicit preconditioned operator, including diagonal scaling and null
  7:     space removal if applicable.

  9:     Collective on KSP

 11:     Input Parameter:
 12: .   ksp - the Krylov subspace context

 14:     Output Parameter:
 15: .   mat - the explict preconditioned operator

 17:     Notes:
 18:     This computation is done by applying the operators to columns of the
 19:     identity matrix.

 21:     Currently, this routine uses a dense matrix format when 1 processor
 22:     is used and a sparse format otherwise.  This routine is costly in general,
 23:     and is recommended for use only with relatively small systems.

 25:     Level: advanced

 27: .keywords: KSP, compute, explicit, operator

 29: .seealso: KSPComputeEigenvaluesExplicitly(), PCComputeExplicitOperator(), KSPSetDiagonalScale(), KSPSetNullSpace()
 30: @*/
 31: PetscErrorCode  KSPComputeExplicitOperator(KSP ksp,Mat *mat)
 32: {
 33:   Vec            in,out,work;
 35:   PetscMPIInt    size;
 36:   PetscInt       i,M,m,*rows,start,end;
 37:   Mat            A;
 38:   MPI_Comm       comm;
 39:   PetscScalar    *array,one = 1.0;

 44:   PetscObjectGetComm((PetscObject)ksp,&comm);
 45:   MPI_Comm_size(comm,&size);

 47:   VecDuplicate(ksp->vec_sol,&in);
 48:   VecDuplicate(ksp->vec_sol,&out);
 49:   VecDuplicate(ksp->vec_sol,&work);
 50:   VecGetSize(in,&M);
 51:   VecGetLocalSize(in,&m);
 52:   VecGetOwnershipRange(in,&start,&end);
 53:   PetscMalloc1(m,&rows);
 54:   for (i=0; i<m; i++) rows[i] = start + i;

 56:   MatCreate(comm,mat);
 57:   MatSetSizes(*mat,m,m,M,M);
 58:   if (size == 1) {
 59:     MatSetType(*mat,MATSEQDENSE);
 60:     MatSeqDenseSetPreallocation(*mat,NULL);
 61:   } else {
 62:     MatSetType(*mat,MATMPIAIJ);
 63:     MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
 64:   }
 65:   MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
 66:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
 67:   PCGetOperators(ksp->pc,&A,NULL);

 69:   for (i=0; i<M; i++) {

 71:     VecSet(in,0.0);
 72:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
 73:     VecAssemblyBegin(in);
 74:     VecAssemblyEnd(in);

 76:     KSP_PCApplyBAorAB(ksp,in,out,work);

 78:     VecGetArray(out,&array);
 79:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
 80:     VecRestoreArray(out,&array);

 82:   }
 83:   PetscFree(rows);
 84:   VecDestroy(&in);
 85:   VecDestroy(&out);
 86:   VecDestroy(&work);
 87:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
 88:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
 89:   return(0);
 90: }

 92: /*@
 93:    KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
 94:    preconditioned operator using LAPACK.

 96:    Collective on KSP

 98:    Input Parameter:
 99: +  ksp - iterative context obtained from KSPCreate()
100: -  n - size of arrays r and c

102:    Output Parameters:
103: +  r - real part of computed eigenvalues, provided by user with a dimension at least of n
104: -  c - complex part of computed eigenvalues, provided by user with a dimension at least of n

106:    Notes:
107:    This approach is very slow but will generally provide accurate eigenvalue
108:    estimates.  This routine explicitly forms a dense matrix representing
109:    the preconditioned operator, and thus will run only for relatively small
110:    problems, say n < 500.

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

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

120:    Level: advanced

122: .keywords: KSP, compute, eigenvalues, explicitly

124: .seealso: KSPComputeEigenvalues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSPSetOperators(), KSPSolve()
125: @*/
126: PetscErrorCode  KSPComputeEigenvaluesExplicitly(KSP ksp,PetscInt nmax,PetscReal r[],PetscReal c[])
127: {
128:   Mat               BA;
129:   PetscErrorCode    ierr;
130:   PetscMPIInt       size,rank;
131:   MPI_Comm          comm;
132:   PetscScalar       *array;
133:   Mat               A;
134:   PetscInt          m,row,nz,i,n,dummy;
135:   const PetscInt    *cols;
136:   const PetscScalar *vals;

139:   PetscObjectGetComm((PetscObject)ksp,&comm);
140:   KSPComputeExplicitOperator(ksp,&BA);
141:   MPI_Comm_size(comm,&size);
142:   MPI_Comm_rank(comm,&rank);

144:   MatGetSize(BA,&n,&n);
145:   if (size > 1) { /* assemble matrix on first processor */
146:     MatCreate(PetscObjectComm((PetscObject)ksp),&A);
147:     if (!rank) {
148:       MatSetSizes(A,n,n,n,n);
149:     } else {
150:       MatSetSizes(A,0,0,n,n);
151:     }
152:     MatSetType(A,MATMPIDENSE);
153:     MatMPIDenseSetPreallocation(A,NULL);
154:     PetscLogObjectParent((PetscObject)BA,(PetscObject)A);

156:     MatGetOwnershipRange(BA,&row,&dummy);
157:     MatGetLocalSize(BA,&m,&dummy);
158:     for (i=0; i<m; i++) {
159:       MatGetRow(BA,row,&nz,&cols,&vals);
160:       MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
161:       MatRestoreRow(BA,row,&nz,&cols,&vals);
162:       row++;
163:     }

165:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
166:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
167:     MatDenseGetArray(A,&array);
168:   } else {
169:     MatDenseGetArray(BA,&array);
170:   }

172: #if defined(PETSC_HAVE_ESSL)
173:   /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */
174:   if (!rank) {
175:     PetscScalar  sdummy,*cwork;
176:     PetscReal    *work,*realpart;
177:     PetscBLASInt clen,idummy,lwork,bn,zero = 0;
178:     PetscInt     *perm;

180: #if !defined(PETSC_USE_COMPLEX)
181:     clen = n;
182: #else
183:     clen = 2*n;
184: #endif
185:     PetscMalloc1(clen,&cwork);
186:     idummy = -1;                /* unused */
187:     PetscBLASIntCast(n,&bn);
188:     lwork  = 5*n;
189:     PetscMalloc1(lwork,&work);
190:     PetscMalloc1(n,&realpart);
191:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
192:     PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_(&zero,array,&bn,cwork,&sdummy,&idummy,&idummy,&bn,work,&lwork));
193:     PetscFPTrapPop();
194:     PetscFree(work);

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

200: #if !defined(PETSC_USE_COMPLEX)
201:     for (i=0; i<n; i++) {
202:       realpart[i] = cwork[2*i];
203:       perm[i]     = i;
204:     }
205:     PetscSortRealWithPermutation(n,realpart,perm);
206:     for (i=0; i<n; i++) {
207:       r[i] = cwork[2*perm[i]];
208:       c[i] = cwork[2*perm[i]+1];
209:     }
210: #else
211:     for (i=0; i<n; i++) {
212:       realpart[i] = PetscRealPart(cwork[i]);
213:       perm[i]     = i;
214:     }
215:     PetscSortRealWithPermutation(n,realpart,perm);
216:     for (i=0; i<n; i++) {
217:       r[i] = PetscRealPart(cwork[perm[i]]);
218:       c[i] = PetscImaginaryPart(cwork[perm[i]]);
219:     }
220: #endif
221:     PetscFree(perm);
222:     PetscFree(realpart);
223:     PetscFree(cwork);
224:   }
225: #elif !defined(PETSC_USE_COMPLEX)
226:   if (!rank) {
227:     PetscScalar  *work;
228:     PetscReal    *realpart,*imagpart;
229:     PetscBLASInt idummy,lwork;
230:     PetscInt     *perm;

232:     idummy   = n;
233:     lwork    = 5*n;
234:     PetscMalloc2(n,&realpart,n,&imagpart);
235:     PetscMalloc1(5*n,&work);
236: #if defined(PETSC_MISSING_LAPACK_GEEV)
237:     SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
238: #else
239:     {
240:       PetscBLASInt lierr;
241:       PetscScalar  sdummy;
242:       PetscBLASInt bn;

244:       PetscBLASIntCast(n,&bn);
245:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
246:       PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,array,&bn,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr));
247:       if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
248:       PetscFPTrapPop();
249:     }
250: #endif
251:     PetscFree(work);
252:     PetscMalloc1(n,&perm);

254:     for (i=0; i<n; i++)  perm[i] = i;
255:     PetscSortRealWithPermutation(n,realpart,perm);
256:     for (i=0; i<n; i++) {
257:       r[i] = realpart[perm[i]];
258:       c[i] = imagpart[perm[i]];
259:     }
260:     PetscFree(perm);
261:     PetscFree2(realpart,imagpart);
262:   }
263: #else
264:   if (!rank) {
265:     PetscScalar  *work,*eigs;
266:     PetscReal    *rwork;
267:     PetscBLASInt idummy,lwork;
268:     PetscInt     *perm;

270:     idummy = n;
271:     lwork  = 5*n;
272:     PetscMalloc1(5*n,&work);
273:     PetscMalloc1(2*n,&rwork);
274:     PetscMalloc1(n,&eigs);
275: #if defined(PETSC_MISSING_LAPACK_GEEV)
276:     SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
277: #else
278:     {
279:       PetscBLASInt lierr;
280:       PetscScalar  sdummy;
281:       PetscBLASInt nb;
282:       PetscBLASIntCast(n,&nb);
283:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
284:       PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&nb,array,&nb,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&lierr));
285:       if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
286:       PetscFPTrapPop();
287:     }
288: #endif
289:     PetscFree(work);
290:     PetscFree(rwork);
291:     PetscMalloc1(n,&perm);
292:     for (i=0; i<n; i++) perm[i] = i;
293:     for (i=0; i<n; i++) r[i]    = PetscRealPart(eigs[i]);
294:     PetscSortRealWithPermutation(n,r,perm);
295:     for (i=0; i<n; i++) {
296:       r[i] = PetscRealPart(eigs[perm[i]]);
297:       c[i] = PetscImaginaryPart(eigs[perm[i]]);
298:     }
299:     PetscFree(perm);
300:     PetscFree(eigs);
301:   }
302: #endif
303:   if (size > 1) {
304:     MatDenseRestoreArray(A,&array);
305:     MatDestroy(&A);
306:   } else {
307:     MatDenseRestoreArray(BA,&array);
308:   }
309:   MatDestroy(&BA);
310:   return(0);
311: }

313: static PetscErrorCode PolyEval(PetscInt nroots,const PetscReal *r,const PetscReal *c,PetscReal x,PetscReal y,PetscReal *px,PetscReal *py)
314: {
315:   PetscInt  i;
316:   PetscReal rprod = 1,iprod = 0;

319:   for (i=0; i<nroots; i++) {
320:     PetscReal rnew = rprod*(x - r[i]) - iprod*(y - c[i]);
321:     PetscReal inew = rprod*(y - c[i]) + iprod*(x - r[i]);
322:     rprod = rnew;
323:     iprod = inew;
324:   }
325:   *px = rprod;
326:   *py = iprod;
327:   return(0);
328: }

330:  #include <petscdraw.h>
331: /* collective on KSP */
332: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp,PetscInt neig,const PetscReal *r,const PetscReal *c)
333: {
335:   PetscReal      xmin,xmax,ymin,ymax,*xloc,*yloc,*value,px0,py0,rscale,iscale;
336:   PetscInt       M,N,i,j;
337:   PetscMPIInt    rank;
338:   PetscViewer    viewer;
339:   PetscDraw      draw;
340:   PetscDrawAxis  drawaxis;

343:   MPI_Comm_rank(PetscObjectComm((PetscObject)ksp),&rank);
344:   if (rank) return(0);
345:   M    = 80;
346:   N    = 80;
347:   xmin = r[0]; xmax = r[0];
348:   ymin = c[0]; ymax = c[0];
349:   for (i=1; i<neig; i++) {
350:     xmin = PetscMin(xmin,r[i]);
351:     xmax = PetscMax(xmax,r[i]);
352:     ymin = PetscMin(ymin,c[i]);
353:     ymax = PetscMax(ymax,c[i]);
354:   }
355:   PetscMalloc3(M,&xloc,N,&yloc,M*N,&value);
356:   for (i=0; i<M; i++) xloc[i] = xmin - 0.1*(xmax-xmin) + 1.2*(xmax-xmin)*i/(M-1);
357:   for (i=0; i<N; i++) yloc[i] = ymin - 0.1*(ymax-ymin) + 1.2*(ymax-ymin)*i/(N-1);
358:   PolyEval(neig,r,c,0,0,&px0,&py0);
359:   rscale = px0/(PetscSqr(px0)+PetscSqr(py0));
360:   iscale = -py0/(PetscSqr(px0)+PetscSqr(py0));
361:   for (j=0; j<N; j++) {
362:     for (i=0; i<M; i++) {
363:       PetscReal px,py,tx,ty,tmod;
364:       PolyEval(neig,r,c,xloc[i],yloc[j],&px,&py);
365:       tx   = px*rscale - py*iscale;
366:       ty   = py*rscale + px*iscale;
367:       tmod = PetscSqr(tx) + PetscSqr(ty); /* modulus of the complex polynomial */
368:       if (tmod > 1) tmod = 1.0;
369:       if (tmod > 0.5 && tmod < 1) tmod = 0.5;
370:       if (tmod > 0.2 && tmod < 0.5) tmod = 0.2;
371:       if (tmod > 0.05 && tmod < 0.2) tmod = 0.05;
372:       if (tmod < 1e-3) tmod = 1e-3;
373:       value[i+j*M] = PetscLogReal(tmod) / PetscLogReal(10.0);
374:     }
375:   }
376:   PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Iteratively Computed Eigen-contours",PETSC_DECIDE,PETSC_DECIDE,450,450,&viewer);
377:   PetscViewerDrawGetDraw(viewer,0,&draw);
378:   PetscDrawTensorContour(draw,M,N,NULL,NULL,value);
379:   if (0) {
380:     PetscDrawAxisCreate(draw,&drawaxis);
381:     PetscDrawAxisSetLimits(drawaxis,xmin,xmax,ymin,ymax);
382:     PetscDrawAxisSetLabels(drawaxis,"Eigen-counters","real","imag");
383:     PetscDrawAxisDraw(drawaxis);
384:     PetscDrawAxisDestroy(&drawaxis);
385:   }
386:   PetscViewerDestroy(&viewer);
387:   PetscFree3(xloc,yloc,value);
388:   return(0);
389: }