1 #include <petsc/private/matimpl.h>
2 #include <petsc/private/pcimpl.h>
3 #include <petsc/private/dmimpl.h>
4 #include <petscksp.h> /*I "petscksp.h" I*/
5 #include <petscdm.h>
6 #include <petscdmda.h>
7 #include <petscdmshell.h>
8
9 #include "../src/ksp/pc/impls/telescope/telescope.h"
10
11 static PetscBool cited = PETSC_FALSE;
12 static const char citation[] = "@inproceedings{MaySananRuppKnepleySmith2016,\n"
13 " title = {Extreme-Scale Multigrid Components within PETSc},\n"
14 " author = {Dave A. May and Patrick Sanan and Karl Rupp and Matthew G. Knepley and Barry F. Smith},\n"
15 " booktitle = {Proceedings of the Platform for Advanced Scientific Computing Conference},\n"
16 " series = {PASC '16},\n"
17 " isbn = {978-1-4503-4126-4},\n"
18 " location = {Lausanne, Switzerland},\n"
19 " pages = {5:1--5:12},\n"
20 " articleno = {5},\n"
21 " numpages = {12},\n"
22 " url = {https://doi.acm.org/10.1145/2929908.2929913},\n"
23 " doi = {10.1145/2929908.2929913},\n"
24 " acmid = {2929913},\n"
25 " publisher = {ACM},\n"
26 " address = {New York, NY, USA},\n"
27 " keywords = {GPU, HPC, agglomeration, coarse-level solver, multigrid, parallel computing, preconditioning},\n"
28 " year = {2016}\n"
29 "}\n";
30
31 typedef struct {
32 DM dm_fine, dm_coarse; /* these DM's should be topologically identical but use different communicators */
33 Mat permutation;
34 Vec xp;
35 PetscErrorCode (*fp_dm_field_scatter)(DM, Vec, ScatterMode, DM, Vec);
36 PetscErrorCode (*fp_dm_state_scatter)(DM, ScatterMode, DM);
37 void *dmksp_context_determined;
38 void *dmksp_context_user;
39 } PC_Telescope_CoarseDMCtx;
40
PCTelescopeSetUp_scatters_CoarseDM(PC pc,PC_Telescope sred,PC_Telescope_CoarseDMCtx * ctx)41 static PetscErrorCode PCTelescopeSetUp_scatters_CoarseDM(PC pc, PC_Telescope sred, PC_Telescope_CoarseDMCtx *ctx)
42 {
43 Vec xred, yred, xtmp, x, xp;
44 VecScatter scatter;
45 IS isin;
46 Mat B;
47 PetscInt m, bs, st, ed;
48 MPI_Comm comm;
49
50 PetscFunctionBegin;
51 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
52 PetscCall(PCGetOperators(pc, NULL, &B));
53 PetscCall(MatCreateVecs(B, &x, NULL));
54 PetscCall(MatGetBlockSize(B, &bs));
55 PetscCall(VecDuplicate(x, &xp));
56 m = 0;
57 xred = NULL;
58 yred = NULL;
59 if (PCTelescope_isActiveRank(sred)) {
60 PetscCall(DMCreateGlobalVector(ctx->dm_coarse, &xred));
61 PetscCall(VecDuplicate(xred, &yred));
62 PetscCall(VecGetOwnershipRange(xred, &st, &ed));
63 PetscCall(ISCreateStride(comm, ed - st, st, 1, &isin));
64 PetscCall(VecGetLocalSize(xred, &m));
65 } else {
66 PetscCall(VecGetOwnershipRange(x, &st, &ed));
67 PetscCall(ISCreateStride(comm, 0, st, 1, &isin));
68 }
69 PetscCall(ISSetBlockSize(isin, bs));
70 PetscCall(VecCreate(comm, &xtmp));
71 PetscCall(VecSetSizes(xtmp, m, PETSC_DECIDE));
72 PetscCall(VecSetBlockSize(xtmp, bs));
73 PetscCall(VecSetType(xtmp, ((PetscObject)x)->type_name));
74 PetscCall(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 PetscCall(VecDestroy(&x));
82 PetscFunctionReturn(PETSC_SUCCESS);
83 }
84
PCTelescopeSetUp_CoarseDM(PC pc,PC_Telescope sred)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;
91
92 PetscFunctionBegin;
93 PetscCall(PetscInfo(pc, "PCTelescope: setup (CoarseDM)\n"));
94 PetscCall(PetscNew(&ctx));
95 sred->dm_ctx = (void *)ctx;
96
97 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
98 PetscCall(PCGetDM(pc, &dm));
99 PetscCall(DMGetCoarseDM(dm, &dm_coarse));
100 ctx->dm_fine = dm;
101 ctx->dm_coarse = dm_coarse;
102
103 /* attach coarse dm to ksp on sub communicator */
104 if (PCTelescope_isActiveRank(sred)) {
105 PetscCall(KSPSetDM(sred->ksp, ctx->dm_coarse));
106 if (sred->ignore_kspcomputeoperators) PetscCall(KSPSetDMActive(sred->ksp, KSP_DMACTIVE_ALL, PETSC_FALSE));
107 }
108
109 /* check if there is a method to provide a permutation */
110 has_perm = PETSC_FALSE;
111 has_kspcomputeoperators = PETSC_FALSE;
112 using_kspcomputeoperators = PETSC_FALSE;
113
114 /* if no permutation is provided, we must rely on KSPSetComputeOperators */
115 {
116 PetscErrorCode (*dmfine_kspfunc)(KSP, Mat, Mat, void *) = NULL;
117 void *dmfine_kspctx = NULL, *dmcoarse_kspctx = NULL;
118 void *dmfine_appctx = NULL, *dmcoarse_appctx = NULL;
119 void *dmfine_shellctx = NULL, *dmcoarse_shellctx = NULL;
120
121 PetscCall(DMKSPGetComputeOperators(dm, &dmfine_kspfunc, &dmfine_kspctx));
122 if (dmfine_kspfunc) has_kspcomputeoperators = PETSC_TRUE;
123
124 PetscCall(DMGetApplicationContext(ctx->dm_fine, &dmfine_appctx));
125 PetscCall(DMShellGetContext(ctx->dm_fine, &dmfine_shellctx));
126
127 /* need to define dmcoarse_kspctx */
128 if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) {
129 PetscCall(PetscInfo(pc, "PCTelescope: KSPSetComputeOperators fetched from parent DM\n"));
130 if (PCTelescope_isActiveRank(sred)) {
131 PetscCall(DMGetApplicationContext(ctx->dm_coarse, &dmcoarse_appctx));
132 PetscCall(DMShellGetContext(ctx->dm_coarse, &dmcoarse_shellctx));
133 }
134
135 /* Assume that if the fine operator didn't require any context, neither will the coarse */
136 if (!dmfine_kspctx) {
137 dmcoarse_kspctx = NULL;
138 PetscCall(PetscInfo(pc, "PCTelescope: KSPSetComputeOperators using NULL context\n"));
139 } else {
140 PetscCall(PetscInfo(pc, "PCTelescope: KSPSetComputeOperators detected non-NULL context from parent DM \n"));
141 if (PCTelescope_isActiveRank(sred)) {
142 if (dmfine_kspctx == dmfine_appctx) {
143 dmcoarse_kspctx = dmcoarse_appctx;
144 PetscCall(PetscInfo(pc, "PCTelescope: KSPSetComputeOperators using context from DM->ApplicationContext\n"));
145 PetscCheck(dmcoarse_kspctx, PETSC_COMM_SELF, PETSC_ERR_USER, "Non NULL dmfine->kspctx == dmfine->appctx. NULL dmcoarse->appctx found. Likely this is an error");
146 } else if (dmfine_kspctx == dmfine_shellctx) {
147 dmcoarse_kspctx = dmcoarse_shellctx;
148 PetscCall(PetscInfo(pc, "PCTelescope: KSPSetComputeOperators using context from DMShell->Context\n"));
149 PetscCheck(dmcoarse_kspctx, PETSC_COMM_SELF, PETSC_ERR_USER, "Non NULL dmfine->kspctx == dmfine.shell->ctx. NULL dmcoarse.shell->ctx found. Likely this is an error");
150 }
151 ctx->dmksp_context_determined = dmcoarse_kspctx;
152
153 /* look for user provided method to fetch the context */
154 {
155 PetscErrorCode (*fp_get_coarsedm_context)(DM, void **) = NULL;
156 void *dmcoarse_context_user = NULL;
157 char dmcoarse_method[PETSC_MAX_PATH_LEN];
158
159 PetscCall(PetscSNPrintf(dmcoarse_method, sizeof(dmcoarse_method), "PCTelescopeGetCoarseDMKSPContext"));
160 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_coarse, dmcoarse_method, &fp_get_coarsedm_context));
161 if (fp_get_coarsedm_context) {
162 PetscCall(PetscInfo(pc, "PCTelescope: Found composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n"));
163 PetscCall(fp_get_coarsedm_context(ctx->dm_coarse, &dmcoarse_context_user));
164 ctx->dmksp_context_user = dmcoarse_context_user;
165 dmcoarse_kspctx = dmcoarse_context_user;
166 } else {
167 PetscCall(PetscInfo(pc, "PCTelescope: Failed to find composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n"));
168 }
169 }
170
171 if (!dmcoarse_kspctx) {
172 PetscCall(PetscInfo(pc, "PCTelescope: KSPSetComputeOperators failed to determine the context to use on sub-communicator\n"));
173 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot determine which context with use for KSPSetComputeOperators() on sub-communicator");
174 }
175 }
176 }
177 }
178
179 if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) {
180 using_kspcomputeoperators = PETSC_TRUE;
181
182 if (PCTelescope_isActiveRank(sred)) {
183 /* sub ksp inherits dmksp_func and context provided by user */
184 PetscCall(KSPSetComputeOperators(sred->ksp, dmfine_kspfunc, dmcoarse_kspctx));
185 /* PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)ctx->dmrepart)); */
186 PetscCall(KSPSetDMActive(sred->ksp, KSP_DMACTIVE_ALL, PETSC_TRUE));
187 }
188 }
189 }
190
191 PetscCheck(has_perm || !has_kspcomputeoperators || using_kspcomputeoperators, 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");
192 PetscCheck(has_perm || has_kspcomputeoperators, 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");
193
194 {
195 char dmfine_method[PETSC_MAX_PATH_LEN];
196
197 PetscCall(PetscSNPrintf(dmfine_method, sizeof(dmfine_method), "PCTelescopeFieldScatter"));
198 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_fine, dmfine_method, &ctx->fp_dm_field_scatter));
199
200 PetscCall(PetscSNPrintf(dmfine_method, sizeof(dmfine_method), "PCTelescopeStateScatter"));
201 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_fine, dmfine_method, &ctx->fp_dm_state_scatter));
202 }
203
204 if (ctx->fp_dm_state_scatter) {
205 PetscCall(PetscInfo(pc, "PCTelescope: Found composed method PCTelescopeStateScatter from parent DM\n"));
206 } else {
207 PetscCall(PetscInfo(pc, "PCTelescope: Failed to find composed method PCTelescopeStateScatter from parent DM\n"));
208 }
209
210 if (ctx->fp_dm_field_scatter) {
211 PetscCall(PetscInfo(pc, "PCTelescope: Found composed method PCTelescopeFieldScatter from parent DM\n"));
212 } else {
213 PetscCall(PetscInfo(pc, "PCTelescope: Failed to find composed method PCTelescopeFieldScatter from parent DM\n"));
214 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");
215 }
216
217 /* PetscCall(PCTelescopeSetUp_permutation_CoarseDM(pc,sred,ctx)); */
218 PetscCall(PCTelescopeSetUp_scatters_CoarseDM(pc, sred, ctx));
219 PetscFunctionReturn(PETSC_SUCCESS);
220 }
221
PCApply_Telescope_CoarseDM(PC pc,Vec x,Vec y)222 PetscErrorCode PCApply_Telescope_CoarseDM(PC pc, Vec x, Vec y)
223 {
224 PC_Telescope sred = (PC_Telescope)pc->data;
225 Vec xred, yred;
226 PC_Telescope_CoarseDMCtx *ctx;
227
228 PetscFunctionBegin;
229 ctx = (PC_Telescope_CoarseDMCtx *)sred->dm_ctx;
230 xred = sred->xred;
231 yred = sred->yred;
232
233 PetscCall(PetscCitationsRegister(citation, &cited));
234
235 if (ctx->fp_dm_state_scatter) PetscCall(ctx->fp_dm_state_scatter(ctx->dm_fine, SCATTER_FORWARD, ctx->dm_coarse));
236
237 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine, x, SCATTER_FORWARD, ctx->dm_coarse, xred));
238
239 /* solve */
240 if (PCTelescope_isActiveRank(sred)) PetscCall(KSPSolve(sred->ksp, xred, yred));
241
242 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine, y, SCATTER_REVERSE, ctx->dm_coarse, yred));
243 PetscFunctionReturn(PETSC_SUCCESS);
244 }
245
PCTelescopeSubNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,MatNullSpace nullspace,MatNullSpace * sub_nullspace)246 static PetscErrorCode PCTelescopeSubNullSpaceCreate_CoarseDM(PC pc, PC_Telescope sred, MatNullSpace nullspace, MatNullSpace *sub_nullspace)
247 {
248 PetscBool has_const;
249 PetscInt k, n = 0;
250 const Vec *vecs;
251 Vec *sub_vecs = NULL;
252 MPI_Comm subcomm;
253 PC_Telescope_CoarseDMCtx *ctx;
254
255 PetscFunctionBegin;
256 ctx = (PC_Telescope_CoarseDMCtx *)sred->dm_ctx;
257 subcomm = sred->subcomm;
258 PetscCall(MatNullSpaceGetVecs(nullspace, &has_const, &n, &vecs));
259
260 if (PCTelescope_isActiveRank(sred)) {
261 /* create new vectors */
262 if (n) PetscCall(VecDuplicateVecs(sred->xred, n, &sub_vecs));
263 }
264
265 /* copy entries */
266 for (k = 0; k < n; k++) PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine, vecs[k], SCATTER_FORWARD, ctx->dm_coarse, sub_vecs[k]));
267
268 if (PCTelescope_isActiveRank(sred)) {
269 /* create new (near) nullspace for redundant object */
270 PetscCall(MatNullSpaceCreate(subcomm, has_const, n, sub_vecs, sub_nullspace));
271 PetscCall(VecDestroyVecs(n, &sub_vecs));
272 }
273 PetscFunctionReturn(PETSC_SUCCESS);
274 }
275
PCTelescopeMatNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,Mat sub_mat)276 PetscErrorCode PCTelescopeMatNullSpaceCreate_CoarseDM(PC pc, PC_Telescope sred, Mat sub_mat)
277 {
278 Mat B;
279 PC_Telescope_CoarseDMCtx *ctx;
280
281 PetscFunctionBegin;
282 ctx = (PC_Telescope_CoarseDMCtx *)sred->dm_ctx;
283 PetscCall(PCGetOperators(pc, NULL, &B));
284 {
285 MatNullSpace nullspace, sub_nullspace;
286 PetscCall(MatGetNullSpace(B, &nullspace));
287 if (nullspace) {
288 PetscCall(PetscInfo(pc, "PCTelescope: generating nullspace (CoarseDM)\n"));
289 PetscCall(PCTelescopeSubNullSpaceCreate_CoarseDM(pc, sred, nullspace, &sub_nullspace));
290
291 /* attach any user nullspace removal methods and contexts */
292 if (PCTelescope_isActiveRank(sred)) {
293 void *context = NULL;
294 if (nullspace->remove && !nullspace->rmctx) {
295 PetscCall(MatNullSpaceSetFunction(sub_nullspace, nullspace->remove, context));
296 } else if (nullspace->remove && nullspace->rmctx) {
297 char dmcoarse_method[PETSC_MAX_PATH_LEN];
298 PetscErrorCode (*fp_get_coarsedm_context)(DM, void **) = NULL;
299
300 PetscCall(PetscSNPrintf(dmcoarse_method, sizeof(dmcoarse_method), "PCTelescopeGetCoarseDMNullSpaceUserContext"));
301 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_coarse, dmcoarse_method, &fp_get_coarsedm_context));
302 PetscCheck(context, 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);
303 PetscCall(MatNullSpaceSetFunction(sub_nullspace, nullspace->remove, context));
304 }
305 }
306
307 if (PCTelescope_isActiveRank(sred)) {
308 PetscCall(MatSetNullSpace(sub_mat, sub_nullspace));
309 PetscCall(MatNullSpaceDestroy(&sub_nullspace));
310 }
311 }
312 }
313 {
314 MatNullSpace nearnullspace, sub_nearnullspace;
315 PetscCall(MatGetNearNullSpace(B, &nearnullspace));
316 if (nearnullspace) {
317 PetscCall(PetscInfo(pc, "PCTelescope: generating near nullspace (CoarseDM)\n"));
318 PetscCall(PCTelescopeSubNullSpaceCreate_CoarseDM(pc, sred, nearnullspace, &sub_nearnullspace));
319
320 /* attach any user nullspace removal methods and contexts */
321 if (PCTelescope_isActiveRank(sred)) {
322 void *context = NULL;
323 if (nearnullspace->remove && !nearnullspace->rmctx) {
324 PetscCall(MatNullSpaceSetFunction(sub_nearnullspace, nearnullspace->remove, context));
325 } else if (nearnullspace->remove && nearnullspace->rmctx) {
326 char dmcoarse_method[PETSC_MAX_PATH_LEN];
327 PetscErrorCode (*fp_get_coarsedm_context)(DM, void **) = NULL;
328
329 PetscCall(PetscSNPrintf(dmcoarse_method, sizeof(dmcoarse_method), "PCTelescopeGetCoarseDMNearNullSpaceUserContext"));
330 PetscCall(PetscObjectQueryFunction((PetscObject)ctx->dm_coarse, dmcoarse_method, &fp_get_coarsedm_context));
331 PetscCheck(context, 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);
332 PetscCall(MatNullSpaceSetFunction(sub_nearnullspace, nearnullspace->remove, context));
333 }
334 }
335
336 if (PCTelescope_isActiveRank(sred)) {
337 PetscCall(MatSetNearNullSpace(sub_mat, sub_nearnullspace));
338 PetscCall(MatNullSpaceDestroy(&sub_nearnullspace));
339 }
340 }
341 }
342 PetscFunctionReturn(PETSC_SUCCESS);
343 }
344
PCReset_Telescope_CoarseDM(PC pc)345 PetscErrorCode PCReset_Telescope_CoarseDM(PC pc)
346 {
347 PC_Telescope sred = (PC_Telescope)pc->data;
348 PC_Telescope_CoarseDMCtx *ctx;
349
350 PetscFunctionBegin;
351 ctx = (PC_Telescope_CoarseDMCtx *)sred->dm_ctx;
352 ctx->dm_fine = NULL; /* since I did not increment the ref counter we set these to NULL */
353 ctx->dm_coarse = NULL; /* since I did not increment the ref counter we set these to NULL */
354 ctx->permutation = NULL; /* this will be fetched from the dm so no need to call destroy */
355 PetscCall(VecDestroy(&ctx->xp));
356 ctx->fp_dm_field_scatter = NULL;
357 ctx->fp_dm_state_scatter = NULL;
358 ctx->dmksp_context_determined = NULL;
359 ctx->dmksp_context_user = NULL;
360 PetscFunctionReturn(PETSC_SUCCESS);
361 }
362
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)363 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)
364 {
365 PC_Telescope sred = (PC_Telescope)pc->data;
366 Vec yred = NULL;
367 PetscBool default_init_guess_value = PETSC_FALSE;
368 PC_Telescope_CoarseDMCtx *ctx;
369
370 PetscFunctionBegin;
371 ctx = (PC_Telescope_CoarseDMCtx *)sred->dm_ctx;
372 yred = sred->yred;
373
374 PetscCheck(its <= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCApplyRichardson_Telescope_CoarseDM only supports max_it = 1");
375 *reason = (PCRichardsonConvergedReason)0;
376
377 if (!zeroguess) {
378 PetscCall(PetscInfo(pc, "PCTelescopeCoarseDM: Scattering y for non-zero-initial guess\n"));
379
380 PetscCall(ctx->fp_dm_field_scatter(ctx->dm_fine, y, SCATTER_FORWARD, ctx->dm_coarse, yred));
381 }
382
383 if (PCTelescope_isActiveRank(sred)) {
384 PetscCall(KSPGetInitialGuessNonzero(sred->ksp, &default_init_guess_value));
385 if (!zeroguess) PetscCall(KSPSetInitialGuessNonzero(sred->ksp, PETSC_TRUE));
386 }
387
388 PetscCall(PCApply_Telescope_CoarseDM(pc, x, y));
389
390 if (PCTelescope_isActiveRank(sred)) PetscCall(KSPSetInitialGuessNonzero(sred->ksp, default_init_guess_value));
391
392 if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
393 *outits = 1;
394 PetscFunctionReturn(PETSC_SUCCESS);
395 }
396