Actual source code: telescope_coarsedm.c


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

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

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

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

 85: PetscErrorCode PCTelescopeSetUp_CoarseDM(PC pc, PC_Telescope sred)
 86: {
 87:   PC_Telescope_CoarseDMCtx *ctx;
 88:   DM                        dm, dm_coarse = NULL;
 89:   MPI_Comm                  comm;
 90:   PetscBool                 has_perm, has_kspcomputeoperators, using_kspcomputeoperators;

 92:   PetscInfo(pc, "PCTelescope: setup (CoarseDM)\n");
 93:   PetscNew(&ctx);
 94:   sred->dm_ctx = (void *)ctx;

 96:   PetscObjectGetComm((PetscObject)pc, &comm);
 97:   PCGetDM(pc, &dm);
 98:   DMGetCoarseDM(dm, &dm_coarse);
 99:   ctx->dm_fine   = dm;
100:   ctx->dm_coarse = dm_coarse;

102:   /* attach coarse dm to ksp on sub communicator */
103:   if (PCTelescope_isActiveRank(sred)) {
104:     KSPSetDM(sred->ksp, ctx->dm_coarse);
105:     if (sred->ignore_kspcomputeoperators) KSPSetDMActive(sred->ksp, PETSC_FALSE);
106:   }

108:   /* check if there is a method to provide a permutation */
109:   has_perm                  = PETSC_FALSE;
110:   has_kspcomputeoperators   = PETSC_FALSE;
111:   using_kspcomputeoperators = PETSC_FALSE;

113:   /* if no permutation is provided, we must rely on KSPSetComputeOperators */
114:   {
115:     PetscErrorCode (*dmfine_kspfunc)(KSP, Mat, Mat, void *) = NULL;
116:     void *dmfine_kspctx = NULL, *dmcoarse_kspctx = NULL;
117:     void *dmfine_appctx = NULL, *dmcoarse_appctx = NULL;
118:     void *dmfine_shellctx = NULL, *dmcoarse_shellctx = NULL;

120:     DMKSPGetComputeOperators(dm, &dmfine_kspfunc, &dmfine_kspctx);
121:     if (dmfine_kspfunc) has_kspcomputeoperators = PETSC_TRUE;

123:     DMGetApplicationContext(ctx->dm_fine, &dmfine_appctx);
124:     DMShellGetContext(ctx->dm_fine, &dmfine_shellctx);

126:     /* need to define dmcoarse_kspctx */
127:     if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) {
128:       PetscInfo(pc, "PCTelescope: KSPSetComputeOperators fetched from parent DM\n");
129:       if (PCTelescope_isActiveRank(sred)) {
130:         DMGetApplicationContext(ctx->dm_coarse, &dmcoarse_appctx);
131:         DMShellGetContext(ctx->dm_coarse, &dmcoarse_shellctx);
132:       }

134:       /* Assume that if the fine operator didn't require any context, neither will the coarse */
135:       if (!dmfine_kspctx) {
136:         dmcoarse_kspctx = NULL;
137:         PetscInfo(pc, "PCTelescope: KSPSetComputeOperators using NULL context\n");
138:       } else {
139:         PetscInfo(pc, "PCTelescope: KSPSetComputeOperators detected non-NULL context from parent DM \n");
140:         if (PCTelescope_isActiveRank(sred)) {
141:           if (dmfine_kspctx == dmfine_appctx) {
142:             dmcoarse_kspctx = dmcoarse_appctx;
143:             PetscInfo(pc, "PCTelescope: KSPSetComputeOperators using context from DM->ApplicationContext\n");
145:           } else if (dmfine_kspctx == dmfine_shellctx) {
146:             dmcoarse_kspctx = dmcoarse_shellctx;
147:             PetscInfo(pc, "PCTelescope: KSPSetComputeOperators using context from DMShell->Context\n");
149:           }
150:           ctx->dmksp_context_determined = dmcoarse_kspctx;

152:           /* look for user provided method to fetch the context */
153:           {
154:             PetscErrorCode (*fp_get_coarsedm_context)(DM, void **) = NULL;
155:             void *dmcoarse_context_user                            = NULL;
156:             char  dmcoarse_method[PETSC_MAX_PATH_LEN];

158:             PetscSNPrintf(dmcoarse_method, sizeof(dmcoarse_method), "PCTelescopeGetCoarseDMKSPContext");
159:             PetscObjectQueryFunction((PetscObject)ctx->dm_coarse, dmcoarse_method, &fp_get_coarsedm_context);
160:             if (fp_get_coarsedm_context) {
161:               PetscInfo(pc, "PCTelescope: Found composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n");
162:               fp_get_coarsedm_context(ctx->dm_coarse, &dmcoarse_context_user);
163:               ctx->dmksp_context_user = dmcoarse_context_user;
164:               dmcoarse_kspctx         = dmcoarse_context_user;
165:             } else {
166:               PetscInfo(pc, "PCTelescope: Failed to find composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n");
167:             }
168:           }

170:           if (!dmcoarse_kspctx) {
171:             PetscInfo(pc, "PCTelescope: KSPSetComputeOperators failed to determine the context to use on sub-communicator\n");
172:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot determine which context with use for KSPSetComputeOperators() on sub-communicator");
173:           }
174:         }
175:       }
176:     }

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

181:       if (PCTelescope_isActiveRank(sred)) {
182:         /* sub ksp inherits dmksp_func and context provided by user */
183:         KSPSetComputeOperators(sred->ksp, dmfine_kspfunc, dmcoarse_kspctx);
184:         /* PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)ctx->dmrepart); */
185:         KSPSetDMActive(sred->ksp, PETSC_TRUE);
186:       }
187:     }
188:   }


