Actual source code: telescope_coarsedm.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

  2: #include <petsc/private/matimpl.h>
  3: #include <petsc/private/pcimpl.h>
  4: #include <petsc/private/dmimpl.h>
  5: #include <petscksp.h>
  6: #include <petscdm.h>
  7: #include <petscdmda.h>
  8: #include <petscdmshell.h>

 10: #include "../src/ksp/pc/impls/telescope/telescope.h"

 12: static PetscBool  cited = PETSC_FALSE;
 13: static const char citation[] =
 14: "@inproceedings{MaySananRuppKnepleySmith2016,\n"
 15: "  title     = {Extreme-Scale Multigrid Components within PETSc},\n"
 16: "  author    = {Dave A. May and Patrick Sanan and Karl Rupp and Matthew G. Knepley and Barry F. Smith},\n"
 17: "  booktitle = {Proceedings of the Platform for Advanced Scientific Computing Conference},\n"
 18: "  series    = {PASC '16},\n"
 19: "  isbn      = {978-1-4503-4126-4},\n"
 20: "  location  = {Lausanne, Switzerland},\n"
 21: "  pages     = {5:1--5:12},\n"
 22: "  articleno = {5},\n"
 23: "  numpages  = {12},\n"
 24: "  url       = {https://doi.acm.org/10.1145/2929908.2929913},\n"
 25: "  doi       = {10.1145/2929908.2929913},\n"
 26: "  acmid     = {2929913},\n"
 27: "  publisher = {ACM},\n"
 28: "  address   = {New York, NY, USA},\n"
 29: "  keywords  = {GPU, HPC, agglomeration, coarse-level solver, multigrid, parallel computing, preconditioning},\n"
 30: "  year      = {2016}\n"
 31: "}\n";

 33: typedef struct {
 34:   DM              dm_fine,dm_coarse; /* these DM's should be topologically identical but use different communicators */
 35:   Mat             permutation;
 36:   Vec             xp;
 37:   PetscErrorCode  (*fp_dm_field_scatter)(DM,Vec,ScatterMode,DM,Vec);
 38:   PetscErrorCode  (*fp_dm_state_scatter)(DM,ScatterMode,DM);
 39:   void            *dmksp_context_determined;
 40:   void            *dmksp_context_user;
 41: } PC_Telescope_CoarseDMCtx;


 44: PetscErrorCode PCTelescopeSetUp_scatters_CoarseDM(PC pc,PC_Telescope sred,PC_Telescope_CoarseDMCtx *ctx)
 45: {
 47:   Vec            xred,yred,xtmp,x,xp;
 48:   VecScatter     scatter;
 49:   IS             isin;
 50:   Mat            B;
 51:   PetscInt       m,bs,st,ed;
 52:   MPI_Comm       comm;

 55:   PetscObjectGetComm((PetscObject)pc,&comm);
 56:   PCGetOperators(pc,NULL,&B);
 57:   MatCreateVecs(B,&x,NULL);
 58:   MatGetBlockSize(B,&bs);
 59:   VecDuplicate(x,&xp);
 60:   m = 0;
 61:   xred = NULL;
 62:   yred = NULL;
 63:   if (PCTelescope_isActiveRank(sred)) {
 64:     DMCreateGlobalVector(ctx->dm_coarse,&xred);
 65:     VecDuplicate(xred,&yred);
 66:     VecGetOwnershipRange(xred,&st,&ed);
 67:     ISCreateStride(comm,ed-st,st,1,&isin);
 68:     VecGetLocalSize(xred,&m);
 69:   } else {
 70:     VecGetOwnershipRange(x,&st,&ed);
 71:     ISCreateStride(comm,0,st,1,&isin);
 72:   }
 73:   ISSetBlockSize(isin,bs);
 74:   VecCreate(comm,&xtmp);
 75:   VecSetSizes(xtmp,m,PETSC_DECIDE);
 76:   VecSetBlockSize(xtmp,bs);
 77:   VecSetType(xtmp,((PetscObject)x)->type_name);
 78:   VecScatterCreate(x,isin,xtmp,NULL,&scatter);
 79:   sred->xred    = xred;
 80:   sred->yred    = yred;
 81:   sred->isin    = isin;
 82:   sred->scatter = scatter;
 83:   sred->xtmp    = xtmp;
 84:   ctx->xp       = xp;
 85:   VecDestroy(&x);
 86:   return(0);
 87: }

 89: PetscErrorCode PCTelescopeSetUp_CoarseDM(PC pc,PC_Telescope sred)
 90: {
 91:   PC_Telescope_CoarseDMCtx *ctx;
 92:   DM                       dm,dm_coarse = NULL;
 93:   MPI_Comm                 comm;
 94:   PetscBool                has_perm,has_kspcomputeoperators,using_kspcomputeoperators;
 95:   PetscErrorCode           ierr;

 98:   PetscInfo(pc,"PCTelescope: setup (CoarseDM)\n");
 99:   PetscNew(&ctx);
100:   sred->dm_ctx = (void*)ctx;

102:   PetscObjectGetComm((PetscObject)pc,&comm);
103:   PCGetDM(pc,&dm);
104:   DMGetCoarseDM(dm,&dm_coarse);
105:   ctx->dm_fine   = dm;
106:   ctx->dm_coarse = dm_coarse;

108:   /* attach coarse dm to ksp on sub communicator */
109:   if (PCTelescope_isActiveRank(sred)) {
110:     KSPSetDM(sred->ksp,ctx->dm_coarse);
111:     if (sred->ignore_kspcomputeoperators) {
112:       KSPSetDMActive(sred->ksp,PETSC_FALSE);
113:     }
114:   }

116:   /* check if there is a method to provide a permutation */
117:   has_perm = PETSC_FALSE;
118:   has_kspcomputeoperators = PETSC_FALSE;
119:   using_kspcomputeoperators = PETSC_FALSE;

121:   /* if no permutation is provided, we must rely on KSPSetComputeOperators */
122:   {
123:     PetscErrorCode (*dmfine_kspfunc)(KSP,Mat,Mat,void*) = NULL;
124:     void           *dmfine_kspctx = NULL,*dmcoarse_kspctx = NULL;
125:     void           *dmfine_appctx = NULL,*dmcoarse_appctx = NULL;
126:     void           *dmfine_shellctx = NULL,*dmcoarse_shellctx = NULL;

128:     DMKSPGetComputeOperators(dm,&dmfine_kspfunc,&dmfine_kspctx);
129:     if (dmfine_kspfunc) { has_kspcomputeoperators = PETSC_TRUE; }

131:     DMGetApplicationContext(ctx->dm_fine,&dmfine_appctx);
132:     DMShellGetContext(ctx->dm_fine,&dmfine_shellctx);

134:     /* need to define dmcoarse_kspctx */
135:     if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) {

137:       PetscInfo(pc,"PCTelescope: KSPSetComputeOperators fetched from parent DM\n");
138:       if (PCTelescope_isActiveRank(sred)) {
139:         DMGetApplicationContext(ctx->dm_coarse,&dmcoarse_appctx);
140:         DMShellGetContext(ctx->dm_coarse,&dmcoarse_shellctx);
141:       }

143:       /* Assume that if the fine operator didn't require any context, neither will the coarse */
144:       if (!dmfine_kspctx) {
145:         dmcoarse_kspctx = NULL;
146:         PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using NULL context\n");
147:       } else {

149:         PetscInfo(pc,"PCTelescope: KSPSetComputeOperators detected non-NULL context from parent DM \n");
150:         if (PCTelescope_isActiveRank(sred)) {

152:           if (dmfine_kspctx == dmfine_appctx) {
153:             dmcoarse_kspctx = dmcoarse_appctx;
154:             PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using context from DM->ApplicationContext\n");
155:             if (!dmcoarse_kspctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Non NULL dmfine->kspctx == dmfine->appctx. NULL dmcoarse->appctx found. Likely this is an error");
156:           } else if (dmfine_kspctx == dmfine_shellctx) {
157:             dmcoarse_kspctx = dmcoarse_shellctx;
158:             PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using context from DMShell->Context\n");
159:             if (!dmcoarse_kspctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Non NULL dmfine->kspctx == dmfine.shell->ctx. NULL dmcoarse.shell->ctx found. Likely this is an error");
160:           }
161:           ctx->dmksp_context_determined = dmcoarse_kspctx;

163:           /* look for user provided method to fetch the context */
164:           {
165:             PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL;
166:             void *dmcoarse_context_user = NULL;
167:             char dmcoarse_method[PETSC_MAX_PATH_LEN];

169:             PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMKSPContext");
170:             PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context);
171:             if (fp_get_coarsedm_context) {
172:               PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n");
173:               fp_get_coarsedm_context(ctx->dm_coarse,&dmcoarse_context_user);
174:               ctx->dmksp_context_user = dmcoarse_context_user;
175:               dmcoarse_kspctx = dmcoarse_context_user;
176:             } else {
177:               PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n");
178:             }
179:           }

181:           if (!dmcoarse_kspctx) {
182:             PetscInfo(pc,"PCTelescope: KSPSetComputeOperators failed to determine the context to use on sub-communicator\n");
183:             SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot determine which context with use for KSPSetComputeOperators() on sub-communicator");
184:           }
185:         }
186:       }
187:     }

189:     if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) {
190:       using_kspcomputeoperators = PETSC_TRUE;

192:       if (PCTelescope_isActiveRank(sred)) {
193:         /* sub ksp inherits dmksp_func and context provided by user */
194:         KSPSetComputeOperators(sred->ksp,dmfine_kspfunc,dmcoarse_kspctx);
195:         /*PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)ctx->dmrepart);*/
196:         KSPSetDMActive(sred->ksp,PETSC_TRUE);
197:       }
198:     }
199:   }

