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