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