201:   if (!has_perm && has_kspcomputeoperators && !using_kspcomputeoperators) SETERRQ(comm,PETSC_ERR_SUP,"No method to permute an operator was found on the parent DM. A method for KSPSetComputeOperators() was provided but it was requested to be ignored. Telescope setup cannot proceed");
202:   if (!has_perm && !has_kspcomputeoperators) SETERRQ(comm,PETSC_ERR_SUP,"No method to permute an operator was found on the parent DM. No method for KSPSetComputeOperators() was provided. Telescope setup cannot proceed");

204:   {
205:     char dmfine_method[PETSC_MAX_PATH_LEN];

207:     PetscSNPrintf(dmfine_method,sizeof(dmfine_method),"PCTelescopeFieldScatter");
208:     PetscObjectQueryFunction((PetscObject)ctx->dm_fine,dmfine_method,&ctx->fp_dm_field_scatter);

210:     PetscSNPrintf(dmfine_method,sizeof(dmfine_method),"PCTelescopeStateScatter");
211:     PetscObjectQueryFunction((PetscObject)ctx->dm_fine,dmfine_method,&ctx->fp_dm_state_scatter);
212:   }

214:   if (ctx->fp_dm_state_scatter) {
215:     PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeStateScatter from parent DM\n");
216:   } else {
217:     PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeStateScatter from parent DM\n");
218:   }

