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