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