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
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: [](chapter_ksp), `KSP`, `KSPSetOperators()`, `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
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 preconditioner operator, rhs vector, solution vectors should be
102: set before this routine is called. i.e use `KSPSetOperators()`, `KSPSolve()`
104: Level: advanced
106: .seealso: [](chapter_ksp), `KSP`, `KSPComputeEigenvalues()`, `KSPMonitorSingularValue()`, `KSPComputeExtremeSingularValues()`, `KSPSetOperators()`, `KSPSolve()`
107: @*/
108: PetscErrorCode KSPComputeEigenvaluesExplicitly(KSP ksp, PetscInt nmax, PetscReal r[], PetscReal c[])
109: {
110: Mat BA;
111: PetscMPIInt size, rank;
112: MPI_Comm comm;
113: PetscScalar *array;
114: Mat A;
115: PetscInt m, row, nz, i, n, dummy;
116: const PetscInt *cols;
117: const PetscScalar *vals;
119: PetscObjectGetComm((PetscObject)ksp, &comm);
120: KSPComputeOperator(ksp, MATDENSE, &BA);
121: MPI_Comm_size(comm, &size);
122: MPI_Comm_rank(comm, &rank);
124: MatGetSize(BA, &n, &n);
125: if (size > 1) { /* assemble matrix on first processor */
126: MatCreate(PetscObjectComm((PetscObject)ksp), &A);
127: if (rank == 0) {
128: MatSetSizes(A, n, n, n, n);
129: } else {
130: MatSetSizes(A, 0, 0, n, n);
131: }
132: MatSetType(A, MATMPIDENSE);
133: MatMPIDenseSetPreallocation(A, NULL);
135: MatGetOwnershipRange(BA, &row, &dummy);
136: MatGetLocalSize(BA, &m, &dummy);
137: for (i = 0; i < m; i++) {
138: MatGetRow(BA, row, &nz, &cols, &vals);
139: MatSetValues(A, 1, &row, nz, cols, vals, INSERT_VALUES);
140: MatRestoreRow(BA, row, &nz, &cols, &vals);
141: row++;
142: }
144: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
145: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
146: MatDenseGetArray(A, &array);
147: } else {
148: MatDenseGetArray(BA, &array);
149: }
151: #if !defined(PETSC_USE_COMPLEX)
152: if (rank == 0) {
153: PetscScalar *work;
154: PetscReal *realpart, *imagpart;
155: PetscBLASInt idummy, lwork;
156: PetscInt *perm;
158: idummy = n;
159: lwork = 5 * n;
160: PetscMalloc2(n, &realpart, n, &imagpart);
161: PetscMalloc1(5 * n, &work);
162: {
163: PetscBLASInt lierr;
164: PetscScalar sdummy;
165: PetscBLASInt bn;
167: PetscBLASIntCast(n, &bn);
168: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
169: PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "N", &bn, array, &bn, realpart, imagpart, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, &lierr));
171: PetscFPTrapPop();
172: }
173: PetscFree(work);
174: PetscMalloc1(n, &perm);
176: for (i = 0; i < n; i++) perm[i] = i;
177: PetscSortRealWithPermutation(n, realpart, perm);
178: for (i = 0; i < n; i++) {
179: r[i] = realpart[perm[i]];
180: c[i] = imagpart[perm[i]];
181: }
182: PetscFree(perm);
183: PetscFree2(realpart, imagpart);
184: }
185: #else
186: if (rank == 0) {
187: PetscScalar *work, *eigs;
188: PetscReal *rwork;
189: PetscBLASInt idummy, lwork;
190: PetscInt *perm;
192: idummy = n;
193: lwork = 5 * n;
194: PetscMalloc1(5 * n, &work);
195: PetscMalloc1(2 * n, &rwork);
196: PetscMalloc1(n, &eigs);
197: {
198: PetscBLASInt lierr;
199: PetscScalar sdummy;
200: PetscBLASInt nb;
201: PetscBLASIntCast(n, &nb);
202: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
203: PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "N", &nb, array, &nb, eigs, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, rwork, &lierr));
205: PetscFPTrapPop();
206: }
207: PetscFree(work);
208: PetscFree(rwork);
209: PetscMalloc1(n, &perm);
210: for (i = 0; i < n; i++) perm[i] = i;
211: for (i = 0; i < n; i++) r[i] = PetscRealPart(eigs[i]);
212: PetscSortRealWithPermutation(n, r, perm);
213: for (i = 0; i < n; i++) {
214: r[i] = PetscRealPart(eigs[perm[i]]);
215: c[i] = PetscImaginaryPart(eigs[perm[i]]);
216: }
217: PetscFree(perm);
218: PetscFree(eigs);
219: }
220: #endif
221: if (size > 1) {
222: MatDenseRestoreArray(A, &array);
223: MatDestroy(&A);
224: } else {
225: MatDenseRestoreArray(BA, &array);
226: }
227: MatDestroy(&BA);
228: return 0;
229: }
231: static PetscErrorCode PolyEval(PetscInt nroots, const PetscReal *r, const PetscReal *c, PetscReal x, PetscReal y, PetscReal *px, PetscReal *py)
232: {
233: PetscInt i;
234: PetscReal rprod = 1, iprod = 0;
236: for (i = 0; i < nroots; i++) {
237: PetscReal rnew = rprod * (x - r[i]) - iprod * (y - c[i]);
238: PetscReal inew = rprod * (y - c[i]) + iprod * (x - r[i]);
239: rprod = rnew;
240: iprod = inew;
241: }
242: *px = rprod;
243: *py = iprod;
244: return 0;
245: }
247: #include <petscdraw.h>
248: /* Collective */
249: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp, PetscInt neig, const PetscReal *r, const PetscReal *c)
250: {
251: PetscReal xmin, xmax, ymin, ymax, *xloc, *yloc, *value, px0, py0, rscale, iscale;
252: PetscInt M, N, i, j;
253: PetscMPIInt rank;
254: PetscViewer viewer;
255: PetscDraw draw;
256: PetscDrawAxis drawaxis;
258: MPI_Comm_rank(PetscObjectComm((PetscObject)ksp), &rank);
259: if (rank) return 0;
260: M = 80;
261: N = 80;
262: xmin = r[0];
263: xmax = r[0];
264: ymin = c[0];
265: 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: }