193:   {
194:     char dmfine_method[PETSC_MAX_PATH_LEN];

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

199:     PetscSNPrintf(dmfine_method, sizeof(dmfine_method), "PCTelescopeStateScatter");
200:     PetscObjectQueryFunction((PetscObject)ctx->dm_fine, dmfine_method, &ctx->fp_dm_state_scatter);
201:   }

203:   if (ctx->fp_dm_state_scatter) {
204:     PetscInfo(pc, "PCTelescope: Found composed method PCTelescopeStateScatter from parent DM\n");
205:   } else {
206:     PetscInfo(pc, "PCTelescope: Failed to find composed method PCTelescopeStateScatter from parent DM\n");
207:   }

209:   if (ctx->fp_dm_field_scatter) {
210:     PetscInfo(pc, "PCTelescope: Found composed method PCTelescopeFieldScatter from parent DM\n");
211:   } else {
212:     PetscInfo(pc, "PCTelescope: Failed to find composed method PCTelescopeFieldScatter from parent DM\n");
213:     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");
214:   }

216:   /* PCTelescopeSetUp_permutation_CoarseDM(pc,sred,ctx); */
217:   PCTelescopeSetUp_scatters_CoarseDM(pc, sred, ctx);
218:   return 0;
219: }

221: PetscErrorCode PCApply_Telescope_CoarseDM(PC pc, Vec x, Vec y)
222: {
223:   PC_Telescope              sred = (PC_Telescope)pc->data;
224:   Vec                       xred, yred;
225:   PC_Telescope_CoarseDMCtx *ctx;

227:   ctx  = (PC_Telescope_CoarseDMCtx *)sred->dm_ctx;
228:   xred = sred->xred;
229:   yred = sred->yred;

231:   PetscCitationsRegister(citation, &cited);

233:   if (ctx->fp_dm_state_scatter) ctx->fp_dm_state_scatter(ctx->dm_fine, SCATTER_FORWARD, ctx->dm_coarse);

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

237:   /* solve */
238:   if (PCTelescope_isActiveRank(sred)) KSPSolve(sred->ksp, xred, yred);

240:   ctx->fp_dm_field_scatter(ctx->dm_fine, y, SCATTER_REVERSE, ctx->dm_coarse, yred);
241:   return 0;
242: }

244: PetscErrorCode PCTelescopeSubNullSpaceCreate_CoarseDM(PC pc, PC_Telescope sred, MatNullSpace nullspace, MatNullSpace *sub_nullspace)
245: {
246:   PetscBool                 has_const;
247:   PetscInt                  k, n = 0;
248:   const Vec                *vecs;
249:   Vec                      *sub_vecs = NULL;
250:   MPI_Comm                  subcomm;
251:   PC_Telescope_CoarseDMCtx *ctx;

253:   ctx     = (PC_Telescope_CoarseDMCtx *)sred->dm_ctx;
254:   subcomm = sred->subcomm;
255:   MatNullSpaceGetVecs(nullspace, &has_const, &n, &vecs);

257:   if (PCTelescope_isActiveRank(sred)) {
258:     /* create new vectors */
259:     if (n) VecDuplicateVecs(sred->xred, n, &sub_vecs);
260:   }

262:   /* copy entries */
263:   for (k = 0; k < n; k++) ctx->fp_dm_field_scatter(ctx->dm_fine, vecs[k], SCATTER_FORWARD, ctx->dm_coarse, sub_vecs[k]);

265:   if (PCTelescope_isActiveRank(sred)) {
266:     /* create new (near) nullspace for redundant object */
267:     MatNullSpaceCreate(subcomm, has_const, n, sub_vecs, sub_nullspace);
268:     VecDestroyVecs(n, &sub_vecs);
269:   }
270:   return 0;
271: }

