Actual source code: eige.c

petsc-3.7.4 2016-10-02
Report Typos and Errors
  2: #include <petsc/private/kspimpl.h>   /*I "petscksp.h" I*/
  3: #include <petscblaslapack.h>

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

 11:     Collective on KSP

 13:     Input Parameter:
 14: .   ksp - the Krylov subspace context

 16:     Output Parameter:
 17: .   mat - the explict preconditioned operator

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

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

 27:     Level: advanced

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

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

 46:   PetscObjectGetComm((PetscObject)ksp,&comm);
 47:   MPI_Comm_size(comm,&size);

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

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

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

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

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

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

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

 96: /*@
 97:    KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
 98:    preconditioned operator using LAPACK.

100:    Collective on KSP

102:    Input Parameter:
103: +  ksp - iterative context obtained from KSPCreate()
104: -  n - size of arrays r and c

106:    Output Parameters:
107: +  r - real part of computed eigenvalues, provided by user with a dimension at least of n
108: -  c - complex part of computed eigenvalues, provided by user with a dimension at least of n

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

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

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

124:    Level: advanced

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

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

143:   PetscObjectGetComm((PetscObject)ksp,&comm);
144:   KSPComputeExplicitOperator(ksp,&BA);
145:   MPI_Comm_size(comm,&size);
146:   MPI_Comm_rank(comm,&rank);

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

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

169:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
170:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
171:     MatDenseGetArray(A,&array);
172:   } else {
173:     MatDenseGetArray(BA,&array);
174:   }

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

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

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

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

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

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

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

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

319: static PetscErrorCode PolyEval(PetscInt nroots,const PetscReal *r,const PetscReal *c,PetscReal x,PetscReal y,PetscReal *px,PetscReal *py)
320: {
321:   PetscInt  i;
322:   PetscReal rprod = 1,iprod = 0;

325:   for (i=0; i<nroots; i++) {
326:     PetscReal rnew = rprod*(x - r[i]) - iprod*(y - c[i]);
327:     PetscReal inew = rprod*(y - c[i]) + iprod*(x - r[i]);
328:     rprod = rnew;
329:     iprod = inew;
330:   }
331:   *px = rprod;
332:   *py = iprod;
333:   return(0);
334: }

336: #include <petscdraw.h>
339: /* collective on KSP */
340: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp,PetscInt neig,const PetscReal *r,const PetscReal *c)
341: {
343:   PetscReal      xmin,xmax,ymin,ymax,*xloc,*yloc,*value,px0,py0,rscale,iscale;
344:   PetscInt       M,N,i,j;
345:   PetscMPIInt    rank;
346:   PetscViewer    viewer;
347:   PetscDraw      draw;
348:   PetscDrawAxis  drawaxis;

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