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