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 diagonal scaling and null
     34:   space removal if applicable.

     36:   Collective

     38:   Input Parameters:
     39: + ksp     - the Krylov subspace context
     40: - mattype - the matrix type to be used

     42:   Output Parameter:
     43: . mat - the explicit preconditioned operator

     45:   Level: advanced

     47:   Notes:
     48:   This computation is done by applying the operators to columns of the
     49:   identity matrix.

     51:   Currently, this routine uses a dense matrix format for the output operator if `mattype` is `NULL`.
     52:   This routine is costly in general, and is recommended for use only with relatively small systems.

     54: .seealso: [](ch_ksp), `KSP`, `KSPSetOperators()`, `KSPComputeEigenvaluesExplicitly()`, `PCComputeOperator()`, `KSPSetDiagonalScale()`, `KSPSetNullSpace()`, `MatType`
     55: @*/
     56: PetscErrorCode KSPComputeOperator(KSP ksp, MatType mattype, Mat *mat)
     57: {
     58:   PetscInt N, M, m, n;
     59:   Mat_KSP  ctx;
     60:   Mat      A, Aksp;

     62:   PetscFunctionBegin;
     64:   PetscAssertPointer(mat, 3);
     65:   PetscCall(KSPGetOperators(ksp, &A, NULL));
     66:   PetscCall(MatGetLocalSize(A, &m, &n));
     67:   PetscCall(MatGetSize(A, &M, &N));
     68:   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)ksp), m, n, M, N, &ctx, &Aksp));
     69:   PetscCall(MatShellSetOperation(Aksp, MATOP_MULT, (PetscErrorCodeFn *)MatMult_KSP));
     70:   PetscCall(MatShellSetOperation(Aksp, MATOP_CREATE_VECS, (PetscErrorCodeFn *)MatCreateVecs_KSP));
     71:   ctx.ksp = ksp;
     72:   PetscCall(MatCreateVecs(A, &ctx.work, NULL));
     73:   PetscCall(MatComputeOperator(Aksp, mattype, mat));
     74:   PetscCall(VecDestroy(&ctx.work));
     75:   PetscCall(MatDestroy(&Aksp));
     76:   PetscFunctionReturn(PETSC_SUCCESS);
     77: }

     79: /*@
     80:   KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
     81:   preconditioned operator using LAPACK.

     83:   Collective

     85:   Input Parameters:
     86: + ksp  - iterative context obtained from `KSPCreate()`
     87: - nmax - size of arrays `r` and `c`

     89:   Output Parameters:
     90: + r - real part of computed eigenvalues, provided by user with a dimension at least of `n`
     91: - c - complex part of computed eigenvalues, provided by user with a dimension at least of `n`

     93:   Level: advanced

     95:   Notes:
     96:   This approach is very slow but will generally provide accurate eigenvalue
     97:   estimates.  This routine explicitly forms a dense matrix representing
     98:   the preconditioned operator, and thus will run only for relatively small
     99:   problems, say `n` < 500.

    101:   Many users may just want to use the monitoring routine
    102:   `KSPMonitorSingularValue()` (which can be set with option -ksp_monitor_singular_value)
    103:   to print the singular values at each iteration of the linear solve.

    105:   The preconditioner operator, rhs vector, and solution vectors should be
    106:   set before this routine is called. i.e use `KSPSetOperators()`, `KSPSolve()`

    108: .seealso: [](ch_ksp), `KSP`, `KSPComputeEigenvalues()`, `KSPMonitorSingularValue()`, `KSPComputeExtremeSingularValues()`, `KSPSetOperators()`, `KSPSolve()`
    109: @*/
    110: PetscErrorCode KSPComputeEigenvaluesExplicitly(KSP ksp, PetscInt nmax, PetscReal r[], PetscReal c[])
    111: {
    112:   Mat                BA;
    113:   PetscMPIInt        size, rank;
    114:   MPI_Comm           comm;
    115:   PetscScalar       *array;
    116:   Mat                A;
    117:   PetscInt           m, row, nz, i, n, dummy;
    118:   const PetscInt    *cols;
    119:   const PetscScalar *vals;

    121:   PetscFunctionBegin;
    122:   PetscCall(PetscObjectGetComm((PetscObject)ksp, &comm));
    123:   PetscCall(KSPComputeOperator(ksp, MATDENSE, &BA));
    124:   PetscCallMPI(MPI_Comm_size(comm, &size));
    125:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    127:   PetscCall(MatGetSize(BA, &n, &n));
    128:   if (size > 1) { /* assemble matrix on first processor */
    129:     PetscCall(MatCreate(PetscObjectComm((PetscObject)ksp), &A));
    130:     if (rank == 0) {
    131:       PetscCall(MatSetSizes(A, n, n, n, n));
    132:     } else {
    133:       PetscCall(MatSetSizes(A, 0, 0, n, n));
    134:     }
    135:     PetscCall(MatSetType(A, MATMPIDENSE));
    136:     PetscCall(MatMPIDenseSetPreallocation(A, NULL));

    138:     PetscCall(MatGetOwnershipRange(BA, &row, &dummy));
    139:     PetscCall(MatGetLocalSize(BA, &m, &dummy));
    140:     for (i = 0; i < m; i++) {
    141:       PetscCall(MatGetRow(BA, row, &nz, &cols, &vals));
    142:       PetscCall(MatSetValues(A, 1, &row, nz, cols, vals, INSERT_VALUES));
    143:       PetscCall(MatRestoreRow(BA, row, &nz, &cols, &vals));
    144:       row++;
    145:     }

    147:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    148:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    149:     PetscCall(MatDenseGetArray(A, &array));
    150:   } else {
    151:     PetscCall(MatDenseGetArray(BA, &array));
    152:   }

    154: #if !defined(PETSC_USE_COMPLEX)
    155:   if (rank == 0) {
    156:     PetscScalar *work;
    157:     PetscReal   *realpart, *imagpart;
    158:     PetscBLASInt idummy, lwork;
    159:     PetscInt    *perm;

    161:     PetscCall(PetscBLASIntCast(n, &idummy));
    162:     PetscCall(PetscBLASIntCast(5 * n, &lwork));
    163:     PetscCall(PetscMalloc2(n, &realpart, n, &imagpart));
    164:     PetscCall(PetscMalloc1(5 * n, &work));
    165:     {
    166:       PetscBLASInt lierr;
    167:       PetscScalar  sdummy;
    168:       PetscBLASInt bn;

    170:       PetscCall(PetscBLASIntCast(n, &bn));
    171:       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    172:       PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "N", &bn, array, &bn, realpart, imagpart, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, &lierr));
    173:       PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %" PetscBLASInt_FMT, lierr);
    174:       PetscCall(PetscFPTrapPop());
    175:     }
    176:     PetscCall(PetscFree(work));
    177:     PetscCall(PetscMalloc1(n, &perm));

    179:     for (i = 0; i < n; i++) perm[i] = i;
    180:     PetscCall(PetscSortRealWithPermutation(n, realpart, perm));
    181:     for (i = 0; i < n; i++) {
    182:       r[i] = realpart[perm[i]];
    183:       c[i] = imagpart[perm[i]];
    184:     }
    185:     PetscCall(PetscFree(perm));
    186:     PetscCall(PetscFree2(realpart, imagpart));
    187:   }
    188: #else
    189:   if (rank == 0) {
    190:     PetscScalar *work, *eigs;
    191:     PetscReal   *rwork;
    192:     PetscBLASInt idummy, lwork;
    193:     PetscInt    *perm;

    195:     PetscCall(PetscBLASIntCast(n, &idummy));
    196:     PetscCall(PetscBLASIntCast(5 * n, &lwork));
    197:     PetscCall(PetscMalloc1(5 * n, &work));
    198:     PetscCall(PetscMalloc1(2 * n, &rwork));
    199:     PetscCall(PetscMalloc1(n, &eigs));
    200:     {
    201:       PetscBLASInt lierr;
    202:       PetscScalar  sdummy;
    203:       PetscBLASInt nb;
    204:       PetscCall(PetscBLASIntCast(n, &nb));
    205:       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    206:       PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "N", &nb, array, &nb, eigs, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, rwork, &lierr));
    207:       PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %" PetscBLASInt_FMT, lierr);
    208:       PetscCall(PetscFPTrapPop());
    209:     }
    210:     PetscCall(PetscFree(work));
    211:     PetscCall(PetscFree(rwork));
    212:     PetscCall(PetscMalloc1(n, &perm));
    213:     for (i = 0; i < n; i++) perm[i] = i;
    214:     for (i = 0; i < n; i++) r[i] = PetscRealPart(eigs[i]);
    215:     PetscCall(PetscSortRealWithPermutation(n, r, perm));
    216:     for (i = 0; i < n; i++) {
    217:       r[i] = PetscRealPart(eigs[perm[i]]);
    218:       c[i] = PetscImaginaryPart(eigs[perm[i]]);
    219:     }
    220:     PetscCall(PetscFree(perm));
    221:     PetscCall(PetscFree(eigs));
    222:   }
    223: #endif
    224:   if (size > 1) {
    225:     PetscCall(MatDenseRestoreArray(A, &array));
    226:     PetscCall(MatDestroy(&A));
    227:   } else {
    228:     PetscCall(MatDenseRestoreArray(BA, &array));
    229:   }
    230:   PetscCall(MatDestroy(&BA));
    231:   PetscFunctionReturn(PETSC_SUCCESS);
    232: }

    234: static PetscErrorCode PolyEval(PetscInt nroots, const PetscReal *r, const PetscReal *c, PetscReal x, PetscReal y, PetscReal *px, PetscReal *py)
    235: {
    236:   PetscInt  i;
    237:   PetscReal rprod = 1, iprod = 0;

    239:   PetscFunctionBegin;
    240:   for (i = 0; i < nroots; i++) {
    241:     PetscReal rnew = rprod * (x - r[i]) - iprod * (y - c[i]);
    242:     PetscReal inew = rprod * (y - c[i]) + iprod * (x - r[i]);
    243:     rprod          = rnew;
    244:     iprod          = inew;
    245:   }
    246:   *px = rprod;
    247:   *py = iprod;
    248:   PetscFunctionReturn(PETSC_SUCCESS);
    249: }

    251: #include <petscdraw.h>
    252: /* Collective */
    253: PetscErrorCode KSPPlotEigenContours_Private(KSP ksp, PetscInt neig, const PetscReal *r, const PetscReal *c)
    254: {
    255:   PetscReal     xmin, xmax, ymin, ymax, *xloc, *yloc, *value, px0, py0, rscale, iscale;
    256:   int           M, N, i, j;
    257:   PetscMPIInt   rank;
    258:   PetscViewer   viewer;
    259:   PetscDraw     draw;
    260:   PetscDrawAxis drawaxis;

    262:   PetscFunctionBegin;
    263:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ksp), &rank));
    264:   if (rank) PetscFunctionReturn(PETSC_SUCCESS);
    265:   M    = 80;
    266:   N    = 80;
    267:   xmin = r[0];
    268:   xmax = r[0];
    269:   ymin = c[0];
    270:   ymax = c[0];
    271:   for (i = 1; i < neig; i++) {
    272:     xmin = PetscMin(xmin, r[i]);
    273:     xmax = PetscMax(xmax, r[i]);
    274:     ymin = PetscMin(ymin, c[i]);
    275:     ymax = PetscMax(ymax, c[i]);
    276:   }
    277:   PetscCall(PetscMalloc3(M, &xloc, N, &yloc, M * N, &value));
    278:   for (i = 0; i < M; i++) xloc[i] = xmin - 0.1 * (xmax - xmin) + 1.2 * (xmax - xmin) * i / (M - 1);
    279:   for (i = 0; i < N; i++) yloc[i] = ymin - 0.1 * (ymax - ymin) + 1.2 * (ymax - ymin) * i / (N - 1);
    280:   PetscCall(PolyEval(neig, r, c, 0, 0, &px0, &py0));
    281:   rscale = px0 / (PetscSqr(px0) + PetscSqr(py0));
    282:   iscale = -py0 / (PetscSqr(px0) + PetscSqr(py0));
    283:   for (j = 0; j < N; j++) {
    284:     for (i = 0; i < M; i++) {
    285:       PetscReal px, py, tx, ty, tmod;
    286:       PetscCall(PolyEval(neig, r, c, xloc[i], yloc[j], &px, &py));
    287:       tx   = px * rscale - py * iscale;
    288:       ty   = py * rscale + px * iscale;
    289:       tmod = PetscSqr(tx) + PetscSqr(ty); /* modulus of the complex polynomial */
    290:       if (tmod > 1) tmod = 1.0;
    291:       if (tmod > 0.5 && tmod < 1) tmod = 0.5;
    292:       if (tmod > 0.2 && tmod < 0.5) tmod = 0.2;
    293:       if (tmod > 0.05 && tmod < 0.2) tmod = 0.05;
    294:       if (tmod < 1e-3) tmod = 1e-3;
    295:       value[i + j * M] = PetscLogReal(tmod) / PetscLogReal(10.0);
    296:     }
    297:   }
    298:   PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, NULL, "Iteratively Computed Eigen-contours", PETSC_DECIDE, PETSC_DECIDE, 450, 450, &viewer));
    299:   PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    300:   PetscCall(PetscDrawTensorContour(draw, M, N, NULL, NULL, value));
    301:   if (0) {
    302:     PetscCall(PetscDrawAxisCreate(draw, &drawaxis));
    303:     PetscCall(PetscDrawAxisSetLimits(drawaxis, xmin, xmax, ymin, ymax));
    304:     PetscCall(PetscDrawAxisSetLabels(drawaxis, "Eigen-counters", "real", "imag"));
    305:     PetscCall(PetscDrawAxisDraw(drawaxis));
    306:     PetscCall(PetscDrawAxisDestroy(&drawaxis));
    307:   }
    308:   PetscCall(PetscViewerDestroy(&viewer));
    309:   PetscCall(PetscFree3(xloc, yloc, value));
    310:   PetscFunctionReturn(PETSC_SUCCESS);
    311: }