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;
18: MatShellGetContext(A,&ctx);
19: KSPGetOperators(ctx->ksp,&M,NULL);
20: MatCreateVecs(M,X,Y);
21: return(0);
22: }
24: static PetscErrorCode MatMult_KSP(Mat A,Vec X,Vec Y)
25: {
26: Mat_KSP *ctx;
30: MatShellGetContext(A,&ctx);
31: KSP_PCApplyBAorAB(ctx->ksp,X,Y,ctx->work);
32: return(0);
33: }
35: /*@
36: KSPComputeOperator - Computes the explicit preconditioned operator, including diagonal scaling and null
37: space removal if applicable.
39: Collective on ksp
41: Input Parameter:
42: + ksp - the Krylov subspace context
43: - mattype - the matrix type to be used
45: Output Parameter:
46: . mat - the explict preconditioned operator
48: Notes:
49: This computation is done by applying the operators to columns of the
50: identity matrix.
52: Currently, this routine uses a dense matrix format for the output operator if mattype == NULL.
53: This routine is costly in general, and is recommended for use only with relatively small systems.
55: Level: advanced
57: .seealso: KSPComputeEigenvaluesExplicitly(), PCComputeOperator(), KSPSetDiagonalScale(), KSPSetNullSpace(), MatType
58: @*/
59: PetscErrorCode KSPComputeOperator(KSP ksp, MatType mattype, Mat *mat)
60: {
62: PetscInt N,M,m,n;
63: Mat_KSP ctx;
64: Mat A,Aksp;
69: KSPGetOperators(ksp,&A,NULL);
70: MatGetLocalSize(A,&m,&n);
71: MatGetSize(A,&M,&N);
72: MatCreateShell(PetscObjectComm((PetscObject)ksp),m,n,M,N,&ctx,&Aksp);
73: MatShellSetOperation(Aksp,MATOP_MULT,(void (*)(void))MatMult_KSP);
74: MatShellSetOperation(Aksp,MATOP_CREATE_VECS,(void (*)(void))MatCreateVecs_KSP);
75: ctx.ksp = ksp;
76: MatCreateVecs(A,&ctx.work,NULL);
77: MatComputeOperator(Aksp,mattype,mat);
78: VecDestroy(&ctx.work);
79: MatDestroy(&Aksp);
80: return(0);
81: }
83: /*@
84: KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
85: preconditioned operator using LAPACK.
87: Collective on ksp
89: Input Parameter:
90: + ksp - iterative context obtained from KSPCreate()
91: - n - size of arrays r and c
93: Output Parameters:
94: + r - real part of computed eigenvalues, provided by user with a dimension at least of n
95: - c - complex part of computed eigenvalues, provided by user with a dimension at least of n
97: Notes:
98: This approach is very slow but will generally provide accurate eigenvalue
99: estimates. This routine explicitly forms a dense matrix representing
100: the preconditioned operator, and thus will run only for relatively small
101: problems, say n < 500.
103: Many users may just want to use the monitoring routine
104: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
105: to print the singular values at each iteration of the linear solve.
107: The preconditoner operator, rhs vector, solution vectors should be
108: set before this routine is called. i.e use KSPSetOperators(),KSPSolve() or
109: KSPSetOperators()
111: Level: advanced
113: .seealso: KSPComputeEigenvalues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSPSetOperators(), KSPSolve()
114: @*/
115: PetscErrorCode KSPComputeEigenvaluesExplicitly(KSP ksp,PetscInt nmax,PetscReal r[],PetscReal c[])
116: {
117: Mat BA;
118: PetscErrorCode ierr;
119: PetscMPIInt size,rank;
120: MPI_Comm comm;
121: PetscScalar *array;
122: Mat A;
123: PetscInt m,row,nz,i,n,dummy;
124: const PetscInt *cols;
125: const PetscScalar *vals;
128: PetscObjectGetComm((PetscObject)ksp,&comm);
129: KSPComputeOperator(ksp,MATDENSE,&BA);
130: MPI_Comm_size(comm,&size);
131: MPI_Comm_rank(comm,&rank);
133: MatGetSize(BA,&n,&n);
134: if (size > 1) { /* assemble matrix on first processor */
135: MatCreate(PetscObjectComm((PetscObject)ksp),&A);
136: if (!rank) {
137: MatSetSizes(A,n,n,n,n);
138: } else {
139: MatSetSizes(A,0,0,n,n);
140: }
141: MatSetType(A,MATMPIDENSE);
142: MatMPIDenseSetPreallocation(A,NULL);
143: PetscLogObjectParent((PetscObject)BA,(PetscObject)A);
145: MatGetOwnershipRange(BA,&row,&dummy);
146: MatGetLocalSize(BA,&m,&dummy);
147: for (i=0; i<m; i++) {
148: MatGetRow(BA,row,&nz,&cols,&vals);
149: MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
150: MatRestoreRow(BA,row,&nz,&cols,&vals);
151: row++;
152: }
154: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
155: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
156: MatDenseGetArray(A,&array);
157: } else {
158: MatDenseGetArray(BA,&array);
159: }
161: #if defined(PETSC_HAVE_ESSL)
162: /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */
163: if (!rank) {
164: PetscScalar sdummy,*cwork;
165: PetscReal *work,*realpart;
166: PetscBLASInt clen,idummy,lwork,bn,zero = 0;
167: PetscInt *perm;
169: #if !defined(PETSC_USE_COMPLEX)
170: clen = n;
171: #else
172: clen = 2*n;
173: #endif
174: PetscMalloc1(clen,&cwork);
175: idummy = -1; /* unused */
176: PetscBLASIntCast(n,&bn);
177: lwork = 5*n;
178: PetscMalloc1(lwork,&work);
179: PetscMalloc1(n,&realpart);
180: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
181: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_(&zero,array,&bn,cwork,&sdummy,&idummy,&idummy,&bn,work,&lwork));
182: PetscFPTrapPop();
183: PetscFree(work);
185: /* For now we stick with the convention of storing the real and imaginary
186: components of evalues separately. But is this what we really want? */
187: PetscMalloc1(n,&perm);
189: #if !defined(PETSC_USE_COMPLEX)
190: for (i=0; i<n; i++) {
191: realpart[i] = cwork[2*i];
192: perm[i] = i;
193: }
194: PetscSortRealWithPermutation(n,realpart,perm);
195: for (i=0; i<n; i++) {
196: r[i] = cwork[2*perm[i]];
197: c[i] = cwork[2*perm[i]+1];
198: }
199: #else
200: for (i=0; i<n; i++) {
201: realpart[i] = PetscRealPart(cwork[i]);
202: perm[i] = i;
203: }
204: PetscSortRealWithPermutation(n,realpart,perm);
205: for (i=0; i<n; i++) {
206: r[i] = PetscRealPart(cwork[perm[i]]);
207: c[i] = PetscImaginaryPart(cwork[perm[i]]);
208: }
209: #endif
210: PetscFree(perm);
211: PetscFree(realpart);
212: PetscFree(cwork);
213: }
214: #elif !defined(PETSC_USE_COMPLEX)
215: if (!rank) {
216: PetscScalar *work;
217: PetscReal *realpart,*imagpart;
218: PetscBLASInt idummy,lwork;
219: PetscInt *perm;
221: idummy = n;
222: lwork = 5*n;
223: PetscMalloc2(n,&realpart,n,&imagpart);
224: PetscMalloc1(5*n,&work);
225: {
226: PetscBLASInt lierr;
227: PetscScalar sdummy;
228: PetscBLASInt bn;
230: PetscBLASIntCast(n,&bn);
231: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
232: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,array,&bn,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr));
233: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
234: PetscFPTrapPop();
235: }
236: PetscFree(work);
237: PetscMalloc1(n,&perm);
239: for (i=0; i<n; i++) perm[i] = i;
240: PetscSortRealWithPermutation(n,realpart,perm);
241: for (i=0; i<n; i++) {
242: r[i] = realpart[perm[i]];
243: c[i] = imagpart[perm[i]];
244: }
245: PetscFree(perm);
246: PetscFree2(realpart,imagpart);
247: }
248: #else
249: if (!rank) {
250: PetscScalar *work,*eigs;
251: PetscReal *rwork;
252: PetscBLASInt idummy,lwork;
253: PetscInt *perm;
255: idummy = n;
256: lwork = 5*n;
257: PetscMalloc1(5*n,&work);
258: PetscMalloc1(2*n,&rwork);
259: PetscMalloc1(n,&eigs);
260: {
261: PetscBLASInt lierr;
262: PetscScalar sdummy;
263: PetscBLASInt nb;
264: PetscBLASIntCast(n,&nb);
265: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
266: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&nb,array,&nb,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&lierr));
267: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
268: PetscFPTrapPop();
269: }
270: PetscFree(work);
271: PetscFree(rwork);
272: PetscMalloc1(n,&perm);
273: for (i=0; i<n; i++) perm[i] = i;
274: for (i=0; i<n; i++) r[i] = PetscRealPart(eigs[i]);
275: PetscSortRealWithPermutation(n,r,perm);
276: for (i=0; i<n; i++) {
277: r[i] = PetscRealPart(eigs[perm[i]]);
278: c[i] = PetscImaginaryPart(eigs[perm[i]]);
279: }
280: PetscFree(perm);
281: PetscFree(eigs);
282: }
283: #endif
284: if (size > 1) {
285: MatDenseRestoreArray(A,&array);
286: MatDestroy(&A);
287: } else {
288: MatDenseRestoreArray(BA,&array);
289: }
290: MatDestroy(&BA);
291: return(0);
292: }
294: static PetscErrorCode PolyEval(PetscInt nroots,const PetscReal *r,const PetscReal *c,PetscReal x,PetscReal y,PetscReal *px,PetscReal *py)
295: {
296: PetscInt i;
297: PetscReal rprod = 1,iprod = 0;
300: for (i=0; i<nroots; i++) {
301: PetscReal rnew = rprod*(x - r[i]) - iprod*(y - c[i]);
302: PetscReal inew = rprod*(y - c[i]) + iprod*(x - r[i]);
303: rprod = rnew;
304: iprod = inew;
305: }
306: *px = rprod;
307: *py = iprod;
308: return(0);
309: }
311: #include <petscdraw.h>
312: /* collective on ksp */
313: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp,PetscInt neig,const PetscReal *r,const PetscReal *c)
314: {
316: PetscReal xmin,xmax,ymin,ymax,*xloc,*yloc,*value,px0,py0,rscale,iscale;
317: PetscInt M,N,i,j;
318: PetscMPIInt rank;
319: PetscViewer viewer;
320: PetscDraw draw;
321: PetscDrawAxis drawaxis;
324: MPI_Comm_rank(PetscObjectComm((PetscObject)ksp),&rank);
325: if (rank) return(0);
326: M = 80;
327: N = 80;
328: xmin = r[0]; xmax = r[0];
329: ymin = c[0]; ymax = c[0];
330: for (i=1; i<neig; i++) {
331: xmin = PetscMin(xmin,r[i]);
332: xmax = PetscMax(xmax,r[i]);
333: ymin = PetscMin(ymin,c[i]);
334: ymax = PetscMax(ymax,c[i]);
335: }
336: PetscMalloc3(M,&xloc,N,&yloc,M*N,&value);
337: for (i=0; i<M; i++) xloc[i] = xmin - 0.1*(xmax-xmin) + 1.2*(xmax-xmin)*i/(M-1);
338: for (i=0; i<N; i++) yloc[i] = ymin - 0.1*(ymax-ymin) + 1.2*(ymax-ymin)*i/(N-1);
339: PolyEval(neig,r,c,0,0,&px0,&py0);
340: rscale = px0/(PetscSqr(px0)+PetscSqr(py0));
341: iscale = -py0/(PetscSqr(px0)+PetscSqr(py0));
342: for (j=0; j<N; j++) {
343: for (i=0; i<M; i++) {
344: PetscReal px,py,tx,ty,tmod;
345: PolyEval(neig,r,c,xloc[i],yloc[j],&px,&py);
346: tx = px*rscale - py*iscale;
347: ty = py*rscale + px*iscale;
348: tmod = PetscSqr(tx) + PetscSqr(ty); /* modulus of the complex polynomial */
349: if (tmod > 1) tmod = 1.0;
350: if (tmod > 0.5 && tmod < 1) tmod = 0.5;
351: if (tmod > 0.2 && tmod < 0.5) tmod = 0.2;
352: if (tmod > 0.05 && tmod < 0.2) tmod = 0.05;
353: if (tmod < 1e-3) tmod = 1e-3;
354: value[i+j*M] = PetscLogReal(tmod) / PetscLogReal(10.0);
355: }
356: }
357: PetscViewerDrawOpen(PETSC_COMM_SELF,NULL,"Iteratively Computed Eigen-contours",PETSC_DECIDE,PETSC_DECIDE,450,450,&viewer);
358: PetscViewerDrawGetDraw(viewer,0,&draw);
359: PetscDrawTensorContour(draw,M,N,NULL,NULL,value);
360: if (0) {
361: PetscDrawAxisCreate(draw,&drawaxis);
362: PetscDrawAxisSetLimits(drawaxis,xmin,xmax,ymin,ymax);
363: PetscDrawAxisSetLabels(drawaxis,"Eigen-counters","real","imag");
364: PetscDrawAxisDraw(drawaxis);
365: PetscDrawAxisDestroy(&drawaxis);
366: }
367: PetscViewerDestroy(&viewer);
368: PetscFree3(xloc,yloc,value);
369: return(0);
370: }