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