220:   if (ctx->fp_dm_field_scatter) {
221:     PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeFieldScatter from parent DM\n");
222:   } else {
223:     PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeFieldScatter from parent DM\n");
224:     SETERRQ(comm,PETSC_ERR_SUP,"No method to scatter fields between the parent DM and coarse DM was found. Must call PetscObjectComposeFunction() with the parent DM. Telescope setup cannot proceed");
225:   }

227:   /*PCTelescopeSetUp_permutation_CoarseDM(pc,sred,ctx);*/
228:   PCTelescopeSetUp_scatters_CoarseDM(pc,sred,ctx);
229:   return(0);
230: }

232: PetscErrorCode PCApply_Telescope_CoarseDM(PC pc,Vec x,Vec y)
233: {
234:   PC_Telescope             sred = (PC_Telescope)pc->data;
235:   PetscErrorCode           ierr;
236:   Vec                      xred,yred;
237:   PC_Telescope_CoarseDMCtx *ctx;

240:   ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
241:   xred = sred->xred;
242:   yred = sred->yred;

244:   PetscCitationsRegister(citation,&cited);

246:   if (ctx->fp_dm_state_scatter) {
247:     ctx->fp_dm_state_scatter(ctx->dm_fine,SCATTER_FORWARD,ctx->dm_coarse);
248:   }

250:   ctx->fp_dm_field_scatter(ctx->dm_fine,x,SCATTER_FORWARD,ctx->dm_coarse,xred);

252:   /* solve */
253:   if (PCTelescope_isActiveRank(sred)) {
254:     KSPSolve(sred->ksp,xred,yred);
255:   }

257:   ctx->fp_dm_field_scatter(ctx->dm_fine,y,SCATTER_REVERSE,ctx->dm_coarse,yred);
258:   return(0);
259: }

