Actual source code: eige.c
1: #include <petsc/private/kspimpl.h>
2: #include <petscdm.h>
3: #include <petscblaslapack.h>
5: typedef struct {
6: KSP ksp;
7: Vec work;
8: } Mat_KSP;
10: static PetscErrorCode MatCreateVecs_KSP(Mat A, Vec *X, Vec *Y)
11: {
12: Mat_KSP *ctx;
13: Mat M;
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: }
22: static PetscErrorCode MatMult_KSP(Mat A, Vec X, Vec Y)
23: {
24: Mat_KSP *ctx;
26: PetscFunctionBegin;
27: PetscCall(MatShellGetContext(A, &ctx));
28: PetscCall(KSP_PCApplyBAorAB(ctx->ksp, X, Y, ctx->work));
29: PetscFunctionReturn(PETSC_SUCCESS);
30: }
32: /*@
33: KSPComputeOperator - Computes the explicit preconditioned operator, including null 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: Level: advanced
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` is `NULL`.
51: This routine is costly in general, and is recommended for use only with relatively small systems.
53: .seealso: [](ch_ksp), `KSP`, `KSPSetOperators()`, `KSPComputeEigenvaluesExplicitly()`, `PCComputeOperator()`, `MatSetNullSpace()`, `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;
61: PetscFunctionBegin;
63: PetscAssertPointer(mat, 3);
64: PetscCall(KSPGetOperators(ksp, &A, NULL));
65: PetscCall(MatGetLocalSize(A, &m, &n));
66: PetscCall(MatGetSize(A, &M, &N));
67: PetscCall(MatCreateShell(PetscObjectComm((PetscObject)ksp), m, n, M, N, &ctx, &Aksp));
68: PetscCall(MatShellSetOperation(Aksp, MATOP_MULT, (PetscErrorCodeFn *)MatMult_KSP));
69: PetscCall(MatShellSetOperation(Aksp, MATOP_CREATE_VECS, (PetscErrorCodeFn *)MatCreateVecs_KSP));
70: ctx.ksp = ksp;
71: PetscCall(MatCreateVecs(A, &ctx.work, NULL));
72: PetscCall(MatComputeOperator(Aksp, mattype, mat));
73: PetscCall(VecDestroy(&ctx.work));
74: PetscCall(MatDestroy(&Aksp));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: /*@
79: KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
80: preconditioned operator using LAPACK.
82: Collective
84: Input Parameters:
85: + ksp - iterative context obtained from `KSPCreate()`
86: - nmax - size of arrays `r` and `c`
88: Output Parameters:
89: + r - real part of computed eigenvalues, provided by user with a dimension at least of `n`
90: - c - complex part of computed eigenvalues, provided by user with a dimension at least of `n`
92: Level: advanced
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, and solution vectors should be
105: set before this routine is called. i.e use `KSPSetOperators()`, `KSPSolve()`
107: .seealso: [](ch_ksp), `KSP`, `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: PetscFunctionBegin;
121: PetscCall(PetscObjectGetComm((PetscObject)ksp, &comm));
122: PetscCall(KSPComputeOperator(ksp, MATDENSE, &BA));
123: PetscCallMPI(MPI_Comm_size(comm, &size));
124: PetscCallMPI(MPI_Comm_rank(comm, &rank));
126: PetscCall(MatGetSize(BA, &n, &n));
127: if (size > 1) { /* assemble matrix on first processor */
128: PetscCall(MatCreate(PetscObjectComm((PetscObject)ksp), &A));
129: if (rank == 0) {
130: PetscCall(MatSetSizes(A, n, n, n, n));
131: } else {
132: PetscCall(MatSetSizes(A, 0, 0, n, n));
133: }
134: PetscCall(MatSetType(A, MATMPIDENSE));
135: PetscCall(MatMPIDenseSetPreallocation(A, NULL));
137: PetscCall(MatGetOwnershipRange(BA, &row, &dummy));
138: PetscCall(MatGetLocalSize(BA, &m, &dummy));
139: for (i = 0; i < m; i++) {
140: PetscCall(MatGetRow(BA, row, &nz, &cols, &vals));
141: PetscCall(MatSetValues(A, 1, &row, nz, cols, vals, INSERT_VALUES));
142: PetscCall(MatRestoreRow(BA, row, &nz, &cols, &vals));
143: row++;
144: }
146: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
147: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
148: PetscCall(MatDenseGetArray(A, &array));
149: } else {
150: PetscCall(MatDenseGetArray(BA, &array));
151: }
153: #if !PetscDefined(USE_COMPLEX)
154: if (rank == 0) {
155: PetscScalar *work;
156: PetscReal *realpart, *imagpart;
157: PetscBLASInt idummy, lwork;
158: PetscInt *perm;
160: PetscCall(PetscBLASIntCast(n, &idummy));
161: PetscCall(PetscBLASIntCast(5 * n, &lwork));
162: PetscCall(PetscMalloc2(n, &realpart, n, &imagpart));
163: PetscCall(PetscMalloc1(5 * n, &work));
164: {
165: PetscScalar sdummy;
166: PetscBLASInt bn;
168: PetscCall(PetscBLASIntCast(n, &bn));
169: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
170: PetscCallLAPACKInfo("LAPACKgeev", LAPACKgeev_("N", "N", &bn, array, &bn, realpart, imagpart, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, &info));
171: PetscCall(PetscFPTrapPop());
172: }
173: PetscCall(PetscFree(work));
174: PetscCall(PetscMalloc1(n, &perm));
176: for (i = 0; i < n; i++) perm[i] = i;
177: PetscCall(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: PetscCall(PetscFree(perm));
183: PetscCall(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: PetscCall(PetscBLASIntCast(n, &idummy));
193: PetscCall(PetscBLASIntCast(5 * n, &lwork));
194: PetscCall(PetscMalloc1(5 * n, &work));
195: PetscCall(PetscMalloc1(2 * n, &rwork));
196: PetscCall(PetscMalloc1(n, &eigs));
197: {
198: PetscScalar sdummy;
199: PetscBLASInt nb;
200: PetscCall(PetscBLASIntCast(n, &nb));
201: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
202: PetscCallLAPACKInfo("LAPACKgeev", LAPACKgeev_("N", "N", &nb, array, &nb, eigs, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, rwork, &info));
203: PetscCall(PetscFPTrapPop());
204: }
205: PetscCall(PetscFree(work));
206: PetscCall(PetscFree(rwork));
207: PetscCall(PetscMalloc1(n, &perm));
208: for (i = 0; i < n; i++) perm[i] = i;
209: for (i = 0; i < n; i++) r[i] = PetscRealPart(eigs[i]);
210: PetscCall(PetscSortRealWithPermutation(n, r, perm));
211: for (i = 0; i < n; i++) {
212: r[i] = PetscRealPart(eigs[perm[i]]);
213: c[i] = PetscImaginaryPart(eigs[perm[i]]);
214: }
215: PetscCall(PetscFree(perm));
216: PetscCall(PetscFree(eigs));
217: }
218: #endif
219: if (size > 1) {
220: PetscCall(MatDenseRestoreArray(A, &array));
221: PetscCall(MatDestroy(&A));
222: } else {
223: PetscCall(MatDenseRestoreArray(BA, &array));
224: }
225: PetscCall(MatDestroy(&BA));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: static PetscErrorCode PolyEval(PetscInt nroots, const PetscReal *r, const PetscReal *c, PetscReal x, PetscReal y, PetscReal *px, PetscReal *py)
230: {
231: PetscInt i;
232: PetscReal rprod = 1, iprod = 0;
234: PetscFunctionBegin;
235: for (i = 0; i < nroots; i++) {
236: PetscReal rnew = rprod * (x - r[i]) - iprod * (y - c[i]);
237: PetscReal inew = rprod * (y - c[i]) + iprod * (x - r[i]);
238: rprod = rnew;
239: iprod = inew;
240: }
241: *px = rprod;
242: *py = iprod;
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: #include <petscdraw.h>
247: /* Collective */
248: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp, PetscInt neig, const PetscReal *r, const PetscReal *c)
249: {
250: PetscReal xmin, xmax, ymin, ymax, *xloc, *yloc, *value, px0, py0, rscale, iscale;
251: int M, N, i, j;
252: PetscMPIInt rank;
253: PetscViewer viewer;
254: PetscDraw draw;
255: PetscDrawAxis drawaxis;
257: PetscFunctionBegin;
258: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ksp), &rank));
259: if (rank) PetscFunctionReturn(PETSC_SUCCESS);
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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, "Iteratively Computed Eigen-contours", PETSC_DECIDE, PETSC_DECIDE, 450, 450, &viewer));
294: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
295: PetscCall(PetscDrawTensorContour(draw, M, N, NULL, NULL, value));
296: if (0) {
297: PetscCall(PetscDrawAxisCreate(draw, &drawaxis));
298: PetscCall(PetscDrawAxisSetLimits(drawaxis, xmin, xmax, ymin, ymax));
299: PetscCall(PetscDrawAxisSetLabels(drawaxis, "Eigen-counters", "real", "imag"));
300: PetscCall(PetscDrawAxisDraw(drawaxis));
301: PetscCall(PetscDrawAxisDestroy(&drawaxis));
302: }
303: PetscCall(PetscViewerDestroy(&viewer));
304: PetscCall(PetscFree3(xloc, yloc, value));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }