Actual source code: eige.c


  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;

 16:   MatShellGetContext(A,&ctx);
 17:   KSPGetOperators(ctx->ksp,&M,NULL);
 18:   MatCreateVecs(M,X,Y);
 19:   return 0;
 20: }

 22: static PetscErrorCode MatMult_KSP(Mat A,Vec X,Vec Y)
 23: {
 24:   Mat_KSP        *ctx;

 26:   MatShellGetContext(A,&ctx);
 27:   KSP_PCApplyBAorAB(ctx->ksp,X,Y,ctx->work);
 28:   return 0;
 29: }

 31: /*@
 32:     KSPComputeOperator - Computes the explicit preconditioned operator, including diagonal scaling and null
 33:     space removal if applicable.

 35:     Collective on ksp

 37:     Input Parameters:
 38: +   ksp - the Krylov subspace context
 39: -   mattype - the matrix type to be used

 41:     Output Parameter:
 42: .   mat - the explicit preconditioned operator

 44:     Notes:
 45:     This computation is done by applying the operators to columns of the
 46:     identity matrix.

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

 51:     Level: advanced

 53: .seealso: KSPComputeEigenvaluesExplicitly(), PCComputeOperator(), KSPSetDiagonalScale(), KSPSetNullSpace(), MatType
 54: @*/
 55: PetscErrorCode  KSPComputeOperator(KSP ksp, MatType mattype, Mat *mat)
 56: {
 57:   PetscInt       N,M,m,n;
 58:   Mat_KSP        ctx;
 59:   Mat            A,Aksp;

 63:   KSPGetOperators(ksp,&A,NULL);
 64:   MatGetLocalSize(A,&m,&n);
 65:   MatGetSize(A,&M,&N);
 66:   MatCreateShell(PetscObjectComm((PetscObject)ksp),m,n,M,N,&ctx,&Aksp);
 67:   MatShellSetOperation(Aksp,MATOP_MULT,(void (*)(void))MatMult_KSP);
 68:   MatShellSetOperation(Aksp,MATOP_CREATE_VECS,(void (*)(void))MatCreateVecs_KSP);
 69:   ctx.ksp = ksp;
 70:   MatCreateVecs(A,&ctx.work,NULL);
 71:   MatComputeOperator(Aksp,mattype,mat);
 72:   VecDestroy(&ctx.work);
 73:   MatDestroy(&Aksp);
 74:   return 0;
 75: }

 77: /*@
 78:    KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
 79:    preconditioned operator using LAPACK.

 81:    Collective on ksp

 83:    Input Parameters:
 84: +  ksp - iterative context obtained from KSPCreate()
 85: -  n - size of arrays r and c

 87:    Output Parameters:
 88: +  r - real part of computed eigenvalues, provided by user with a dimension at least of n
 89: -  c - complex part of computed eigenvalues, provided by user with a dimension at least of n

 91:    Notes:
 92:    This approach is very slow but will generally provide accurate eigenvalue
 93:    estimates.  This routine explicitly forms a dense matrix representing
 94:    the preconditioned operator, and thus will run only for relatively small
 95:    problems, say n < 500.

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

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

105:    Level: advanced

107: .seealso: KSPComputeEigenvalues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSPSetOperators(), KSPSolve()
108: @*/
109: PetscErrorCode  KSPComputeEigenvaluesExplicitly(KSP ksp,PetscInt nmax,PetscReal r[],PetscReal c[])
110: {
111:   Mat               BA;
112:   PetscMPIInt       size,rank;
113:   MPI_Comm          comm;
114:   PetscScalar       *array;
115:   Mat               A;
116:   PetscInt          m,row,nz,i,n,dummy;
117:   const PetscInt    *cols;
118:   const PetscScalar *vals;

120:   PetscObjectGetComm((PetscObject)ksp,&comm);
121:   KSPComputeOperator(ksp,MATDENSE,&BA);
122:   MPI_Comm_size(comm,&size);
123:   MPI_Comm_rank(comm,&rank);

125:   MatGetSize(BA,&n,&n);
126:   if (size > 1) { /* assemble matrix on first processor */
127:     MatCreate(PetscObjectComm((PetscObject)ksp),&A);
128:     if (rank == 0) {
129:       MatSetSizes(A,n,n,n,n);
130:     } else {
131:       MatSetSizes(A,0,0,n,n);
132:     }
133:     MatSetType(A,MATMPIDENSE);
134:     MatMPIDenseSetPreallocation(A,NULL);
135:     PetscLogObjectParent((PetscObject)BA,(PetscObject)A);

137:     MatGetOwnershipRange(BA,&row,&dummy);
138:     MatGetLocalSize(BA,&m,&dummy);
139:     for (i=0; i<m; i++) {
140:       MatGetRow(BA,row,&nz,&cols,&vals);
141:       MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
142:       MatRestoreRow(BA,row,&nz,&cols,&vals);
143:       row++;
144:     }

146:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
147:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
148:     MatDenseGetArray(A,&array);
149:   } else {
150:     MatDenseGetArray(BA,&array);
151:   }

153: #if !defined(PETSC_USE_COMPLEX)
154:   if (rank == 0) {
155:     PetscScalar  *work;
156:     PetscReal    *realpart,*imagpart;
157:     PetscBLASInt idummy,lwork;
158:     PetscInt     *perm;

160:     idummy   = n;
161:     lwork    = 5*n;
162:     PetscMalloc2(n,&realpart,n,&imagpart);
163:     PetscMalloc1(5*n,&work);
164:     {
165:       PetscBLASInt lierr;
166:       PetscScalar  sdummy;
167:       PetscBLASInt bn;

169:       PetscBLASIntCast(n,&bn);
170:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
171:       PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,array,&bn,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr));
173:       PetscFPTrapPop();
174:     }
175:     PetscFree(work);
176:     PetscMalloc1(n,&perm);

