Actual source code: eige.c
2: #include <petsc/private/kspimpl.h>
3: #include <petscdm.h>
4: #include <petscblaslapack.h>
6: typedef struct {
7: KSP ksp;
8: Vec work;
9: } Mat_KSP;
11: static PetscErrorCode MatCreateVecs_KSP(Mat A,Vec *X,Vec *Y)
12: {
13: Mat_KSP *ctx;
14: Mat M;
16: MatShellGetContext(A,&ctx);
17: KSPGetOperators(ctx->ksp,&M,NULL);
18: MatCreateVecs(M,X,Y);
19: return 0;
20: }
22: static PetscErrorCode MatMult_KSP(Mat A,Vec X,Vec Y)
23: {
24: Mat_KSP *ctx;
26: MatShellGetContext(A,&ctx);
27: KSP_PCApplyBAorAB(ctx->ksp,X,Y,ctx->work);
28: return 0;
29: }
31: /*@
32: KSPComputeOperator - Computes the explicit preconditioned operator, including diagonal scaling and null
33: space removal if applicable.
35: Collective on ksp
37: Input Parameters:
38: + ksp - the Krylov subspace context
39: - mattype - the matrix type to be used
41: Output Parameter:
42: . mat - the explicit preconditioned operator
44: Notes:
45: This computation is done by applying the operators to columns of the
46: identity matrix.
48: Currently, this routine uses a dense matrix format for the output operator if mattype == NULL.
49: This routine is costly in general, and is recommended for use only with relatively small systems.
51: Level: advanced
53: .seealso: KSPComputeEigenvaluesExplicitly(), PCComputeOperator(), KSPSetDiagonalScale(), KSPSetNullSpace(), MatType
54: @*/
55: PetscErrorCode KSPComputeOperator(KSP ksp, MatType mattype, Mat *mat)
56: {
57: PetscInt N,M,m,n;
58: Mat_KSP ctx;
59: Mat A,Aksp;
63: KSPGetOperators(ksp,&A,NULL);
64: MatGetLocalSize(A,&m,&n);
65: MatGetSize(A,&M,&N);
66: MatCreateShell(PetscObjectComm((PetscObject)ksp),m,n,M,N,&ctx,&Aksp);
67: MatShellSetOperation(Aksp,MATOP_MULT,(void (*)(void))MatMult_KSP);
68: MatShellSetOperation(Aksp,MATOP_CREATE_VECS,(void (*)(void))MatCreateVecs_KSP);
69: ctx.ksp = ksp;
70: MatCreateVecs(A,&ctx.work,NULL);
71: MatComputeOperator(Aksp,mattype,mat);
72: VecDestroy(&ctx.work);
73: MatDestroy(&Aksp);
74: return 0;
75: }
77: /*@
78: KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
79: preconditioned operator using LAPACK.
81: Collective on ksp
83: Input Parameters:
84: + ksp - iterative context obtained from KSPCreate()
85: - n - size of arrays r and c
87: Output Parameters:
88: + r - real part of computed eigenvalues, provided by user with a dimension at least of n
89: - c - complex part of computed eigenvalues, provided by user with a dimension at least of n
91: Notes:
92: This approach is very slow but will generally provide accurate eigenvalue
93: estimates. This routine explicitly forms a dense matrix representing
94: the preconditioned operator, and thus will run only for relatively small
95: problems, say n < 500.
97: Many users may just want to use the monitoring routine
98: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
99: to print the singular values at each iteration of the linear solve.
101: The preconditoner operator, rhs vector, solution vectors should be
102: set before this routine is called. i.e use KSPSetOperators(),KSPSolve() or
103: KSPSetOperators()
105: Level: advanced
107: .seealso: KSPComputeEigenvalues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSPSetOperators(), KSPSolve()
108: @*/
109: PetscErrorCode KSPComputeEigenvaluesExplicitly(KSP ksp,PetscInt nmax,PetscReal r[],PetscReal c[])
110: {
111: Mat BA;
112: PetscMPIInt size,rank;
113: MPI_Comm comm;
114: PetscScalar *array;
115: Mat A;
116: PetscInt m,row,nz,i,n,dummy;
117: const PetscInt *cols;
118: const PetscScalar *vals;
120: PetscObjectGetComm((PetscObject)ksp,&comm);
121: KSPComputeOperator(ksp,MATDENSE,&BA);
122: MPI_Comm_size(comm,&size);
123: MPI_Comm_rank(comm,&rank);
125: MatGetSize(BA,&n,&n);
126: if (size > 1) { /* assemble matrix on first processor */
127: MatCreate(PetscObjectComm((PetscObject)ksp),&A);
128: if (rank == 0) {
129: MatSetSizes(A,n,n,n,n);
130: } else {
131: MatSetSizes(A,0,0,n,n);
132: }
133: MatSetType(A,MATMPIDENSE);
134: MatMPIDenseSetPreallocation(A,NULL);
135: PetscLogObjectParent((PetscObject)BA,(PetscObject)A);
137: MatGetOwnershipRange(BA,&row,&dummy);
138: MatGetLocalSize(BA,&m,&dummy);
139: for (i=0; i<m; i++) {
140: MatGetRow(BA,row,&nz,&cols,&vals);
141: MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
142: MatRestoreRow(BA,row,&nz,&cols,&vals);
143: row++;
144: }
146: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
147: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
148: MatDenseGetArray(A,&array);
149: } else {
150: MatDenseGetArray(BA,&array);
151: }
153: #if !defined(PETSC_USE_COMPLEX)
154: if (rank == 0) {
155: PetscScalar *work;
156: PetscReal *realpart,*imagpart;
157: PetscBLASInt idummy,lwork;
158: PetscInt *perm;
160: idummy = n;
161: lwork = 5*n;
162: PetscMalloc2(n,&realpart,n,&imagpart);
163: PetscMalloc1(5*n,&work);
164: {
165: PetscBLASInt lierr;
166: PetscScalar sdummy;
167: PetscBLASInt bn;
169: PetscBLASIntCast(n,&bn);
170: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
171: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&bn,array,&bn,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr));
173: PetscFPTrapPop();
174: }
175: PetscFree(work);
176: PetscMalloc1(n,&perm);
178: for (i=0; i<n; i++) perm[i] = i;
179: PetscSortRealWithPermutation(n,realpart,perm);
180: for (i=0; i<n; i++) {
181: r[i] = realpart[perm[i]];
182: c[i] = imagpart[perm[i]];
183: }
184: PetscFree(perm);
185: PetscFree2(realpart,imagpart);
186: }
187: #else
188: if (rank == 0) {
189: PetscScalar *work,*eigs;
190: PetscReal *rwork;
191: PetscBLASInt idummy,lwork;
192: PetscInt *perm;
194: idummy = n;
195: lwork = 5*n;
196: PetscMalloc1(5*n,&work);
197: PetscMalloc1(2*n,&rwork);
198: PetscMalloc1(n,&eigs);
199: {
200: PetscBLASInt lierr;
201: PetscScalar sdummy;
202: PetscBLASInt nb;
203: PetscBLASIntCast(n,&nb);
204: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
205: PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&nb,array,&nb,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&lierr));
207: PetscFPTrapPop();
208: }
209: PetscFree(work);
210: PetscFree(rwork);
211: PetscMalloc1(n,&perm);
212: for (i=0; i<n; i++) perm[i] = i;
213: for (i=0; i<n; i++) r[i] = PetscRealPart(eigs[i]);
214: PetscSortRealWithPermutation(n,r,perm);
215: for (i=0; i<n; i++) {
216: r[i] = PetscRealPart(eigs[perm[i]]);
217: c[i] = PetscImaginaryPart(eigs[perm[i]]);
218: }
219: PetscFree(perm);
220: PetscFree(eigs);
221: }
222: #endif
223: if (size > 1) {
224: MatDenseRestoreArray(A,&array);
225: MatDestroy(&A);
226: } else {
227: MatDenseRestoreArray(BA,&array);
228: }
229: MatDestroy(&BA);
230: return 0;
231: }
233: static PetscErrorCode PolyEval(PetscInt nroots,const PetscReal *r,const PetscReal *c,PetscReal x,PetscReal y,PetscReal *px,PetscReal *py)
234: {
235: PetscInt i;
236: PetscReal rprod = 1,iprod = 0;
238: for (i=0; i<nroots; i++) {
239: PetscReal rnew = rprod*(x - r[i]) - iprod*(y - c[i]);
240: PetscReal inew = rprod*(y - c[i]) + iprod*(x - r[i]);
241: rprod = rnew;
242: iprod = inew;
243: }
244: *px = rprod;
245: *py = iprod;
246: return 0;
247: }
249: #include <petscdraw.h>
250: /* collective on ksp */
251: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp,PetscInt neig,const PetscReal *r,const PetscReal *c)
252: {
253: PetscReal xmin,xmax,ymin,ymax,*xloc,*yloc,*value,px0,py0,rscale,iscale;
254: PetscInt M,N,i,j;
255: PetscMPIInt rank;
256: PetscViewer viewer;
257: PetscDraw draw;
258: PetscDrawAxis drawaxis;
260: MPI_Comm_rank(PetscObjectComm((PetscObject)ksp),&rank);
261: if (rank) return 0;
262: M = 80;
263: N = 80;
264: xmin = r[0]; xmax = r[0];
265: ymin = c[0]; ymax = c[0];
266: for (i=1; i<neig; i++) {
267: xmin = PetscMin(xmin,r[i]);
268: xmax = PetscMax(xmax,r[i]);
269: ymin = PetscMin(ymin,c[i]);
270: ymax = PetscMax(ymax,c[i]);
271: }
272: PetscMalloc3(M,&xloc,N,&yloc,M*N,&value);
273: for (i=0; i<M; i++) xloc[i] = xmin - 0.1*(xmax-xmin) + 1.2*(xmax-xmin)*i/(M-1);
274: for (i=0; i<N; i++) yloc[i] = ymin - 0.1*(ymax-ymin) + 1.2*(ymax-ymin)*i/(N-1);
275: PolyEval(neig,r,c,0,0,&px0,&py0);
276: rscale = px0/(PetscSqr(px0)+PetscSqr(py0));
277: iscale = -py0/(PetscSqr(px0)+PetscSqr(py0));
278: for (j=0; j<N; j++) {
279: for (i=0; i<M; i++) {
280: PetscReal px,py,tx,ty,tmod;
281: PolyEval(neig,r,c,xloc[i],yloc[j],&px,&py);
282: tx = px*rscale - py*iscale;
283: ty = py*rscale + px*iscale;
284: tmod = PetscSqr(tx) + PetscSqr(ty); /* modulus of the complex polynomial */
285: if (tmod > 1) tmod = 1.0;
286: if (tmod > 0.5 && tmod < 1) tmod = 0.5;
287: if (tmod > 0.2 && tmod < 0.5) tmod = 0.2;
288: if (tmod > 0.05 && tmod < 0.2) tmod = 0.05;
289: if (tmod < 1e-3) tmod = 1e-3;
290: value[i+j*M] = PetscLogReal(tmod) / PetscLogReal(10.0);
291: }
292: }
293: PetscViewerDrawOpen(PETSC_COMM_SELF,NULL,"Iteratively Computed Eigen-contours",PETSC_DECIDE,PETSC_DECIDE,450,450,&viewer);
294: PetscViewerDrawGetDraw(viewer,0,&draw);
295: PetscDrawTensorContour(draw,M,N,NULL,NULL,value);
296: if (0) {
297: PetscDrawAxisCreate(draw,&drawaxis);
298: PetscDrawAxisSetLimits(drawaxis,xmin,xmax,ymin,ymax);
299: PetscDrawAxisSetLabels(drawaxis,"Eigen-counters","real","imag");
300: PetscDrawAxisDraw(drawaxis);
301: PetscDrawAxisDestroy(&drawaxis);
302: }
303: PetscViewerDestroy(&viewer);
304: PetscFree3(xloc,yloc,value);
305: return 0;
306: }