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