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