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: }