Actual source code: pcmat.c

      1: #include <petsc/private/pcimpl.h>

      3: typedef enum {
      4:   PCMATOP_UNSPECIFIED,
      5:   PCMATOP_MULT,
      6:   PCMATOP_MULT_TRANSPOSE,
      7:   PCMATOP_MULT_HERMITIAN_TRANSPOSE,
      8:   PCMATOP_SOLVE,
      9:   PCMATOP_SOLVE_TRANSPOSE,
     10: } PCMatOperation;

     12: const char *const PCMatOpTypes[] = {"Unspecified", "Mult", "MultTranspose", "MultHermitianTranspose", "Solve", "SolveTranspose", NULL};

     14: typedef struct _PCMAT {
     15:   PCMatOperation apply;
     16: } PC_Mat;

     18: static PetscErrorCode PCApply_Mat(PC pc, Vec x, Vec y)
     19: {
     20:   PC_Mat *pcmat = (PC_Mat *)pc->data;

     22:   PetscFunctionBegin;
     23:   switch (pcmat->apply) {
     24:   case PCMATOP_MULT:
     25:     PetscCall(MatMult(pc->pmat, x, y));
     26:     break;
     27:   case PCMATOP_MULT_TRANSPOSE:
     28:     PetscCall(MatMultTranspose(pc->pmat, x, y));
     29:     break;
     30:   case PCMATOP_SOLVE:
     31:     PetscCall(MatSolve(pc->pmat, x, y));
     32:     break;
     33:   case PCMATOP_SOLVE_TRANSPOSE:
     34:     PetscCall(MatSolveTranspose(pc->pmat, x, y));
     35:     break;
     36:   case PCMATOP_MULT_HERMITIAN_TRANSPOSE:
     37:     PetscCall(MatMultHermitianTranspose(pc->pmat, x, y));
     38:     break;
     39:   default:
     40:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported %s case", PCMatOpTypes[pcmat->apply]);
     41:   }
     42:   PetscFunctionReturn(PETSC_SUCCESS);
     43: }

     45: static PetscErrorCode PCSetUp_Mat(PC pc)
     46: {
     47:   PC_Mat *pcmat = (PC_Mat *)pc->data;

     49:   PetscFunctionBegin;
     50:   if (pcmat->apply == PCMATOP_UNSPECIFIED) {
     51:     PetscBool hassolve;
     52:     PetscCall(MatHasOperation(pc->pmat, MATOP_SOLVE, &hassolve));
     53:     if (hassolve) pcmat->apply = PCMATOP_SOLVE;
     54:     else pcmat->apply = PCMATOP_MULT;
     55:   }
     56:   PetscFunctionReturn(PETSC_SUCCESS);
     57: }

     59: static PetscErrorCode PCMatApply_Mat(PC pc, Mat X, Mat Y)
     60: {
     61:   PC_Mat *pcmat = (PC_Mat *)pc->data;
     62:   Mat     W;

     64:   PetscFunctionBegin;
     65:   switch (pcmat->apply) {
     66:   case PCMATOP_MULT:
     67:     PetscCall(MatMatMult(pc->pmat, X, MAT_REUSE_MATRIX, PETSC_CURRENT, &Y));
     68:     break;
     69:   case PCMATOP_MULT_TRANSPOSE:
     70:     PetscCall(MatTransposeMatMult(pc->pmat, X, MAT_REUSE_MATRIX, PETSC_CURRENT, &Y));
     71:     break;
     72:   case PCMATOP_SOLVE:
     73:     PetscCall(MatMatSolve(pc->pmat, X, Y));
     74:     break;
     75:   case PCMATOP_SOLVE_TRANSPOSE:
     76:     PetscCall(MatMatSolveTranspose(pc->pmat, X, Y));
     77:     break;
     78:   case PCMATOP_MULT_HERMITIAN_TRANSPOSE:
     79:     PetscCall(MatDuplicate(X, MAT_COPY_VALUES, &W));
     80:     PetscCall(MatConjugate(W));
     81:     PetscCall(MatTransposeMatMult(pc->pmat, W, MAT_REUSE_MATRIX, PETSC_CURRENT, &Y));
     82:     PetscCall(MatConjugate(Y));
     83:     PetscCall(MatDestroy(&W));
     84:     break;
     85:   default:
     86:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported %s case", PCMatOpTypes[pcmat->apply]);
     87:   }
     88:   PetscFunctionReturn(PETSC_SUCCESS);
     89: }

     91: static PetscErrorCode PCApplyTranspose_Mat(PC pc, Vec x, Vec y)
     92: {
     93:   PC_Mat *pcmat = (PC_Mat *)pc->data;
     94:   Vec     w;

     96:   PetscFunctionBegin;
     97:   switch (pcmat->apply) {
     98:   case PCMATOP_MULT:
     99:     PetscCall(MatMultTranspose(pc->pmat, x, y));
    100:     break;
    101:   case PCMATOP_MULT_TRANSPOSE:
    102:     PetscCall(MatMult(pc->pmat, x, y));
    103:     break;
    104:   case PCMATOP_SOLVE:
    105:     PetscCall(MatSolveTranspose(pc->pmat, x, y));
    106:     break;
    107:   case PCMATOP_SOLVE_TRANSPOSE:
    108:     PetscCall(MatSolve(pc->pmat, x, y));
    109:     break;
    110:   case PCMATOP_MULT_HERMITIAN_TRANSPOSE:
    111:     PetscCall(VecDuplicate(x, &w));
    112:     PetscCall(VecCopy(x, w));
    113:     PetscCall(VecConjugate(w));
    114:     PetscCall(MatMult(pc->pmat, w, y));
    115:     PetscCall(VecConjugate(y));
    116:     PetscCall(VecDestroy(&w));
    117:     break;
    118:   default:
    119:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported %s case", PCMatOpTypes[pcmat->apply]);
    120:   }
    121:   PetscFunctionReturn(PETSC_SUCCESS);
    122: }

    124: static PetscErrorCode PCDestroy_Mat(PC pc)
    125: {
    126:   PetscFunctionBegin;
    127:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatSetApplyOperation_C", NULL));
    128:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatGetApplyOperation_C", NULL));
    129:   PetscCall(PetscFree(pc->data));
    130:   PetscFunctionReturn(PETSC_SUCCESS);
    131: }

    133: /*@
    134:   PCMatSetApplyOperation - Set which matrix operation of the preconditioning matrix implements `PCApply()` for `PCMAT`.

    136:   Logically collective

    138:   Input Parameters:
    139: + pc    - An instance of `PCMAT`
    140: - matop - The selected `MatOperation`

    142:   Level: intermediate

    144:   Note:
    145:   If you have a matrix type that implements an exact inverse that isn't a factorization,
    146:   you can use `PCMatSetApplyOperation(pc, MATOP_SOLVE)`.

    148: .seealso: [](ch_ksp), `PCMAT`, `PCMatGetApplyOperation()`, `PCApply()`, `MatOperation`
    149: @*/
    150: PetscErrorCode PCMatSetApplyOperation(PC pc, MatOperation matop)
    151: {
    152:   PetscFunctionBegin;
    154:   PetscTryMethod((PetscObject)pc, "PCMatSetApplyOperation_C", (PC, MatOperation), (pc, matop));
    155:   PetscFunctionReturn(PETSC_SUCCESS);
    156: }

    158: /*@
    159:   PCMatGetApplyOperation - Get which matrix operation of the preconditioning matrix implements `PCApply()` for `PCMAT`.

    161:   Logically collective

    163:   Input Parameter:
    164: . pc - An instance of `PCMAT`

    166:   Output Parameter:
    167: . matop - The `MatOperation`

    169:   Level: intermediate

    171: .seealso: [](ch_ksp), `PCMAT`, `PCMatSetApplyOperation()`, `PCApply()`, `MatOperation`
    172: @*/
    173: PetscErrorCode PCMatGetApplyOperation(PC pc, MatOperation *matop)
    174: {
    175:   PetscFunctionBegin;
    177:   PetscAssertPointer(matop, 2);
    178:   PetscUseMethod((PetscObject)pc, "PCMatGetApplyOperation_C", (PC, MatOperation *), (pc, matop));
    179:   PetscFunctionReturn(PETSC_SUCCESS);
    180: }

    182: static PetscErrorCode PCMatSetApplyOperation_Mat(PC pc, MatOperation matop)
    183: {
    184:   PC_Mat        *pcmat = (PC_Mat *)pc->data;
    185:   PCMatOperation pcmatop;

    187:   PetscFunctionBegin;
    188:   // clang-format off
    189: #define MATOP_TO_PCMATOP_CASE(var, OP) case MATOP_##OP: (var) = PCMATOP_##OP; break
    190:   switch (matop) {
    191:   MATOP_TO_PCMATOP_CASE(pcmatop, MULT);
    192:   MATOP_TO_PCMATOP_CASE(pcmatop, MULT_TRANSPOSE);
    193:   MATOP_TO_PCMATOP_CASE(pcmatop, MULT_HERMITIAN_TRANSPOSE);
    194:   MATOP_TO_PCMATOP_CASE(pcmatop, SOLVE);
    195:   MATOP_TO_PCMATOP_CASE(pcmatop, SOLVE_TRANSPOSE);
    196:   default:
    197:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported MatOperation %d for PCMatSetApplyOperation()", (int)matop);
    198:   }
    199: #undef MATOP_TO_PCMATOP_CASE
    200:   // clang-format on

    202:   pcmat->apply = pcmatop;
    203:   PetscFunctionReturn(PETSC_SUCCESS);
    204: }

    206: static PetscErrorCode PCMatGetApplyOperation_Mat(PC pc, MatOperation *matop_p)
    207: {
    208:   PC_Mat      *pcmat = (PC_Mat *)pc->data;
    209:   MatOperation matop = MATOP_MULT;

    211:   PetscFunctionBegin;
    212:   if (!pc->setupcalled) PetscCall(PCSetUp(pc));

    214:   // clang-format off
    215: #define PCMATOP_TO_MATOP_CASE(var, OP) case PCMATOP_##OP: (var) = MATOP_##OP; break
    216:   switch (pcmat->apply) {
    217:   PCMATOP_TO_MATOP_CASE(matop, MULT);
    218:   PCMATOP_TO_MATOP_CASE(matop, MULT_TRANSPOSE);
    219:   PCMATOP_TO_MATOP_CASE(matop, MULT_HERMITIAN_TRANSPOSE);
    220:   PCMATOP_TO_MATOP_CASE(matop, SOLVE);
    221:   PCMATOP_TO_MATOP_CASE(matop, SOLVE_TRANSPOSE);
    222:   default:
    223:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported %s case", PCMatOpTypes[pcmat->apply]);
    224:   }
    225: #undef PCMATOP_TO_MATOP_CASE
    226:   // clang-format on

    228:   *matop_p = matop;
    229:   PetscFunctionReturn(PETSC_SUCCESS);
    230: }

    232: static PetscErrorCode PCView_Mat(PC pc, PetscViewer viewer)
    233: {
    234:   PC_Mat   *pcmat = (PC_Mat *)pc->data;
    235:   PetscBool iascii;

    237:   PetscFunctionBegin;
    238:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    239:   if (iascii) { PetscCall(PetscViewerASCIIPrintf(viewer, "PCApply() == Mat%s()\n", PCMatOpTypes[pcmat->apply])); }
    240:   PetscFunctionReturn(PETSC_SUCCESS);
    241: }

    243: /*MC
    244:      PCMAT - A preconditioner obtained by applying an operation of the `pmat` provided in
    245:              in `PCSetOperators(pc,amat,pmat)` or `KSPSetOperators(ksp,amat,pmat)`.  By default, the operation is `MATOP_MULT`,
    246:              meaning that the `pmat` provides an approximate inverse of `amat`.
    247:              If some other operation of `pmat` implements the approximate inverse,
    248:              use `PCMatSetApplyOperation()` to select that operation.

    250:    Level: intermediate

    252: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSHELL`, `MatOperation`, `PCMatSetApplyOperation()`, `PCMatGetApplyOperation()`
    253: M*/

    255: PETSC_EXTERN PetscErrorCode PCCreate_Mat(PC pc)
    256: {
    257:   PC_Mat *data;

    259:   PetscFunctionBegin;
    260:   PetscCall(PetscNew(&data));
    261:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatSetApplyOperation_C", PCMatSetApplyOperation_Mat));
    262:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatGetApplyOperation_C", PCMatGetApplyOperation_Mat));
    263:   pc->data                     = data;
    264:   pc->ops->apply               = PCApply_Mat;
    265:   pc->ops->matapply            = PCMatApply_Mat;
    266:   pc->ops->applytranspose      = PCApplyTranspose_Mat;
    267:   pc->ops->setup               = PCSetUp_Mat;
    268:   pc->ops->destroy             = PCDestroy_Mat;
    269:   pc->ops->setfromoptions      = NULL;
    270:   pc->ops->view                = PCView_Mat;
    271:   pc->ops->applyrichardson     = NULL;
    272:   pc->ops->applysymmetricleft  = NULL;
    273:   pc->ops->applysymmetricright = NULL;
    274:   PetscFunctionReturn(PETSC_SUCCESS);
    275: }