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:   PetscBool allocated;
 14: } PC_LMVM;

 16: /*@
 17:   PCLMVMSetMatLMVM - Replaces the `MATLMVM` matrix inside the preconditioner with
 18:   the one provided by the user.

 20:   Input Parameters:
 21: + pc - An `PCLMVM` preconditioner
 22: - B  - An  LMVM-type matrix (`MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`)

 24:   Level: intermediate

 26: .seealso: [](ch_ksp), `PCLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`, `PCLMVMGetMatLMVM()`
 27: @*/
 28: PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B)
 29: {
 30:   PC_LMVM  *ctx = (PC_LMVM *)pc->data;
 31:   PetscBool same;

 33:   PetscFunctionBegin;
 36:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
 37:   PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
 38:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
 39:   PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
 40:   PetscCall(MatDestroy(&ctx->B));
 41:   PetscCall(PetscObjectReference((PetscObject)B));
 42:   ctx->B = B;
 43:   PetscFunctionReturn(PETSC_SUCCESS);
 44: }

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

 49:   Input Parameter:
 50: . pc - An `PCLMVM` preconditioner

 52:   Output Parameter:
 53: . B - matrix inside the preconditioner, one of type `MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`

 55:   Level: intermediate

 57: .seealso: [](ch_ksp), `PCLMVM`, `MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`, `PCLMVMSetMatLMVM()`
 58: @*/
 59: PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B)
 60: {
 61:   PC_LMVM  *ctx = (PC_LMVM *)pc->data;
 62:   PetscBool same;

 64:   PetscFunctionBegin;
 66:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
 67:   PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
 68:   *B = ctx->B;
 69:   PetscFunctionReturn(PETSC_SUCCESS);
 70: }

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

 75:   Input Parameters:
 76: + pc       - An `PCLMVM` preconditioner
 77: - inactive - Index set defining the variables removed from the problem

 79:   Level: intermediate

 81:   Developer Notes:
 82:   Need to explain the purpose of this `IS`

 84: .seealso: [](ch_ksp), `PCLMVM`, `MatLMVMUpdate()`
 85: @*/
 86: PetscErrorCode PCLMVMSetIS(PC pc, IS inactive)
 87: {
 88:   PC_LMVM  *ctx = (PC_LMVM *)pc->data;
 89:   PetscBool same;

 91:   PetscFunctionBegin;
 94:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
 95:   PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
 96:   PetscCall(PCLMVMClearIS(pc));
 97:   PetscCall(PetscObjectReference((PetscObject)inactive));
 98:   ctx->inactive = inactive;
 99:   PetscFunctionReturn(PETSC_SUCCESS);
100: }

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

105:   Input Parameter:
106: . pc - An `PCLMVM` preconditioner

108:   Level: intermediate

110: .seealso: [](ch_ksp), `PCLMVM`, `MatLMVMUpdate()`
111: @*/
112: PetscErrorCode PCLMVMClearIS(PC pc)
113: {
114:   PC_LMVM  *ctx = (PC_LMVM *)pc->data;
115:   PetscBool same;

117:   PetscFunctionBegin;
119:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
120:   PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
121:   if (ctx->inactive) PetscCall(ISDestroy(&ctx->inactive));
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: static PetscErrorCode PCApply_LMVM(PC pc, Vec x, Vec y)
126: {
127:   PC_LMVM *ctx = (PC_LMVM *)pc->data;
128:   Vec      xsub, ysub;

130:   PetscFunctionBegin;
131:   if (ctx->inactive) {
132:     PetscCall(VecZeroEntries(ctx->xwork));
133:     PetscCall(VecGetSubVector(ctx->xwork, ctx->inactive, &xsub));
134:     PetscCall(VecCopy(x, xsub));
135:     PetscCall(VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub));
136:   } else {
137:     PetscCall(VecCopy(x, ctx->xwork));
138:   }
139:   PetscCall(MatSolve(ctx->B, ctx->xwork, ctx->ywork));
140:   if (ctx->inactive) {
141:     PetscCall(VecGetSubVector(ctx->ywork, ctx->inactive, &ysub));
142:     PetscCall(VecCopy(ysub, y));
143:     PetscCall(VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub));
144:   } else {
145:     PetscCall(VecCopy(ctx->ywork, y));
146:   }
147:   PetscFunctionReturn(PETSC_SUCCESS);
148: }

