Actual source code: eige.c
petsc-3.11.3 2019-06-26
2: #include <petsc/private/kspimpl.h>
3: #include <petscblaslapack.h>
5: /*@
6: KSPComputeExplicitOperator - Computes the explicit preconditioned operator, including diagonal scaling and null
7: space removal if applicable.
9: Collective on KSP
11: Input Parameter:
12: . ksp - the Krylov subspace context
14: Output Parameter:
15: . mat - the explict preconditioned operator
17: Notes:
18: This computation is done by applying the operators to columns of the
19: identity matrix.
21: Currently, this routine uses a dense matrix format when 1 processor
22: is used and a sparse format otherwise. This routine is costly in general,
23: and is recommended for use only with relatively small systems.
25: Level: advanced
27: .keywords: KSP, compute, explicit, operator
29: .seealso: KSPComputeEigenvaluesExplicitly(), PCComputeExplicitOperator(), KSPSetDiagonalScale(), KSPSetNullSpace()
30: @*/
31: PetscErrorCode KSPComputeExplicitOperator(KSP ksp,Mat *mat)
32: {
33: Vec in,out,work;
35: PetscMPIInt size;
36: PetscInt i,M,m,*rows,start,end;
37: Mat A;
38: MPI_Comm comm;
39: PetscScalar *array,one = 1.0;
44: PetscObjectGetComm((PetscObject)ksp,&comm);
45: MPI_Comm_size(comm,&size);
47: VecDuplicate(ksp->vec_sol,&in);
48: VecDuplicate(ksp->vec_sol,&out);
49: VecDuplicate(ksp->vec_sol,&work);
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_PCApplyBAorAB(ksp,in,out,work);
78: VecGetArray(out,&array);
79: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
80: VecRestoreArray(out,&array);
82: }
83: PetscFree(rows);
84: VecDestroy(&in);
85: VecDestroy(&out);
86: VecDestroy(&work);
87: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
89: return(0);
90: }
92: /*@
93: KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
94: preconditioned operator using LAPACK.
96: Collective on KSP
98: Input Parameter:
99: + ksp - iterative context obtained from KSPCreate()
100: - n - size of arrays r and c
102: Output Parameters:
103: + r - real part of computed eigenvalues, provided by user with a dimension at least of n
104: - c - complex part of computed eigenvalues, provided by user with a dimension at least of n
106: Notes:
107: This approach is very slow but will generally provide accurate eigenvalue
108: estimates. This routine explicitly forms a dense matrix representing
109: the preconditioned operator, and thus will run only for relatively small
110: problems, say n < 500.
112: Many users may just want to use the monitoring routine
113: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
114: to print the singular values at each iteration of the linear solve.
116: The preconditoner operator, rhs vector, solution vectors should be
117: set before this routine is called. i.e use KSPSetOperators(),KSPSolve() or
118: KSPSetOperators()
120: Level: advanced
122: .keywords: KSP, compute, eigenvalues, explicitly
124: .seealso: KSPComputeEigenvalues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSPSetOperators(), KSPSolve()
125: @*/
126: PetscErrorCode KSPComputeEigenvaluesExplicitly(KSP ksp,PetscInt nmax,PetscReal r[],PetscReal c[])
127: {
128: Mat BA;
129: PetscErrorCode ierr;
130: PetscMPIInt size,rank;
131: MPI_Comm comm;
132: PetscScalar *array;
133: Mat A;
134: PetscInt m,row,nz,i,n,dummy;
135: const PetscInt *cols;
136: const PetscScalar *vals;
139: PetscObjectGetComm((PetscObject)ksp,&comm);
140: KSPComputeExplicitOperator(ksp,&BA);
141: MPI_Comm_size(comm,&size);
142: MPI_Comm_rank(comm,&rank);
144: MatGetSize(BA,&n,&n);
145: if (size > 1) { /* assemble matrix on first processor */
146: MatCreate(PetscObjectComm((PetscObject)ksp),&A);
147: if (!rank) {
148: MatSetSizes(A,n,n,n,n);
149: } else {
150: MatSetSizes(A,0,0,n,n);
151: }
152: MatSetType(A,MATMPIDENSE);
153: MatMPIDenseSetPreallocation(A,NULL);
154: PetscLogObjectParent((PetscObject)BA,(PetscObject)A);
156: MatGetOwnershipRange(BA,&row,&dummy);
157: MatGetLocalSize(BA,&m,&dummy);
158: for (i=0; i<m; i++) {
159: MatGetRow(BA,row,&nz,&cols,&vals);
160: MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
161: MatRestoreRow(BA,row,&nz,&cols,&vals);
162: row++;
163: }
165: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
166: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
167: MatDenseGetArray(A,&array);
168: } else {
169: MatDenseGetArray(BA,&array);
170: }
172: #if defined(PETSC_HAVE_ESSL)
173: /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */
174: if (!rank) {
175: PetscScalar sdummy,*cwork;
176: PetscReal *work,*realpart;
177: PetscBLASInt clen,idummy,lwork,bn,zero = 0;
178: PetscInt *perm;
180: #if !defined(PETSC_USE_COMPLEX)
181: clen = n;
182: #else
183: clen = 2*n;
184: #endif
185: PetscMalloc1(clen,&cwork);
186: idummy = -1; /* unused */
187: PetscBLASIntCast(n,&bn);
188: lwork = 5*n;
189: PetscMalloc1(lwork,&work);
190: PetscMalloc1(n,&realpart);
191: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
192: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_(&zero,array,&bn,cwork,&sdummy,&idummy,&idummy,&bn,work,&lwork));
193: PetscFPTrapPop();
194: PetscFree(work);
196: /* For now we stick with the convention of storing the real and imaginary
197: components of evalues separately. But is this what we really want? */
198: PetscMalloc1(n,&perm);
200: #if !defined(PETSC_USE_COMPLEX)
201: for (i=0; i<n; i++) {
202: realpart[i] = cwork[2*i];
203: perm[i] = i;
204: }
205: PetscSortRealWithPermutation(n,realpart,perm);
206: for (i=0; i<n; i++) {
207: r[i] = cwork[2*perm[i]];
208: c[i] = cwork[2*perm[i]+1];
209: }
210: #else
211: for (i=0; i<n; i++) {
212: realpart[i] = PetscRealPart(cwork[i]);
213: perm[i] = i;
214: }
215: PetscSortRealWithPermutation(n,realpart,perm);
216: for (i=0; i<n; i++) {
217: r[i] = PetscRealPart(cwork[perm[i]]);
218: c[i] = PetscImaginaryPart(cwork[perm[i]]);
219: }
220: #endif
221: PetscFree(perm);
222: PetscFree(realpart);
223: PetscFree(cwork);
224: }
225: #elif !defined(PETSC_USE_COMPLEX)
226: if (!rank) {
227: PetscScalar *work;
228: PetscReal *realpart,*imagpart;
229: PetscBLASInt idummy,lwork;
230: PetscInt *perm;
232: idummy = n;
233: lwork = 5*n;
234: PetscMalloc2(n,&realpart,n,&imagpart);
235: PetscMalloc1(5*n,&work);
236: #if defined(PETSC_MISSING_LAPACK_GEEV)
237: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
238: #else
239: {
240: PetscBLASInt lierr;
241: PetscScalar sdummy;
242: PetscBLASInt bn;
244: PetscBLASIntCast(n,&bn);
245: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
246: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,array,&bn,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr));
247: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
248: PetscFPTrapPop();
249: }
250: #endif
251: PetscFree(work);
252: PetscMalloc1(n,&perm);
254: for (i=0; i<n; i++) perm[i] = i;
255: PetscSortRealWithPermutation(n,realpart,perm);
256: for (i=0; i<n; i++) {
257: r[i] = realpart[perm[i]];
258: c[i] = imagpart[perm[i]];
259: }
260: PetscFree(perm);
261: PetscFree2(realpart,imagpart);
262: }
263: #else
264: if (!rank) {
265: PetscScalar *work,*eigs;
266: PetscReal *rwork;
267: PetscBLASInt idummy,lwork;
268: PetscInt *perm;
270: idummy = n;
271: lwork = 5*n;
272: PetscMalloc1(5*n,&work);
273: PetscMalloc1(2*n,&rwork);
274: PetscMalloc1(n,&eigs);
275: #if defined(PETSC_MISSING_LAPACK_GEEV)
276: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
277: #else
278: {
279: PetscBLASInt lierr;
280: PetscScalar sdummy;
281: PetscBLASInt nb;
282: PetscBLASIntCast(n,&nb);
283: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
284: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&nb,array,&nb,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&lierr));
285: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
286: PetscFPTrapPop();
287: }
288: #endif
289: PetscFree(work);
290: PetscFree(rwork);
291: PetscMalloc1(n,&perm);
292: for (i=0; i<n; i++) perm[i] = i;
293: for (i=0; i<n; i++) r[i] = PetscRealPart(eigs[i]);
294: PetscSortRealWithPermutation(n,r,perm);
295: for (i=0; i<n; i++) {
296: r[i] = PetscRealPart(eigs[perm[i]]);
297: c[i] = PetscImaginaryPart(eigs[perm[i]]);
298: }
299: PetscFree(perm);
300: PetscFree(eigs);
301: }
302: #endif
303: if (size > 1) {
304: MatDenseRestoreArray(A,&array);
305: MatDestroy(&A);
306: } else {
307: MatDenseRestoreArray(BA,&array);
308: }
309: MatDestroy(&BA);
310: return(0);
311: }
313: static PetscErrorCode PolyEval(PetscInt nroots,const PetscReal *r,const PetscReal *c,PetscReal x,PetscReal y,PetscReal *px,PetscReal *py)
314: {
315: PetscInt i;
316: PetscReal rprod = 1,iprod = 0;
319: for (i=0; i<nroots; i++) {
320: PetscReal rnew = rprod*(x - r[i]) - iprod*(y - c[i]);
321: PetscReal inew = rprod*(y - c[i]) + iprod*(x - r[i]);
322: rprod = rnew;
323: iprod = inew;
324: }
325: *px = rprod;
326: *py = iprod;
327: return(0);
328: }
330: #include <petscdraw.h>
331: /* collective on KSP */
332: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp,PetscInt neig,const PetscReal *r,const PetscReal *c)
333: {
335: PetscReal xmin,xmax,ymin,ymax,*xloc,*yloc,*value,px0,py0,rscale,iscale;
336: PetscInt M,N,i,j;
337: PetscMPIInt rank;
338: PetscViewer viewer;
339: PetscDraw draw;
340: PetscDrawAxis drawaxis;
343: MPI_Comm_rank(PetscObjectComm((PetscObject)ksp),&rank);
344: if (rank) return(0);
345: M = 80;
346: N = 80;
347: xmin = r[0]; xmax = r[0];
348: ymin = c[0]; ymax = c[0];
349: for (i=1; i<neig; i++) {
350: xmin = PetscMin(xmin,r[i]);
351: xmax = PetscMax(xmax,r[i]);
352: ymin = PetscMin(ymin,c[i]);
353: ymax = PetscMax(ymax,c[i]);
354: }
355: PetscMalloc3(M,&xloc,N,&yloc,M*N,&value);
356: for (i=0; i<M; i++) xloc[i] = xmin - 0.1*(xmax-xmin) + 1.2*(xmax-xmin)*i/(M-1);
357: for (i=0; i<N; i++) yloc[i] = ymin - 0.1*(ymax-ymin) + 1.2*(ymax-ymin)*i/(N-1);
358: PolyEval(neig,r,c,0,0,&px0,&py0);
359: rscale = px0/(PetscSqr(px0)+PetscSqr(py0));
360: iscale = -py0/(PetscSqr(px0)+PetscSqr(py0));
361: for (j=0; j<N; j++) {
362: for (i=0; i<M; i++) {
363: PetscReal px,py,tx,ty,tmod;
364: PolyEval(neig,r,c,xloc[i],yloc[j],&px,&py);
365: tx = px*rscale - py*iscale;
366: ty = py*rscale + px*iscale;
367: tmod = PetscSqr(tx) + PetscSqr(ty); /* modulus of the complex polynomial */
368: if (tmod > 1) tmod = 1.0;
369: if (tmod > 0.5 && tmod < 1) tmod = 0.5;
370: if (tmod > 0.2 && tmod < 0.5) tmod = 0.2;
371: if (tmod > 0.05 && tmod < 0.2) tmod = 0.05;
372: if (tmod < 1e-3) tmod = 1e-3;
373: value[i+j*M] = PetscLogReal(tmod) / PetscLogReal(10.0);
374: }
375: }
376: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Iteratively Computed Eigen-contours",PETSC_DECIDE,PETSC_DECIDE,450,450,&viewer);
377: PetscViewerDrawGetDraw(viewer,0,&draw);
378: PetscDrawTensorContour(draw,M,N,NULL,NULL,value);
379: if (0) {
380: PetscDrawAxisCreate(draw,&drawaxis);
381: PetscDrawAxisSetLimits(drawaxis,xmin,xmax,ymin,ymax);
382: PetscDrawAxisSetLabels(drawaxis,"Eigen-counters","real","imag");
383: PetscDrawAxisDraw(drawaxis);
384: PetscDrawAxisDestroy(&drawaxis);
385: }
386: PetscViewerDestroy(&viewer);
387: PetscFree3(xloc,yloc,value);
388: return(0);
389: }