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