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