Actual source code: lmvmpc.c

      1: /*
      2:    This provides a thin wrapper around LMVM matrices in order to use their MatLMVMSolve
      3:    methods as preconditioner applications in KSP solves.
      4: */

      6: #include <petsc/private/pcimpl.h>
      7: #include <petsc/private/matimpl.h>

      9: typedef struct {
     10:   Vec              xwork, ywork;
     11:   IS               inactive;
     12:   Mat              B;
     13:   Vec              X;
     14:   PetscObjectState Xstate;
     15:   PetscBool        setfromoptionscalled;
     16: } PC_LMVM;

     18: /*@
     19:   PCLMVMSetUpdateVec - Set the vector to be used as solution update for the internal LMVM matrix.

     21:   Input Parameters:
     22: + pc - The preconditioner
     23: - X  - Solution vector

     25:   Level: intermediate

     27:   Notes:
     28:   This is only needed if you want the preconditioner to automatically update the internal matrix.
     29:   It is called in some `SNES` implementations to update the preconditioner.
     30:   The right-hand side of the linear system is used as function vector.

     32: .seealso: `MatLMVMUpdate()`, `PCLMVMSetMatLMVM()`
     33: @*/
     34: PetscErrorCode PCLMVMSetUpdateVec(PC pc, Vec X)
     35: {
     36:   PC_LMVM  *ctx;
     37:   PetscBool same;

     39:   PetscFunctionBegin;
     42:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
     43:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
     44:   ctx = (PC_LMVM *)pc->data;
     45:   PetscCall(PetscObjectReference((PetscObject)X));
     46:   PetscCall(VecDestroy(&ctx->X));
     47:   ctx->X      = X;
     48:   ctx->Xstate = -1;
     49:   PetscFunctionReturn(PETSC_SUCCESS);
     50: }

     52: /*@
     53:   PCLMVMSetMatLMVM - Replaces the `MATLMVM` matrix inside the preconditioner with the one provided by the user.

     55:   Input Parameters:
     56: + pc - An `PCLMVM` preconditioner
     57: - B  - An `MATLMVM` type matrix

     59:   Level: intermediate

     61: .seealso: [](ch_ksp), `PCLMVMGetMatLMVM()`
     62: @*/
     63: PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B)
     64: {
     65:   PC_LMVM  *ctx;
     66:   PetscBool same;

     68:   PetscFunctionBegin;
     71:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
     72:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
     73:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
     74:   PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an MATLMVM.");
     75:   ctx = (PC_LMVM *)pc->data;
     76:   PetscCall(PetscObjectReference((PetscObject)B));
     77:   PetscCall(MatDestroy(&ctx->B));
     78:   ctx->B = B;
     79:   PetscFunctionReturn(PETSC_SUCCESS);
     80: }

     82: /*@
     83:   PCLMVMGetMatLMVM - Returns a pointer to the underlying `MATLMVM` matrix.

     85:   Input Parameter:
     86: . pc - An `PCLMVM` preconditioner

     88:   Output Parameter:
     89: . B - `MATLMVM` matrix used by the preconditioner

     91:   Level: intermediate

     93: .seealso: [](ch_ksp), `PCLMVMSetMatLMVM()`
     94: @*/
     95: PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B)
     96: {
     97:   PC_LMVM  *ctx;
     98:   PetscBool same;

    100:   PetscFunctionBegin;
    102:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
    103:   PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
    104:   ctx = (PC_LMVM *)pc->data;
    105:   if (!ctx->B) {
    106:     Mat J;

    108:     if (pc->useAmat) J = pc->mat;
    109:     else J = pc->pmat;
    110:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)J, MATLMVM, &same));
    111:     if (same) *B = J;
    112:     else {
    113:       const char *prefix;

    115:       PetscCall(PCGetOptionsPrefix(pc, &prefix));
    116:       PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B));
    117:       PetscCall(MatSetOptionsPrefix(ctx->B, prefix));
    118:       PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_"));
    119:       PetscCall(MatSetType(ctx->B, MATLMVMBFGS));
    120:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1));
    121:       *B = ctx->B;
    122:     }
    123:   } else *B = ctx->B;
    124:   PetscFunctionReturn(PETSC_SUCCESS);
    125: }

    127: /*@
    128:   PCLMVMSetIS - Sets the index sets that reduce the `PC` application.

    130:   Input Parameters:
    131: + pc       - An `PCLMVM` preconditioner
    132: - inactive - Index set defining the variables removed from the problem

    134:   Level: intermediate

    136:   Developer Notes:
    137:   Need to explain the purpose of this `IS`

    139: .seealso: [](ch_ksp), `PCLMVMClearIS()`
    140: @*/
    141: PetscErrorCode PCLMVMSetIS(PC pc, IS inactive)
    142: {
    143:   PC_LMVM  *ctx;
    144:   PetscBool same;

    146:   PetscFunctionBegin;
    149:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
    150:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    151:   ctx = (PC_LMVM *)pc->data;
    152:   PetscCall(PCLMVMClearIS(pc));
    153:   PetscCall(PetscObjectReference((PetscObject)inactive));
    154:   ctx->inactive = inactive;
    155:   PetscFunctionReturn(PETSC_SUCCESS);
    156: }

    158: /*@
    159:   PCLMVMClearIS - Removes the inactive variable index set from a `PCLMVM`

    161:   Input Parameter:
    162: . pc - An `PCLMVM` preconditioner

    164:   Level: intermediate

    166: .seealso: [](ch_ksp), `PCLMVMSetIS()`
    167: @*/
    168: PetscErrorCode PCLMVMClearIS(PC pc)
    169: {
    170:   PC_LMVM  *ctx;
    171:   PetscBool same;

    173:   PetscFunctionBegin;
    175:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
    176:   if (!same) PetscFunctionReturn(PETSC_SUCCESS);
    177:   ctx = (PC_LMVM *)pc->data;
    178:   PetscCall(ISDestroy(&ctx->inactive));
    179:   PetscFunctionReturn(PETSC_SUCCESS);
    180: }

    182: static PetscErrorCode PCApply_LMVM(PC pc, Vec x, Vec y)
    183: {
    184:   PC_LMVM *ctx = (PC_LMVM *)pc->data;
    185:   Vec      xsub, ysub, Bx = x, By = y;
    186:   Mat      B = ctx->B ? ctx->B : (pc->useAmat ? pc->mat : pc->pmat);

    188:   PetscFunctionBegin;
    189:   if (ctx->inactive) {
    190:     if (!ctx->xwork) PetscCall(MatCreateVecs(B, &ctx->xwork, &ctx->ywork));
    191:     PetscCall(VecZeroEntries(ctx->xwork));
    192:     PetscCall(VecGetSubVector(ctx->xwork, ctx->inactive, &xsub));
    193:     PetscCall(VecCopy(x, xsub));
    194:     PetscCall(VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub));
    195:     Bx = ctx->xwork;
    196:     By = ctx->ywork;
    197:   }
    198:   PetscCall(MatSolve(B, Bx, By));
    199:   if (ctx->inactive) {
    200:     PetscCall(VecGetSubVector(ctx->ywork, ctx->inactive, &ysub));
    201:     PetscCall(VecCopy(ysub, y));
    202:     PetscCall(VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub));
    203:   }
    204:   PetscFunctionReturn(PETSC_SUCCESS);
    205: }

    207: static PetscErrorCode PCReset_LMVM(PC pc)
    208: {
    209:   PC_LMVM *ctx = (PC_LMVM *)pc->data;

    211:   PetscFunctionBegin;
    212:   PetscCall(VecDestroy(&ctx->xwork));
    213:   PetscCall(VecDestroy(&ctx->ywork));
    214:   PetscFunctionReturn(PETSC_SUCCESS);
    215: }

    217: static PetscErrorCode PCSetUp_LMVM(PC pc)
    218: {
    219:   PetscInt  n, N;
    220:   PetscBool allocated;
    221:   Mat       B;
    222:   PC_LMVM  *ctx = (PC_LMVM *)pc->data;

    224:   PetscFunctionBegin;
    225:   PetscCall(PCLMVMGetMatLMVM(pc, &B));
    226:   PetscCall(MatLMVMIsAllocated(B, &allocated));
    227:   if (!allocated) {
    228:     Vec t1, t2;

    230:     PetscCall(MatCreateVecs(pc->mat, &t1, &t2));
    231:     PetscCall(VecGetLocalSize(t1, &n));
    232:     PetscCall(VecGetSize(t1, &N));
    233:     PetscCall(MatSetSizes(B, n, n, N, N));
    234:     PetscCall(MatLMVMAllocate(B, t1, t2));
    235:     PetscCall(VecDestroy(&t1));
    236:     PetscCall(VecDestroy(&t2));
    237:   }
    238:   /* Only call SetFromOptions if we internally handle the LMVM matrix */
    239:   if (B == ctx->B && ctx->setfromoptionscalled) PetscCall(MatSetFromOptions(ctx->B));
    240:   ctx->setfromoptionscalled = PETSC_FALSE;
    241:   PetscFunctionReturn(PETSC_SUCCESS);
    242: }

    244: static PetscErrorCode PCView_LMVM(PC pc, PetscViewer viewer)
    245: {
    246:   PC_LMVM  *ctx = (PC_LMVM *)pc->data;
    247:   PetscBool isascii;

    249:   PetscFunctionBegin;
    250:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    251:   if (isascii && ctx->B && ctx->B->assembled) {
    252:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
    253:     PetscCall(MatView(ctx->B, viewer));
    254:     PetscCall(PetscViewerPopFormat(viewer));
    255:   }
    256:   PetscFunctionReturn(PETSC_SUCCESS);
    257: }

    259: static PetscErrorCode PCSetFromOptions_LMVM(PC pc, PetscOptionItems PetscOptionsObject)
    260: {
    261:   PC_LMVM *ctx = (PC_LMVM *)pc->data;

    263:   PetscFunctionBegin;
    264:   /* defer SetFromOptions calls to PCSetUp_LMVM */
    265:   ctx->setfromoptionscalled = PETSC_TRUE;
    266:   PetscFunctionReturn(PETSC_SUCCESS);
    267: }

    269: static PetscErrorCode PCPreSolve_LMVM(PC pc, KSP ksp, Vec F, Vec X)
    270: {
    271:   PC_LMVM *ctx = (PC_LMVM *)pc->data;

    273:   PetscFunctionBegin;
    274:   if (ctx->X && ctx->B) { /* Perform update only if requested. Otherwise we assume the user, e.g. TAO, has already taken care of it */
    275:     PetscObjectState Xstate;

    277:     PetscCall(PetscObjectStateGet((PetscObject)ctx->X, &Xstate));
    278:     if (ctx->Xstate != Xstate) PetscCall(MatLMVMUpdate(ctx->B, ctx->X, F));
    279:     ctx->Xstate = Xstate;
    280:   }
    281:   PetscFunctionReturn(PETSC_SUCCESS);
    282: }

    284: static PetscErrorCode PCDestroy_LMVM(PC pc)
    285: {
    286:   PC_LMVM *ctx = (PC_LMVM *)pc->data;

    288:   PetscFunctionBegin;
    289:   PetscCall(ISDestroy(&ctx->inactive));
    290:   PetscCall(VecDestroy(&ctx->xwork));
    291:   PetscCall(VecDestroy(&ctx->ywork));
    292:   PetscCall(VecDestroy(&ctx->X));
    293:   PetscCall(MatDestroy(&ctx->B));
    294:   PetscCall(PetscFree(pc->data));
    295:   PetscFunctionReturn(PETSC_SUCCESS);
    296: }

    298: /*MC
    299:    PCLMVM - A preconditioner constructed from a `MATLMVM` matrix.
    300:             If the matrix used to construct the preconditioner is not of type `MATLMVM`, an internal matrix is used.
    301:             Options for the internal `MATLMVM` matrix can be accessed with the `-pc_lmvm_` prefix.
    302:             Alternatively, the user can pass a suitable matrix with `PCLMVMSetMatLMVM()`.

    304:    Level: intermediate

    306: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCLMVMSetUpdateVec()`, `PCLMVMSetMatLMVM()`, `PCLMVMGetMatLMVM()`
    307: M*/
    308: PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
    309: {
    310:   PC_LMVM *ctx;

    312:   PetscFunctionBegin;
    313:   PetscCall(PetscNew(&ctx));
    314:   pc->data = (void *)ctx;

    316:   pc->ops->reset               = PCReset_LMVM;
    317:   pc->ops->setup               = PCSetUp_LMVM;
    318:   pc->ops->destroy             = PCDestroy_LMVM;
    319:   pc->ops->view                = PCView_LMVM;
    320:   pc->ops->apply               = PCApply_LMVM;
    321:   pc->ops->setfromoptions      = PCSetFromOptions_LMVM;
    322:   pc->ops->applysymmetricleft  = NULL;
    323:   pc->ops->applysymmetricright = NULL;
    324:   pc->ops->applytranspose      = NULL;
    325:   pc->ops->applyrichardson     = NULL;
    326:   pc->ops->presolve            = PCPreSolve_LMVM;
    327:   pc->ops->postsolve           = NULL;
    328:   PetscFunctionReturn(PETSC_SUCCESS);
    329: }