178:     for (i=0; i<n; i++)  perm[i] = i;
179:     PetscSortRealWithPermutation(n,realpart,perm);
180:     for (i=0; i<n; i++) {
181:       r[i] = realpart[perm[i]];
182:       c[i] = imagpart[perm[i]];
183:     }
184:     PetscFree(perm);
185:     PetscFree2(realpart,imagpart);
186:   }
187: #else
188:   if (rank == 0) {
189:     PetscScalar  *work,*eigs;
190:     PetscReal    *rwork;
191:     PetscBLASInt idummy,lwork;
192:     PetscInt     *perm;

194:     idummy = n;
195:     lwork  = 5*n;
196:     PetscMalloc1(5*n,&work);
197:     PetscMalloc1(2*n,&rwork);
198:     PetscMalloc1(n,&eigs);
199:     {
200:       PetscBLASInt lierr;
201:       PetscScalar  sdummy;
202:       PetscBLASInt nb;
203:       PetscBLASIntCast(n,&nb);
204:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
205:       PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&nb,array,&nb,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&lierr));
207:       PetscFPTrapPop();
208:     }
209:     PetscFree(work);
210:     PetscFree(rwork);
211:     PetscMalloc1(n,&perm);
212:     for (i=0; i<n; i++) perm[i] = i;
213:     for (i=0; i<n; i++) r[i]    = PetscRealPart(eigs[i]);
214:     PetscSortRealWithPermutation(n,r,perm);
215:     for (i=0; i<n; i++) {
216:       r[i] = PetscRealPart(eigs[perm[i]]);
217:       c[i] = PetscImaginaryPart(eigs[perm[i]]);
218:     }
219:     PetscFree(perm);
220:     PetscFree(eigs);
221:   }
222: #endif
223:   if (size > 1) {
224:     MatDenseRestoreArray(A,&array);
225:     MatDestroy(&A);
226:   } else {
227:     MatDenseRestoreArray(BA,&array);
228:   }
229:   MatDestroy(&BA);
230:   return 0;
231: }