273: PetscErrorCode PCTelescopeMatNullSpaceCreate_CoarseDM(PC pc, PC_Telescope sred, Mat sub_mat)
274: {
275:   Mat                       B;
276:   PC_Telescope_CoarseDMCtx *ctx;

278:   ctx = (PC_Telescope_CoarseDMCtx *)sred->dm_ctx;
279:   PCGetOperators(pc, NULL, &B);
280:   {
281:     MatNullSpace nullspace, sub_nullspace;
282:     MatGetNullSpace(B, &nullspace);
283:     if (nullspace) {
284:       PetscInfo(pc, "PCTelescope: generating nullspace (CoarseDM)\n");
285:       PCTelescopeSubNullSpaceCreate_CoarseDM(pc, sred, nullspace, &sub_nullspace);

287:       /* attach any user nullspace removal methods and contexts */
288:       if (PCTelescope_isActiveRank(sred)) {
289:         void *context = NULL;
290:         if (nullspace->remove && !nullspace->rmctx) {
291:           MatNullSpaceSetFunction(sub_nullspace, nullspace->remove, context);
292:         } else if (nullspace->remove && nullspace->rmctx) {
293:           char dmcoarse_method[PETSC_MAX_PATH_LEN];
294:           PetscErrorCode (*fp_get_coarsedm_context)(DM, void **) = NULL;

296:           PetscSNPrintf(dmcoarse_method, sizeof(dmcoarse_method), "PCTelescopeGetCoarseDMNullSpaceUserContext");
297:           PetscObjectQueryFunction((PetscObject)ctx->dm_coarse, dmcoarse_method, &fp_get_coarsedm_context);
299:           MatNullSpaceSetFunction(sub_nullspace, nullspace->remove, context);
300:         }
301:       }

303:       if (PCTelescope_isActiveRank(sred)) {
304:         MatSetNullSpace(sub_mat, sub_nullspace);
305:         MatNullSpaceDestroy(&sub_nullspace);
306:       }
307:     }
308:   }
309:   {
310:     MatNullSpace nearnullspace, sub_nearnullspace;
311:     MatGetNearNullSpace(B, &nearnullspace);
312:     if (nearnullspace) {
313:       PetscInfo(pc, "PCTelescope: generating near nullspace (CoarseDM)\n");
314:       PCTelescopeSubNullSpaceCreate_CoarseDM(pc, sred, nearnullspace, &sub_nearnullspace);

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

325:           PetscSNPrintf(dmcoarse_method, sizeof(dmcoarse_method), "PCTelescopeGetCoarseDMNearNullSpaceUserContext");
326:           PetscObjectQueryFunction((PetscObject)ctx->dm_coarse, dmcoarse_method, &fp_get_coarsedm_context);
328:           MatNullSpaceSetFunction(sub_nearnullspace, nearnullspace->remove, context);
329:         }
330:       }

332:       if (PCTelescope_isActiveRank(sred)) {
333:         MatSetNearNullSpace(sub_mat, sub_nearnullspace);
334:         MatNullSpaceDestroy(&sub_nearnullspace);
335:       }
336:     }
337:   }
338:   return 0;
339: }

341: PetscErrorCode PCReset_Telescope_CoarseDM(PC pc)
342: {
343:   PC_Telescope              sred = (PC_Telescope)pc->data;
344:   PC_Telescope_CoarseDMCtx *ctx;

346:   ctx              = (PC_Telescope_CoarseDMCtx *)sred->dm_ctx;
347:   ctx->dm_fine     = NULL; /* since I did not increment the ref counter we set these to NULL */
348:   ctx->dm_coarse   = NULL; /* since I did not increment the ref counter we set these to NULL */
349:   ctx->permutation = NULL; /* this will be fetched from the dm so no need to call destroy */
350:   VecDestroy(&ctx->xp);
351:   ctx->fp_dm_field_scatter      = NULL;
352:   ctx->fp_dm_state_scatter      = NULL;
353:   ctx->dmksp_context_determined = NULL;
354:   ctx->dmksp_context_user       = NULL;
355:   return 0;
356: }

358: 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)
359: {
360:   PC_Telescope              sred                     = (PC_Telescope)pc->data;
361:   Vec                       yred                     = NULL;
362:   PetscBool                 default_init_guess_value = PETSC_FALSE;
363:   PC_Telescope_CoarseDMCtx *ctx;

365:   ctx  = (PC_Telescope_CoarseDMCtx *)sred->dm_ctx;
366:   yred = sred->yred;

369:   *reason = (PCRichardsonConvergedReason)0;

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

374:     ctx->fp_dm_field_scatter(ctx->dm_fine, y, SCATTER_FORWARD, ctx->dm_coarse, yred);
375:   }

377:   if (PCTelescope_isActiveRank(sred)) {
378:     KSPGetInitialGuessNonzero(sred->ksp, &default_init_guess_value);
379:     if (!zeroguess) KSPSetInitialGuessNonzero(sred->ksp, PETSC_TRUE);
380:   }

382:   PCApply_Telescope_CoarseDM(pc, x, y);

384:   if (PCTelescope_isActiveRank(sred)) KSPSetInitialGuessNonzero(sred->ksp, default_init_guess_value);

386:   if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
387:   *outits = 1;
388:   return 0;
389: }