Actual source code: precon.c

      1: /*
      2:     The PC (preconditioner) interface routines, callable by users.
      3: */
      4: #include <petsc/private/pcimpl.h>
      5: #include <petscdm.h>

      7: /* Logging support */
      8: PetscClassId  PC_CLASSID;
      9: PetscLogEvent PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_MatApply, PC_ApplyCoarse, PC_ApplySymmetricLeft;
     10: PetscLogEvent PC_ApplySymmetricRight, PC_ModifySubMatrices, PC_ApplyOnBlocks, PC_ApplyTransposeOnBlocks;
     11: PetscInt      PetscMGLevelId;
     12: PetscLogStage PCMPIStage;

     14: PETSC_INTERN PetscErrorCode PCGetDefaultType_Private(PC pc, const char *type[])
     15: {
     16:   PetscMPIInt size;
     17:   PetscBool   hasopblock, hasopsolve, flg1, flg2, set, flg3, isnormal;

     19:   PetscFunctionBegin;
     20:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
     21:   if (pc->pmat) {
     22:     PetscCall(MatHasOperation(pc->pmat, MATOP_GET_DIAGONAL_BLOCK, &hasopblock));
     23:     PetscCall(MatHasOperation(pc->pmat, MATOP_SOLVE, &hasopsolve));
     24:     if (size == 1) {
     25:       PetscCall(MatGetFactorAvailable(pc->pmat, "petsc", MAT_FACTOR_ICC, &flg1));
     26:       PetscCall(MatGetFactorAvailable(pc->pmat, "petsc", MAT_FACTOR_ILU, &flg2));
     27:       PetscCall(MatIsSymmetricKnown(pc->pmat, &set, &flg3));
     28:       PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &isnormal, MATNORMAL, MATNORMALHERMITIAN, NULL));
     29:       if (flg1 && (!flg2 || (set && flg3))) {
     30:         *type = PCICC;
     31:       } else if (flg2) {
     32:         *type = PCILU;
     33:       } else if (isnormal) {
     34:         *type = PCNONE;
     35:       } else if (hasopblock) { /* likely is a parallel matrix run on one processor */
     36:         if (pc->kspnestlevel > 0) {
     37:           Mat D;

     39:           PetscCall(MatGetDiagonalBlock(pc->pmat, &D));
     40:           PetscCall(PetscObjectTypeCompare((PetscObject)D, ((PetscObject)pc->pmat)->type_name, &flg1)); /* make sure there is no recursive call to PCGetDefaultType_Private() */
     41:         } else flg1 = PETSC_FALSE;
     42:         if (!flg1) *type = PCBJACOBI;
     43:         else *type = PCNONE;
     44:       } else if (hasopsolve) {
     45:         *type = PCMAT;
     46:       } else {
     47:         *type = PCNONE;
     48:       }
     49:     } else {
     50:       if (hasopblock) {
     51:         *type = PCBJACOBI;
     52:       } else if (hasopsolve) {
     53:         *type = PCMAT;
     54:       } else {
     55:         *type = PCNONE;
     56:       }
     57:     }
     58:   } else *type = NULL;
     59:   PetscFunctionReturn(PETSC_SUCCESS);
     60: }

     62: /* do not log solves, setup and applications of preconditioners while constructing preconditioners; perhaps they should be logged separately from the regular solves */
     63: PETSC_EXTERN PetscLogEvent KSP_Solve, KSP_SetUp;

     65: static PetscErrorCode PCLogEventsDeactivatePush(void)
     66: {
     67:   PetscFunctionBegin;
     68:   PetscCall(KSPInitializePackage());
     69:   PetscCall(PetscLogEventDeactivatePush(KSP_Solve));
     70:   PetscCall(PetscLogEventDeactivatePush(KSP_SetUp));
     71:   PetscCall(PetscLogEventDeactivatePush(PC_Apply));
     72:   PetscCall(PetscLogEventDeactivatePush(PC_SetUp));
     73:   PetscCall(PetscLogEventDeactivatePush(PC_SetUpOnBlocks));
     74:   PetscFunctionReturn(PETSC_SUCCESS);
     75: }

     77: static PetscErrorCode PCLogEventsDeactivatePop(void)
     78: {
     79:   PetscFunctionBegin;
     80:   PetscCall(KSPInitializePackage());
     81:   PetscCall(PetscLogEventDeactivatePop(KSP_Solve));
     82:   PetscCall(PetscLogEventDeactivatePop(KSP_SetUp));
     83:   PetscCall(PetscLogEventDeactivatePop(PC_Apply));
     84:   PetscCall(PetscLogEventDeactivatePop(PC_SetUp));
     85:   PetscCall(PetscLogEventDeactivatePop(PC_SetUpOnBlocks));
     86:   PetscFunctionReturn(PETSC_SUCCESS);
     87: }

     89: /*@
     90:   PCReset - Resets a `PC` context to the state it was in before `PCSetUp()` was called, and removes any allocated `Vec` and `Mat` from its data structure

     92:   Collective

     94:   Input Parameter:
     95: . pc - the `PC` preconditioner context

     97:   Level: developer

     99:   Notes:
    100:   Any options set, including those set with `KSPSetFromOptions()` remain.

    102:   This allows a `PC` to be reused for a different sized linear system but using the same options that have been previously set in `pc`

    104: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`
    105: @*/
    106: PetscErrorCode PCReset(PC pc)
    107: {
    108:   PetscFunctionBegin;
    110:   PetscTryTypeMethod(pc, reset);
    111:   PetscCall(VecDestroy(&pc->diagonalscaleright));
    112:   PetscCall(VecDestroy(&pc->diagonalscaleleft));
    113:   PetscCall(MatDestroy(&pc->pmat));
    114:   PetscCall(MatDestroy(&pc->mat));

    116:   pc->setupcalled = PETSC_FALSE;
    117:   PetscFunctionReturn(PETSC_SUCCESS);
    118: }

    120: /*@
    121:   PCDestroy - Destroys `PC` context that was created with `PCCreate()`.

    123:   Collective

    125:   Input Parameter:
    126: . pc - the `PC` preconditioner context

    128:   Level: developer

    130: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`
    131: @*/
    132: PetscErrorCode PCDestroy(PC *pc)
    133: {
    134:   PetscFunctionBegin;
    135:   if (!*pc) PetscFunctionReturn(PETSC_SUCCESS);
    137:   if (--((PetscObject)*pc)->refct > 0) {
    138:     *pc = NULL;
    139:     PetscFunctionReturn(PETSC_SUCCESS);
    140:   }

    142:   PetscCall(PCReset(*pc));

    144:   /* if memory was published with SAWs then destroy it */
    145:   PetscCall(PetscObjectSAWsViewOff((PetscObject)*pc));
    146:   PetscTryTypeMethod(*pc, destroy);
    147:   PetscCall(DMDestroy(&(*pc)->dm));
    148:   PetscCall(PetscHeaderDestroy(pc));
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: /*@
    153:   PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
    154:   scaling as needed by certain time-stepping codes.

    156:   Logically Collective

    158:   Input Parameter:
    159: . pc - the `PC` preconditioner context

    161:   Output Parameter:
    162: . flag - `PETSC_TRUE` if it applies the scaling

    164:   Level: developer

    166:   Note:
    167:   If this returns `PETSC_TRUE` then the system solved via the Krylov method is, for left and right preconditioning,

    169:   $$
    170:   \begin{align*}
    171:   D M A D^{-1} y = D M b  \\
    172:   D A M D^{-1} z = D b.
    173:   \end{align*}
    174:   $$

    176: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCSetDiagonalScale()`
    177: @*/
    178: PetscErrorCode PCGetDiagonalScale(PC pc, PetscBool *flag)
    179: {
    180:   PetscFunctionBegin;
    182:   PetscAssertPointer(flag, 2);
    183:   *flag = pc->diagonalscale;
    184:   PetscFunctionReturn(PETSC_SUCCESS);
    185: }

    187: /*@
    188:   PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
    189:   scaling as needed by certain time-stepping codes.

    191:   Logically Collective

    193:   Input Parameters:
    194: + pc - the `PC` preconditioner context
    195: - s  - scaling vector

    197:   Level: intermediate

    199:   Notes:
    200:   The system solved via the Krylov method is, for left and right preconditioning,
    201:   $$
    202:   \begin{align*}
    203:   D M A D^{-1} y = D M b \\
    204:   D A M D^{-1} z = D b.
    205:   \end{align*}
    206:   $$

    208:   `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.

    210: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCDiagonalScaleRight()`, `PCGetDiagonalScale()`
    211: @*/
    212: PetscErrorCode PCSetDiagonalScale(PC pc, Vec s)
    213: {
    214:   PetscFunctionBegin;
    217:   pc->diagonalscale = PETSC_TRUE;

    219:   PetscCall(PetscObjectReference((PetscObject)s));
    220:   PetscCall(VecDestroy(&pc->diagonalscaleleft));

    222:   pc->diagonalscaleleft = s;

    224:   PetscCall(VecDuplicate(s, &pc->diagonalscaleright));
    225:   PetscCall(VecCopy(s, pc->diagonalscaleright));
    226:   PetscCall(VecReciprocal(pc->diagonalscaleright));
    227:   PetscFunctionReturn(PETSC_SUCCESS);
    228: }

    230: /*@
    231:   PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.

    233:   Logically Collective

    235:   Input Parameters:
    236: + pc  - the `PC` preconditioner context
    237: . in  - input vector
    238: - out - scaled vector (maybe the same as in)

    240:   Level: intermediate

    242:   Notes:
    243:   The system solved via the Krylov method is, for left and right preconditioning,

    245:   $$
    246:   \begin{align*}
    247:   D M A D^{-1} y = D M b  \\
    248:   D A M D^{-1} z = D b.
    249:   \end{align*}
    250:   $$

    252:   `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.

    254:   If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`

    256: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCSetDiagonalScale()`, `PCDiagonalScaleRight()`, `MatDiagonalScale()`
    257: @*/
    258: PetscErrorCode PCDiagonalScaleLeft(PC pc, Vec in, Vec out)
    259: {
    260:   PetscFunctionBegin;
    264:   if (pc->diagonalscale) PetscCall(VecPointwiseMult(out, pc->diagonalscaleleft, in));
    265:   else if (in != out) PetscCall(VecCopy(in, out));
    266:   PetscFunctionReturn(PETSC_SUCCESS);
    267: }

    269: /*@
    270:   PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.

    272:   Logically Collective

    274:   Input Parameters:
    275: + pc  - the `PC` preconditioner context
    276: . in  - input vector
    277: - out - scaled vector (maybe the same as in)

    279:   Level: intermediate

    281:   Notes:
    282:   The system solved via the Krylov method is, for left and right preconditioning,

    284:   $$
    285:   \begin{align*}
    286:   D M A D^{-1} y = D M b  \\
    287:   D A M D^{-1} z = D b.
    288:   \end{align*}
    289:   $$

    291:   `PCDiagonalScaleLeft()` scales a vector by $D$. `PCDiagonalScaleRight()` scales a vector by $D^{-1}$.

    293:   If diagonal scaling is turned off and `in` is not `out` then `in` is copied to `out`

    295: .seealso: [](ch_ksp), `PCCreate()`, `PCSetUp()`, `PCDiagonalScaleLeft()`, `PCSetDiagonalScale()`, `MatDiagonalScale()`
    296: @*/
    297: PetscErrorCode PCDiagonalScaleRight(PC pc, Vec in, Vec out)
    298: {
    299:   PetscFunctionBegin;
    303:   if (pc->diagonalscale) {
    304:     PetscCall(VecPointwiseMult(out, pc->diagonalscaleright, in));
    305:   } else if (in != out) {
    306:     PetscCall(VecCopy(in, out));
    307:   }
    308:   PetscFunctionReturn(PETSC_SUCCESS);
    309: }

    311: /*@
    312:   PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
    313:   operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
    314:   `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.

    316:   Logically Collective

    318:   Input Parameters:
    319: + pc  - the `PC` preconditioner context
    320: - flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)

    322:   Options Database Key:
    323: . -pc_use_amat (true|false) - use the amat argument to `KSPSetOperators()` or `PCSetOperators()` to apply the operator

    325:   Level: intermediate

    327:   Note:
    328:   For the common case in which the linear system matrix and the matrix used to construct the
    329:   preconditioner are identical, this routine has no affect.

    331: .seealso: [](ch_ksp), `PC`, `PCGetUseAmat()`, `PCBJACOBI`, `PCMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`,
    332:           `KSPSetOperators()`, `PCSetOperators()`
    333: @*/
    334: PetscErrorCode PCSetUseAmat(PC pc, PetscBool flg)
    335: {
    336:   PetscFunctionBegin;
    338:   pc->useAmat = flg;
    339:   PetscFunctionReturn(PETSC_SUCCESS);
    340: }

    342: /*@
    343:   PCSetErrorIfFailure - Causes `PC` to generate an error if a floating point exception, for example a zero pivot, is detected.

    345:   Logically Collective

    347:   Input Parameters:
    348: + pc  - iterative context obtained from `PCCreate()`
    349: - flg - `PETSC_TRUE` indicates you want the error generated

    351:   Level: advanced

    353:   Notes:
    354:   Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call `KSPGetConvergedReason()` after a `KSPSolve()`
    355:   to determine if it has converged or failed. Or use -ksp_error_if_not_converged to cause the program to terminate as soon as lack of convergence is
    356:   detected.

    358:   This is propagated into `KSP`s used by this `PC`, which then propagate it into `PC`s used by those `KSP`s

    360: .seealso: [](ch_ksp), `PC`, `KSPSetErrorIfNotConverged()`, `PCGetInitialGuessNonzero()`, `PCSetInitialGuessKnoll()`, `PCGetInitialGuessKnoll()`
    361: @*/
    362: PetscErrorCode PCSetErrorIfFailure(PC pc, PetscBool flg)
    363: {
    364:   PetscFunctionBegin;
    367:   pc->erroriffailure = flg;
    368:   PetscFunctionReturn(PETSC_SUCCESS);
    369: }

    371: /*@
    372:   PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
    373:   operator during the preconditioning process it applies the Amat provided to `TSSetRHSJacobian()`,
    374:   `TSSetIJacobian()`, `SNESSetJacobian()`, `KSPSetOperators()` or `PCSetOperators()` not the Pmat.

    376:   Logically Collective

    378:   Input Parameter:
    379: . pc - the `PC` preconditioner context

    381:   Output Parameter:
    382: . flg - `PETSC_TRUE` to use the Amat, `PETSC_FALSE` to use the Pmat (default is false)

    384:   Level: intermediate

    386:   Note:
    387:   For the common case in which the linear system matrix and the matrix used to construct the
    388:   preconditioner are identical, this routine is does nothing.

    390: .seealso: [](ch_ksp), `PC`, `PCSetUseAmat()`, `PCBJACOBI`, `PCMG`, `PCFIELDSPLIT`, `PCCOMPOSITE`
    391: @*/
    392: PetscErrorCode PCGetUseAmat(PC pc, PetscBool *flg)
    393: {
    394:   PetscFunctionBegin;
    396:   *flg = pc->useAmat;
    397:   PetscFunctionReturn(PETSC_SUCCESS);
    398: }

    400: /*@
    401:   PCSetKSPNestLevel - sets the amount of nesting the `KSP` that contains this `PC` has

    403:   Collective

    405:   Input Parameters:
    406: + pc    - the `PC`
    407: - level - the nest level

    409:   Level: developer

    411: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPGetNestLevel()`, `PCGetKSPNestLevel()`, `KSPSetNestLevel()`
    412: @*/
    413: PetscErrorCode PCSetKSPNestLevel(PC pc, PetscInt level)
    414: {
    415:   PetscFunctionBegin;
    418:   pc->kspnestlevel = level;
    419:   PetscFunctionReturn(PETSC_SUCCESS);
    420: }

    422: /*@
    423:   PCGetKSPNestLevel - gets the amount of nesting the `KSP` that contains this `PC` has

    425:   Not Collective

    427:   Input Parameter:
    428: . pc - the `PC`

    430:   Output Parameter:
    431: . level - the nest level

    433:   Level: developer

    435: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPSetNestLevel()`, `PCSetKSPNestLevel()`, `KSPGetNestLevel()`
    436: @*/
    437: PetscErrorCode PCGetKSPNestLevel(PC pc, PetscInt *level)
    438: {
    439:   PetscFunctionBegin;
    441:   PetscAssertPointer(level, 2);
    442:   *level = pc->kspnestlevel;
    443:   PetscFunctionReturn(PETSC_SUCCESS);
    444: }

    446: /*@
    447:   PCCreate - Creates a preconditioner context, `PC`

    449:   Collective

    451:   Input Parameter:
    452: . comm - MPI communicator

    454:   Output Parameter:
    455: . newpc - location to put the `PC` preconditioner context

    457:   Level: developer

    459:   Notes:
    460:   This is rarely called directly by users since `KSP` manages the `PC` objects it uses. Use `KSPGetPC()` to access the `PC` used by a `KSP`.

    462:   Use `PCSetType()` or `PCSetFromOptions()` with the option `-pc_type pctype` to set the `PCType` for this `PC`

    464:   The default preconditioner type `PCType` for sparse matrices is `PCILU` or `PCICC` with 0 fill on one process and block Jacobi (`PCBJACOBI`) with `PCILU` or `PCICC`
    465:   in parallel. For dense matrices it is always `PCNONE`.

    467: .seealso: [](ch_ksp), `PC`, `PCType`, `PCSetType`, `PCSetUp()`, `PCApply()`, `PCDestroy()`, `KSP`, `KSPGetPC()`
    468: @*/
    469: PetscErrorCode PCCreate(MPI_Comm comm, PC *newpc)
    470: {
    471:   PC pc;

    473:   PetscFunctionBegin;
    474:   PetscAssertPointer(newpc, 2);
    475:   PetscCall(PCInitializePackage());

    477:   PetscCall(PetscHeaderCreate(pc, PC_CLASSID, "PC", "Preconditioner", "PC", comm, PCDestroy, PCView));
    478:   pc->mat                  = NULL;
    479:   pc->pmat                 = NULL;
    480:   pc->setupcalled          = PETSC_FALSE;
    481:   pc->setfromoptionscalled = 0;
    482:   pc->data                 = NULL;
    483:   pc->diagonalscale        = PETSC_FALSE;
    484:   pc->diagonalscaleleft    = NULL;
    485:   pc->diagonalscaleright   = NULL;

    487:   pc->modifysubmatrices  = NULL;
    488:   pc->modifysubmatricesP = NULL;

    490:   *newpc = pc;
    491:   PetscFunctionReturn(PETSC_SUCCESS);
    492: }

    494: /*@
    495:   PCApply - Applies the preconditioner to a vector.

    497:   Collective

    499:   Input Parameters:
    500: + pc - the `PC` preconditioner context
    501: - x  - input vector

    503:   Output Parameter:
    504: . y - output vector

    506:   Level: developer

    508: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`, `PCApplyBAorAB()`
    509: @*/
    510: PetscErrorCode PCApply(PC pc, Vec x, Vec y)
    511: {
    512:   PetscInt m, n, mv, nv;

    514:   PetscFunctionBegin;
    518:   PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
    519:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
    520:   /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
    521:   PetscCall(MatGetLocalSize(pc->pmat, &m, &n));
    522:   PetscCall(VecGetLocalSize(x, &mv));
    523:   PetscCall(VecGetLocalSize(y, &nv));
    524:   /* check pmat * y = x is feasible */
    525:   PetscCheck(mv == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Preconditioner number of local rows %" PetscInt_FMT " does not equal input vector size %" PetscInt_FMT, m, mv);
    526:   PetscCheck(nv == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Preconditioner number of local columns %" PetscInt_FMT " does not equal output vector size %" PetscInt_FMT, n, nv);
    527:   PetscCall(VecSetErrorIfLocked(y, 3));

    529:   PetscCall(PCSetUp(pc));
    530:   PetscCall(VecLockReadPush(x));
    531:   PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
    532:   PetscUseTypeMethod(pc, apply, x, y);
    533:   PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
    534:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
    535:   PetscCall(VecLockReadPop(x));
    536:   PetscFunctionReturn(PETSC_SUCCESS);
    537: }

    539: static PetscErrorCode PCMatApplyTranspose_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
    540: {
    541:   Mat       A;
    542:   Vec       cy, cx;
    543:   PetscInt  m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
    544:   PetscBool match;

    546:   PetscFunctionBegin;
    550:   PetscCheckSameComm(pc, 1, X, 2);
    551:   PetscCheckSameComm(pc, 1, Y, 3);
    552:   PetscCheck(Y != X, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "Y and X must be different matrices");
    553:   PetscCall(PCGetOperators(pc, NULL, &A));
    554:   PetscCall(MatGetLocalSize(A, &m3, &n3));
    555:   PetscCall(MatGetLocalSize(X, &m2, &n2));
    556:   PetscCall(MatGetLocalSize(Y, &m1, &n1));
    557:   PetscCall(MatGetSize(A, &M3, &N3));
    558:   PetscCall(MatGetSize(X, &M2, &N2));
    559:   PetscCall(MatGetSize(Y, &M1, &N1));
    560:   PetscCheck(n1 == n2 && N1 == N2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible number of columns between block of input vectors (n,N) = (%" PetscInt_FMT ",%" PetscInt_FMT ") and block of output vectors (n,N) = (%" PetscInt_FMT ",%" PetscInt_FMT ")", n2, N2, n1, N1);
    561:   PetscCheck(m2 == m3 && M2 == M3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout between block of input vectors (m,M) = (%" PetscInt_FMT ",%" PetscInt_FMT ") and Pmat (m,M)x(n,N) = (%" PetscInt_FMT ",%" PetscInt_FMT ")x(%" PetscInt_FMT ",%" PetscInt_FMT ")", m2, M2, m3, M3, n3, N3);
    562:   PetscCheck(m1 == n3 && M1 == N3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible layout between block of output vectors (m,M) = (%" PetscInt_FMT ",%" PetscInt_FMT ") and Pmat (m,M)x(n,N) = (%" PetscInt_FMT ",%" PetscInt_FMT ")x(%" PetscInt_FMT ",%" PetscInt_FMT ")", m1, M1, m3, M3, n3, N3);
    563:   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, ""));
    564:   PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of output vectors not stored in a dense Mat");
    565:   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, ""));
    566:   PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of input vectors not stored in a dense Mat");
    567:   PetscCall(PCSetUp(pc));
    568:   if (!transpose && pc->ops->matapply) {
    569:     PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
    570:     PetscUseTypeMethod(pc, matapply, X, Y);
    571:     PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
    572:   } else if (transpose && pc->ops->matapplytranspose) {
    573:     PetscCall(PetscLogEventBegin(PC_MatApply, pc, X, Y, 0));
    574:     PetscUseTypeMethod(pc, matapplytranspose, X, Y);
    575:     PetscCall(PetscLogEventEnd(PC_MatApply, pc, X, Y, 0));
    576:   } else {
    577:     PetscCall(PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name));
    578:     for (n1 = 0; n1 < N1; ++n1) {
    579:       PetscCall(MatDenseGetColumnVecRead(X, n1, &cx));
    580:       PetscCall(MatDenseGetColumnVecWrite(Y, n1, &cy));
    581:       if (!transpose) PetscCall(PCApply(pc, cx, cy));
    582:       else PetscCall(PCApplyTranspose(pc, cx, cy));
    583:       PetscCall(MatDenseRestoreColumnVecWrite(Y, n1, &cy));
    584:       PetscCall(MatDenseRestoreColumnVecRead(X, n1, &cx));
    585:     }
    586:   }
    587:   PetscFunctionReturn(PETSC_SUCCESS);
    588: }

    590: /*@
    591:   PCMatApply - Applies the preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApply()`, `Y` and `X` must be different matrices.

    593:   Collective

    595:   Input Parameters:
    596: + pc - the `PC` preconditioner context
    597: - X  - block of input vectors

    599:   Output Parameter:
    600: . Y - block of output vectors

    602:   Level: developer

    604: .seealso: [](ch_ksp), `PC`, `PCApply()`, `KSPMatSolve()`
    605: @*/
    606: PetscErrorCode PCMatApply(PC pc, Mat X, Mat Y)
    607: {
    608:   PetscFunctionBegin;
    609:   PetscCall(PCMatApplyTranspose_Private(pc, X, Y, PETSC_FALSE));
    610:   PetscFunctionReturn(PETSC_SUCCESS);
    611: }

    613: /*@
    614:   PCMatApplyTranspose - Applies the transpose of preconditioner to multiple vectors stored as a `MATDENSE`. Like `PCApplyTranspose()`, `Y` and `X` must be different matrices.

    616:   Collective

    618:   Input Parameters:
    619: + pc - the `PC` preconditioner context
    620: - X  - block of input vectors

    622:   Output Parameter:
    623: . Y - block of output vectors

    625:   Level: developer

    627: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`, `KSPMatSolveTranspose()`
    628: @*/
    629: PetscErrorCode PCMatApplyTranspose(PC pc, Mat X, Mat Y)
    630: {
    631:   PetscFunctionBegin;
    632:   PetscCall(PCMatApplyTranspose_Private(pc, X, Y, PETSC_TRUE));
    633:   PetscFunctionReturn(PETSC_SUCCESS);
    634: }

    636: /*@
    637:   PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.

    639:   Collective

    641:   Input Parameters:
    642: + pc - the `PC` preconditioner context
    643: - x  - input vector

    645:   Output Parameter:
    646: . y - output vector

    648:   Level: developer

    650:   Note:
    651:   Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.

    653: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricRight()`
    654: @*/
    655: PetscErrorCode PCApplySymmetricLeft(PC pc, Vec x, Vec y)
    656: {
    657:   PetscFunctionBegin;
    661:   PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
    662:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
    663:   PetscCall(PCSetUp(pc));
    664:   PetscCall(VecLockReadPush(x));
    665:   PetscCall(PetscLogEventBegin(PC_ApplySymmetricLeft, pc, x, y, 0));
    666:   PetscUseTypeMethod(pc, applysymmetricleft, x, y);
    667:   PetscCall(PetscLogEventEnd(PC_ApplySymmetricLeft, pc, x, y, 0));
    668:   PetscCall(VecLockReadPop(x));
    669:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
    670:   PetscFunctionReturn(PETSC_SUCCESS);
    671: }

    673: /*@
    674:   PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.

    676:   Collective

    678:   Input Parameters:
    679: + pc - the `PC` preconditioner context
    680: - x  - input vector

    682:   Output Parameter:
    683: . y - output vector

    685:   Level: developer

    687:   Note:
    688:   Currently, this routine is implemented only for `PCICC` and `PCJACOBI` preconditioners.

    690: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplySymmetricLeft()`
    691: @*/
    692: PetscErrorCode PCApplySymmetricRight(PC pc, Vec x, Vec y)
    693: {
    694:   PetscFunctionBegin;
    698:   PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
    699:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
    700:   PetscCall(PCSetUp(pc));
    701:   PetscCall(VecLockReadPush(x));
    702:   PetscCall(PetscLogEventBegin(PC_ApplySymmetricRight, pc, x, y, 0));
    703:   PetscUseTypeMethod(pc, applysymmetricright, x, y);
    704:   PetscCall(PetscLogEventEnd(PC_ApplySymmetricRight, pc, x, y, 0));
    705:   PetscCall(VecLockReadPop(x));
    706:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
    707:   PetscFunctionReturn(PETSC_SUCCESS);
    708: }

    710: /*@
    711:   PCApplyTranspose - Applies the transpose of preconditioner to a vector.

    713:   Collective

    715:   Input Parameters:
    716: + pc - the `PC` preconditioner context
    717: - x  - input vector

    719:   Output Parameter:
    720: . y - output vector

    722:   Level: developer

    724:   Note:
    725:   For complex numbers this applies the non-Hermitian transpose.

    727:   Developer Note:
    728:   We need to implement a `PCApplyHermitianTranspose()`

    730: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyBAorAB()`, `PCApplyBAorABTranspose()`, `PCApplyTransposeExists()`
    731: @*/
    732: PetscErrorCode PCApplyTranspose(PC pc, Vec x, Vec y)
    733: {
    734:   PetscFunctionBegin;
    738:   PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
    739:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 2, PETSC_TRUE));
    740:   PetscCall(PCSetUp(pc));
    741:   PetscCall(VecLockReadPush(x));
    742:   PetscCall(PetscLogEventBegin(PC_Apply, pc, x, y, 0));
    743:   PetscUseTypeMethod(pc, applytranspose, x, y);
    744:   PetscCall(PetscLogEventEnd(PC_Apply, pc, x, y, 0));
    745:   PetscCall(VecLockReadPop(x));
    746:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 3, PETSC_FALSE));
    747:   PetscFunctionReturn(PETSC_SUCCESS);
    748: }

    750: /*@
    751:   PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation

    753:   Collective

    755:   Input Parameter:
    756: . pc - the `PC` preconditioner context

    758:   Output Parameter:
    759: . flg - `PETSC_TRUE` if a transpose operation is defined

    761:   Level: developer

    763: .seealso: [](ch_ksp), `PC`, `PCApplyTranspose()`
    764: @*/
    765: PetscErrorCode PCApplyTransposeExists(PC pc, PetscBool *flg)
    766: {
    767:   PetscFunctionBegin;
    769:   PetscAssertPointer(flg, 2);
    770:   if (pc->ops->applytranspose) *flg = PETSC_TRUE;
    771:   else *flg = PETSC_FALSE;
    772:   PetscFunctionReturn(PETSC_SUCCESS);
    773: }

    775: /*@
    776:   PCApplyBAorAB - Applies the preconditioner and operator to a vector. $y = B*A*x $ or $ y = A*B*x$.

    778:   Collective

    780:   Input Parameters:
    781: + pc   - the `PC` preconditioner context
    782: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
    783: . x    - input vector
    784: - work - work vector

    786:   Output Parameter:
    787: . y - output vector

    789:   Level: developer

    791:   Note:
    792:   If the `PC` has had `PCSetDiagonalScale()` set then $ D M A D^{-1} $ for left preconditioning or $ D A M D^{-1} $ is actually applied.
    793:   The specific `KSPSolve()` method must also be written to handle the post-solve "correction" for the diagonal scaling.

    795: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorABTranspose()`
    796: @*/
    797: PetscErrorCode PCApplyBAorAB(PC pc, PCSide side, Vec x, Vec y, Vec work)
    798: {
    799:   PetscFunctionBegin;
    805:   PetscCheckSameComm(pc, 1, x, 3);
    806:   PetscCheckSameComm(pc, 1, y, 4);
    807:   PetscCheckSameComm(pc, 1, work, 5);
    808:   PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
    809:   PetscCheck(side == PC_LEFT || side == PC_SYMMETRIC || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right, left, or symmetric");
    810:   PetscCheck(!pc->diagonalscale || side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot include diagonal scaling with symmetric preconditioner application");
    811:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));

    813:   PetscCall(PCSetUp(pc));
    814:   if (pc->diagonalscale) {
    815:     if (pc->ops->applyBA) {
    816:       Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
    817:       PetscCall(VecDuplicate(x, &work2));
    818:       PetscCall(PCDiagonalScaleRight(pc, x, work2));
    819:       PetscUseTypeMethod(pc, applyBA, side, work2, y, work);
    820:       PetscCall(PCDiagonalScaleLeft(pc, y, y));
    821:       PetscCall(VecDestroy(&work2));
    822:     } else if (side == PC_RIGHT) {
    823:       PetscCall(PCDiagonalScaleRight(pc, x, y));
    824:       PetscCall(PCApply(pc, y, work));
    825:       PetscCall(MatMult(pc->mat, work, y));
    826:       PetscCall(PCDiagonalScaleLeft(pc, y, y));
    827:     } else if (side == PC_LEFT) {
    828:       PetscCall(PCDiagonalScaleRight(pc, x, y));
    829:       PetscCall(MatMult(pc->mat, y, work));
    830:       PetscCall(PCApply(pc, work, y));
    831:       PetscCall(PCDiagonalScaleLeft(pc, y, y));
    832:     } else PetscCheck(side != PC_SYMMETRIC, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot provide diagonal scaling with symmetric application of preconditioner");
    833:   } else {
    834:     if (pc->ops->applyBA) {
    835:       PetscUseTypeMethod(pc, applyBA, side, x, y, work);
    836:     } else if (side == PC_RIGHT) {
    837:       PetscCall(PCApply(pc, x, work));
    838:       PetscCall(MatMult(pc->mat, work, y));
    839:     } else if (side == PC_LEFT) {
    840:       PetscCall(MatMult(pc->mat, x, work));
    841:       PetscCall(PCApply(pc, work, y));
    842:     } else if (side == PC_SYMMETRIC) {
    843:       /* There's an extra copy here; maybe should provide 2 work vectors instead? */
    844:       PetscCall(PCApplySymmetricRight(pc, x, work));
    845:       PetscCall(MatMult(pc->mat, work, y));
    846:       PetscCall(VecCopy(y, work));
    847:       PetscCall(PCApplySymmetricLeft(pc, work, y));
    848:     }
    849:   }
    850:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
    851:   PetscFunctionReturn(PETSC_SUCCESS);
    852: }

    854: /*@
    855:   PCApplyBAorABTranspose - Applies the transpose of the preconditioner
    856:   and operator to a vector. That is, applies $B^T * A^T$ with left preconditioning,
    857:   NOT $(B*A)^T = A^T*B^T$.

    859:   Collective

    861:   Input Parameters:
    862: + pc   - the `PC` preconditioner context
    863: . side - indicates the preconditioner side, one of `PC_LEFT`, `PC_RIGHT`, or `PC_SYMMETRIC`
    864: . x    - input vector
    865: - work - work vector

    867:   Output Parameter:
    868: . y - output vector

    870:   Level: developer

    872:   Note:
    873:   This routine is used internally so that the same Krylov code can be used to solve $A x = b$ and $A^T x = b$, with a preconditioner
    874:   defined by $B^T$. This is why this has the funny form that it computes $B^T * A^T$

    876: .seealso: [](ch_ksp), `PC`, `PCApply()`, `PCApplyTranspose()`, `PCApplyBAorAB()`
    877: @*/
    878: PetscErrorCode PCApplyBAorABTranspose(PC pc, PCSide side, Vec x, Vec y, Vec work)
    879: {
    880:   PetscFunctionBegin;
    885:   PetscCheck(x != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "x and y must be different vectors");
    886:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(x, 3, PETSC_TRUE));
    887:   if (pc->ops->applyBAtranspose) {
    888:     PetscUseTypeMethod(pc, applyBAtranspose, side, x, y, work);
    889:     if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
    890:     PetscFunctionReturn(PETSC_SUCCESS);
    891:   }
    892:   PetscCheck(side == PC_LEFT || side == PC_RIGHT, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Side must be right or left");

    894:   PetscCall(PCSetUp(pc));
    895:   if (side == PC_RIGHT) {
    896:     PetscCall(PCApplyTranspose(pc, x, work));
    897:     PetscCall(MatMultTranspose(pc->mat, work, y));
    898:   } else if (side == PC_LEFT) {
    899:     PetscCall(MatMultTranspose(pc->mat, x, work));
    900:     PetscCall(PCApplyTranspose(pc, work, y));
    901:   }
    902:   /* add support for PC_SYMMETRIC */
    903:   if (pc->erroriffailure) PetscCall(VecValidValues_Internal(y, 4, PETSC_FALSE));
    904:   PetscFunctionReturn(PETSC_SUCCESS);
    905: }

    907: /*@
    908:   PCApplyRichardsonExists - Determines whether a particular preconditioner has a
    909:   built-in fast application of Richardson's method.

    911:   Not Collective

    913:   Input Parameter:
    914: . pc - the preconditioner

    916:   Output Parameter:
    917: . exists - `PETSC_TRUE` or `PETSC_FALSE`

    919:   Level: developer

    921: .seealso: [](ch_ksp), `PC`, `KSPRICHARDSON`, `PCApplyRichardson()`
    922: @*/
    923: PetscErrorCode PCApplyRichardsonExists(PC pc, PetscBool *exists)
    924: {
    925:   PetscFunctionBegin;
    927:   PetscAssertPointer(exists, 2);
    928:   if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
    929:   else *exists = PETSC_FALSE;
    930:   PetscFunctionReturn(PETSC_SUCCESS);
    931: }

    933: /*@
    934:   PCApplyRichardson - Applies several steps of Richardson iteration with
    935:   the particular preconditioner. This routine is usually used by the
    936:   Krylov solvers and not the application code directly.

    938:   Collective

    940:   Input Parameters:
    941: + pc        - the `PC` preconditioner context
    942: . b         - the right-hand side
    943: . w         - one work vector
    944: . rtol      - relative decrease in residual norm convergence criteria
    945: . abstol    - absolute residual norm convergence criteria
    946: . dtol      - divergence residual norm increase criteria
    947: . its       - the number of iterations to apply.
    948: - guesszero - if the input x contains nonzero initial guess

    950:   Output Parameters:
    951: + outits - number of iterations actually used (for SOR this always equals its)
    952: . reason - the reason the apply terminated
    953: - y      - the solution (also contains initial guess if guesszero is `PETSC_FALSE`

    955:   Level: developer

    957:   Notes:
    958:   Most preconditioners do not support this function. Use the command
    959:   `PCApplyRichardsonExists()` to determine if one does.

    961:   Except for the `PCMG` this routine ignores the convergence tolerances
    962:   and always runs for the number of iterations

    964: .seealso: [](ch_ksp), `PC`, `PCApplyRichardsonExists()`
    965: @*/
    966: PetscErrorCode PCApplyRichardson(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason)
    967: {
    968:   PetscFunctionBegin;
    973:   PetscCheck(b != y, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "b and y must be different vectors");
    974:   PetscCall(PCSetUp(pc));
    975:   PetscUseTypeMethod(pc, applyrichardson, b, y, w, rtol, abstol, dtol, its, guesszero, outits, reason);
    976:   PetscFunctionReturn(PETSC_SUCCESS);
    977: }

    979: /*@
    980:   PCSetFailedReason - Sets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail

    982:   Logically Collective

    984:   Input Parameters:
    985: + pc     - the `PC` preconditioner context
    986: - reason - the reason it failed

    988:   Level: advanced

    990: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCFailedReason`
    991: @*/
    992: PetscErrorCode PCSetFailedReason(PC pc, PCFailedReason reason)
    993: {
    994:   PetscFunctionBegin;
    996:   pc->failedreason = reason;
    997:   PetscFunctionReturn(PETSC_SUCCESS);
    998: }

    1000: /*@
    1001:   PCGetFailedReason - Gets the reason a `PCSetUp()` failed or `PC_NOERROR` if it did not fail

    1003:   Not Collective

    1005:   Input Parameter:
    1006: . pc - the `PC` preconditioner context

    1008:   Output Parameter:
    1009: . reason - the reason it failed

    1011:   Level: advanced

    1013:   Note:
    1014:   After a call to `KSPCheckDot()` or  `KSPCheckNorm()` inside a `KSPSolve()` or a call to `PCReduceFailedReason()`
    1015:   this is the maximum reason over all MPI processes in the `PC` communicator and hence logically collective.
    1016:   Otherwise it returns the local value.

    1018: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCSetFailedReason()`, `PCFailedReason`
    1019: @*/
    1020: PetscErrorCode PCGetFailedReason(PC pc, PCFailedReason *reason)
    1021: {
    1022:   PetscFunctionBegin;
    1024:   *reason = pc->failedreason;
    1025:   PetscFunctionReturn(PETSC_SUCCESS);
    1026: }

    1028: /*@
    1029:   PCReduceFailedReason - Reduce the failed reason among the MPI processes that share the `PC`

    1031:   Collective

    1033:   Input Parameter:
    1034: . pc - the `PC` preconditioner context

    1036:   Level: advanced

    1038:   Note:
    1039:   Different MPI processes may have different reasons or no reason, see `PCGetFailedReason()`. This routine
    1040:   makes them have a common value (failure if any MPI process had a failure).

    1042: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `PCGetFailedReason()`, `PCSetFailedReason()`, `PCFailedReason`
    1043: @*/
    1044: PetscErrorCode PCReduceFailedReason(PC pc)
    1045: {
    1046:   PetscInt buf;

    1048:   PetscFunctionBegin;
    1050:   buf = (PetscInt)pc->failedreason;
    1051:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &buf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
    1052:   pc->failedreason = (PCFailedReason)buf;
    1053:   PetscFunctionReturn(PETSC_SUCCESS);
    1054: }

    1056: /*
    1057:       a setupcall of 0 indicates never setup,
    1058:                      1 indicates has been previously setup
    1059:                     -1 indicates a PCSetUp() was attempted and failed
    1060: */
    1061: /*@
    1062:   PCSetUp - Prepares for the use of a preconditioner. Performs all the one-time operations needed before the preconditioner
    1063:   can be used with `PCApply()`

    1065:   Collective

    1067:   Input Parameter:
    1068: . pc - the `PC` preconditioner context

    1070:   Level: developer

    1072:   Notes:
    1073:   For example, for `PCLU` this will compute the factorization.

    1075:   This is called automatically by `KSPSetUp()` or `PCApply()` so rarely needs to be called directly.

    1077:   For nested preconditioners, such as `PCFIELDSPLIT` or `PCBJACOBI` this may not finish the construction of the preconditioner
    1078:   on the inner levels, the routine `PCSetUpOnBlocks()` may compute more of the preconditioner in those situations.

    1080: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCApply()`, `PCDestroy()`, `KSPSetUp()`, `PCSetUpOnBlocks()`
    1081: @*/
    1082: PetscErrorCode PCSetUp(PC pc)
    1083: {
    1084:   const char      *def;
    1085:   PetscObjectState matstate, matnonzerostate;

    1087:   PetscFunctionBegin;
    1089:   PetscCheck(pc->mat, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Matrix must be set first");

    1091:   if (pc->setupcalled && pc->reusepreconditioner) {
    1092:     PetscCall(PetscInfo(pc, "Leaving PC with identical preconditioner since reuse preconditioner is set\n"));
    1093:     PetscFunctionReturn(PETSC_SUCCESS);
    1094:   }

    1096:   PetscCall(PetscObjectStateGet((PetscObject)pc->pmat, &matstate));
    1097:   PetscCall(MatGetNonzeroState(pc->pmat, &matnonzerostate));
    1098:   if (!pc->setupcalled) {
    1099:     //PetscCall(PetscInfo(pc, "Setting up PC for first time\n"));
    1100:     pc->flag = DIFFERENT_NONZERO_PATTERN;
    1101:   } else if (matstate == pc->matstate) PetscFunctionReturn(PETSC_SUCCESS);
    1102:   else {
    1103:     if (matnonzerostate != pc->matnonzerostate) {
    1104:       PetscCall(PetscInfo(pc, "Setting up PC with different nonzero pattern\n"));
    1105:       pc->flag = DIFFERENT_NONZERO_PATTERN;
    1106:     } else {
    1107:       //PetscCall(PetscInfo(pc, "Setting up PC with same nonzero pattern\n"));
    1108:       pc->flag = SAME_NONZERO_PATTERN;
    1109:     }
    1110:   }
    1111:   pc->matstate        = matstate;
    1112:   pc->matnonzerostate = matnonzerostate;

    1114:   if (!((PetscObject)pc)->type_name) {
    1115:     PetscCall(PCGetDefaultType_Private(pc, &def));
    1116:     PetscCall(PCSetType(pc, def));
    1117:   }

    1119:   PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
    1120:   PetscCall(MatSetErrorIfFailure(pc->mat, pc->erroriffailure));
    1121:   PetscCall(PetscLogEventBegin(PC_SetUp, pc, 0, 0, 0));
    1122:   if (pc->ops->setup) {
    1123:     PetscCall(PCLogEventsDeactivatePush());
    1124:     PetscUseTypeMethod(pc, setup);
    1125:     PetscCall(PCLogEventsDeactivatePop());
    1126:   }
    1127:   PetscCall(PetscLogEventEnd(PC_SetUp, pc, 0, 0, 0));
    1128:   if (pc->postsetup) PetscCall((*pc->postsetup)(pc));
    1129:   if (!pc->setupcalled) pc->setupcalled = PETSC_TRUE;
    1130:   PetscFunctionReturn(PETSC_SUCCESS);
    1131: }

    1133: /*@
    1134:   PCSetUpOnBlocks - Sets up the preconditioner for each block in
    1135:   the block Jacobi, overlapping Schwarz, and fieldsplit methods.

    1137:   Collective

    1139:   Input Parameter:
    1140: . pc - the `PC` preconditioner context

    1142:   Level: developer

    1144:   Notes:
    1145:   For nested preconditioners such as `PCBJACOBI`, `PCSetUp()` is not called on each sub-`KSP` when `PCSetUp()` is
    1146:   called on the outer `PC`, this routine ensures it is called.

    1148:   It calls `PCSetUp()` if not yet called.

    1150: .seealso: [](ch_ksp), `PC`, `PCSetUp()`, `PCCreate()`, `PCApply()`, `PCDestroy()`
    1151: @*/
    1152: PetscErrorCode PCSetUpOnBlocks(PC pc)
    1153: {
    1154:   PetscFunctionBegin;
    1156:   if (!pc->setupcalled) PetscCall(PCSetUp(pc)); /* "if" to prevent -info extra prints */
    1157:   if (!pc->ops->setuponblocks) PetscFunctionReturn(PETSC_SUCCESS);
    1158:   PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
    1159:   PetscCall(PetscLogEventBegin(PC_SetUpOnBlocks, pc, 0, 0, 0));
    1160:   PetscCall(PCLogEventsDeactivatePush());
    1161:   PetscUseTypeMethod(pc, setuponblocks);
    1162:   PetscCall(PCLogEventsDeactivatePop());
    1163:   PetscCall(PetscLogEventEnd(PC_SetUpOnBlocks, pc, 0, 0, 0));
    1164:   PetscFunctionReturn(PETSC_SUCCESS);
    1165: }

    1167: /*@C
    1168:   PCSetModifySubMatrices - Sets a user-defined routine for modifying the
    1169:   submatrices that arise within certain subdomain-based preconditioners such as `PCASM`

    1171:   Logically Collective

    1173:   Input Parameters:
    1174: + pc   - the `PC` preconditioner context
    1175: . func - routine for modifying the submatrices, see `PCModifySubMatricesFn`
    1176: - ctx  - optional user-defined context (may be `NULL`)

    1178:   Level: advanced

    1180:   Notes:
    1181:   The basic submatrices are extracted from the matrix used to construct the preconditioner as
    1182:   usual; the user can then alter these (for example, to set different boundary
    1183:   conditions for each submatrix) before they are used for the local solves.

    1185:   `PCSetModifySubMatrices()` MUST be called before `KSPSetUp()` and
    1186:   `KSPSolve()`.

    1188:   A routine set by `PCSetModifySubMatrices()` is currently called within
    1189:   `PCBJACOBI`, `PCASM`, `PCGASM`, and `PCHPDDM`.
    1190:   All other preconditioners ignore this routine.

    1192: .seealso: [](ch_ksp), `PC`, `PCModifySubMatricesFn`, `PCBJACOBI`, `PCASM`, `PCModifySubMatrices()`
    1193: @*/
    1194: PetscErrorCode PCSetModifySubMatrices(PC pc, PCModifySubMatricesFn *func, PetscCtx ctx)
    1195: {
    1196:   PetscFunctionBegin;
    1198:   pc->modifysubmatrices  = func;
    1199:   pc->modifysubmatricesP = ctx;
    1200:   PetscFunctionReturn(PETSC_SUCCESS);
    1201: }

    1203: /*@C
    1204:   PCModifySubMatrices - Calls an optional user-defined routine within
    1205:   certain preconditioners if one has been set with `PCSetModifySubMatrices()`.

    1207:   Collective

    1209:   Input Parameters:
    1210: + pc     - the `PC` preconditioner context
    1211: . nsub   - the number of local submatrices
    1212: . row    - an array of index sets that contain the global row numbers
    1213:          that comprise each local submatrix
    1214: . col    - an array of index sets that contain the global column numbers
    1215:          that comprise each local submatrix
    1216: . submat - array of local submatrices
    1217: - ctx    - optional user-defined context for private data for the
    1218:          user-defined routine (may be `NULL`)

    1220:   Output Parameter:
    1221: . submat - array of local submatrices (the entries of which may
    1222:             have been modified)

    1224:   Level: developer

    1226:   Note:
    1227:   The user should NOT generally call this routine, as it will
    1228:   automatically be called within certain preconditioners.

    1230: .seealso: [](ch_ksp), `PC`, `PCModifySubMatricesFn`, `PCSetModifySubMatrices()`
    1231: @*/
    1232: PetscErrorCode PCModifySubMatrices(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], PetscCtx ctx)
    1233: {
    1234:   PetscFunctionBegin;
    1236:   if (!pc->modifysubmatrices) PetscFunctionReturn(PETSC_SUCCESS);
    1237:   PetscCall(PetscLogEventBegin(PC_ModifySubMatrices, pc, 0, 0, 0));
    1238:   PetscCall((*pc->modifysubmatrices)(pc, nsub, row, col, submat, ctx));
    1239:   PetscCall(PetscLogEventEnd(PC_ModifySubMatrices, pc, 0, 0, 0));
    1240:   PetscFunctionReturn(PETSC_SUCCESS);
    1241: }

    1243: /*@
    1244:   PCSetOperators - Sets the matrix associated with the linear system and
    1245:   a (possibly) different one from which the preconditioner will be constructed.

    1247:   Logically Collective

    1249:   Input Parameters:
    1250: + pc   - the `PC` preconditioner context
    1251: . Amat - the matrix that defines the linear system
    1252: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

    1254:   Level: advanced

    1256:   Notes:
    1257:   Using this routine directly is rarely needed, the preferred, and equivalent, usage is `KSPSetOperators()`.

    1259:   Passing a `NULL` for `Amat` or `Pmat` removes the matrix that is currently used.

    1261:   If you wish to replace either `Amat` or `Pmat` but leave the other one untouched then
    1262:   first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
    1263:   on it and then pass it back in in your call to `KSPSetOperators()`.

    1265:   More Notes about Repeated Solution of Linear Systems:
    1266:   PETSc does NOT reset the matrix entries of either `Amat` or `Pmat`
    1267:   to zero after a linear solve; the user is completely responsible for
    1268:   matrix assembly.  See the routine `MatZeroEntries()` if desiring to
    1269:   zero all elements of a matrix.

    1271: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`
    1272:  @*/
    1273: PetscErrorCode PCSetOperators(PC pc, Mat Amat, Mat Pmat)
    1274: {
    1275:   PetscInt m1, n1, m2, n2;

    1277:   PetscFunctionBegin;
    1281:   if (Amat) PetscCheckSameComm(pc, 1, Amat, 2);
    1282:   if (Pmat) PetscCheckSameComm(pc, 1, Pmat, 3);
    1283:   if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
    1284:     PetscCall(MatGetLocalSize(Amat, &m1, &n1));
    1285:     PetscCall(MatGetLocalSize(pc->mat, &m2, &n2));
    1286:     PetscCheck(m1 == m2 && n1 == n2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot change local size of Amat after use old sizes %" PetscInt_FMT " %" PetscInt_FMT " new sizes %" PetscInt_FMT " %" PetscInt_FMT, m2, n2, m1, n1);
    1287:     PetscCall(MatGetLocalSize(Pmat, &m1, &n1));
    1288:     PetscCall(MatGetLocalSize(pc->pmat, &m2, &n2));
    1289:     PetscCheck(m1 == m2 && n1 == n2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot change local size of Pmat after use old sizes %" PetscInt_FMT " %" PetscInt_FMT " new sizes %" PetscInt_FMT " %" PetscInt_FMT, m2, n2, m1, n1);
    1290:   }

    1292:   if (Pmat != pc->pmat) {
    1293:     /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
    1294:     pc->matnonzerostate = -1;
    1295:     pc->matstate        = -1;
    1296:   }

    1298:   /* reference first in case the matrices are the same */
    1299:   PetscCall(PetscObjectReference((PetscObject)Amat));
    1300:   PetscCall(MatDestroy(&pc->mat));
    1301:   PetscCall(PetscObjectReference((PetscObject)Pmat));
    1302:   PetscCall(MatDestroy(&pc->pmat));
    1303:   pc->mat  = Amat;
    1304:   pc->pmat = Pmat;
    1305:   PetscFunctionReturn(PETSC_SUCCESS);
    1306: }

    1308: /*@
    1309:   PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner `PC` has changed.

    1311:   Logically Collective

    1313:   Input Parameters:
    1314: + pc   - the `PC` preconditioner context
    1315: - flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner

    1317:   Level: intermediate

    1319:   Note:
    1320:   Normally if a matrix inside a `PC` changes the `PC` automatically updates itself using information from the changed matrix. This option
    1321:   prevents this.

    1323: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCGetReusePreconditioner()`, `KSPSetReusePreconditioner()`
    1324:  @*/
    1325: PetscErrorCode PCSetReusePreconditioner(PC pc, PetscBool flag)
    1326: {
    1327:   PetscFunctionBegin;
    1330:   pc->reusepreconditioner = flag;
    1331:   PetscTryMethod(pc, "PCSetReusePreconditioner_C", (PC, PetscBool), (pc, flag));
    1332:   PetscFunctionReturn(PETSC_SUCCESS);
    1333: }

    1335: /*@
    1336:   PCGetReusePreconditioner - Determines if the `PC` reuses the current preconditioner even if the operator in the preconditioner has changed.

    1338:   Not Collective

    1340:   Input Parameter:
    1341: . pc - the `PC` preconditioner context

    1343:   Output Parameter:
    1344: . flag - `PETSC_TRUE` do not compute a new preconditioner, `PETSC_FALSE` do compute a new preconditioner

    1346:   Level: intermediate

    1348: .seealso: [](ch_ksp), `PC`, `PCGetOperators()`, `MatZeroEntries()`, `PCSetReusePreconditioner()`
    1349:  @*/
    1350: PetscErrorCode PCGetReusePreconditioner(PC pc, PetscBool *flag)
    1351: {
    1352:   PetscFunctionBegin;
    1354:   PetscAssertPointer(flag, 2);
    1355:   *flag = pc->reusepreconditioner;
    1356:   PetscFunctionReturn(PETSC_SUCCESS);
    1357: }

    1359: /*@
    1360:   PCGetOperators - Gets the matrix associated with the linear system and
    1361:   possibly a different one which is used to construct the preconditioner.

    1363:   Not Collective, though parallel `Mat`s are returned if `pc` is parallel

    1365:   Input Parameter:
    1366: . pc - the `PC` preconditioner context

    1368:   Output Parameters:
    1369: + Amat - the matrix defining the linear system
    1370: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.

    1372:   Level: intermediate

    1374:   Note:
    1375:   Does not increase the reference count of the matrices, so you should not destroy them

    1377:   Alternative usage: If the operators have NOT been set with `KSPSetOperators()`/`PCSetOperators()` then the operators
    1378:   are created in `PC` and returned to the user. In this case, if both operators
    1379:   mat and pmat are requested, two DIFFERENT operators will be returned. If
    1380:   only one is requested both operators in the PC will be the same (i.e. as
    1381:   if one had called `KSPSetOperators()`/`PCSetOperators()` with the same argument for both Mats).
    1382:   The user must set the sizes of the returned matrices and their type etc just
    1383:   as if the user created them with `MatCreate()`. For example,

    1385: .vb
    1386:          KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
    1387:            set size, type, etc of Amat

    1389:          MatCreate(comm,&mat);
    1390:          KSP/PCSetOperators(ksp/pc,Amat,Amat);
    1391:          PetscObjectDereference((PetscObject)mat);
    1392:            set size, type, etc of Amat
    1393: .ve

    1395:   and

    1397: .vb
    1398:          KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
    1399:            set size, type, etc of Amat and Pmat

    1401:          MatCreate(comm,&Amat);
    1402:          MatCreate(comm,&Pmat);
    1403:          KSP/PCSetOperators(ksp/pc,Amat,Pmat);
    1404:          PetscObjectDereference((PetscObject)Amat);
    1405:          PetscObjectDereference((PetscObject)Pmat);
    1406:            set size, type, etc of Amat and Pmat
    1407: .ve

    1409:   The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
    1410:   of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
    1411:   managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
    1412:   at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
    1413:   the `SNES` object (the `SNES` object manages it for you). Similarly when you create a KSP
    1414:   you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
    1415:   Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
    1416:   it can be created for you?

    1418: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperatorsSet()`
    1419: @*/
    1420: PetscErrorCode PCGetOperators(PC pc, Mat *Amat, Mat *Pmat)
    1421: {
    1422:   PetscFunctionBegin;
    1424:   if (Amat) {
    1425:     if (!pc->mat) {
    1426:       if (pc->pmat && !Pmat) { /* Pmat has been set, but user did not request it, so use for Amat */
    1427:         pc->mat = pc->pmat;
    1428:         PetscCall(PetscObjectReference((PetscObject)pc->mat));
    1429:       } else { /* both Amat and Pmat are empty */
    1430:         PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->mat));
    1431:         if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
    1432:           pc->pmat = pc->mat;
    1433:           PetscCall(PetscObjectReference((PetscObject)pc->pmat));
    1434:         }
    1435:       }
    1436:     }
    1437:     *Amat = pc->mat;
    1438:   }
    1439:   if (Pmat) {
    1440:     if (!pc->pmat) {
    1441:       if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
    1442:         pc->pmat = pc->mat;
    1443:         PetscCall(PetscObjectReference((PetscObject)pc->pmat));
    1444:       } else {
    1445:         PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &pc->pmat));
    1446:         if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
    1447:           pc->mat = pc->pmat;
    1448:           PetscCall(PetscObjectReference((PetscObject)pc->mat));
    1449:         }
    1450:       }
    1451:     }
    1452:     *Pmat = pc->pmat;
    1453:   }
    1454:   PetscFunctionReturn(PETSC_SUCCESS);
    1455: }

    1457: /*@
    1458:   PCGetOperatorsSet - Determines if the matrix associated with the linear system and
    1459:   possibly a different one associated with the preconditioner have been set in the `PC`.

    1461:   Not Collective, though the results on all processes should be the same

    1463:   Input Parameter:
    1464: . pc - the `PC` preconditioner context

    1466:   Output Parameters:
    1467: + mat  - the matrix associated with the linear system was set
    1468: - pmat - matrix associated with the preconditioner was set, usually the same

    1470:   Level: intermediate

    1472: .seealso: [](ch_ksp), `PC`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`
    1473: @*/
    1474: PetscErrorCode PCGetOperatorsSet(PC pc, PetscBool *mat, PetscBool *pmat)
    1475: {
    1476:   PetscFunctionBegin;
    1478:   if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
    1479:   if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
    1480:   PetscFunctionReturn(PETSC_SUCCESS);
    1481: }

    1483: /*@
    1484:   PCFactorGetMatrix - Gets the factored matrix from the
    1485:   preconditioner context.  This routine is valid only for the `PCLU`,
    1486:   `PCILU`, `PCCHOLESKY`, and `PCICC` methods.

    1488:   Not Collective though `mat` is parallel if `pc` is parallel

    1490:   Input Parameter:
    1491: . pc - the `PC` preconditioner context

    1493:   Output Parameters:
    1494: . mat - the factored matrix

    1496:   Level: advanced

    1498:   Note:
    1499:   Does not increase the reference count for `mat` so DO NOT destroy it

    1501: .seealso: [](ch_ksp), `PC`, `PCLU`, `PCILU`, `PCCHOLESKY`, `PCICC`
    1502: @*/
    1503: PetscErrorCode PCFactorGetMatrix(PC pc, Mat *mat)
    1504: {
    1505:   PetscFunctionBegin;
    1507:   PetscAssertPointer(mat, 2);
    1508:   PetscCall(PCFactorSetUpMatSolverType(pc));
    1509:   PetscUseTypeMethod(pc, getfactoredmatrix, mat);
    1510:   PetscFunctionReturn(PETSC_SUCCESS);
    1511: }

    1513: /*@
    1514:   PCSetOptionsPrefix - Sets the prefix used for searching for all
    1515:   `PC` options in the database.

    1517:   Logically Collective

    1519:   Input Parameters:
    1520: + pc     - the `PC` preconditioner context
    1521: - prefix - the prefix string to prepend to all `PC` option requests

    1523:   Level: advanced

    1525:   Note:
    1526:   A hyphen (-) must NOT be given at the beginning of the prefix name.
    1527:   The first character of all runtime options is AUTOMATICALLY the
    1528:   hyphen.

    1530: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions()`, `PCAppendOptionsPrefix()`, `PCGetOptionsPrefix()`
    1531: @*/
    1532: PetscErrorCode PCSetOptionsPrefix(PC pc, const char prefix[])
    1533: {
    1534:   PetscFunctionBegin;
    1536:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc, prefix));
    1537:   PetscFunctionReturn(PETSC_SUCCESS);
    1538: }

    1540: /*@
    1541:   PCAppendOptionsPrefix - Appends to the prefix used for searching for all
    1542:   `PC` options in the database.

    1544:   Logically Collective

    1546:   Input Parameters:
    1547: + pc     - the `PC` preconditioner context
    1548: - prefix - the prefix string to prepend to all `PC` option requests

    1550:   Level: advanced

    1552:   Note:
    1553:   A hyphen (-) must NOT be given at the beginning of the prefix name.
    1554:   The first character of all runtime options is AUTOMATICALLY the
    1555:   hyphen.

    1557: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions()`, `PCSetOptionsPrefix()`, `PCGetOptionsPrefix()`
    1558: @*/
    1559: PetscErrorCode PCAppendOptionsPrefix(PC pc, const char prefix[])
    1560: {
    1561:   PetscFunctionBegin;
    1563:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc, prefix));
    1564:   PetscFunctionReturn(PETSC_SUCCESS);
    1565: }

    1567: /*@
    1568:   PCGetOptionsPrefix - Gets the prefix used for searching for all
    1569:   `PC` options in the database.

    1571:   Not Collective

    1573:   Input Parameter:
    1574: . pc - the `PC` preconditioner context

    1576:   Output Parameter:
    1577: . prefix - pointer to the prefix string used, is returned

    1579:   Level: advanced

    1581: .seealso: [](ch_ksp), `PC`, `PCSetFromOptions()`, `PCSetOptionsPrefix()`, `PCAppendOptionsPrefix()`
    1582: @*/
    1583: PetscErrorCode PCGetOptionsPrefix(PC pc, const char *prefix[])
    1584: {
    1585:   PetscFunctionBegin;
    1587:   PetscAssertPointer(prefix, 2);
    1588:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, prefix));
    1589:   PetscFunctionReturn(PETSC_SUCCESS);
    1590: }

    1592: /*
    1593:    Indicates the right-hand side will be changed by KSPSolve(), this occurs for a few
    1594:   preconditioners including BDDC and Eisentat that transform the equations before applying
    1595:   the Krylov methods
    1596: */
    1597: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc, PetscBool *change)
    1598: {
    1599:   PetscFunctionBegin;
    1601:   PetscAssertPointer(change, 2);
    1602:   *change = PETSC_FALSE;
    1603:   PetscTryMethod(pc, "PCPreSolveChangeRHS_C", (PC, PetscBool *), (pc, change));
    1604:   PetscFunctionReturn(PETSC_SUCCESS);
    1605: }

    1607: /*@
    1608:   PCPreSolve - Optional pre-solve phase, intended for any preconditioner-specific actions that must be performed before
    1609:   the iterative solve itself. Used in conjunction with `PCPostSolve()`

    1611:   Collective

    1613:   Input Parameters:
    1614: + pc  - the `PC` preconditioner context
    1615: - ksp - the Krylov subspace context

    1617:   Level: developer

    1619:   Notes:
    1620:   `KSPSolve()` calls this directly, so is rarely called by the user.

    1622:   Certain preconditioners, such as the `PCType` of `PCEISENSTAT`, change the formulation of the linear system to be solved iteratively.
    1623:   This function performs that transformation. `PCPostSolve()` then transforms the system back to its original form after the solve.
    1624:   `PCPostSolve()` also transforms the resulting solution of the transformed system to the solution of the original problem.

    1626:   `KSPSetPostSolve()` provides an alternative way to provide such transformations.

    1628: .seealso: [](ch_ksp), `PC`, `PCPostSolve()`, `KSP`, `PCSetPostSetUp()`, `KSPSetPreSolve()`, `KSPSetPostSolve()`
    1629: @*/
    1630: PetscErrorCode PCPreSolve(PC pc, KSP ksp)
    1631: {
    1632:   Vec x, rhs;

    1634:   PetscFunctionBegin;
    1637:   pc->presolvedone++;
    1638:   PetscCheck(pc->presolvedone <= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot embed PCPreSolve() more than twice");
    1639:   PetscCall(KSPGetSolution(ksp, &x));
    1640:   PetscCall(KSPGetRhs(ksp, &rhs));
    1641:   PetscTryTypeMethod(pc, presolve, ksp, rhs, x);
    1642:   PetscFunctionReturn(PETSC_SUCCESS);
    1643: }

    1645: /*@C
    1646:   PCSetPostSetUp - Sets function called at the end of `PCSetUp()` to adjust the computed preconditioner

    1648:   Logically Collective

    1650:   Input Parameters:
    1651: + pc        - the preconditioner object
    1652: - postsetup - the function to call after `PCSetUp()`

    1654:   Calling sequence of `postsetup`:
    1655: . pc - the `PC` context

    1657:   Level: developer

    1659: .seealso: [](ch_ksp), `PC`, `PCSetUp()`
    1660: @*/
    1661: PetscErrorCode PCSetPostSetUp(PC pc, PetscErrorCode (*postsetup)(PC pc))
    1662: {
    1663:   PetscFunctionBegin;
    1665:   pc->postsetup = postsetup;
    1666:   PetscFunctionReturn(PETSC_SUCCESS);
    1667: }

    1669: /*@
    1670:   PCPostSolve - Optional post-solve phase, intended for any
    1671:   preconditioner-specific actions that must be performed after
    1672:   the iterative solve itself.

    1674:   Collective

    1676:   Input Parameters:
    1677: + pc  - the `PC` preconditioner context
    1678: - ksp - the `KSP` Krylov subspace context

    1680:   Example Usage:
    1681: .vb
    1682:     PCPreSolve(pc,ksp);
    1683:     KSPSolve(ksp,b,x);
    1684:     PCPostSolve(pc,ksp);
    1685: .ve

    1687:   Level: developer

    1689:   Note:
    1690:   `KSPSolve()` calls this routine directly, so it is rarely called by the user.

    1692: .seealso: [](ch_ksp), `PC`, `KSPSetPostSolve()`, `KSPSetPreSolve()`, `PCPreSolve()`, `KSPSolve()`
    1693: @*/
    1694: PetscErrorCode PCPostSolve(PC pc, KSP ksp)
    1695: {
    1696:   Vec x, rhs;

    1698:   PetscFunctionBegin;
    1701:   pc->presolvedone--;
    1702:   PetscCall(KSPGetSolution(ksp, &x));
    1703:   PetscCall(KSPGetRhs(ksp, &rhs));
    1704:   PetscTryTypeMethod(pc, postsolve, ksp, rhs, x);
    1705:   PetscFunctionReturn(PETSC_SUCCESS);
    1706: }

    1708: /*@
    1709:   PCLoad - Loads a `PC` that has been stored in binary  with `PCView()`.

    1711:   Collective

    1713:   Input Parameters:
    1714: + newdm  - the newly loaded `PC`, this needs to have been created with `PCCreate()` or
    1715:            some related function before a call to `PCLoad()`.
    1716: - viewer - binary file viewer `PETSCVIEWERBINARY`, obtained from `PetscViewerBinaryOpen()`

    1718:   Level: intermediate

    1720:   Note:
    1721:   The type is determined by the data in the file, any `PCType` set into the `PC` before this call is ignored.

    1723: .seealso: [](ch_ksp), `PC`, `PetscViewerBinaryOpen()`, `PCView()`, `MatLoad()`, `VecLoad()`, `PETSCVIEWERBINARY`
    1724: @*/
    1725: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
    1726: {
    1727:   PetscBool isbinary;
    1728:   PetscInt  classid;
    1729:   char      type[256];

    1731:   PetscFunctionBegin;
    1734:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
    1735:   PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");

    1737:   PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
    1738:   PetscCheck(classid == PC_FILE_CLASSID, PetscObjectComm((PetscObject)newdm), PETSC_ERR_ARG_WRONG, "Not PC next in file");
    1739:   PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
    1740:   PetscCall(PCSetType(newdm, type));
    1741:   PetscTryTypeMethod(newdm, load, viewer);
    1742:   PetscFunctionReturn(PETSC_SUCCESS);
    1743: }

    1745: #include <petscdraw.h>
    1746: #if defined(PETSC_HAVE_SAWS)
    1747: #include <petscviewersaws.h>
    1748: #endif

    1750: /*@
    1751:   PCViewFromOptions - View (print or provide information about) the `PC`, based on options in the options database

    1753:   Collective

    1755:   Input Parameters:
    1756: + A    - the `PC` context
    1757: . obj  - Optional object that provides the options prefix
    1758: - name - command line option name

    1760:   Options Database Key:
    1761: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments

    1763:   Level: developer

    1765: .seealso: [](ch_ksp), `PC`, `PCView`, `PetscObjectViewFromOptions()`, `PCCreate()`
    1766: @*/
    1767: PetscErrorCode PCViewFromOptions(PC A, PetscObject obj, const char name[])
    1768: {
    1769:   PetscFunctionBegin;
    1771:   PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
    1772:   PetscFunctionReturn(PETSC_SUCCESS);
    1773: }

    1775: /*@
    1776:   PCView - Prints information about the `PC`

    1778:   Collective

    1780:   Input Parameters:
    1781: + pc     - the `PC` preconditioner context
    1782: - viewer - optional `PetscViewer` visualization context

    1784:   Level: intermediate

    1786:   Notes:
    1787:   The available visualization contexts include
    1788: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
    1789: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
    1790:   output where only the first processor opens
    1791:   the file. All other processors send their
    1792:   data to the first processor to print.

    1794:   The user can open an alternative visualization contexts with
    1795:   `PetscViewerASCIIOpen()` (output to a specified file).

    1797: .seealso: [](ch_ksp), `PC`, `PetscViewer`, `PetscViewerType`, `KSPView()`, `PetscViewerASCIIOpen()`
    1798: @*/
    1799: PetscErrorCode PCView(PC pc, PetscViewer viewer)
    1800: {
    1801:   PCType            cstr;
    1802:   PetscViewerFormat format;
    1803:   PetscBool         isascii, isstring, isbinary, isdraw, pop = PETSC_FALSE;
    1804: #if defined(PETSC_HAVE_SAWS)
    1805:   PetscBool issaws;
    1806: #endif

    1808:   PetscFunctionBegin;
    1810:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));
    1812:   PetscCheckSameComm(pc, 1, viewer, 2);

    1814:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    1815:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
    1816:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
    1817:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
    1818: #if defined(PETSC_HAVE_SAWS)
    1819:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
    1820: #endif

    1822:   if (isascii) {
    1823:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)pc, viewer));
    1824:     if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, "  PC has not been set up so information may be incomplete\n"));
    1825:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1826:     PetscTryTypeMethod(pc, view, viewer);
    1827:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1828:     if (pc->mat) {
    1829:       PetscCall(PetscViewerGetFormat(viewer, &format));
    1830:       if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
    1831:         PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
    1832:         pop = PETSC_TRUE;
    1833:       }
    1834:       if (pc->pmat == pc->mat) {
    1835:         PetscCall(PetscViewerASCIIPrintf(viewer, "  linear system matrix, which is also used to construct the preconditioner:\n"));
    1836:         PetscCall(PetscViewerASCIIPushTab(viewer));
    1837:         PetscCall(MatView(pc->mat, viewer));
    1838:         PetscCall(PetscViewerASCIIPopTab(viewer));
    1839:       } else {
    1840:         if (pc->pmat) {
    1841:           PetscCall(PetscViewerASCIIPrintf(viewer, "  linear system matrix, followed by the matrix used to construct the preconditioner:\n"));
    1842:         } else {
    1843:           PetscCall(PetscViewerASCIIPrintf(viewer, "  linear system matrix:\n"));
    1844:         }
    1845:         PetscCall(PetscViewerASCIIPushTab(viewer));
    1846:         PetscCall(MatView(pc->mat, viewer));
    1847:         if (pc->pmat) PetscCall(MatView(pc->pmat, viewer));
    1848:         PetscCall(PetscViewerASCIIPopTab(viewer));
    1849:       }
    1850:       if (pop) PetscCall(PetscViewerPopFormat(viewer));
    1851:     }
    1852:   } else if (isstring) {
    1853:     PetscCall(PCGetType(pc, &cstr));
    1854:     PetscCall(PetscViewerStringSPrintf(viewer, " PCType: %-7.7s", cstr));
    1855:     PetscTryTypeMethod(pc, view, viewer);
    1856:     if (pc->mat) PetscCall(MatView(pc->mat, viewer));
    1857:     if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
    1858:   } else if (isbinary) {
    1859:     PetscInt    classid = PC_FILE_CLASSID;
    1860:     MPI_Comm    comm;
    1861:     PetscMPIInt rank;
    1862:     char        type[256];

    1864:     PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
    1865:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1866:     if (rank == 0) {
    1867:       PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
    1868:       PetscCall(PetscStrncpy(type, ((PetscObject)pc)->type_name, 256));
    1869:       PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
    1870:     }
    1871:     PetscTryTypeMethod(pc, view, viewer);
    1872:   } else if (isdraw) {
    1873:     PetscDraw draw;
    1874:     char      str[25];
    1875:     PetscReal x, y, bottom, h;
    1876:     PetscInt  n;

    1878:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    1879:     PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
    1880:     if (pc->mat) {
    1881:       PetscCall(MatGetSize(pc->mat, &n, NULL));
    1882:       PetscCall(PetscSNPrintf(str, 25, "PC: %s (%" PetscInt_FMT ")", ((PetscObject)pc)->type_name, n));
    1883:     } else {
    1884:       PetscCall(PetscSNPrintf(str, 25, "PC: %s", ((PetscObject)pc)->type_name));
    1885:     }
    1886:     PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
    1887:     bottom = y - h;
    1888:     PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
    1889:     PetscTryTypeMethod(pc, view, viewer);
    1890:     PetscCall(PetscDrawPopCurrentPoint(draw));
    1891: #if defined(PETSC_HAVE_SAWS)
    1892:   } else if (issaws) {
    1893:     PetscMPIInt rank;

    1895:     PetscCall(PetscObjectName((PetscObject)pc));
    1896:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
    1897:     if (!((PetscObject)pc)->amsmem && rank == 0) PetscCall(PetscObjectViewSAWs((PetscObject)pc, viewer));
    1898:     if (pc->mat) PetscCall(MatView(pc->mat, viewer));
    1899:     if (pc->pmat && pc->pmat != pc->mat) PetscCall(MatView(pc->pmat, viewer));
    1900: #endif
    1901:   }
    1902:   PetscFunctionReturn(PETSC_SUCCESS);
    1903: }

    1905: /*@C
    1906:   PCRegister -  Adds a method (`PCType`) to the PETSc preconditioner package.

    1908:   Not collective. No Fortran Support

    1910:   Input Parameters:
    1911: + sname    - name of a new user-defined solver
    1912: - function - routine to create the method context which will be stored in a `PC` when `PCSetType()` is called

    1914:   Example Usage:
    1915: .vb
    1916:    PCRegister("my_solver", MySolverCreate);
    1917: .ve

    1919:   Then, your solver can be chosen with the procedural interface via
    1920: .vb
    1921:   PCSetType(pc, "my_solver")
    1922: .ve
    1923:   or at runtime via the option
    1924: .vb
    1925:   -pc_type my_solver
    1926: .ve

    1928:   Level: advanced

    1930:   Note:
    1931:   A simpler alternative to using `PCRegister()` for an application specific preconditioner is to use a `PC` of `PCType` `PCSHELL` and
    1932:   provide your customizations with `PCShellSetContext()` and `PCShellSetApply()`

    1934:   `PCRegister()` may be called multiple times to add several user-defined preconditioners.

    1936: .seealso: [](ch_ksp), `PC`, `PCType`, `PCRegisterAll()`, `PCSetType()`, `PCShellSetContext()`, `PCShellSetApply()`, `PCSHELL`
    1937: @*/
    1938: PetscErrorCode PCRegister(const char sname[], PetscErrorCode (*function)(PC))
    1939: {
    1940:   PetscFunctionBegin;
    1941:   PetscCall(PCInitializePackage());
    1942:   PetscCall(PetscFunctionListAdd(&PCList, sname, function));
    1943:   PetscFunctionReturn(PETSC_SUCCESS);
    1944: }

    1946: static PetscErrorCode MatMult_PC(Mat A, Vec X, Vec Y)
    1947: {
    1948:   PC pc;

    1950:   PetscFunctionBegin;
    1951:   PetscCall(MatShellGetContext(A, &pc));
    1952:   PetscCall(PCApply(pc, X, Y));
    1953:   PetscFunctionReturn(PETSC_SUCCESS);
    1954: }

    1956: /*@
    1957:   PCComputeOperator - Computes the explicit preconditioned operator as a matrix `Mat`.

    1959:   Collective

    1961:   Input Parameters:
    1962: + pc      - the `PC` preconditioner object
    1963: - mattype - the `MatType` to be used for the operator

    1965:   Output Parameter:
    1966: . mat - the explicit preconditioned operator

    1968:   Level: advanced

    1970:   Note:
    1971:   This computation is done by applying the operators to columns of the identity matrix.
    1972:   This routine is costly in general, and is recommended for use only with relatively small systems.
    1973:   Currently, this routine uses a dense matrix format when `mattype` == `NULL`

    1975:   Developer Note:
    1976:   This should be called `PCCreateExplicitOperator()`

    1978: .seealso: [](ch_ksp), `PC`, `KSPComputeOperator()`, `MatType`
    1979: @*/
    1980: PetscErrorCode PCComputeOperator(PC pc, MatType mattype, Mat *mat)
    1981: {
    1982:   PetscInt N, M, m, n;
    1983:   Mat      A, Apc;

    1985:   PetscFunctionBegin;
    1987:   PetscAssertPointer(mat, 3);
    1988:   PetscCall(PCGetOperators(pc, &A, NULL));
    1989:   PetscCall(MatGetLocalSize(A, &m, &n));
    1990:   PetscCall(MatGetSize(A, &M, &N));
    1991:   PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pc), m, n, M, N, pc, &Apc));
    1992:   PetscCall(MatShellSetOperation(Apc, MATOP_MULT, (PetscErrorCodeFn *)MatMult_PC));
    1993:   PetscCall(MatComputeOperator(Apc, mattype, mat));
    1994:   PetscCall(MatDestroy(&Apc));
    1995:   PetscFunctionReturn(PETSC_SUCCESS);
    1996: }

    1998: /*@
    1999:   PCSetCoordinates - sets the coordinates of all the nodes (degrees of freedom in the vector) on the local process

    2001:   Collective

    2003:   Input Parameters:
    2004: + pc     - the `PC` preconditioner context
    2005: . dim    - the dimension of the coordinates 1, 2, or 3
    2006: . nloc   - the blocked size of the coordinates array
    2007: - coords - the coordinates array

    2009:   Level: intermediate

    2011:   Notes:
    2012:   `coords` is an array of the dim coordinates for the nodes on
    2013:   the local processor, of size `dim`*`nloc`.
    2014:   If there are 108 equations (dofs) on a processor
    2015:   for a 3d displacement finite element discretization of elasticity (so
    2016:   that there are nloc = 36 = 108/3 nodes) then the array must have 108
    2017:   double precision values (ie, 3 * 36).  These x y z coordinates
    2018:   should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
    2019:   ... , N-1.z ].

    2021:   The information provided here can be used by some preconditioners, such as `PCGAMG`, to produce a better preconditioner.
    2022:   See also  `MatSetNearNullSpace()`.

    2024: .seealso: [](ch_ksp), `PC`, `MatSetNearNullSpace()`
    2025: @*/
    2026: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
    2027: {
    2028:   PetscFunctionBegin;
    2031:   PetscTryMethod(pc, "PCSetCoordinates_C", (PC, PetscInt, PetscInt, PetscReal[]), (pc, dim, nloc, coords));
    2032:   PetscFunctionReturn(PETSC_SUCCESS);
    2033: }

    2035: /*@
    2036:   PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)

    2038:   Logically Collective

    2040:   Input Parameter:
    2041: . pc - the precondition context

    2043:   Output Parameters:
    2044: + num_levels     - the number of levels
    2045: - interpolations - the interpolation matrices (size of `num_levels`-1)

    2047:   Level: advanced

    2049:   Developer Note:
    2050:   Why is this here instead of in `PCMG` etc?

    2052: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetInterpolation()`, `PCGetCoarseOperators()`
    2053: @*/
    2054: PetscErrorCode PCGetInterpolations(PC pc, PetscInt *num_levels, Mat *interpolations[])
    2055: {
    2056:   PetscFunctionBegin;
    2058:   PetscAssertPointer(num_levels, 2);
    2059:   PetscAssertPointer(interpolations, 3);
    2060:   PetscUseMethod(pc, "PCGetInterpolations_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, interpolations));
    2061:   PetscFunctionReturn(PETSC_SUCCESS);
    2062: }

    2064: /*@
    2065:   PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)

    2067:   Logically Collective

    2069:   Input Parameter:
    2070: . pc - the precondition context

    2072:   Output Parameters:
    2073: + num_levels      - the number of levels
    2074: - coarseOperators - the coarse operator matrices (size of `num_levels`-1)

    2076:   Level: advanced

    2078:   Developer Note:
    2079:   Why is this here instead of in `PCMG` etc?

    2081: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
    2082: @*/
    2083: PetscErrorCode PCGetCoarseOperators(PC pc, PetscInt *num_levels, Mat *coarseOperators[])
    2084: {
    2085:   PetscFunctionBegin;
    2087:   PetscAssertPointer(num_levels, 2);
    2088:   PetscAssertPointer(coarseOperators, 3);
    2089:   PetscUseMethod(pc, "PCGetCoarseOperators_C", (PC, PetscInt *, Mat *[]), (pc, num_levels, coarseOperators));
    2090:   PetscFunctionReturn(PETSC_SUCCESS);
    2091: }