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[] =
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;
43: PetscErrorCode PCTelescopeSetUp_scatters_CoarseDM(PC pc,PC_Telescope sred,PC_Telescope_CoarseDMCtx *ctx)
44: {
45: Vec xred,yred,xtmp,x,xp;
46: VecScatter scatter;
47: IS isin;
48: Mat B;
49: PetscInt m,bs,st,ed;
50: MPI_Comm comm;
52: PetscObjectGetComm((PetscObject)pc,&comm);
53: PCGetOperators(pc,NULL,&B);
54: MatCreateVecs(B,&x,NULL);
55: MatGetBlockSize(B,&bs);
56: VecDuplicate(x,&xp);
57: m = 0;
58: xred = NULL;
59: yred = NULL;
60: if (PCTelescope_isActiveRank(sred)) {
61: DMCreateGlobalVector(ctx->dm_coarse,&xred);
62: VecDuplicate(xred,&yred);
63: VecGetOwnershipRange(xred,&st,&ed);
64: ISCreateStride(comm,ed-st,st,1,&isin);
65: VecGetLocalSize(xred,&m);
66: } else {
67: VecGetOwnershipRange(x,&st,&ed);
68: ISCreateStride(comm,0,st,1,&isin);
69: }
70: ISSetBlockSize(isin,bs);
71: VecCreate(comm,&xtmp);
72: VecSetSizes(xtmp,m,PETSC_DECIDE);
73: VecSetBlockSize(xtmp,bs);
74: VecSetType(xtmp,((PetscObject)x)->type_name);
75: VecScatterCreate(x,isin,xtmp,NULL,&scatter);
76: sred->xred = xred;
77: sred->yred = yred;
78: sred->isin = isin;
79: sred->scatter = scatter;
80: sred->xtmp = xtmp;
81: ctx->xp = xp;
82: VecDestroy(&x);
83: return 0;
84: }
86: PetscErrorCode PCTelescopeSetUp_CoarseDM(PC pc,PC_Telescope sred)
87: {
88: PC_Telescope_CoarseDMCtx *ctx;
89: DM dm,dm_coarse = NULL;
90: MPI_Comm comm;
91: PetscBool has_perm,has_kspcomputeoperators,using_kspcomputeoperators;
93: PetscInfo(pc,"PCTelescope: setup (CoarseDM)\n");
94: PetscNew(&ctx);
95: sred->dm_ctx = (void*)ctx;
97: PetscObjectGetComm((PetscObject)pc,&comm);
98: PCGetDM(pc,&dm);
99: DMGetCoarseDM(dm,&dm_coarse);
100: ctx->dm_fine = dm;
101: ctx->dm_coarse = dm_coarse;
103: /* attach coarse dm to ksp on sub communicator */
104: if (PCTelescope_isActiveRank(sred)) {
105: KSPSetDM(sred->ksp,ctx->dm_coarse);
106: if (sred->ignore_kspcomputeoperators) {
107: KSPSetDMActive(sred->ksp,PETSC_FALSE);
108: }
109: }
111: /* check if there is a method to provide a permutation */
112: has_perm = PETSC_FALSE;
113: has_kspcomputeoperators = PETSC_FALSE;
114: using_kspcomputeoperators = PETSC_FALSE;
116: /* if no permutation is provided, we must rely on KSPSetComputeOperators */
117: {
118: PetscErrorCode (*dmfine_kspfunc)(KSP,Mat,Mat,void*) = NULL;
119: void *dmfine_kspctx = NULL,*dmcoarse_kspctx = NULL;
120: void *dmfine_appctx = NULL,*dmcoarse_appctx = NULL;
121: void *dmfine_shellctx = NULL,*dmcoarse_shellctx = NULL;
123: DMKSPGetComputeOperators(dm,&dmfine_kspfunc,&dmfine_kspctx);
124: if (dmfine_kspfunc) { has_kspcomputeoperators = PETSC_TRUE; }
126: DMGetApplicationContext(ctx->dm_fine,&dmfine_appctx);
127: DMShellGetContext(ctx->dm_fine,&dmfine_shellctx);
129: /* need to define dmcoarse_kspctx */
130: if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) {
132: PetscInfo(pc,"PCTelescope: KSPSetComputeOperators fetched from parent DM\n");
133: if (PCTelescope_isActiveRank(sred)) {
134: DMGetApplicationContext(ctx->dm_coarse,&dmcoarse_appctx);
135: DMShellGetContext(ctx->dm_coarse,&dmcoarse_shellctx);
136: }
138: /* Assume that if the fine operator didn't require any context, neither will the coarse */
139: if (!dmfine_kspctx) {
140: dmcoarse_kspctx = NULL;
141: PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using NULL context\n");
142: } else {
144: PetscInfo(pc,"PCTelescope: KSPSetComputeOperators detected non-NULL context from parent DM \n");
145: if (PCTelescope_isActiveRank(sred)) {
147: if (dmfine_kspctx == dmfine_appctx) {
148: dmcoarse_kspctx = dmcoarse_appctx;
149: PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using context from DM->ApplicationContext\n");
151: } else if (dmfine_kspctx == dmfine_shellctx) {
152: dmcoarse_kspctx = dmcoarse_shellctx;
153: PetscInfo(pc,"PCTelescope: KSPSetComputeOperators using context from DMShell->Context\n");
155: }
156: ctx->dmksp_context_determined = dmcoarse_kspctx;
158: /* look for user provided method to fetch the context */
159: {
160: PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL;
161: void *dmcoarse_context_user = NULL;
162: char dmcoarse_method[PETSC_MAX_PATH_LEN];
164: PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMKSPContext");
165: PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context);
166: if (fp_get_coarsedm_context) {
167: PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n");
168: fp_get_coarsedm_context(ctx->dm_coarse,&dmcoarse_context_user);
169: ctx->dmksp_context_user = dmcoarse_context_user;
170: dmcoarse_kspctx = dmcoarse_context_user;
171: } else {
172: PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeGetCoarseDMKSPContext from coarse DM\n");
173: }
174: }
176: if (!dmcoarse_kspctx) {
177: PetscInfo(pc,"PCTelescope: KSPSetComputeOperators failed to determine the context to use on sub-communicator\n");
178: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot determine which context with use for KSPSetComputeOperators() on sub-communicator");
179: }
180: }
181: }
182: }
184: if (dmfine_kspfunc && !sred->ignore_kspcomputeoperators) {
185: using_kspcomputeoperators = PETSC_TRUE;
187: if (PCTelescope_isActiveRank(sred)) {
188: /* sub ksp inherits dmksp_func and context provided by user */
189: KSPSetComputeOperators(sred->ksp,dmfine_kspfunc,dmcoarse_kspctx);
190: /* PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)ctx->dmrepart); */
191: KSPSetDMActive(sred->ksp,PETSC_TRUE);
192: }
193: }
194: }
199: {
200: char dmfine_method[PETSC_MAX_PATH_LEN];
202: PetscSNPrintf(dmfine_method,sizeof(dmfine_method),"PCTelescopeFieldScatter");
203: PetscObjectQueryFunction((PetscObject)ctx->dm_fine,dmfine_method,&ctx->fp_dm_field_scatter);
205: PetscSNPrintf(dmfine_method,sizeof(dmfine_method),"PCTelescopeStateScatter");
206: PetscObjectQueryFunction((PetscObject)ctx->dm_fine,dmfine_method,&ctx->fp_dm_state_scatter);
207: }
209: if (ctx->fp_dm_state_scatter) {
210: PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeStateScatter from parent DM\n");
211: } else {
212: PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeStateScatter from parent DM\n");
213: }
215: if (ctx->fp_dm_field_scatter) {
216: PetscInfo(pc,"PCTelescope: Found composed method PCTelescopeFieldScatter from parent DM\n");
217: } else {
218: PetscInfo(pc,"PCTelescope: Failed to find composed method PCTelescopeFieldScatter from parent DM\n");
219: 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");
220: }
222: /* PCTelescopeSetUp_permutation_CoarseDM(pc,sred,ctx); */
223: PCTelescopeSetUp_scatters_CoarseDM(pc,sred,ctx);
224: return 0;
225: }
227: PetscErrorCode PCApply_Telescope_CoarseDM(PC pc,Vec x,Vec y)
228: {
229: PC_Telescope sred = (PC_Telescope)pc->data;
230: Vec xred,yred;
231: PC_Telescope_CoarseDMCtx *ctx;
233: ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
234: xred = sred->xred;
235: yred = sred->yred;
237: PetscCitationsRegister(citation,&cited);
239: if (ctx->fp_dm_state_scatter) {
240: ctx->fp_dm_state_scatter(ctx->dm_fine,SCATTER_FORWARD,ctx->dm_coarse);
241: }
243: ctx->fp_dm_field_scatter(ctx->dm_fine,x,SCATTER_FORWARD,ctx->dm_coarse,xred);
245: /* solve */
246: if (PCTelescope_isActiveRank(sred)) {
247: KSPSolve(sred->ksp,xred,yred);
248: }
250: ctx->fp_dm_field_scatter(ctx->dm_fine,y,SCATTER_REVERSE,ctx->dm_coarse,yred);
251: return 0;
252: }
254: PetscErrorCode PCTelescopeSubNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,MatNullSpace nullspace,MatNullSpace *sub_nullspace)
255: {
256: PetscBool has_const;
257: PetscInt k,n = 0;
258: const Vec *vecs;
259: Vec *sub_vecs = NULL;
260: MPI_Comm subcomm;
261: PC_Telescope_CoarseDMCtx *ctx;
263: ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
264: subcomm = sred->subcomm;
265: MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs);
267: if (PCTelescope_isActiveRank(sred)) {
268: /* create new vectors */
269: if (n) {
270: VecDuplicateVecs(sred->xred,n,&sub_vecs);
271: }
272: }
274: /* copy entries */
275: for (k=0; k<n; k++) {
276: ctx->fp_dm_field_scatter(ctx->dm_fine,vecs[k],SCATTER_FORWARD,ctx->dm_coarse,sub_vecs[k]);
277: }
279: if (PCTelescope_isActiveRank(sred)) {
280: /* create new (near) nullspace for redundant object */
281: MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,sub_nullspace);
282: VecDestroyVecs(n,&sub_vecs);
283: }
284: return 0;
285: }
287: PetscErrorCode PCTelescopeMatNullSpaceCreate_CoarseDM(PC pc,PC_Telescope sred,Mat sub_mat)
288: {
289: Mat B;
290: PC_Telescope_CoarseDMCtx *ctx;
292: ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
293: PCGetOperators(pc,NULL,&B);
294: {
295: MatNullSpace nullspace,sub_nullspace;
296: MatGetNullSpace(B,&nullspace);
297: if (nullspace) {
298: PetscInfo(pc,"PCTelescope: generating nullspace (CoarseDM)\n");
299: PCTelescopeSubNullSpaceCreate_CoarseDM(pc,sred,nullspace,&sub_nullspace);
301: /* attach any user nullspace removal methods and contexts */
302: if (PCTelescope_isActiveRank(sred)) {
303: void *context = NULL;
304: if (nullspace->remove && !nullspace->rmctx) {
305: MatNullSpaceSetFunction(sub_nullspace,nullspace->remove,context);
306: } else if (nullspace->remove && nullspace->rmctx) {
307: char dmcoarse_method[PETSC_MAX_PATH_LEN];
308: PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL;
310: PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMNullSpaceUserContext");
311: PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context);
313: MatNullSpaceSetFunction(sub_nullspace,nullspace->remove,context);
314: }
315: }
317: if (PCTelescope_isActiveRank(sred)) {
318: MatSetNullSpace(sub_mat,sub_nullspace);
319: MatNullSpaceDestroy(&sub_nullspace);
320: }
321: }
322: }
323: {
324: MatNullSpace nearnullspace,sub_nearnullspace;
325: MatGetNearNullSpace(B,&nearnullspace);
326: if (nearnullspace) {
327: PetscInfo(pc,"PCTelescope: generating near nullspace (CoarseDM)\n");
328: PCTelescopeSubNullSpaceCreate_CoarseDM(pc,sred,nearnullspace,&sub_nearnullspace);
330: /* attach any user nullspace removal methods and contexts */
331: if (PCTelescope_isActiveRank(sred)) {
332: void *context = NULL;
333: if (nearnullspace->remove && !nearnullspace->rmctx) {
334: MatNullSpaceSetFunction(sub_nearnullspace,nearnullspace->remove,context);
335: } else if (nearnullspace->remove && nearnullspace->rmctx) {
336: char dmcoarse_method[PETSC_MAX_PATH_LEN];
337: PetscErrorCode (*fp_get_coarsedm_context)(DM,void**) = NULL;
339: PetscSNPrintf(dmcoarse_method,sizeof(dmcoarse_method),"PCTelescopeGetCoarseDMNearNullSpaceUserContext");
340: PetscObjectQueryFunction((PetscObject)ctx->dm_coarse,dmcoarse_method,&fp_get_coarsedm_context);
342: MatNullSpaceSetFunction(sub_nearnullspace,nearnullspace->remove,context);
343: }
344: }
346: if (PCTelescope_isActiveRank(sred)) {
347: MatSetNearNullSpace(sub_mat,sub_nearnullspace);
348: MatNullSpaceDestroy(&sub_nearnullspace);
349: }
350: }
351: }
352: return 0;
353: }
355: PetscErrorCode PCReset_Telescope_CoarseDM(PC pc)
356: {
357: PC_Telescope sred = (PC_Telescope)pc->data;
358: PC_Telescope_CoarseDMCtx *ctx;
360: ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
361: ctx->dm_fine = NULL; /* since I did not increment the ref counter we set these to NULL */
362: ctx->dm_coarse = NULL; /* since I did not increment the ref counter we set these to NULL */
363: ctx->permutation = NULL; /* this will be fetched from the dm so no need to call destroy */
364: VecDestroy(&ctx->xp);
365: ctx->fp_dm_field_scatter = NULL;
366: ctx->fp_dm_state_scatter = NULL;
367: ctx->dmksp_context_determined = NULL;
368: ctx->dmksp_context_user = NULL;
369: return 0;
370: }
372: 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)
373: {
374: PC_Telescope sred = (PC_Telescope)pc->data;
375: Vec yred = NULL;
376: PetscBool default_init_guess_value = PETSC_FALSE;
377: PC_Telescope_CoarseDMCtx *ctx;
379: ctx = (PC_Telescope_CoarseDMCtx*)sred->dm_ctx;
380: yred = sred->yred;
383: *reason = (PCRichardsonConvergedReason)0;
385: if (!zeroguess) {
386: PetscInfo(pc,"PCTelescopeCoarseDM: Scattering y for non-zero-initial guess\n");
388: ctx->fp_dm_field_scatter(ctx->dm_fine,y,SCATTER_FORWARD,ctx->dm_coarse,yred);
389: }
391: if (PCTelescope_isActiveRank(sred)) {
392: KSPGetInitialGuessNonzero(sred->ksp,&default_init_guess_value);
393: if (!zeroguess) KSPSetInitialGuessNonzero(sred->ksp,PETSC_TRUE);
394: }
396: PCApply_Telescope_CoarseDM(pc,x,y);
398: if (PCTelescope_isActiveRank(sred)) {
399: KSPSetInitialGuessNonzero(sred->ksp,default_init_guess_value);
400: }
402: if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
403: *outits = 1;
404: return 0;
405: }