261: PetscErrorCode PCTelescopeSubNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,MatNullSpace nullspace,MatNullSpace *sub_nullspace)
262: {
263:   PetscErrorCode           ierr;
264:   PetscBool                has_const;
265:   PetscInt                 k,n = 0;
266:   const Vec                *vecs;
267:   Vec                      *sub_vecs = NULL;
268:   MPI_Comm                 subcomm;
269:   PC_Telescope_CoarseDMCtx *ctx;

272:   ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
273:   subcomm = sred->subcomm;
274:   MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs);

276:   if (PCTelescope_isActiveRank(sred)) {
277:     /* create new vectors */
278:     if (n) {
279:       VecDuplicateVecs(sred->xred,n,&sub_vecs);
280:     }
281:   }

283:   /* copy entries */
284:   for (k=0; k<n; k++) {
285:     ctx->fp_dm_field_scatter(ctx->dm_fine,vecs[k],SCATTER_FORWARD,ctx->dm_coarse,sub_vecs[k]);
286:   }

288:   if (PCTelescope_isActiveRank(sred)) {
289:     /* create new (near) nullspace for redundant object */
290:     MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,sub_nullspace);
291:     VecDestroyVecs(n,&sub_vecs);
292:   }
293:   return(0);
294: }

296: PetscErrorCode PCTelescopeMatNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,Mat sub_mat)
297: {
298:   PetscErrorCode           ierr;
299:   Mat                      B;
300:   PC_Telescope_CoarseDMCtx *ctx;

303:   ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
304:   PCGetOperators(pc,NULL,&B);
305:   {
306:     MatNullSpace nullspace,sub_nullspace;
307:     MatGetNullSpace(B,&nullspace);
308:     if (nullspace) {
309:       PetscInfo(pc,"PCTelescope: generating nullspace (CoarseDM)\n");
310:       PCTelescopeSubNullSpaceCreate_CoarseDM(pc,sred,nullspace,&sub_nullspace);

312:       /* attach any user nullspace removal methods and contexts */
313:       if (PCTelescope_isActiveRank(sred)) {
314:         void *context = NULL;
315:         if (nullspace->remove && !nullspace->rmctx){
316:           MatNullSpaceSetFunction(sub_nullspace,nullspace->remove,context);
317:         } else if (nullspace->remove && nullspace->rmctx) {
318:           char           dmcoarse_method[PETSC_MAX_PATH_LEN];
319:           PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL;

321:           PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMNullSpaceUserContext");
322:           PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context);
323:           if (!context) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Propagation of user null-space removal method with non-NULL context requires the coarse DM be composed with a function named \"%s\"",dmcoarse_method);
324:           MatNullSpaceSetFunction(sub_nullspace,nullspace->remove,context);
325:         }
326:       }

328:       if (PCTelescope_isActiveRank(sred)) {
329:         MatSetNullSpace(sub_mat,sub_nullspace);
330:         MatNullSpaceDestroy(&sub_nullspace);
331:       }
332:     }
333:   }
334:   {
335:     MatNullSpace nearnullspace,sub_nearnullspace;
336:     MatGetNearNullSpace(B,&nearnullspace);
337:     if (nearnullspace) {
338:       PetscInfo(pc,"PCTelescope: generating near nullspace (CoarseDM)\n");
339:       PCTelescopeSubNullSpaceCreate_CoarseDM(pc,sred,nearnullspace,&sub_nearnullspace);

341:       /* attach any user nullspace removal methods and contexts */
342:       if (PCTelescope_isActiveRank(sred)) {
343:         void *context = NULL;
344:         if (nearnullspace->remove && !nearnullspace->rmctx){
345:           MatNullSpaceSetFunction(sub_nearnullspace,nearnullspace->remove,context);
346:         } else if (nearnullspace->remove && nearnullspace->rmctx) {
347:           char           dmcoarse_method[PETSC_MAX_PATH_LEN];
348:           PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL;

350:           PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMNearNullSpaceUserContext");
351:           PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context);
352:           if (!context) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Propagation of user near null-space removal method with non-NULL context requires the coarse DM be composed with a function named \"%s\"",dmcoarse_method);
353:           MatNullSpaceSetFunction(sub_nearnullspace,nearnullspace->remove,context);
354:         }
355:       }