150: static PetscErrorCode PCReset_LMVM(PC pc)
151: {
152:   PC_LMVM *ctx = (PC_LMVM *)pc->data;

154:   PetscFunctionBegin;
155:   if (ctx->xwork) PetscCall(VecDestroy(&ctx->xwork));
156:   if (ctx->ywork) PetscCall(VecDestroy(&ctx->ywork));
157:   PetscFunctionReturn(PETSC_SUCCESS);
158: }

160: static PetscErrorCode PCSetUp_LMVM(PC pc)
161: {
162:   PC_LMVM    *ctx = (PC_LMVM *)pc->data;
163:   PetscInt    n, N;
164:   PetscBool   allocated;
165:   const char *prefix;

167:   PetscFunctionBegin;
168:   if (pc->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
169:   PetscCall(MatLMVMIsAllocated(ctx->B, &allocated));
170:   if (!allocated) {
171:     PetscCall(MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork));
172:     PetscCall(VecGetLocalSize(ctx->xwork, &n));
173:     PetscCall(VecGetSize(ctx->xwork, &N));
174:     PetscCall(MatSetSizes(ctx->B, n, n, N, N));
175:     PetscCall(MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork));
176:   } else {
177:     PetscCall(MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork));
178:   }
179:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
180:   PetscCall(MatSetOptionsPrefix(ctx->B, prefix));
181:   PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_"));
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: static PetscErrorCode PCView_LMVM(PC pc, PetscViewer viewer)
186: {
187:   PC_LMVM  *ctx = (PC_LMVM *)pc->data;
188:   PetscBool iascii;

190:   PetscFunctionBegin;
191:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
192:   if (iascii && ctx->B->assembled) {
193:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
194:     PetscCall(MatView(ctx->B, viewer));
195:     PetscCall(PetscViewerPopFormat(viewer));
196:   }
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }

200: static PetscErrorCode PCSetFromOptions_LMVM(PC pc, PetscOptionItems *PetscOptionsObject)
201: {
202:   PC_LMVM    *ctx = (PC_LMVM *)pc->data;
203:   const char *prefix;

205:   PetscFunctionBegin;
206:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
207:   PetscCall(MatSetOptionsPrefix(ctx->B, prefix));
208:   PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_"));
209:   PetscCall(MatSetFromOptions(ctx->B));
210:   PetscFunctionReturn(PETSC_SUCCESS);
211: }

213: static PetscErrorCode PCDestroy_LMVM(PC pc)
214: {
215:   PC_LMVM *ctx = (PC_LMVM *)pc->data;

217:   PetscFunctionBegin;
218:   if (ctx->inactive) PetscCall(ISDestroy(&ctx->inactive));
219:   if (pc->setupcalled) {
220:     PetscCall(VecDestroy(&ctx->xwork));
221:     PetscCall(VecDestroy(&ctx->ywork));
222:   }
223:   PetscCall(MatDestroy(&ctx->B));
224:   PetscCall(PetscFree(pc->data));
225:   PetscFunctionReturn(PETSC_SUCCESS);
226: }

228: /*MC
229:    PCLMVM - A a preconditioner constructed from a `MATLMVM` matrix. Options for the
230:             underlying `MATLMVM` matrix can be access with the -pc_lmvm_ prefix.

232:    Level: intermediate

234: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PCLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`,
235:           `PC`, `MATLMVM`, `PCLMVMUpdate()`, `PCLMVMSetMatLMVM()`, `PCLMVMGetMatLMVM()`
236: M*/
237: PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
238: {
239:   PC_LMVM *ctx;

241:   PetscFunctionBegin;
242:   PetscCall(PetscNew(&ctx));
243:   pc->data = (void *)ctx;

245:   pc->ops->reset               = PCReset_LMVM;
246:   pc->ops->setup               = PCSetUp_LMVM;
247:   pc->ops->destroy             = PCDestroy_LMVM;
248:   pc->ops->view                = PCView_LMVM;
249:   pc->ops->apply               = PCApply_LMVM;
250:   pc->ops->setfromoptions      = PCSetFromOptions_LMVM;
251:   pc->ops->applysymmetricleft  = NULL;
252:   pc->ops->applysymmetricright = NULL;
253:   pc->ops->applytranspose      = NULL;
254:   pc->ops->applyrichardson     = NULL;
255:   pc->ops->presolve            = NULL;
256:   pc->ops->postsolve           = NULL;

258:   PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B));
259:   PetscCall(MatSetType(ctx->B, MATLMVMBFGS));
260:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1));
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }