Actual source code: eige.c

petsc-3.5.1 2014-07-24
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.

 10:     Collective on KSP

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

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

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

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

 26:     Level: advanced

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

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

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

 48:   VecDuplicate(ksp->vec_sol,&in);
 49:   VecDuplicate(ksp->vec_sol,&out);
 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_MatMult(ksp,A,in,out);
 77:     KSP_PCApply(ksp,out,in);

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

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

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

 98:    Collective on KSP

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

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

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

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

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

122:    Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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