357:       if (PCTelescope_isActiveRank(sred)) {
358:         MatSetNearNullSpace(sub_mat,sub_nearnullspace);
359:         MatNullSpaceDestroy(&sub_nearnullspace);
360:       }
361:     }
362:   }
363:   return(0);
364: }

366: PetscErrorCode PCReset_Telescope_CoarseDM(PC pc)
367: {
368:   PetscErrorCode           ierr;
369:   PC_Telescope             sred = (PC_Telescope)pc->data;
370:   PC_Telescope_CoarseDMCtx *ctx;

373:   ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
374:   ctx->dm_fine = NULL; /* since I did not increment the ref counter we set these to NULL */
375:   ctx->dm_coarse = NULL; /* since I did not increment the ref counter we set these to NULL */
376:   ctx->permutation = NULL; /* this will be fetched from the dm so no need to call destroy */
377:   VecDestroy(&ctx->xp);
378:   ctx->fp_dm_field_scatter = NULL;
379:   ctx->fp_dm_state_scatter = NULL;
380:   ctx->dmksp_context_determined = NULL;
381:   ctx->dmksp_context_user = NULL;
382:   return(0);
383: }

385: PetscErrorCode PCApplyRichardson_Telescope_CoarseDM(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt its,PetscBool zeroguess,PetscInt *outits,PCRichardsonConvergedReason *reason)
386: {
387:   PC_Telescope             sred = (PC_Telescope)pc->data;
388:   PetscErrorCode           ierr;
389:   Vec                      yred = NULL;
390:   PetscBool                default_init_guess_value = PETSC_FALSE;
391:   PC_Telescope_CoarseDMCtx *ctx;

394:   ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
395:   yred = sred->yred;

397:   if (its > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCApplyRichardson_Telescope_CoarseDM only supports max_it = 1");
398:   *reason = (PCRichardsonConvergedReason)0;

400:   if (!zeroguess) {
401:     PetscInfo(pc,"PCTelescopeCoarseDM: Scattering y for non-zero-initial guess\n");

403:     ctx->fp_dm_field_scatter(ctx->dm_fine,y,SCATTER_FORWARD,ctx->dm_coarse,yred);
404:   }

406:   if (PCTelescope_isActiveRank(sred)) {
407:     KSPGetInitialGuessNonzero(sred->ksp,&default_init_guess_value);
408:     if (!zeroguess) KSPSetInitialGuessNonzero(sred->ksp,PETSC_TRUE);
409:   }

411:   PCApply_Telescope_CoarseDM(pc,x,y);

413:   if (PCTelescope_isActiveRank(sred)) {
414:     KSPSetInitialGuessNonzero(sred->ksp,default_init_guess_value);
415:   }

417:   if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
418:   *outits = 1;
419:   return(0);
420: }