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 LMVM matrix inside the preconditioner with
 18:    the one provided by the user.

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

 24:    Level: intermediate
 25: @*/
 26: PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B)
 27: {
 28:   PC_LMVM          *ctx = (PC_LMVM*)pc->data;
 29:   PetscBool        same;

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

 43: /*@
 44:    PCLMVMGetMatLMVM - Returns a pointer to the underlying LMVM matrix.

 46:    Input Parameters:
 47: .  pc - An LMVM preconditioner

 49:    Output Parameters:
 50: .  B - LMVM matrix inside the preconditioner

 52:    Level: intermediate
 53: @*/
 54: PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B)
 55: {
 56:   PC_LMVM          *ctx = (PC_LMVM*)pc->data;
 57:   PetscBool        same;

 60:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
 62:   *B = ctx->B;
 63:   return 0;
 64: }

 66: /*@
 67:    PCLMVMSetIS - Sets the index sets that reduce the PC application.

 69:    Input Parameters:
 70: +  pc - An LMVM preconditioner
 71: -  inactive - Index set defining the variables removed from the problem

 73:    Level: intermediate

 75: .seealso:  MatLMVMUpdate()
 76: @*/
 77: PetscErrorCode PCLMVMSetIS(PC pc, IS inactive)
 78: {
 79:   PC_LMVM          *ctx = (PC_LMVM*)pc->data;
 80:   PetscBool        same;

 84:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
 86:   PCLMVMClearIS(pc);
 87:   PetscObjectReference((PetscObject)inactive);
 88:   ctx->inactive = inactive;
 89:   return 0;
 90: }

 92: /*@
 93:    PCLMVMClearIS - Removes the inactive variable index set.

 95:    Input Parameters:
 96: .  pc - An LMVM preconditioner

 98:    Level: intermediate

100: .seealso:  MatLMVMUpdate()
101: @*/
102: PetscErrorCode PCLMVMClearIS(PC pc)
103: {
104:   PC_LMVM          *ctx = (PC_LMVM*)pc->data;
105:   PetscBool        same;

108:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
110:   if (ctx->inactive) {
111:     ISDestroy(&ctx->inactive);
112:   }
113:   return 0;
114: }

116: static PetscErrorCode PCApply_LMVM(PC pc,Vec x,Vec y)
117: {
118:   PC_LMVM          *ctx = (PC_LMVM*)pc->data;
119:   Vec              xsub, ysub;

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

140: static PetscErrorCode PCReset_LMVM(PC pc)
141: {
142:   PC_LMVM        *ctx = (PC_LMVM*)pc->data;

144:   if (ctx->xwork) {
145:     VecDestroy(&ctx->xwork);
146:   }
147:   if (ctx->ywork) {
148:     VecDestroy(&ctx->ywork);
149:   }
150:   return 0;
151: }

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

159:   MatLMVMIsAllocated(ctx->B, &allocated);
160:   if (!allocated) {
161:     MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork);
162:     VecGetLocalSize(ctx->xwork, &n);
163:     VecGetSize(ctx->xwork, &N);
164:     MatSetSizes(ctx->B, n, n, N, N);
165:     MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork);
166:   } else {
167:     MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork);
168:   }
169:   return 0;
170: }

172: static PetscErrorCode PCView_LMVM(PC pc,PetscViewer viewer)
173: {
174:   PC_LMVM        *ctx = (PC_LMVM*)pc->data;
175:   PetscBool      iascii;

177:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
178:   if (iascii && ctx->B->assembled) {
179:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
180:     MatView(ctx->B, viewer);
181:     PetscViewerPopFormat(viewer);
182:   }
183:   return 0;
184: }

186: static PetscErrorCode PCSetFromOptions_LMVM(PetscOptionItems* PetscOptionsObject, PC pc)
187: {
188:   PC_LMVM        *ctx = (PC_LMVM*)pc->data;

190:   MatSetFromOptions(ctx->B);
191:   return 0;
192: }

194: static PetscErrorCode PCDestroy_LMVM(PC pc)
195: {
196:   PC_LMVM        *ctx = (PC_LMVM*)pc->data;

198:   if (ctx->inactive) {
199:     ISDestroy(&ctx->inactive);
200:   }
201:   if (pc->setupcalled) {
202:     VecDestroy(&ctx->xwork);
203:     VecDestroy(&ctx->ywork);
204:   }
205:   MatDestroy(&ctx->B);
206:   PetscFree(pc->data);
207:   return 0;
208: }

210: /*MC
211:    PCLMVM - Creates a preconditioner around an LMVM matrix. Options for the
212:             underlying LMVM matrix can be access with the "-pc_lmvm_" prefix.

214:    Level: intermediate

216: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types),
217:            PC, MATLMVM, PCLMVMUpdate(), PCLMVMSetMatLMVM(), PCLMVMGetMatLMVM()
218: M*/
219: PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
220: {
221:   PC_LMVM        *ctx;

223:   PetscNewLog(pc,&ctx);
224:   pc->data = (void*)ctx;

226:   pc->ops->reset           = PCReset_LMVM;
227:   pc->ops->setup           = PCSetUp_LMVM;
228:   pc->ops->destroy         = PCDestroy_LMVM;
229:   pc->ops->view            = PCView_LMVM;
230:   pc->ops->apply           = PCApply_LMVM;
231:   pc->ops->setfromoptions  = PCSetFromOptions_LMVM;
232:   pc->ops->applysymmetricleft  = NULL;
233:   pc->ops->applysymmetricright = NULL;
234:   pc->ops->applytranspose  = NULL;
235:   pc->ops->applyrichardson = NULL;
236:   pc->ops->presolve        = NULL;
237:   pc->ops->postsolve       = NULL;

239:   PCSetReusePreconditioner(pc, PETSC_TRUE);

241:   MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B);
242:   MatSetType(ctx->B, MATLMVMBFGS);
243:   PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1);
244:   MatSetOptionsPrefix(ctx->B, "pc_lmvm_");
245:   return 0;
246: }