Actual source code: eige.c
petsc-3.4.5 2014-06-29
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: PetscMalloc(m*sizeof(PetscInt),&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,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
106: - c - complex part of computed eigenvalues
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(BA,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: PetscMalloc(clen*sizeof(PetscScalar),&cwork);
188: idummy = -1; /* unused */
189: PetscBLASIntCast(n,&bn);
190: lwork = 5*n;
191: PetscMalloc(lwork*sizeof(PetscReal),&work);
192: PetscMalloc(n*sizeof(PetscReal),&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: PetscMalloc(n*sizeof(PetscInt),&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: PetscMalloc(2*n*sizeof(PetscReal),&realpart);
237: imagpart = realpart + n;
238: PetscMalloc(5*n*sizeof(PetscReal),&work);
239: #if defined(PETSC_MISSING_LAPACK_GEEV)
240: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
241: #else
242: {
243: PetscBLASInt lierr;
244: PetscScalar sdummy;
245: PetscBLASInt bn;
247: PetscBLASIntCast(n,&bn);
248: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
249: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,array,&bn,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr));
250: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
251: PetscFPTrapPop();
252: }
253: #endif
254: PetscFree(work);
255: PetscMalloc(n*sizeof(PetscInt),&perm);
257: for (i=0; i<n; i++) perm[i] = i;
258: PetscSortRealWithPermutation(n,realpart,perm);
259: for (i=0; i<n; i++) {
260: r[i] = realpart[perm[i]];
261: c[i] = imagpart[perm[i]];
262: }
263: PetscFree(perm);
264: PetscFree(realpart);
265: }
266: #else
267: if (!rank) {
268: PetscScalar *work,*eigs;
269: PetscReal *rwork;
270: PetscBLASInt idummy,lwork;
271: PetscInt *perm;
273: idummy = n;
274: lwork = 5*n;
275: PetscMalloc(5*n*sizeof(PetscScalar),&work);
276: PetscMalloc(2*n*sizeof(PetscReal),&rwork);
277: PetscMalloc(n*sizeof(PetscScalar),&eigs);
278: #if defined(PETSC_MISSING_LAPACK_GEEV)
279: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
280: #else
281: {
282: PetscBLASInt lierr;
283: PetscScalar sdummy;
284: PetscBLASInt nb;
285: PetscBLASIntCast(n,&nb);
286: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
287: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&nb,array,&nb,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&lierr));
288: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
289: PetscFPTrapPop();
290: }
291: #endif
292: PetscFree(work);
293: PetscFree(rwork);
294: PetscMalloc(n*sizeof(PetscInt),&perm);
295: for (i=0; i<n; i++) perm[i] = i;
296: for (i=0; i<n; i++) r[i] = PetscRealPart(eigs[i]);
297: PetscSortRealWithPermutation(n,r,perm);
298: for (i=0; i<n; i++) {
299: r[i] = PetscRealPart(eigs[perm[i]]);
300: c[i] = PetscImaginaryPart(eigs[perm[i]]);
301: }
302: PetscFree(perm);
303: PetscFree(eigs);
304: }
305: #endif
306: if (size > 1) {
307: MatDenseRestoreArray(A,&array);
308: MatDestroy(&A);
309: } else {
310: MatDenseRestoreArray(BA,&array);
311: }
312: MatDestroy(&BA);
313: return(0);
314: }
318: static PetscErrorCode PolyEval(PetscInt nroots,const PetscReal *r,const PetscReal *c,PetscReal x,PetscReal y,PetscReal *px,PetscReal *py)
319: {
320: PetscInt i;
321: PetscReal rprod = 1,iprod = 0;
324: for (i=0; i<nroots; i++) {
325: PetscReal rnew = rprod*(x - r[i]) - iprod*(y - c[i]);
326: PetscReal inew = rprod*(y - c[i]) + iprod*(x - r[i]);
327: rprod = rnew;
328: iprod = inew;
329: }
330: *px = rprod;
331: *py = iprod;
332: return(0);
333: }
335: #include <petscdraw.h>
338: /* collective on KSP */
339: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp,PetscInt neig,const PetscReal *r,const PetscReal *c)
340: {
342: PetscReal xmin,xmax,ymin,ymax,*xloc,*yloc,*value,px0,py0,rscale,iscale;
343: PetscInt M,N,i,j;
344: PetscMPIInt rank;
345: PetscViewer viewer;
346: PetscDraw draw;
347: PetscDrawAxis drawaxis;
350: MPI_Comm_rank(PetscObjectComm((PetscObject)ksp),&rank);
351: if (rank) return(0);
352: M = 80;
353: N = 80;
354: xmin = r[0]; xmax = r[0];
355: ymin = c[0]; ymax = c[0];
356: for (i=1; i<neig; i++) {
357: xmin = PetscMin(xmin,r[i]);
358: xmax = PetscMax(xmax,r[i]);
359: ymin = PetscMin(ymin,c[i]);
360: ymax = PetscMax(ymax,c[i]);
361: }
362: PetscMalloc3(M,PetscReal,&xloc,N,PetscReal,&yloc,M*N,PetscReal,&value);
363: for (i=0; i<M; i++) xloc[i] = xmin - 0.1*(xmax-xmin) + 1.2*(xmax-xmin)*i/(M-1);
364: for (i=0; i<N; i++) yloc[i] = ymin - 0.1*(ymax-ymin) + 1.2*(ymax-ymin)*i/(N-1);
365: PolyEval(neig,r,c,0,0,&px0,&py0);
366: rscale = px0/(PetscSqr(px0)+PetscSqr(py0));
367: iscale = -py0/(PetscSqr(px0)+PetscSqr(py0));
368: for (j=0; j<N; j++) {
369: for (i=0; i<M; i++) {
370: PetscReal px,py,tx,ty,tmod;
371: PolyEval(neig,r,c,xloc[i],yloc[j],&px,&py);
372: tx = px*rscale - py*iscale;
373: ty = py*rscale + px*iscale;
374: tmod = PetscSqr(tx) + PetscSqr(ty); /* modulus of the complex polynomial */
375: if (tmod > 1) tmod = 1.0;
376: if (tmod > 0.5 && tmod < 1) tmod = 0.5;
377: if (tmod > 0.2 && tmod < 0.5) tmod = 0.2;
378: if (tmod > 0.05 && tmod < 0.2) tmod = 0.05;
379: if (tmod < 1e-3) tmod = 1e-3;
380: value[i+j*M] = PetscLogReal(tmod) / PetscLogReal(10.0);
381: }
382: }
383: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Iteratively Computed Eigen-contours",PETSC_DECIDE,PETSC_DECIDE,450,450,&viewer);
384: PetscViewerDrawGetDraw(viewer,0,&draw);
385: PetscDrawTensorContour(draw,M,N,NULL,NULL,value);
386: if (0) {
387: PetscDrawAxisCreate(draw,&drawaxis);
388: PetscDrawAxisSetLimits(drawaxis,xmin,xmax,ymin,ymax);
389: PetscDrawAxisSetLabels(drawaxis,"Eigen-counters","real","imag");
390: PetscDrawAxisDraw(drawaxis);
391: PetscDrawAxisDestroy(&drawaxis);
392: }
393: PetscViewerDestroy(&viewer);
394: PetscFree3(xloc,yloc,value);
395: return(0);
396: }