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