233: static PetscErrorCode PolyEval(PetscInt nroots,const PetscReal *r,const PetscReal *c,PetscReal x,PetscReal y,PetscReal *px,PetscReal *py)
234: {
235:   PetscInt  i;
236:   PetscReal rprod = 1,iprod = 0;

238:   for (i=0; i<nroots; i++) {
239:     PetscReal rnew = rprod*(x - r[i]) - iprod*(y - c[i]);
240:     PetscReal inew = rprod*(y - c[i]) + iprod*(x - r[i]);
241:     rprod = rnew;
242:     iprod = inew;
243:   }
244:   *px = rprod;
245:   *py = iprod;
246:   return 0;
247: }

249: #include <petscdraw.h>
250: /* collective on ksp */
251: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp,PetscInt neig,const PetscReal *r,const PetscReal *c)
252: {
253:   PetscReal      xmin,xmax,ymin,ymax,*xloc,*yloc,*value,px0,py0,rscale,iscale;
254:   PetscInt       M,N,i,j;
255:   PetscMPIInt    rank;
256:   PetscViewer    viewer;
257:   PetscDraw      draw;
258:   PetscDrawAxis  drawaxis;

260:   MPI_Comm_rank(PetscObjectComm((PetscObject)ksp),&rank);
261:   if (rank) return 0;
262:   M    = 80;
263:   N    = 80;
264:   xmin = r[0]; xmax = r[0];
265:   ymin = c[0]; ymax = c[0];
266:   for (i=1; i<neig; i++) {
267:     xmin = PetscMin(xmin,r[i]);
268:     xmax = PetscMax(xmax,r[i]);
269:     ymin = PetscMin(ymin,c[i]);
270:     ymax = PetscMax(ymax,c[i]);
271:   }
272:   PetscMalloc3(M,&xloc,N,&yloc,M*N,&value);
273:   for (i=0; i<M; i++) xloc[i] = xmin - 0.1*(xmax-xmin) + 1.2*(xmax-xmin)*i/(M-1);
274:   for (i=0; i<N; i++) yloc[i] = ymin - 0.1*(ymax-ymin) + 1.2*(ymax-ymin)*i/(N-1);
275:   PolyEval(neig,r,c,0,0,&px0,&py0);
276:   rscale = px0/(PetscSqr(px0)+PetscSqr(py0));
277:   iscale = -py0/(PetscSqr(px0)+PetscSqr(py0));
278:   for (j=0; j<N; j++) {
279:     for (i=0; i<M; i++) {
280:       PetscReal px,py,tx,ty,tmod;
281:       PolyEval(neig,r,c,xloc[i],yloc[j],&px,&py);
282:       tx   = px*rscale - py*iscale;
283:       ty   = py*rscale + px*iscale;
284:       tmod = PetscSqr(tx) + PetscSqr(ty); /* modulus of the complex polynomial */
285:       if (tmod > 1) tmod = 1.0;
286:       if (tmod > 0.5 && tmod < 1) tmod = 0.5;
287:       if (tmod > 0.2 && tmod < 0.5) tmod = 0.2;
288:       if (tmod > 0.05 && tmod < 0.2) tmod = 0.05;
289:       if (tmod < 1e-3) tmod = 1e-3;
290:       value[i+j*M] = PetscLogReal(tmod) / PetscLogReal(10.0);
291:     }
292:   }
293:   PetscViewerDrawOpen(PETSC_COMM_SELF,NULL,"Iteratively Computed Eigen-contours",PETSC_DECIDE,PETSC_DECIDE,450,450,&viewer);
294:   PetscViewerDrawGetDraw(viewer,0,&draw);
295:   PetscDrawTensorContour(draw,M,N,NULL,NULL,value);
296:   if (0) {
297:     PetscDrawAxisCreate(draw,&drawaxis);
298:     PetscDrawAxisSetLimits(drawaxis,xmin,xmax,ymin,ymax);
299:     PetscDrawAxisSetLabels(drawaxis,"Eigen-counters","real","imag");
300:     PetscDrawAxisDraw(drawaxis);
301:     PetscDrawAxisDestroy(&drawaxis);
302:   }
303:   PetscViewerDestroy(&viewer);
304:   PetscFree3(xloc,yloc,value);
305:   return 0;
306: }