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: }