Actual source code: lmvmpc.c

petsc-3.11.4 2019-09-28
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>

  8: typedef struct {
  9:   Vec  xwork, ywork;
 10:   IS   inactive;
 11:   Mat  B;
 12:   PetscBool allocated;
 13: } PC_LMVM;

 15: /*@
 16:    PCLMVMSetMatLMVM - Replaces the LMVM matrix inside the preconditioner with 
 17:    the one provided by the user.
 18:    
 19:    Input Parameters:
 20: +  pc - An LMVM preconditioner
 21: -  B  - An LMVM-type matrix (MATLDFP, MATLBFGS, MATLSR1, MATLBRDN, MATLMBRDN, MATLSBRDN)
 22:   
 23:    Level: intermediate
 24: @*/
 25: PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B)
 26: {
 27:   PC_LMVM          *ctx = (PC_LMVM*)pc->data;
 28:   PetscErrorCode   ierr;
 29:   PetscBool        same;
 30: 
 34:   PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
 35:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
 36:   PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
 37:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
 38:   MatDestroy(&ctx->B);
 39:   PetscObjectReference((PetscObject)B);
 40:   ctx->B = B;
 41:   return(0);
 42: }

 44: /*@
 45:    PCLMVMGetMatLMVM - Returns a pointer to the underlying LMVM matrix.
 46:    
 47:    Input Parameters:
 48: .  pc - An LMVM preconditioner

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

233:    Level: intermediate

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

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