Actual source code: lmvmpc.c

petsc-master 2020-06-03
Report Typos and Errors
  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.
 19:    
 20:    Input Parameters:
 21: +  pc - An LMVM preconditioner
 22: -  B  - An LMVM-type matrix (MATLDFP, MATLBFGS, MATLSR1, MATLBRDN, MATLMBRDN, MATLSBRDN)
 23:   
 24:    Level: intermediate
 25: @*/
 26: PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B)
 27: {
 28:   PC_LMVM          *ctx = (PC_LMVM*)pc->data;
 29:   PetscErrorCode   ierr;
 30:   PetscBool        same;
 31:   
 35:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
 36:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
 37:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
 38:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
 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 LMVM matrix.
 47:    
 48:    Input Parameters:
 49: .  pc - An LMVM preconditioner

 51:    Output Parameters:
 52: .  B - LMVM matrix inside the preconditioner
 53:   
 54:    Level: intermediate
 55: @*/
 56: PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B)
 57: {
 58:   PC_LMVM          *ctx = (PC_LMVM*)pc->data;
 59:   PetscErrorCode   ierr;
 60:   PetscBool        same;
 61:   
 64:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
 65:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
 66:   *B = ctx->B;
 67:   return(0);
 68: }

 70: /*@
 71:    PCLMVMSetIS - Sets the index sets that reduce the PC application.
 72:    
 73:    Input Parameters:
 74: +  pc - An LMVM preconditioner
 75: -  inactive - Index set defining the variables removed from the problem
 76:   
 77:    Level: intermediate

 79: .seealso:  MatLMVMUpdate()
 80: @*/
 81: PetscErrorCode PCLMVMSetIS(PC pc, IS inactive)
 82: {
 83:   PC_LMVM          *ctx = (PC_LMVM*)pc->data;
 84:   PetscErrorCode   ierr;
 85:   PetscBool        same;
 86:   
 90:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
 91:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
 92:   PCLMVMClearIS(pc);
 93:   PetscObjectReference((PetscObject)inactive);
 94:   ctx->inactive = inactive;
 95:   return(0);
 96: }

 98: /*@
 99:    PCLMVMClearIS - Removes the inactive variable index set.
100:    
101:    Input Parameters:
102: .  pc - An LMVM preconditioner
103:   
104:    Level: intermediate

106: .seealso:  MatLMVMUpdate()
107: @*/
108: PetscErrorCode PCLMVMClearIS(PC pc)
109: {
110:   PC_LMVM          *ctx = (PC_LMVM*)pc->data;
111:   PetscErrorCode   ierr;
112:   PetscBool        same;
113:   
116:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
117:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
118:   if (ctx->inactive) {
119:     ISDestroy(&ctx->inactive);
120:   }
121:   return(0);
122: }

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

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

150: static PetscErrorCode PCReset_LMVM(PC pc)
151: {
152:   PC_LMVM        *ctx = (PC_LMVM*)pc->data;
154:   
156:   if (ctx->xwork) {
157:     VecDestroy(&ctx->xwork);
158:   }
159:   if (ctx->ywork) {
160:     VecDestroy(&ctx->ywork);
161:   }
162:   return(0);
163: }

165: static PetscErrorCode PCSetUp_LMVM(PC pc)
166: {
167:   PC_LMVM        *ctx = (PC_LMVM*)pc->data;
169:   PetscInt       n, N;
170:   PetscBool      allocated;

173:   MatLMVMIsAllocated(ctx->B, &allocated);
174:   if (!allocated) {
175:     MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork);
176:     VecGetLocalSize(ctx->xwork, &n);
177:     VecGetSize(ctx->xwork, &N);
178:     MatSetSizes(ctx->B, n, n, N, N);
179:     MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork);
180:   } else {
181:     MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork);
182:   }
183:   return(0);
184: }

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

193:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
194:   if (iascii && ctx->B->assembled) {
195:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
196:     MatView(ctx->B, viewer);
197:     PetscViewerPopFormat(viewer);
198:   }
199:   return(0);
200: }

202: static PetscErrorCode PCSetFromOptions_LMVM(PetscOptionItems* PetscOptionsObject, PC pc)
203: {
204:   PC_LMVM        *ctx = (PC_LMVM*)pc->data;
206:   
208:   MatSetFromOptions(ctx->B);
209:   return(0);
210: }

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

218:   if (ctx->inactive) {
219:     ISDestroy(&ctx->inactive);
220:   }
221:   if (pc->setupcalled) {
222:     VecDestroy(&ctx->xwork);
223:     VecDestroy(&ctx->ywork);
224:   }
225:   MatDestroy(&ctx->B);
226:   PetscFree(pc->data);
227:   return(0);
228: }

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

234:    Level: intermediate

236: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), 
237:            PC, MATLMVM, PCLMVMUpdate(), PCLMVMSetMatLMVM(), PCLMVMGetMatLMVM()
238: M*/
239: PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
240: {
242:   PC_LMVM        *ctx;

245:   PetscNewLog(pc,&ctx);
246:   pc->data = (void*)ctx;
247:   
248:   pc->ops->reset           = PCReset_LMVM;
249:   pc->ops->setup           = PCSetUp_LMVM;
250:   pc->ops->destroy         = PCDestroy_LMVM;
251:   pc->ops->view            = PCView_LMVM;
252:   pc->ops->apply           = PCApply_LMVM;
253:   pc->ops->setfromoptions  = PCSetFromOptions_LMVM;
254:   pc->ops->applysymmetricleft  = NULL;
255:   pc->ops->applysymmetricright = NULL;
256:   pc->ops->applytranspose  = NULL;
257:   pc->ops->applyrichardson = NULL;
258:   pc->ops->presolve        = NULL;
259:   pc->ops->postsolve       = NULL;
260:   
261:   PCSetReusePreconditioner(pc, PETSC_TRUE);
262:   
263:   MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B);
264:   MatSetType(ctx->B, MATLMVMBFGS);
265:   PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1);
266:   MatSetOptionsPrefix(ctx->B, "pc_lmvm_");
267:   return(0);
268: }