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: `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;

 35:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
 37:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
 39:   MatDestroy(&ctx->B);
 40:   PetscObjectReference((PetscObject)B);
 41:   ctx->B = B;
 42:   return 0;
 43: }

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

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

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

 54:    Level: intermediate

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

 64:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
 66:   *B = ctx->B;
 67:   return 0;
 68: }

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

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

 77:    Level: intermediate

 79:    Developer Note:
 80:    Need to explain the purpose of this `IS`

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

 91:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
 93:   PCLMVMClearIS(pc);
 94:   PetscObjectReference((PetscObject)inactive);
 95:   ctx->inactive = inactive;
 96:   return 0;
 97: }

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

102:    Input Parameter:
103: .  pc - An `PCLMVM` preconditioner

105:    Level: intermediate

107: .seealso: `PCLMVM`, `MatLMVMUpdate()`
108: @*/
109: PetscErrorCode PCLMVMClearIS(PC pc)
110: {
111:   PC_LMVM  *ctx = (PC_LMVM *)pc->data;
112:   PetscBool same;

115:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
117:   if (ctx->inactive) ISDestroy(&ctx->inactive);
118:   return 0;
119: }

121: static PetscErrorCode PCApply_LMVM(PC pc, Vec x, Vec y)
122: {
123:   PC_LMVM *ctx = (PC_LMVM *)pc->data;
124:   Vec      xsub, ysub;

126:   if (ctx->inactive) {
127:     VecZeroEntries(ctx->xwork);
128:     VecGetSubVector(ctx->xwork, ctx->inactive, &xsub);
129:     VecCopy(x, xsub);
130:     VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub);
131:   } else {
132:     VecCopy(x, ctx->xwork);
133:   }
134:   MatSolve(ctx->B, ctx->xwork, ctx->ywork);
135:   if (ctx->inactive) {
136:     VecGetSubVector(ctx->ywork, ctx->inactive, &ysub);
137:     VecCopy(ysub, y);
138:     VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub);
139:   } else {
140:     VecCopy(ctx->ywork, y);
141:   }
142:   return 0;
143: }

145: static PetscErrorCode PCReset_LMVM(PC pc)
146: {
147:   PC_LMVM *ctx = (PC_LMVM *)pc->data;

149:   if (ctx->xwork) VecDestroy(&ctx->xwork);
150:   if (ctx->ywork) VecDestroy(&ctx->ywork);
151:   return 0;
152: }

154: static PetscErrorCode PCSetUp_LMVM(PC pc)
155: {
156:   PC_LMVM    *ctx = (PC_LMVM *)pc->data;
157:   PetscInt    n, N;
158:   PetscBool   allocated;
159:   const char *prefix;

161:   if (pc->setupcalled) return 0;
162:   MatLMVMIsAllocated(ctx->B, &allocated);
163:   if (!allocated) {
164:     MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork);
165:     VecGetLocalSize(ctx->xwork, &n);
166:     VecGetSize(ctx->xwork, &N);
167:     MatSetSizes(ctx->B, n, n, N, N);
168:     MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork);
169:   } else {
170:     MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork);
171:   }
172:   PCGetOptionsPrefix(pc, &prefix);
173:   MatSetOptionsPrefix(ctx->B, prefix);
174:   MatAppendOptionsPrefix(ctx->B, "pc_lmvm_");
175:   return 0;
176: }

178: static PetscErrorCode PCView_LMVM(PC pc, PetscViewer viewer)
179: {
180:   PC_LMVM  *ctx = (PC_LMVM *)pc->data;
181:   PetscBool iascii;

183:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
184:   if (iascii && ctx->B->assembled) {
185:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
186:     MatView(ctx->B, viewer);
187:     PetscViewerPopFormat(viewer);
188:   }
189:   return 0;
190: }

192: static PetscErrorCode PCSetFromOptions_LMVM(PC pc, PetscOptionItems *PetscOptionsObject)
193: {
194:   PC_LMVM    *ctx = (PC_LMVM *)pc->data;
195:   const char *prefix;

197:   PCGetOptionsPrefix(pc, &prefix);
198:   MatSetOptionsPrefix(ctx->B, prefix);
199:   MatAppendOptionsPrefix(ctx->B, "pc_lmvm_");
200:   MatSetFromOptions(ctx->B);
201:   return 0;
202: }

204: static PetscErrorCode PCDestroy_LMVM(PC pc)
205: {
206:   PC_LMVM *ctx = (PC_LMVM *)pc->data;

208:   if (ctx->inactive) ISDestroy(&ctx->inactive);
209:   if (pc->setupcalled) {
210:     VecDestroy(&ctx->xwork);
211:     VecDestroy(&ctx->ywork);
212:   }
213:   MatDestroy(&ctx->B);
214:   PetscFree(pc->data);
215:   return 0;
216: }

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

222:    Level: intermediate

224: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`,
225:           `PC`, `MATLMVM`, `PCLMVMUpdate()`, `PCLMVMSetMatLMVM()`, `PCLMVMGetMatLMVM()`
226: M*/
227: PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
228: {
229:   PC_LMVM *ctx;

231:   PetscNew(&ctx);
232:   pc->data = (void *)ctx;

234:   pc->ops->reset               = PCReset_LMVM;
235:   pc->ops->setup               = PCSetUp_LMVM;
236:   pc->ops->destroy             = PCDestroy_LMVM;
237:   pc->ops->view                = PCView_LMVM;
238:   pc->ops->apply               = PCApply_LMVM;
239:   pc->ops->setfromoptions      = PCSetFromOptions_LMVM;
240:   pc->ops->applysymmetricleft  = NULL;
241:   pc->ops->applysymmetricright = NULL;
242:   pc->ops->applytranspose      = NULL;
243:   pc->ops->applyrichardson     = NULL;
244:   pc->ops->presolve            = NULL;
245:   pc->ops->postsolve           = NULL;

247:   MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B);
248:   MatSetType(ctx->B, MATLMVMBFGS);
249:   PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1);
250:   return 0;
251: }