Actual source code: telescope.c
petsc-3.8.4 2018-03-24
2: #include <petsc/private/matimpl.h>
3: #include <petsc/private/pcimpl.h>
4: #include <petscksp.h>
5: #include <petscdm.h>
6: #include "../src/ksp/pc/impls/telescope/telescope.h"
8: static PetscBool cited = PETSC_FALSE;
9: static const char citation[] =
10: "@inproceedings{MaySananRuppKnepleySmith2016,\n"
11: " title = {Extreme-Scale Multigrid Components within PETSc},\n"
12: " author = {Dave A. May and Patrick Sanan and Karl Rupp and Matthew G. Knepley and Barry F. Smith},\n"
13: " booktitle = {Proceedings of the Platform for Advanced Scientific Computing Conference},\n"
14: " series = {PASC '16},\n"
15: " isbn = {978-1-4503-4126-4},\n"
16: " location = {Lausanne, Switzerland},\n"
17: " pages = {5:1--5:12},\n"
18: " articleno = {5},\n"
19: " numpages = {12},\n"
20: " url = {http://doi.acm.org/10.1145/2929908.2929913},\n"
21: " doi = {10.1145/2929908.2929913},\n"
22: " acmid = {2929913},\n"
23: " publisher = {ACM},\n"
24: " address = {New York, NY, USA},\n"
25: " keywords = {GPU, HPC, agglomeration, coarse-level solver, multigrid, parallel computing, preconditioning},\n"
26: " year = {2016}\n"
27: "}\n";
29: /*
30: PCTelescopeSetUp_default()
31: PCTelescopeMatCreate_default()
33: default
35: // scatter in
36: x(comm) -> xtmp(comm)
38: xred(subcomm) <- xtmp
39: yred(subcomm)
41: yred(subcomm) --> xtmp
43: // scatter out
44: xtmp(comm) -> y(comm)
45: */
47: PetscBool isActiveRank(PetscSubcomm scomm)
48: {
49: if (scomm->color == 0) { return PETSC_TRUE; }
50: else { return PETSC_FALSE; }
51: }
53: DM private_PCTelescopeGetSubDM(PC_Telescope sred)
54: {
55: DM subdm = NULL;
57: if (!isActiveRank(sred->psubcomm)) { subdm = NULL; }
58: else {
59: switch (sred->sr_type) {
60: case TELESCOPE_DEFAULT: subdm = NULL;
61: break;
62: case TELESCOPE_DMDA: subdm = ((PC_Telescope_DMDACtx*)sred->dm_ctx)->dmrepart;
63: break;
64: case TELESCOPE_DMPLEX: subdm = NULL;
65: break;
66: }
67: }
68: return(subdm);
69: }
71: PetscErrorCode PCTelescopeSetUp_default(PC pc,PC_Telescope sred)
72: {
74: PetscInt m,M,bs,st,ed;
75: Vec x,xred,yred,xtmp;
76: Mat B;
77: MPI_Comm comm,subcomm;
78: VecScatter scatter;
79: IS isin;
82: PetscInfo(pc,"PCTelescope: setup (default)\n");
83: comm = PetscSubcommParent(sred->psubcomm);
84: subcomm = PetscSubcommChild(sred->psubcomm);
86: PCGetOperators(pc,NULL,&B);
87: MatGetSize(B,&M,NULL);
88: MatGetBlockSize(B,&bs);
89: MatCreateVecs(B,&x,NULL);
91: xred = NULL;
92: m = 0;
93: if (isActiveRank(sred->psubcomm)) {
94: VecCreate(subcomm,&xred);
95: VecSetSizes(xred,PETSC_DECIDE,M);
96: VecSetBlockSize(xred,bs);
97: VecSetFromOptions(xred);
98: VecGetLocalSize(xred,&m);
99: }
101: yred = NULL;
102: if (isActiveRank(sred->psubcomm)) {
103: VecDuplicate(xred,&yred);
104: }
106: VecCreate(comm,&xtmp);
107: VecSetSizes(xtmp,m,PETSC_DECIDE);
108: VecSetBlockSize(xtmp,bs);
109: VecSetType(xtmp,((PetscObject)x)->type_name);
111: if (isActiveRank(sred->psubcomm)) {
112: VecGetOwnershipRange(xred,&st,&ed);
113: ISCreateStride(comm,(ed-st),st,1,&isin);
114: } else {
115: VecGetOwnershipRange(x,&st,&ed);
116: ISCreateStride(comm,0,st,1,&isin);
117: }
118: ISSetBlockSize(isin,bs);
120: VecScatterCreate(x,isin,xtmp,NULL,&scatter);
122: sred->isin = isin;
123: sred->scatter = scatter;
124: sred->xred = xred;
125: sred->yred = yred;
126: sred->xtmp = xtmp;
127: VecDestroy(&x);
128: return(0);
129: }
131: PetscErrorCode PCTelescopeMatCreate_default(PC pc,PC_Telescope sred,MatReuse reuse,Mat *A)
132: {
134: MPI_Comm comm,subcomm;
135: Mat Bred,B;
136: PetscInt nr,nc;
137: IS isrow,iscol;
138: Mat Blocal,*_Blocal;
141: PetscInfo(pc,"PCTelescope: updating the redundant preconditioned operator (default)\n");
142: PetscObjectGetComm((PetscObject)pc,&comm);
143: subcomm = PetscSubcommChild(sred->psubcomm);
144: PCGetOperators(pc,NULL,&B);
145: MatGetSize(B,&nr,&nc);
146: isrow = sred->isin;
147: ISCreateStride(comm,nc,0,1,&iscol);
148: MatCreateSubMatrices(B,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&_Blocal);
149: Blocal = *_Blocal;
150: PetscFree(_Blocal);
151: Bred = NULL;
152: if (isActiveRank(sred->psubcomm)) {
153: PetscInt mm;
155: if (reuse != MAT_INITIAL_MATRIX) { Bred = *A; }
157: MatGetSize(Blocal,&mm,NULL);
158: MatCreateMPIMatConcatenateSeqMat(subcomm,Blocal,mm,reuse,&Bred);
159: }
160: *A = Bred;
161: ISDestroy(&iscol);
162: MatDestroy(&Blocal);
163: return(0);
164: }
166: static PetscErrorCode PCTelescopeSubNullSpaceCreate_Telescope(PC pc,PC_Telescope sred,MatNullSpace nullspace,MatNullSpace *sub_nullspace)
167: {
168: PetscErrorCode ierr;
169: PetscBool has_const;
170: const Vec *vecs;
171: Vec *sub_vecs = NULL;
172: PetscInt i,k,n = 0;
173: MPI_Comm subcomm;
176: subcomm = PetscSubcommChild(sred->psubcomm);
177: MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs);
179: if (isActiveRank(sred->psubcomm)) {
180: if (n) {
181: VecDuplicateVecs(sred->xred,n,&sub_vecs);
182: }
183: }
185: /* copy entries */
186: for (k=0; k<n; k++) {
187: const PetscScalar *x_array;
188: PetscScalar *LA_sub_vec;
189: PetscInt st,ed;
191: /* pull in vector x->xtmp */
192: VecScatterBegin(sred->scatter,vecs[k],sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);
193: VecScatterEnd(sred->scatter,vecs[k],sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);
194: if (sub_vecs) {
195: /* copy vector entries into xred */
196: VecGetArrayRead(sred->xtmp,&x_array);
197: if (sub_vecs[k]) {
198: VecGetOwnershipRange(sub_vecs[k],&st,&ed);
199: VecGetArray(sub_vecs[k],&LA_sub_vec);
200: for (i=0; i<ed-st; i++) {
201: LA_sub_vec[i] = x_array[i];
202: }
203: VecRestoreArray(sub_vecs[k],&LA_sub_vec);
204: }
205: VecRestoreArrayRead(sred->xtmp,&x_array);
206: }
207: }
209: if (isActiveRank(sred->psubcomm)) {
210: /* create new (near) nullspace for redundant object */
211: MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,sub_nullspace);
212: VecDestroyVecs(n,&sub_vecs);
213: if (nullspace->remove) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Propagation of custom remove callbacks not supported when propagating (near) nullspaces with PCTelescope");
214: if (nullspace->rmctx) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Propagation of custom remove callback context not supported when propagating (near) nullspaces with PCTelescope");
215: }
216: return(0);
217: }
219: static PetscErrorCode PCTelescopeMatNullSpaceCreate_default(PC pc,PC_Telescope sred,Mat sub_mat)
220: {
221: PetscErrorCode ierr;
222: Mat B;
225: PCGetOperators(pc,NULL,&B);
227: /* Propagate the nullspace if it exists */
228: {
229: MatNullSpace nullspace,sub_nullspace;
230: MatGetNullSpace(B,&nullspace);
231: if (nullspace) {
232: PetscInfo(pc,"PCTelescope: generating nullspace (default)\n");
233: PCTelescopeSubNullSpaceCreate_Telescope(pc,sred,nullspace,&sub_nullspace);
234: if (isActiveRank(sred->psubcomm)) {
235: MatSetNullSpace(sub_mat,sub_nullspace);
236: MatNullSpaceDestroy(&sub_nullspace);
237: }
238: }
239: }
241: /* Propagate the near nullspace if it exists */
242: {
243: MatNullSpace nearnullspace,sub_nearnullspace;
244: MatGetNearNullSpace(B,&nearnullspace);
245: if (nearnullspace) {
246: PetscInfo(pc,"PCTelescope: generating near nullspace (default)\n");
247: PCTelescopeSubNullSpaceCreate_Telescope(pc,sred,nearnullspace,&sub_nearnullspace);
248: if (isActiveRank(sred->psubcomm)) {
249: MatSetNearNullSpace(sub_mat,sub_nearnullspace);
250: MatNullSpaceDestroy(&sub_nearnullspace);
251: }
252: }
253: }
254: return(0);
255: }
257: static PetscErrorCode PCView_Telescope(PC pc,PetscViewer viewer)
258: {
259: PC_Telescope sred = (PC_Telescope)pc->data;
260: PetscErrorCode ierr;
261: PetscBool iascii,isstring;
262: PetscViewer subviewer;
265: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
266: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
267: if (iascii) {
268: if (!sred->psubcomm) {
269: PetscViewerASCIIPrintf(viewer," preconditioner not yet setup\n");
270: } else {
271: MPI_Comm comm,subcomm;
272: PetscMPIInt comm_size,subcomm_size;
273: DM dm,subdm;
275: PCGetDM(pc,&dm);
276: subdm = private_PCTelescopeGetSubDM(sred);
277: comm = PetscSubcommParent(sred->psubcomm);
278: subcomm = PetscSubcommChild(sred->psubcomm);
279: MPI_Comm_size(comm,&comm_size);
280: MPI_Comm_size(subcomm,&subcomm_size);
282: PetscViewerASCIIPrintf(viewer," parent comm size reduction factor = %D\n",sred->redfactor);
283: PetscViewerASCIIPrintf(viewer," comm_size = %d , subcomm_size = %d\n",(int)comm_size,(int)subcomm_size);
284: switch (sred->subcommtype) {
285: case PETSC_SUBCOMM_INTERLACED :
286: PetscViewerASCIIPrintf(viewer," subcomm type: interlaced\n",sred->subcommtype);
287: break;
288: case PETSC_SUBCOMM_CONTIGUOUS :
289: PetscViewerASCIIPrintf(viewer," subcomm type: contiguous\n",sred->subcommtype);
290: break;
291: default :
292: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"General subcomm type not supported by PCTelescope");
293: }
294: PetscViewerGetSubViewer(viewer,subcomm,&subviewer);
295: if (isActiveRank(sred->psubcomm)) {
296: PetscViewerASCIIPushTab(subviewer);
298: if (dm && sred->ignore_dm) {
299: PetscViewerASCIIPrintf(subviewer," ignoring DM\n");
300: }
301: if (sred->ignore_kspcomputeoperators) {
302: PetscViewerASCIIPrintf(subviewer," ignoring KSPComputeOperators\n");
303: }
304: switch (sred->sr_type) {
305: case TELESCOPE_DEFAULT:
306: PetscViewerASCIIPrintf(subviewer," using default setup\n");
307: break;
308: case TELESCOPE_DMDA:
309: PetscViewerASCIIPrintf(subviewer," DMDA detected\n");
310: DMView_DMDAShort(subdm,subviewer);
311: break;
312: case TELESCOPE_DMPLEX:
313: PetscViewerASCIIPrintf(subviewer," DMPLEX detected\n");
314: break;
315: }
317: KSPView(sred->ksp,subviewer);
318: PetscViewerASCIIPopTab(subviewer);
319: }
320: PetscViewerRestoreSubViewer(viewer,subcomm,&subviewer);
321: }
322: }
323: return(0);
324: }
326: static PetscErrorCode PCSetUp_Telescope(PC pc)
327: {
328: PC_Telescope sred = (PC_Telescope)pc->data;
329: PetscErrorCode ierr;
330: MPI_Comm comm,subcomm=0;
331: PCTelescopeType sr_type;
334: PetscObjectGetComm((PetscObject)pc,&comm);
336: /* subcomm definition */
337: if (!pc->setupcalled) {
338: if (!sred->psubcomm) {
339: PetscSubcommCreate(comm,&sred->psubcomm);
340: PetscSubcommSetNumber(sred->psubcomm,sred->redfactor);
341: PetscSubcommSetType(sred->psubcomm,sred->subcommtype);
342: PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
343: }
344: }
345: subcomm = PetscSubcommChild(sred->psubcomm);
347: /* internal KSP */
348: if (!pc->setupcalled) {
349: const char *prefix;
351: if (isActiveRank(sred->psubcomm)) {
352: KSPCreate(subcomm,&sred->ksp);
353: KSPSetErrorIfNotConverged(sred->ksp,pc->erroriffailure);
354: PetscObjectIncrementTabLevel((PetscObject)sred->ksp,(PetscObject)pc,1);
355: PetscLogObjectParent((PetscObject)pc,(PetscObject)sred->ksp);
356: KSPSetType(sred->ksp,KSPPREONLY);
357: PCGetOptionsPrefix(pc,&prefix);
358: KSPSetOptionsPrefix(sred->ksp,prefix);
359: KSPAppendOptionsPrefix(sred->ksp,"telescope_");
360: }
361: }
362: /* Determine type of setup/update */
363: if (!pc->setupcalled) {
364: PetscBool has_dm,same;
365: DM dm;
367: sr_type = TELESCOPE_DEFAULT;
368: has_dm = PETSC_FALSE;
369: PCGetDM(pc,&dm);
370: if (dm) { has_dm = PETSC_TRUE; }
371: if (has_dm) {
372: /* check for dmda */
373: PetscObjectTypeCompare((PetscObject)dm,DMDA,&same);
374: if (same) {
375: PetscInfo(pc,"PCTelescope: found DMDA\n");
376: sr_type = TELESCOPE_DMDA;
377: }
378: /* check for dmplex */
379: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&same);
380: if (same) {
381: PetscInfo(pc,"PCTelescope: found DMPLEX\n");
382: sr_type = TELESCOPE_DMPLEX;
383: }
384: }
386: if (sred->ignore_dm) {
387: PetscInfo(pc,"PCTelescope: ignore DM\n");
388: sr_type = TELESCOPE_DEFAULT;
389: }
390: sred->sr_type = sr_type;
391: } else {
392: sr_type = sred->sr_type;
393: }
395: /* set function pointers for repartition setup, matrix creation/update, matrix (near) nullspace, and reset functionality */
396: switch (sr_type) {
397: case TELESCOPE_DEFAULT:
398: sred->pctelescope_setup_type = PCTelescopeSetUp_default;
399: sred->pctelescope_matcreate_type = PCTelescopeMatCreate_default;
400: sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_default;
401: sred->pctelescope_reset_type = NULL;
402: break;
403: case TELESCOPE_DMDA:
404: pc->ops->apply = PCApply_Telescope_dmda;
405: pc->ops->applyrichardson = PCApplyRichardson_Telescope_dmda;
406: sred->pctelescope_setup_type = PCTelescopeSetUp_dmda;
407: sred->pctelescope_matcreate_type = PCTelescopeMatCreate_dmda;
408: sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_dmda;
409: sred->pctelescope_reset_type = PCReset_Telescope_dmda;
410: break;
411: case TELESCOPE_DMPLEX: SETERRQ(comm,PETSC_ERR_SUP,"Support for DMPLEX is currently not available");
412: break;
413: default: SETERRQ(comm,PETSC_ERR_SUP,"Only support for repartitioning DMDA is provided");
414: break;
415: }
417: /* setup */
418: if (sred->pctelescope_setup_type) {
419: sred->pctelescope_setup_type(pc,sred);
420: }
421: /* update */
422: if (!pc->setupcalled) {
423: if (sred->pctelescope_matcreate_type) {
424: sred->pctelescope_matcreate_type(pc,sred,MAT_INITIAL_MATRIX,&sred->Bred);
425: }
426: if (sred->pctelescope_matnullspacecreate_type) {
427: sred->pctelescope_matnullspacecreate_type(pc,sred,sred->Bred);
428: }
429: } else {
430: if (sred->pctelescope_matcreate_type) {
431: sred->pctelescope_matcreate_type(pc,sred,MAT_REUSE_MATRIX,&sred->Bred);
432: }
433: }
435: /* common - no construction */
436: if (isActiveRank(sred->psubcomm)) {
437: KSPSetOperators(sred->ksp,sred->Bred,sred->Bred);
438: if (pc->setfromoptionscalled && !pc->setupcalled){
439: KSPSetFromOptions(sred->ksp);
440: }
441: }
442: return(0);
443: }
445: static PetscErrorCode PCApply_Telescope(PC pc,Vec x,Vec y)
446: {
447: PC_Telescope sred = (PC_Telescope)pc->data;
448: PetscErrorCode ierr;
449: Vec xtmp,xred,yred;
450: PetscInt i,st,ed;
451: VecScatter scatter;
452: PetscScalar *array;
453: const PetscScalar *x_array;
456: PetscCitationsRegister(citation,&cited);
458: xtmp = sred->xtmp;
459: scatter = sred->scatter;
460: xred = sred->xred;
461: yred = sred->yred;
463: /* pull in vector x->xtmp */
464: VecScatterBegin(scatter,x,xtmp,INSERT_VALUES,SCATTER_FORWARD);
465: VecScatterEnd(scatter,x,xtmp,INSERT_VALUES,SCATTER_FORWARD);
467: /* copy vector entries into xred */
468: VecGetArrayRead(xtmp,&x_array);
469: if (xred) {
470: PetscScalar *LA_xred;
471: VecGetOwnershipRange(xred,&st,&ed);
472: VecGetArray(xred,&LA_xred);
473: for (i=0; i<ed-st; i++) {
474: LA_xred[i] = x_array[i];
475: }
476: VecRestoreArray(xred,&LA_xred);
477: }
478: VecRestoreArrayRead(xtmp,&x_array);
479: /* solve */
480: if (isActiveRank(sred->psubcomm)) {
481: KSPSolve(sred->ksp,xred,yred);
482: }
483: /* return vector */
484: VecGetArray(xtmp,&array);
485: if (yred) {
486: const PetscScalar *LA_yred;
487: VecGetOwnershipRange(yred,&st,&ed);
488: VecGetArrayRead(yred,&LA_yred);
489: for (i=0; i<ed-st; i++) {
490: array[i] = LA_yred[i];
491: }
492: VecRestoreArrayRead(yred,&LA_yred);
493: }
494: VecRestoreArray(xtmp,&array);
495: VecScatterBegin(scatter,xtmp,y,INSERT_VALUES,SCATTER_REVERSE);
496: VecScatterEnd(scatter,xtmp,y,INSERT_VALUES,SCATTER_REVERSE);
497: return(0);
498: }
500: static PetscErrorCode PCApplyRichardson_Telescope(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool zeroguess,PetscInt *outits,PCRichardsonConvergedReason *reason)
501: {
502: PC_Telescope sred = (PC_Telescope)pc->data;
503: PetscErrorCode ierr;
504: Vec xtmp,yred;
505: PetscInt i,st,ed;
506: VecScatter scatter;
507: const PetscScalar *x_array;
508: PetscBool default_init_guess_value;
511: xtmp = sred->xtmp;
512: scatter = sred->scatter;
513: yred = sred->yred;
515: if (its > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCApplyRichardson_Telescope only supports max_it = 1");
516: *reason = (PCRichardsonConvergedReason)0;
518: if (!zeroguess) {
519: PetscInfo(pc,"PCTelescope: Scattering y for non-zero initial guess\n");
520: /* pull in vector y->xtmp */
521: VecScatterBegin(scatter,y,xtmp,INSERT_VALUES,SCATTER_FORWARD);
522: VecScatterEnd(scatter,y,xtmp,INSERT_VALUES,SCATTER_FORWARD);
524: /* copy vector entries into xred */
525: VecGetArrayRead(xtmp,&x_array);
526: if (yred) {
527: PetscScalar *LA_yred;
528: VecGetOwnershipRange(yred,&st,&ed);
529: VecGetArray(yred,&LA_yred);
530: for (i=0; i<ed-st; i++) {
531: LA_yred[i] = x_array[i];
532: }
533: VecRestoreArray(yred,&LA_yred);
534: }
535: VecRestoreArrayRead(xtmp,&x_array);
536: }
538: if (isActiveRank(sred->psubcomm)) {
539: KSPGetInitialGuessNonzero(sred->ksp,&default_init_guess_value);
540: if (!zeroguess) KSPSetInitialGuessNonzero(sred->ksp,PETSC_TRUE);
541: }
543: PCApply_Telescope(pc,x,y);
545: if (isActiveRank(sred->psubcomm)) {
546: KSPSetInitialGuessNonzero(sred->ksp,default_init_guess_value);
547: }
549: if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
550: *outits = 1;
551: return(0);
552: }
554: static PetscErrorCode PCReset_Telescope(PC pc)
555: {
556: PC_Telescope sred = (PC_Telescope)pc->data;
559: ISDestroy(&sred->isin);
560: VecScatterDestroy(&sred->scatter);
561: VecDestroy(&sred->xred);
562: VecDestroy(&sred->yred);
563: VecDestroy(&sred->xtmp);
564: MatDestroy(&sred->Bred);
565: KSPReset(sred->ksp);
566: if (sred->pctelescope_reset_type) {
567: sred->pctelescope_reset_type(pc);
568: }
569: return(0);
570: }
572: static PetscErrorCode PCDestroy_Telescope(PC pc)
573: {
574: PC_Telescope sred = (PC_Telescope)pc->data;
575: PetscErrorCode ierr;
578: PCReset_Telescope(pc);
579: KSPDestroy(&sred->ksp);
580: PetscSubcommDestroy(&sred->psubcomm);
581: PetscFree(sred->dm_ctx);
582: PetscFree(pc->data);
583: return(0);
584: }
586: static PetscErrorCode PCSetFromOptions_Telescope(PetscOptionItems *PetscOptionsObject,PC pc)
587: {
588: PC_Telescope sred = (PC_Telescope)pc->data;
589: PetscErrorCode ierr;
590: MPI_Comm comm;
591: PetscMPIInt size;
592: PetscBool flg;
593: PetscSubcommType subcommtype;
596: PetscObjectGetComm((PetscObject)pc,&comm);
597: MPI_Comm_size(comm,&size);
598: PetscOptionsHead(PetscOptionsObject,"Telescope options");
599: PetscOptionsEnum("-pc_telescope_subcomm_type","Subcomm type (interlaced or contiguous)","PCTelescopeSetSubcommType",PetscSubcommTypes,(PetscEnum)sred->subcommtype,(PetscEnum*)&subcommtype,&flg);
600: if (flg) {
601: PCTelescopeSetSubcommType(pc,subcommtype);
602: }
603: PetscOptionsInt("-pc_telescope_reduction_factor","Factor to reduce comm size by","PCTelescopeSetReductionFactor",sred->redfactor,&sred->redfactor,0);
604: if (sred->redfactor > size) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"-pc_telescope_reduction_factor <= comm size");
605: PetscOptionsBool("-pc_telescope_ignore_dm","Ignore any DM attached to the PC","PCTelescopeSetIgnoreDM",sred->ignore_dm,&sred->ignore_dm,0);
606: PetscOptionsBool("-pc_telescope_ignore_kspcomputeoperators","Ignore method used to compute A","PCTelescopeSetIgnoreKSPComputeOperators",sred->ignore_kspcomputeoperators,&sred->ignore_kspcomputeoperators,0);
607: PetscOptionsTail();
608: return(0);
609: }
611: /* PC simplementation specific API's */
613: static PetscErrorCode PCTelescopeGetKSP_Telescope(PC pc,KSP *ksp)
614: {
615: PC_Telescope red = (PC_Telescope)pc->data;
617: if (ksp) *ksp = red->ksp;
618: return(0);
619: }
621: static PetscErrorCode PCTelescopeGetSubcommType_Telescope(PC pc,PetscSubcommType *subcommtype)
622: {
623: PC_Telescope red = (PC_Telescope)pc->data;
625: if (subcommtype) *subcommtype = red->subcommtype;
626: return(0);
627: }
629: static PetscErrorCode PCTelescopeSetSubcommType_Telescope(PC pc,PetscSubcommType subcommtype)
630: {
631: PC_Telescope red = (PC_Telescope)pc->data;
634: if (pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"You cannot change the subcommunicator type for PCTelescope after it has been set up.");
635: red->subcommtype = subcommtype;
636: return(0);
637: }
639: static PetscErrorCode PCTelescopeGetReductionFactor_Telescope(PC pc,PetscInt *fact)
640: {
641: PC_Telescope red = (PC_Telescope)pc->data;
643: if (fact) *fact = red->redfactor;
644: return(0);
645: }
647: static PetscErrorCode PCTelescopeSetReductionFactor_Telescope(PC pc,PetscInt fact)
648: {
649: PC_Telescope red = (PC_Telescope)pc->data;
650: PetscMPIInt size;
651: PetscErrorCode ierr;
654: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
655: if (fact <= 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Reduction factor of telescoping PC %D must be positive",fact);
656: if (fact > size) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Reduction factor of telescoping PC %D must be <= comm.size",fact);
657: red->redfactor = fact;
658: return(0);
659: }
661: static PetscErrorCode PCTelescopeGetIgnoreDM_Telescope(PC pc,PetscBool *v)
662: {
663: PC_Telescope red = (PC_Telescope)pc->data;
665: if (v) *v = red->ignore_dm;
666: return(0);
667: }
669: static PetscErrorCode PCTelescopeSetIgnoreDM_Telescope(PC pc,PetscBool v)
670: {
671: PC_Telescope red = (PC_Telescope)pc->data;
673: red->ignore_dm = v;
674: return(0);
675: }
677: static PetscErrorCode PCTelescopeGetIgnoreKSPComputeOperators_Telescope(PC pc,PetscBool *v)
678: {
679: PC_Telescope red = (PC_Telescope)pc->data;
681: if (v) *v = red->ignore_kspcomputeoperators;
682: return(0);
683: }
685: static PetscErrorCode PCTelescopeSetIgnoreKSPComputeOperators_Telescope(PC pc,PetscBool v)
686: {
687: PC_Telescope red = (PC_Telescope)pc->data;
689: red->ignore_kspcomputeoperators = v;
690: return(0);
691: }
693: static PetscErrorCode PCTelescopeGetDM_Telescope(PC pc,DM *dm)
694: {
695: PC_Telescope red = (PC_Telescope)pc->data;
697: *dm = private_PCTelescopeGetSubDM(red);
698: return(0);
699: }
701: /*@
702: PCTelescopeGetKSP - Gets the KSP created by the telescoping PC.
704: Not Collective
706: Input Parameter:
707: . pc - the preconditioner context
709: Output Parameter:
710: . subksp - the KSP defined the smaller set of processes
712: Level: advanced
714: .keywords: PC, telescopting solve
715: @*/
716: PetscErrorCode PCTelescopeGetKSP(PC pc,KSP *subksp)
717: {
720: PetscUseMethod(pc,"PCTelescopeGetKSP_C",(PC,KSP*),(pc,subksp));
721: return(0);
722: }
724: /*@
725: PCTelescopeGetReductionFactor - Gets the factor by which the original number of processes has been reduced by.
727: Not Collective
729: Input Parameter:
730: . pc - the preconditioner context
732: Output Parameter:
733: . fact - the reduction factor
735: Level: advanced
737: .keywords: PC, telescoping solve
738: @*/
739: PetscErrorCode PCTelescopeGetReductionFactor(PC pc,PetscInt *fact)
740: {
743: PetscUseMethod(pc,"PCTelescopeGetReductionFactor_C",(PC,PetscInt*),(pc,fact));
744: return(0);
745: }
747: /*@
748: PCTelescopeSetReductionFactor - Sets the factor by which the original number of processes has been reduced by.
750: Not Collective
752: Input Parameter:
753: . pc - the preconditioner context
755: Output Parameter:
756: . fact - the reduction factor
758: Level: advanced
760: .keywords: PC, telescoping solve
761: @*/
762: PetscErrorCode PCTelescopeSetReductionFactor(PC pc,PetscInt fact)
763: {
766: PetscTryMethod(pc,"PCTelescopeSetReductionFactor_C",(PC,PetscInt),(pc,fact));
767: return(0);
768: }
770: /*@
771: PCTelescopeGetIgnoreDM - Get the flag indicating if any DM attached to the PC will be used.
773: Not Collective
775: Input Parameter:
776: . pc - the preconditioner context
778: Output Parameter:
779: . v - the flag
781: Level: advanced
783: .keywords: PC, telescoping solve
784: @*/
785: PetscErrorCode PCTelescopeGetIgnoreDM(PC pc,PetscBool *v)
786: {
789: PetscUseMethod(pc,"PCTelescopeGetIgnoreDM_C",(PC,PetscBool*),(pc,v));
790: return(0);
791: }
793: /*@
794: PCTelescopeSetIgnoreDM - Set a flag to ignore any DM attached to the PC.
796: Not Collective
798: Input Parameter:
799: . pc - the preconditioner context
801: Output Parameter:
802: . v - Use PETSC_TRUE to ignore any DM
804: Level: advanced
806: .keywords: PC, telescoping solve
807: @*/
808: PetscErrorCode PCTelescopeSetIgnoreDM(PC pc,PetscBool v)
809: {
812: PetscTryMethod(pc,"PCTelescopeSetIgnoreDM_C",(PC,PetscBool),(pc,v));
813: return(0);
814: }
816: /*@
817: PCTelescopeGetIgnoreKSPComputeOperators - Get the flag indicating if KSPComputeOperators will be used.
819: Not Collective
821: Input Parameter:
822: . pc - the preconditioner context
824: Output Parameter:
825: . v - the flag
827: Level: advanced
829: .keywords: PC, telescoping solve
830: @*/
831: PetscErrorCode PCTelescopeGetIgnoreKSPComputeOperators(PC pc,PetscBool *v)
832: {
835: PetscUseMethod(pc,"PCTelescopeGetIgnoreKSPComputeOperators_C",(PC,PetscBool*),(pc,v));
836: return(0);
837: }
839: /*@
840: PCTelescopeSetIgnoreKSPComputeOperators - Set a flag to ignore KSPComputeOperators.
842: Not Collective
844: Input Parameter:
845: . pc - the preconditioner context
847: Output Parameter:
848: . v - Use PETSC_TRUE to ignore the method (if defined) set via KSPSetComputeOperators on pc
850: Level: advanced
852: .keywords: PC, telescoping solve
853: @*/
854: PetscErrorCode PCTelescopeSetIgnoreKSPComputeOperators(PC pc,PetscBool v)
855: {
858: PetscTryMethod(pc,"PCTelescopeSetIgnoreKSPComputeOperators_C",(PC,PetscBool),(pc,v));
859: return(0);
860: }
862: /*@
863: PCTelescopeGetDM - Get the re-partitioned DM attached to the sub KSP.
865: Not Collective
867: Input Parameter:
868: . pc - the preconditioner context
870: Output Parameter:
871: . subdm - The re-partitioned DM
873: Level: advanced
875: .keywords: PC, telescoping solve
876: @*/
877: PetscErrorCode PCTelescopeGetDM(PC pc,DM *subdm)
878: {
881: PetscUseMethod(pc,"PCTelescopeGetDM_C",(PC,DM*),(pc,subdm));
882: return(0);
883: }
885: /*@
886: PCTelescopeSetSubcommType - set subcommunicator type (interlaced or contiguous)
888: Logically Collective
890: Input Parameter:
891: + pc - the preconditioner context
892: - subcommtype - the subcommunicator type (see PetscSubcommType)
894: Level: advanced
896: .keywords: PC, telescoping solve
898: .seealso: PetscSubcommType, PetscSubcomm, PCTELESCOPE
899: @*/
900: PetscErrorCode PCTelescopeSetSubcommType(PC pc, PetscSubcommType subcommtype)
901: {
904: PetscTryMethod(pc,"PCTelescopeSetSubcommType_C",(PC,PetscSubcommType),(pc,subcommtype));
905: return(0);
906: }
908: /*@
909: PCTelescopeGetSubcommType - Get the subcommunicator type (interlaced or contiguous)
911: Not Collective
913: Input Parameter:
914: . pc - the preconditioner context
916: Output Parameter:
917: . subcommtype - the subcommunicator type (see PetscSubcommType)
919: Level: advanced
921: .keywords: PC, telescoping solve
923: .seealso: PetscSubcomm, PetscSubcommType, PCTELESCOPE
924: @*/
925: PetscErrorCode PCTelescopeGetSubcommType(PC pc, PetscSubcommType *subcommtype)
926: {
929: PetscUseMethod(pc,"PCTelescopeGetSubcommType_C",(PC,PetscSubcommType*),(pc,subcommtype));
930: return(0);
931: }
933: /* -------------------------------------------------------------------------------------*/
934: /*MC
935: PCTELESCOPE - Runs a KSP solver on a sub-group of processors. MPI processes not in the sub-communicator are idle during the solve.
937: Options Database:
938: + -pc_telescope_reduction_factor <r> - factor to use communicator size by. e.g. with 64 MPI processes and r=4, the new sub-communicator will have 64/4 = 16 ranks.
939: - -pc_telescope_ignore_dm - flag to indicate whether an attached DM should be ignored
940: - -pc_telescope_subcomm_type <interlaced,contiguous> - how to define the reduced communicator. see PetscSubcomm for more.
942: Level: advanced
944: Notes:
945: The preconditioner is deemed telescopic as it only calls KSPSolve() on a single
946: sub-communicator, in contrast with PCREDUNDANT which calls KSPSolve() on N sub-communicators.
947: This means there will be MPI processes which will be idle during the application of this preconditioner.
949: The default KSP is PREONLY. If a DM is attached to the PC, it is re-partitioned on the sub-communicator.
950: Both the Bmat operator and the right hand side vector are permuted into the new DOF ordering defined by the re-partitioned DM.
951: Currently only support for re-partitioning a DMDA is provided.
952: Any nullspace attached to the original Bmat operator is extracted, re-partitioned and set on the repartitioned Bmat operator.
953: KSPSetComputeOperators() is not propagated to the sub KSP.
954: Currently there is no support for the flag -pc_use_amat
956: Assuming that the parent preconditioner (PC) is defined on a communicator c, this implementation
957: creates a child sub-communicator (c') containing fewer MPI processes than the original parent preconditioner (PC).
959: Developer Notes:
960: During PCSetup, the B operator is scattered onto c'.
961: Within PCApply, the RHS vector (x) is scattered into a redundant vector, xred (defined on c').
962: Then, KSPSolve() is executed on the c' communicator.
964: The communicator used within the telescoping preconditioner is defined by a PetscSubcomm using the INTERLACED
965: creation routine by default (this can be changed with -pc_telescope_subcomm_type). We run the sub KSP on only the ranks within the communicator which have a color equal to zero.
967: The telescoping preconditioner is aware of nullspaces and near nullspaces which are attached to the B operator.
968: In the case where B has a (near) nullspace attached, the (near) nullspace vectors are extracted from B and mapped into
969: a new (near) nullspace, defined on the sub-communicator, which is attached to B' (the B operator which was scattered to c')
971: The telescoping preconditioner is aware of an attached DM. In the event that the DM is of type DMDA (2D or 3D -
972: 1D support for 1D DMDAs is not provided), a new DMDA is created on c' (e.g. it is re-partitioned), and this new DM
973: is attached the sub KSPSolve(). The design of telescope is such that it should be possible to extend support
974: for re-partitioning other to DM's (e.g. DMPLEX). The user can supply a flag to ignore attached DMs.
976: By default, B' is defined by simply fusing rows from different MPI processes
978: When a DMDA is attached to the parent preconditioner, B' is defined by: (i) performing a symmetric permutation of B
979: into the ordering defined by the DMDA on c', (ii) extracting the local chunks via MatCreateSubMatrices(), (iii) fusing the
980: locally (sequential) matrices defined on the ranks common to c and c' into B' using MatCreateMPIMatConcatenateSeqMat()
982: Limitations/improvements include the following.
983: VecPlaceArray() could be used within PCApply() to improve efficiency and reduce memory usage.
985: The symmetric permutation used when a DMDA is encountered is performed via explicitly assmbleming a permutation matrix P,
986: and performing P^T.A.P. Possibly it might be more efficient to use MatPermute(). We opted to use P^T.A.P as it appears
987: VecPermute() does not supported for the use case required here. By computing P, one can permute both the operator and RHS in a
988: consistent manner.
990: Mapping of vectors is performed in the following way.
991: Suppose the parent comm size was 4, and we set a reduction factor of 2; this would give a comm size on c' of 2.
992: Using the interlaced creation routine, the ranks in c with color = 0 will be rank 0 and 2.
993: We perform the scatter to the sub-comm in the following way.
994: [1] Given a vector x defined on comm c
996: rank(c) : _________ 0 ______ ________ 1 _______ ________ 2 _____________ ___________ 3 __________
997: x : [0, 1, 2, 3, 4, 5] [6, 7, 8, 9, 10, 11] [12, 13, 14, 15, 16, 17] [18, 19, 20, 21, 22, 23]
999: scatter to xtmp defined also on comm c so that we have the following values
1001: rank(c) : ___________________ 0 ________________ _1_ ______________________ 2 _______________________ __3_
1002: xtmp : [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] [ ] [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] [ ]
1004: The entries on rank 1 and 3 (ranks which do not have a color = 0 in c') have no values
1007: [2] Copy the value from rank 0, 2 (indices with respect to comm c) into the vector xred which is defined on communicator c'.
1008: Ranks 0 and 2 are the only ranks in the subcomm which have a color = 0.
1010: rank(c') : ___________________ 0 _______________ ______________________ 1 _____________________
1011: xred : [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]
1014: Contributed by Dave May
1016: Reference:
1017: Dave A. May, Patrick Sanan, Karl Rupp, Matthew G. Knepley, and Barry F. Smith, "Extreme-Scale Multigrid Components within PETSc". 2016. In Proceedings of the Platform for Advanced Scientific Computing Conference (PASC '16). DOI: 10.1145/2929908.2929913
1019: .seealso: PCTelescopeGetKSP(), PCTelescopeGetDM(), PCTelescopeGetReductionFactor(), PCTelescopeSetReductionFactor(), PCTelescopeGetIgnoreDM(), PCTelescopeSetIgnoreDM(), PCREDUNDANT
1020: M*/
1021: PETSC_EXTERN PetscErrorCode PCCreate_Telescope(PC pc)
1022: {
1023: PetscErrorCode ierr;
1024: struct _PC_Telescope *sred;
1027: PetscNewLog(pc,&sred);
1028: sred->subcommtype = PETSC_SUBCOMM_INTERLACED;
1029: sred->redfactor = 1;
1030: sred->ignore_dm = PETSC_FALSE;
1031: sred->ignore_kspcomputeoperators = PETSC_FALSE;
1032: pc->data = (void*)sred;
1034: pc->ops->apply = PCApply_Telescope;
1035: pc->ops->applytranspose = NULL;
1036: pc->ops->applyrichardson = PCApplyRichardson_Telescope;
1037: pc->ops->setup = PCSetUp_Telescope;
1038: pc->ops->destroy = PCDestroy_Telescope;
1039: pc->ops->reset = PCReset_Telescope;
1040: pc->ops->setfromoptions = PCSetFromOptions_Telescope;
1041: pc->ops->view = PCView_Telescope;
1043: sred->pctelescope_setup_type = PCTelescopeSetUp_default;
1044: sred->pctelescope_matcreate_type = PCTelescopeMatCreate_default;
1045: sred->pctelescope_matnullspacecreate_type = PCTelescopeMatNullSpaceCreate_default;
1046: sred->pctelescope_reset_type = NULL;
1048: PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetKSP_C",PCTelescopeGetKSP_Telescope);
1049: PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetSubcommType_C",PCTelescopeGetSubcommType_Telescope);
1050: PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetSubcommType_C",PCTelescopeSetSubcommType_Telescope);
1051: PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetReductionFactor_C",PCTelescopeGetReductionFactor_Telescope);
1052: PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetReductionFactor_C",PCTelescopeSetReductionFactor_Telescope);
1053: PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetIgnoreDM_C",PCTelescopeGetIgnoreDM_Telescope);
1054: PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetIgnoreDM_C",PCTelescopeSetIgnoreDM_Telescope);
1055: PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetIgnoreKSPComputeOperators_C",PCTelescopeGetIgnoreKSPComputeOperators_Telescope);
1056: PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeSetIgnoreKSPComputeOperators_C",PCTelescopeSetIgnoreKSPComputeOperators_Telescope);
1057: PetscObjectComposeFunction((PetscObject)pc,"PCTelescopeGetDM_C",PCTelescopeGetDM_Telescope);
1058: return(0);
1059: }