xref: /petsc/src/ksp/pc/impls/telescope/telescope_dmda.c (revision e5a36eccef3d6b83a2c625c30d0dfd5adb4001f2) 
1 
2 
3 #include <petsc/private/matimpl.h>
4 #include <petsc/private/pcimpl.h>
5 #include <petsc/private/dmimpl.h>
6 #include <petscksp.h>           /*I "petscksp.h" I*/
7 #include <petscdm.h>
8 #include <petscdmda.h>
9 
10 #include "../src/ksp/pc/impls/telescope/telescope.h"
11 
12 static PetscBool  cited = PETSC_FALSE;
13 static const char citation[] =
14 "@inproceedings{MaySananRuppKnepleySmith2016,\n"
15 "  title     = {Extreme-Scale Multigrid Components within PETSc},\n"
16 "  author    = {Dave A. May and Patrick Sanan and Karl Rupp and Matthew G. Knepley and Barry F. Smith},\n"
17 "  booktitle = {Proceedings of the Platform for Advanced Scientific Computing Conference},\n"
18 "  series    = {PASC '16},\n"
19 "  isbn      = {978-1-4503-4126-4},\n"
20 "  location  = {Lausanne, Switzerland},\n"
21 "  pages     = {5:1--5:12},\n"
22 "  articleno = {5},\n"
23 "  numpages  = {12},\n"
24 "  url       = {http://doi.acm.org/10.1145/2929908.2929913},\n"
25 "  doi       = {10.1145/2929908.2929913},\n"
26 "  acmid     = {2929913},\n"
27 "  publisher = {ACM},\n"
28 "  address   = {New York, NY, USA},\n"
29 "  keywords  = {GPU, HPC, agglomeration, coarse-level solver, multigrid, parallel computing, preconditioning},\n"
30 "  year      = {2016}\n"
31 "}\n";
32 
33 static PetscErrorCode _DMDADetermineRankFromGlobalIJK(PetscInt dim,PetscInt i,PetscInt j,PetscInt k,
34                                                       PetscInt Mp,PetscInt Np,PetscInt Pp,
35                                                       PetscInt start_i[],PetscInt start_j[],PetscInt start_k[],
36                                                       PetscInt span_i[],PetscInt span_j[],PetscInt span_k[],
37                                                       PetscMPIInt *_pi,PetscMPIInt *_pj,PetscMPIInt *_pk,PetscMPIInt *rank_re)
38 {
39   PetscInt pi,pj,pk,n;
40 
41   PetscFunctionBegin;
42   *rank_re = -1;
43   if (_pi) *_pi = -1;
44   if (_pj) *_pj = -1;
45   if (_pk) *_pk = -1;
46   pi = pj = pk = -1;
47   if (_pi) {
48     for (n=0; n<Mp; n++) {
49       if ( (i >= start_i[n]) && (i < start_i[n]+span_i[n]) ) {
50         pi = n;
51         break;
52       }
53     }
54     if (pi == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmda-ijk] pi cannot be determined : range %D, val %D",Mp,i);
55     *_pi = pi;
56   }
57 
58   if (_pj) {
59     for (n=0; n<Np; n++) {
60       if ( (j >= start_j[n]) && (j < start_j[n]+span_j[n]) ) {
61         pj = n;
62         break;
63       }
64     }
65     if (pj == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmda-ijk] pj cannot be determined : range %D, val %D",Np,j);
66     *_pj = pj;
67   }
68 
69   if (_pk) {
70     for (n=0; n<Pp; n++) {
71       if ( (k >= start_k[n]) && (k < start_k[n]+span_k[n]) ) {
72         pk = n;
73         break;
74       }
75     }
76     if (pk == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmda-ijk] pk cannot be determined : range %D, val %D",Pp,k);
77     *_pk = pk;
78   }
79 
80   switch (dim) {
81   case 1:
82     *rank_re = pi;
83     break;
84   case 2:
85     *rank_re = pi + pj * Mp;
86     break;
87   case 3:
88     *rank_re = pi + pj * Mp + pk * (Mp*Np);
89     break;
90   }
91   PetscFunctionReturn(0);
92 }
93 
94 static PetscErrorCode _DMDADetermineGlobalS0(PetscInt dim,PetscMPIInt rank_re,PetscInt Mp_re,PetscInt Np_re,PetscInt Pp_re,
95                                       PetscInt range_i_re[],PetscInt range_j_re[],PetscInt range_k_re[],PetscInt *s0)
96 {
97   PetscInt i,j,k,start_IJK = 0;
98   PetscInt rank_ijk;
99 
100   PetscFunctionBegin;
101   switch (dim) {
102   case 1:
103     for (i=0; i<Mp_re; i++) {
104       rank_ijk = i;
105       if (rank_ijk < rank_re) {
106         start_IJK += range_i_re[i];
107       }
108     }
109     break;
110     case 2:
111     for (j=0; j<Np_re; j++) {
112       for (i=0; i<Mp_re; i++) {
113         rank_ijk = i + j*Mp_re;
114         if (rank_ijk < rank_re) {
115           start_IJK += range_i_re[i]*range_j_re[j];
116         }
117       }
118     }
119     break;
120     case 3:
121     for (k=0; k<Pp_re; k++) {
122       for (j=0; j<Np_re; j++) {
123         for (i=0; i<Mp_re; i++) {
124           rank_ijk = i + j*Mp_re + k*Mp_re*Np_re;
125           if (rank_ijk < rank_re) {
126             start_IJK += range_i_re[i]*range_j_re[j]*range_k_re[k];
127           }
128         }
129       }
130     }
131     break;
132   }
133   *s0 = start_IJK;
134   PetscFunctionReturn(0);
135 }
136 
137 static PetscErrorCode PCTelescopeSetUp_dmda_repart_coors2d(PC_Telescope sred,DM dm,DM subdm)
138 {
139   PetscErrorCode ierr;
140   DM             cdm;
141   Vec            coor,coor_natural,perm_coors;
142   PetscInt       i,j,si,sj,ni,nj,M,N,Ml,Nl,c,nidx;
143   PetscInt       *fine_indices;
144   IS             is_fine,is_local;
145   VecScatter     sctx;
146 
147   PetscFunctionBegin;
148   ierr = DMGetCoordinates(dm,&coor);CHKERRQ(ierr);
149   if (!coor) return(0);
150   if (PCTelescope_isActiveRank(sred)) {
151     ierr = DMDASetUniformCoordinates(subdm,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);
152   }
153   /* Get the coordinate vector from the distributed array */
154   ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr);
155   ierr = DMDACreateNaturalVector(cdm,&coor_natural);CHKERRQ(ierr);
156 
157   ierr = DMDAGlobalToNaturalBegin(cdm,coor,INSERT_VALUES,coor_natural);CHKERRQ(ierr);
158   ierr = DMDAGlobalToNaturalEnd(cdm,coor,INSERT_VALUES,coor_natural);CHKERRQ(ierr);
159 
160   /* get indices of the guys I want to grab */
161   ierr = DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
162   if (PCTelescope_isActiveRank(sred)) {
163     ierr = DMDAGetCorners(subdm,&si,&sj,NULL,&ni,&nj,NULL);CHKERRQ(ierr);
164     Ml = ni;
165     Nl = nj;
166   } else {
167     si = sj = 0;
168     ni = nj = 0;
169     Ml = Nl = 0;
170   }
171 
172   ierr = PetscMalloc1(Ml*Nl*2,&fine_indices);CHKERRQ(ierr);
173   c = 0;
174   if (PCTelescope_isActiveRank(sred)) {
175     for (j=sj; j<sj+nj; j++) {
176       for (i=si; i<si+ni; i++) {
177         nidx = (i) + (j)*M;
178         fine_indices[c  ] = 2 * nidx     ;
179         fine_indices[c+1] = 2 * nidx + 1 ;
180         c = c + 2;
181       }
182     }
183     if (c != Ml*Nl*2) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"c %D should equal 2 * Ml %D * Nl %D",c,Ml,Nl);
184   }
185 
186   /* generate scatter */
187   ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),Ml*Nl*2,fine_indices,PETSC_USE_POINTER,&is_fine);CHKERRQ(ierr);
188   ierr = ISCreateStride(PETSC_COMM_SELF,Ml*Nl*2,0,1,&is_local);CHKERRQ(ierr);
189 
190   /* scatter */
191   ierr = VecCreate(PETSC_COMM_SELF,&perm_coors);CHKERRQ(ierr);
192   ierr = VecSetSizes(perm_coors,PETSC_DECIDE,Ml*Nl*2);CHKERRQ(ierr);
193   ierr = VecSetType(perm_coors,VECSEQ);CHKERRQ(ierr);
194 
195   ierr = VecScatterCreate(coor_natural,is_fine,perm_coors,is_local,&sctx);CHKERRQ(ierr);
196   ierr = VecScatterBegin(sctx,coor_natural,perm_coors,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
197   ierr = VecScatterEnd(  sctx,coor_natural,perm_coors,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
198   /* access */
199   if (PCTelescope_isActiveRank(sred)) {
200     Vec               _coors;
201     const PetscScalar *LA_perm;
202     PetscScalar       *LA_coors;
203 
204     ierr = DMGetCoordinates(subdm,&_coors);CHKERRQ(ierr);
205     ierr = VecGetArrayRead(perm_coors,&LA_perm);CHKERRQ(ierr);
206     ierr = VecGetArray(_coors,&LA_coors);CHKERRQ(ierr);
207     for (i=0; i<Ml*Nl*2; i++) {
208       LA_coors[i] = LA_perm[i];
209     }
210     ierr = VecRestoreArray(_coors,&LA_coors);CHKERRQ(ierr);
211     ierr = VecRestoreArrayRead(perm_coors,&LA_perm);CHKERRQ(ierr);
212   }
213 
214   /* update local coords */
215   if (PCTelescope_isActiveRank(sred)) {
216     DM  _dmc;
217     Vec _coors,_coors_local;
218     ierr = DMGetCoordinateDM(subdm,&_dmc);CHKERRQ(ierr);
219     ierr = DMGetCoordinates(subdm,&_coors);CHKERRQ(ierr);
220     ierr = DMGetCoordinatesLocal(subdm,&_coors_local);CHKERRQ(ierr);
221     ierr = DMGlobalToLocalBegin(_dmc,_coors,INSERT_VALUES,_coors_local);CHKERRQ(ierr);
222     ierr = DMGlobalToLocalEnd(_dmc,_coors,INSERT_VALUES,_coors_local);CHKERRQ(ierr);
223   }
224   ierr = VecScatterDestroy(&sctx);CHKERRQ(ierr);
225   ierr = ISDestroy(&is_fine);CHKERRQ(ierr);
226   ierr = PetscFree(fine_indices);CHKERRQ(ierr);
227   ierr = ISDestroy(&is_local);CHKERRQ(ierr);
228   ierr = VecDestroy(&perm_coors);CHKERRQ(ierr);
229   ierr = VecDestroy(&coor_natural);CHKERRQ(ierr);
230   PetscFunctionReturn(0);
231 }
232 
233 static PetscErrorCode PCTelescopeSetUp_dmda_repart_coors3d(PC_Telescope sred,DM dm,DM subdm)
234 {
235   PetscErrorCode ierr;
236   DM             cdm;
237   Vec            coor,coor_natural,perm_coors;
238   PetscInt       i,j,k,si,sj,sk,ni,nj,nk,M,N,P,Ml,Nl,Pl,c,nidx;
239   PetscInt       *fine_indices;
240   IS             is_fine,is_local;
241   VecScatter     sctx;
242 
243   PetscFunctionBegin;
244   ierr = DMGetCoordinates(dm,&coor);CHKERRQ(ierr);
245   if (!coor) PetscFunctionReturn(0);
246 
247   if (PCTelescope_isActiveRank(sred)) {
248     ierr = DMDASetUniformCoordinates(subdm,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);
249   }
250 
251   /* Get the coordinate vector from the distributed array */
252   ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr);
253   ierr = DMDACreateNaturalVector(cdm,&coor_natural);CHKERRQ(ierr);
254   ierr = DMDAGlobalToNaturalBegin(cdm,coor,INSERT_VALUES,coor_natural);CHKERRQ(ierr);
255   ierr = DMDAGlobalToNaturalEnd(cdm,coor,INSERT_VALUES,coor_natural);CHKERRQ(ierr);
256 
257   /* get indices of the guys I want to grab */
258   ierr = DMDAGetInfo(dm,NULL,&M,&N,&P,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
259 
260   if (PCTelescope_isActiveRank(sred)) {
261     ierr = DMDAGetCorners(subdm,&si,&sj,&sk,&ni,&nj,&nk);CHKERRQ(ierr);
262     Ml = ni;
263     Nl = nj;
264     Pl = nk;
265   } else {
266     si = sj = sk = 0;
267     ni = nj = nk = 0;
268     Ml = Nl = Pl = 0;
269   }
270 
271   ierr = PetscMalloc1(Ml*Nl*Pl*3,&fine_indices);CHKERRQ(ierr);
272 
273   c = 0;
274   if (PCTelescope_isActiveRank(sred)) {
275     for (k=sk; k<sk+nk; k++) {
276       for (j=sj; j<sj+nj; j++) {
277         for (i=si; i<si+ni; i++) {
278           nidx = (i) + (j)*M + (k)*M*N;
279           fine_indices[c  ] = 3 * nidx     ;
280           fine_indices[c+1] = 3 * nidx + 1 ;
281           fine_indices[c+2] = 3 * nidx + 2 ;
282           c = c + 3;
283         }
284       }
285     }
286   }
287 
288   /* generate scatter */
289   ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),Ml*Nl*Pl*3,fine_indices,PETSC_USE_POINTER,&is_fine);CHKERRQ(ierr);
290   ierr = ISCreateStride(PETSC_COMM_SELF,Ml*Nl*Pl*3,0,1,&is_local);CHKERRQ(ierr);
291 
292   /* scatter */
293   ierr = VecCreate(PETSC_COMM_SELF,&perm_coors);CHKERRQ(ierr);
294   ierr = VecSetSizes(perm_coors,PETSC_DECIDE,Ml*Nl*Pl*3);CHKERRQ(ierr);
295   ierr = VecSetType(perm_coors,VECSEQ);CHKERRQ(ierr);
296   ierr = VecScatterCreate(coor_natural,is_fine,perm_coors,is_local,&sctx);CHKERRQ(ierr);
297   ierr = VecScatterBegin(sctx,coor_natural,perm_coors,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
298   ierr = VecScatterEnd(  sctx,coor_natural,perm_coors,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
299 
300   /* access */
301   if (PCTelescope_isActiveRank(sred)) {
302     Vec               _coors;
303     const PetscScalar *LA_perm;
304     PetscScalar       *LA_coors;
305 
306     ierr = DMGetCoordinates(subdm,&_coors);CHKERRQ(ierr);
307     ierr = VecGetArrayRead(perm_coors,&LA_perm);CHKERRQ(ierr);
308     ierr = VecGetArray(_coors,&LA_coors);CHKERRQ(ierr);
309     for (i=0; i<Ml*Nl*Pl*3; i++) {
310       LA_coors[i] = LA_perm[i];
311     }
312     ierr = VecRestoreArray(_coors,&LA_coors);CHKERRQ(ierr);
313     ierr = VecRestoreArrayRead(perm_coors,&LA_perm);CHKERRQ(ierr);
314   }
315 
316   /* update local coords */
317   if (PCTelescope_isActiveRank(sred)) {
318     DM  _dmc;
319     Vec _coors,_coors_local;
320 
321     ierr = DMGetCoordinateDM(subdm,&_dmc);CHKERRQ(ierr);
322     ierr = DMGetCoordinates(subdm,&_coors);CHKERRQ(ierr);
323     ierr = DMGetCoordinatesLocal(subdm,&_coors_local);CHKERRQ(ierr);
324     ierr = DMGlobalToLocalBegin(_dmc,_coors,INSERT_VALUES,_coors_local);CHKERRQ(ierr);
325     ierr = DMGlobalToLocalEnd(_dmc,_coors,INSERT_VALUES,_coors_local);CHKERRQ(ierr);
326   }
327 
328   ierr = VecScatterDestroy(&sctx);CHKERRQ(ierr);
329   ierr = ISDestroy(&is_fine);CHKERRQ(ierr);
330   ierr = PetscFree(fine_indices);CHKERRQ(ierr);
331   ierr = ISDestroy(&is_local);CHKERRQ(ierr);
332   ierr = VecDestroy(&perm_coors);CHKERRQ(ierr);
333   ierr = VecDestroy(&coor_natural);CHKERRQ(ierr);
334   PetscFunctionReturn(0);
335 }
336 
337 static PetscErrorCode PCTelescopeSetUp_dmda_repart_coors(PC pc,PC_Telescope sred,PC_Telescope_DMDACtx *ctx)
338 {
339   PetscInt       dim;
340   DM             dm,subdm;
341   PetscSubcomm   psubcomm;
342   MPI_Comm       comm;
343   Vec            coor;
344   PetscErrorCode ierr;
345 
346   PetscFunctionBegin;
347   ierr = PCGetDM(pc,&dm);CHKERRQ(ierr);
348   ierr = DMGetCoordinates(dm,&coor);CHKERRQ(ierr);
349   if (!coor) PetscFunctionReturn(0);
350   psubcomm = sred->psubcomm;
351   comm = PetscSubcommParent(psubcomm);
352   subdm = ctx->dmrepart;
353 
354   ierr = PetscInfo(pc,"PCTelescope: setting up the coordinates (DMDA)\n");CHKERRQ(ierr);
355   ierr = DMDAGetInfo(dm,&dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
356   switch (dim) {
357   case 1: SETERRQ(comm,PETSC_ERR_SUP,"Telescope: DMDA (1D) repartitioning not provided");
358     break;
359   case 2: ierr = PCTelescopeSetUp_dmda_repart_coors2d(sred,dm,subdm);CHKERRQ(ierr);
360     break;
361   case 3: ierr = PCTelescopeSetUp_dmda_repart_coors3d(sred,dm,subdm);CHKERRQ(ierr);
362     break;
363   }
364   PetscFunctionReturn(0);
365 }
366 
367 /* setup repartitioned dm */
368 PetscErrorCode PCTelescopeSetUp_dmda_repart(PC pc,PC_Telescope sred,PC_Telescope_DMDACtx *ctx)
369 {
370   PetscErrorCode        ierr;
371   DM                    dm;
372   PetscInt              dim,nx,ny,nz,ndof,nsw,sum,k;
373   DMBoundaryType        bx,by,bz;
374   DMDAStencilType       stencil;
375   const PetscInt        *_range_i_re;
376   const PetscInt        *_range_j_re;
377   const PetscInt        *_range_k_re;
378   DMDAInterpolationType itype;
379   PetscInt              refine_x,refine_y,refine_z;
380   MPI_Comm              comm,subcomm;
381   const char            *prefix;
382 
383   PetscFunctionBegin;
384   comm = PetscSubcommParent(sred->psubcomm);
385   subcomm = PetscSubcommChild(sred->psubcomm);
386   ierr = PCGetDM(pc,&dm);CHKERRQ(ierr);
387 
388   ierr = DMDAGetInfo(dm,&dim,&nx,&ny,&nz,NULL,NULL,NULL,&ndof,&nsw,&bx,&by,&bz,&stencil);CHKERRQ(ierr);
389   ierr = DMDAGetInterpolationType(dm,&itype);CHKERRQ(ierr);
390   ierr = DMDAGetRefinementFactor(dm,&refine_x,&refine_y,&refine_z);CHKERRQ(ierr);
391 
392   ctx->dmrepart = NULL;
393   _range_i_re = _range_j_re = _range_k_re = NULL;
394   /* Create DMDA on the child communicator */
395   if (PCTelescope_isActiveRank(sred)) {
396     switch (dim) {
397     case 1:
398       ierr = PetscInfo(pc,"PCTelescope: setting up the DMDA on comm subset (1D)\n");CHKERRQ(ierr);
399       /*ierr = DMDACreate1d(subcomm,bx,nx,ndof,nsw,NULL,&ctx->dmrepart);CHKERRQ(ierr);*/
400       ny = nz = 1;
401       by = bz = DM_BOUNDARY_NONE;
402       break;
403     case 2:
404       ierr = PetscInfo(pc,"PCTelescope: setting up the DMDA on comm subset (2D)\n");CHKERRQ(ierr);
405       /*ierr = DMDACreate2d(subcomm,bx,by,stencil,nx,ny, PETSC_DECIDE,PETSC_DECIDE, ndof,nsw, NULL,NULL,&ctx->dmrepart);CHKERRQ(ierr);*/
406       nz = 1;
407       bz = DM_BOUNDARY_NONE;
408       break;
409     case 3:
410       ierr = PetscInfo(pc,"PCTelescope: setting up the DMDA on comm subset (3D)\n");CHKERRQ(ierr);
411       /*ierr = DMDACreate3d(subcomm,bx,by,bz,stencil,nx,ny,nz, PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, ndof,nsw, NULL,NULL,NULL,&ctx->dmrepart);CHKERRQ(ierr);*/
412       break;
413     }
414     /*
415      The API DMDACreate1d(), DMDACreate2d(), DMDACreate3d() does not allow us to set/append
416      a unique option prefix for the DM, thus I prefer to expose the contents of these API's here.
417      This allows users to control the partitioning of the subDM.
418     */
419     ierr = DMDACreate(subcomm,&ctx->dmrepart);CHKERRQ(ierr);
420     /* Set unique option prefix name */
421     ierr = KSPGetOptionsPrefix(sred->ksp,&prefix);CHKERRQ(ierr);
422     ierr = DMSetOptionsPrefix(ctx->dmrepart,prefix);CHKERRQ(ierr);
423     ierr = DMAppendOptionsPrefix(ctx->dmrepart,"repart_");CHKERRQ(ierr);
424     /* standard setup from DMDACreate{1,2,3}d() */
425     ierr = DMSetDimension(ctx->dmrepart,dim);CHKERRQ(ierr);
426     ierr = DMDASetSizes(ctx->dmrepart,nx,ny,nz);CHKERRQ(ierr);
427     ierr = DMDASetNumProcs(ctx->dmrepart,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
428     ierr = DMDASetBoundaryType(ctx->dmrepart,bx,by,bz);CHKERRQ(ierr);
429     ierr = DMDASetDof(ctx->dmrepart,ndof);CHKERRQ(ierr);
430     ierr = DMDASetStencilType(ctx->dmrepart,stencil);CHKERRQ(ierr);
431     ierr = DMDASetStencilWidth(ctx->dmrepart,nsw);CHKERRQ(ierr);
432     ierr = DMDASetOwnershipRanges(ctx->dmrepart,NULL,NULL,NULL);CHKERRQ(ierr);
433     ierr = DMSetFromOptions(ctx->dmrepart);CHKERRQ(ierr);
434     ierr = DMSetUp(ctx->dmrepart);CHKERRQ(ierr);
435     /* Set refinement factors and interpolation type from the partent */
436     ierr = DMDASetRefinementFactor(ctx->dmrepart,refine_x,refine_y,refine_z);CHKERRQ(ierr);
437     ierr = DMDASetInterpolationType(ctx->dmrepart,itype);CHKERRQ(ierr);
438 
439     ierr = DMDAGetInfo(ctx->dmrepart,NULL,NULL,NULL,NULL,&ctx->Mp_re,&ctx->Np_re,&ctx->Pp_re,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
440     ierr = DMDAGetOwnershipRanges(ctx->dmrepart,&_range_i_re,&_range_j_re,&_range_k_re);CHKERRQ(ierr);
441 
442     ctx->dmrepart->ops->creatematrix = dm->ops->creatematrix;
443     ctx->dmrepart->ops->createdomaindecomposition = dm->ops->createdomaindecomposition;
444   }
445 
446   /* generate ranges for repartitioned dm */
447   /* note - assume rank 0 always participates */
448   ierr = MPI_Bcast(&ctx->Mp_re,1,MPIU_INT,0,comm);CHKERRQ(ierr);
449   ierr = MPI_Bcast(&ctx->Np_re,1,MPIU_INT,0,comm);CHKERRQ(ierr);
450   ierr = MPI_Bcast(&ctx->Pp_re,1,MPIU_INT,0,comm);CHKERRQ(ierr);
451 
452   ierr = PetscCalloc1(ctx->Mp_re,&ctx->range_i_re);CHKERRQ(ierr);
453   ierr = PetscCalloc1(ctx->Np_re,&ctx->range_j_re);CHKERRQ(ierr);
454   ierr = PetscCalloc1(ctx->Pp_re,&ctx->range_k_re);CHKERRQ(ierr);
455 
456   if (_range_i_re) {ierr = PetscMemcpy(ctx->range_i_re,_range_i_re,sizeof(PetscInt)*ctx->Mp_re);CHKERRQ(ierr);}
457   if (_range_j_re) {ierr = PetscMemcpy(ctx->range_j_re,_range_j_re,sizeof(PetscInt)*ctx->Np_re);CHKERRQ(ierr);}
458   if (_range_k_re) {ierr = PetscMemcpy(ctx->range_k_re,_range_k_re,sizeof(PetscInt)*ctx->Pp_re);CHKERRQ(ierr);}
459 
460   ierr = MPI_Bcast(ctx->range_i_re,ctx->Mp_re,MPIU_INT,0,comm);CHKERRQ(ierr);
461   ierr = MPI_Bcast(ctx->range_j_re,ctx->Np_re,MPIU_INT,0,comm);CHKERRQ(ierr);
462   ierr = MPI_Bcast(ctx->range_k_re,ctx->Pp_re,MPIU_INT,0,comm);CHKERRQ(ierr);
463 
464   ierr = PetscMalloc1(ctx->Mp_re,&ctx->start_i_re);CHKERRQ(ierr);
465   ierr = PetscMalloc1(ctx->Np_re,&ctx->start_j_re);CHKERRQ(ierr);
466   ierr = PetscMalloc1(ctx->Pp_re,&ctx->start_k_re);CHKERRQ(ierr);
467 
468   sum = 0;
469   for (k=0; k<ctx->Mp_re; k++) {
470     ctx->start_i_re[k] = sum;
471     sum += ctx->range_i_re[k];
472   }
473 
474   sum = 0;
475   for (k=0; k<ctx->Np_re; k++) {
476     ctx->start_j_re[k] = sum;
477     sum += ctx->range_j_re[k];
478   }
479 
480   sum = 0;
481   for (k=0; k<ctx->Pp_re; k++) {
482     ctx->start_k_re[k] = sum;
483     sum += ctx->range_k_re[k];
484   }
485 
486   /* attach repartitioned dm to child ksp */
487   {
488     PetscErrorCode (*dmksp_func)(KSP,Mat,Mat,void*);
489     void           *dmksp_ctx;
490 
491     ierr = DMKSPGetComputeOperators(dm,&dmksp_func,&dmksp_ctx);CHKERRQ(ierr);
492 
493     /* attach dm to ksp on sub communicator */
494     if (PCTelescope_isActiveRank(sred)) {
495       ierr = KSPSetDM(sred->ksp,ctx->dmrepart);CHKERRQ(ierr);
496 
497       if (!dmksp_func || sred->ignore_kspcomputeoperators) {
498         ierr = KSPSetDMActive(sred->ksp,PETSC_FALSE);CHKERRQ(ierr);
499       } else {
500         /* sub ksp inherits dmksp_func and context provided by user */
501         ierr = KSPSetComputeOperators(sred->ksp,dmksp_func,dmksp_ctx);CHKERRQ(ierr);
502         ierr = KSPSetDMActive(sred->ksp,PETSC_TRUE);CHKERRQ(ierr);
503       }
504     }
505   }
506   PetscFunctionReturn(0);
507 }
508 
509 PetscErrorCode PCTelescopeSetUp_dmda_permutation_3d(PC pc,PC_Telescope sred,PC_Telescope_DMDACtx *ctx)
510 {
511   PetscErrorCode ierr;
512   DM             dm;
513   MPI_Comm       comm;
514   Mat            Pscalar,P;
515   PetscInt       ndof;
516   PetscInt       i,j,k,location,startI[3],endI[3],lenI[3],nx,ny,nz;
517   PetscInt       sr,er,Mr;
518   Vec            V;
519 
520   PetscFunctionBegin;
521   ierr = PetscInfo(pc,"PCTelescope: setting up the permutation matrix (DMDA-3D)\n");CHKERRQ(ierr);
522   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
523 
524   ierr = PCGetDM(pc,&dm);CHKERRQ(ierr);
525   ierr = DMDAGetInfo(dm,NULL,&nx,&ny,&nz,NULL,NULL,NULL,&ndof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
526 
527   ierr = DMGetGlobalVector(dm,&V);CHKERRQ(ierr);
528   ierr = VecGetSize(V,&Mr);CHKERRQ(ierr);
529   ierr = VecGetOwnershipRange(V,&sr,&er);CHKERRQ(ierr);
530   ierr = DMRestoreGlobalVector(dm,&V);CHKERRQ(ierr);
531   sr = sr / ndof;
532   er = er / ndof;
533   Mr = Mr / ndof;
534 
535   ierr = MatCreate(comm,&Pscalar);CHKERRQ(ierr);
536   ierr = MatSetSizes(Pscalar,(er-sr),(er-sr),Mr,Mr);CHKERRQ(ierr);
537   ierr = MatSetType(Pscalar,MATAIJ);CHKERRQ(ierr);
538   ierr = MatSeqAIJSetPreallocation(Pscalar,1,NULL);CHKERRQ(ierr);
539   ierr = MatMPIAIJSetPreallocation(Pscalar,1,NULL,1,NULL);CHKERRQ(ierr);
540 
541   ierr = DMDAGetCorners(dm,NULL,NULL,NULL,&lenI[0],&lenI[1],&lenI[2]);CHKERRQ(ierr);
542   ierr = DMDAGetCorners(dm,&startI[0],&startI[1],&startI[2],&endI[0],&endI[1],&endI[2]);CHKERRQ(ierr);
543   endI[0] += startI[0];
544   endI[1] += startI[1];
545   endI[2] += startI[2];
546 
547   for (k=startI[2]; k<endI[2]; k++) {
548     for (j=startI[1]; j<endI[1]; j++) {
549       for (i=startI[0]; i<endI[0]; i++) {
550         PetscMPIInt rank_ijk_re,rank_reI[3];
551         PetscInt    s0_re;
552         PetscInt    ii,jj,kk,local_ijk_re,mapped_ijk;
553         PetscInt    lenI_re[3];
554 
555         location = (i - startI[0]) + (j - startI[1])*lenI[0] + (k - startI[2])*lenI[0]*lenI[1];
556         ierr = _DMDADetermineRankFromGlobalIJK(3,i,j,k,   ctx->Mp_re,ctx->Np_re,ctx->Pp_re,
557                                                ctx->start_i_re,ctx->start_j_re,ctx->start_k_re,
558                                                ctx->range_i_re,ctx->range_j_re,ctx->range_k_re,
559                                                &rank_reI[0],&rank_reI[1],&rank_reI[2],&rank_ijk_re);CHKERRQ(ierr);
560         ierr = _DMDADetermineGlobalS0(3,rank_ijk_re, ctx->Mp_re,ctx->Np_re,ctx->Pp_re, ctx->range_i_re,ctx->range_j_re,ctx->range_k_re, &s0_re);CHKERRQ(ierr);
561         ii = i - ctx->start_i_re[ rank_reI[0] ];
562         if (ii < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmdarepart-perm3d] index error ii");
563         jj = j - ctx->start_j_re[ rank_reI[1] ];
564         if (jj < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmdarepart-perm3d] index error jj");
565         kk = k - ctx->start_k_re[ rank_reI[2] ];
566         if (kk < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmdarepart-perm3d] index error kk");
567         lenI_re[0] = ctx->range_i_re[ rank_reI[0] ];
568         lenI_re[1] = ctx->range_j_re[ rank_reI[1] ];
569         lenI_re[2] = ctx->range_k_re[ rank_reI[2] ];
570         local_ijk_re = ii + jj * lenI_re[0] + kk * lenI_re[0] * lenI_re[1];
571         mapped_ijk = s0_re + local_ijk_re;
572         ierr = MatSetValue(Pscalar,sr+location,mapped_ijk,1.0,INSERT_VALUES);CHKERRQ(ierr);
573       }
574     }
575   }
576   ierr = MatAssemblyBegin(Pscalar,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
577   ierr = MatAssemblyEnd(Pscalar,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
578   ierr = MatCreateMAIJ(Pscalar,ndof,&P);CHKERRQ(ierr);
579   ierr = MatDestroy(&Pscalar);CHKERRQ(ierr);
580   ctx->permutation = P;
581   PetscFunctionReturn(0);
582 }
583 
584 PetscErrorCode PCTelescopeSetUp_dmda_permutation_2d(PC pc,PC_Telescope sred,PC_Telescope_DMDACtx *ctx)
585 {
586   PetscErrorCode ierr;
587   DM             dm;
588   MPI_Comm       comm;
589   Mat            Pscalar,P;
590   PetscInt       ndof;
591   PetscInt       i,j,location,startI[2],endI[2],lenI[2],nx,ny,nz;
592   PetscInt       sr,er,Mr;
593   Vec            V;
594 
595   PetscFunctionBegin;
596   ierr = PetscInfo(pc,"PCTelescope: setting up the permutation matrix (DMDA-2D)\n");CHKERRQ(ierr);
597   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
598   ierr = PCGetDM(pc,&dm);CHKERRQ(ierr);
599   ierr = DMDAGetInfo(dm,NULL,&nx,&ny,&nz,NULL,NULL,NULL,&ndof,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
600   ierr = DMGetGlobalVector(dm,&V);CHKERRQ(ierr);
601   ierr = VecGetSize(V,&Mr);CHKERRQ(ierr);
602   ierr = VecGetOwnershipRange(V,&sr,&er);CHKERRQ(ierr);
603   ierr = DMRestoreGlobalVector(dm,&V);CHKERRQ(ierr);
604   sr = sr / ndof;
605   er = er / ndof;
606   Mr = Mr / ndof;
607 
608   ierr = MatCreate(comm,&Pscalar);CHKERRQ(ierr);
609   ierr = MatSetSizes(Pscalar,(er-sr),(er-sr),Mr,Mr);CHKERRQ(ierr);
610   ierr = MatSetType(Pscalar,MATAIJ);CHKERRQ(ierr);
611   ierr = MatSeqAIJSetPreallocation(Pscalar,1,NULL);CHKERRQ(ierr);
612   ierr = MatMPIAIJSetPreallocation(Pscalar,1,NULL,1,NULL);CHKERRQ(ierr);
613 
614   ierr = DMDAGetCorners(dm,NULL,NULL,NULL,&lenI[0],&lenI[1],NULL);CHKERRQ(ierr);
615   ierr = DMDAGetCorners(dm,&startI[0],&startI[1],NULL,&endI[0],&endI[1],NULL);CHKERRQ(ierr);
616   endI[0] += startI[0];
617   endI[1] += startI[1];
618 
619   for (j=startI[1]; j<endI[1]; j++) {
620     for (i=startI[0]; i<endI[0]; i++) {
621       PetscMPIInt rank_ijk_re,rank_reI[3];
622       PetscInt    s0_re;
623       PetscInt    ii,jj,local_ijk_re,mapped_ijk;
624       PetscInt    lenI_re[3];
625 
626       location = (i - startI[0]) + (j - startI[1])*lenI[0];
627       ierr = _DMDADetermineRankFromGlobalIJK(2,i,j,0,   ctx->Mp_re,ctx->Np_re,ctx->Pp_re,
628                                              ctx->start_i_re,ctx->start_j_re,ctx->start_k_re,
629                                              ctx->range_i_re,ctx->range_j_re,ctx->range_k_re,
630                                              &rank_reI[0],&rank_reI[1],NULL,&rank_ijk_re);CHKERRQ(ierr);
631 
632       ierr = _DMDADetermineGlobalS0(2,rank_ijk_re, ctx->Mp_re,ctx->Np_re,ctx->Pp_re, ctx->range_i_re,ctx->range_j_re,ctx->range_k_re, &s0_re);CHKERRQ(ierr);
633 
634       ii = i - ctx->start_i_re[ rank_reI[0] ];
635       if (ii < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmdarepart-perm2d] index error ii");
636       jj = j - ctx->start_j_re[ rank_reI[1] ];
637       if (jj < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"[dmdarepart-perm2d] index error jj");
638 
639       lenI_re[0] = ctx->range_i_re[ rank_reI[0] ];
640       lenI_re[1] = ctx->range_j_re[ rank_reI[1] ];
641       local_ijk_re = ii + jj * lenI_re[0];
642       mapped_ijk = s0_re + local_ijk_re;
643       ierr = MatSetValue(Pscalar,sr+location,mapped_ijk,1.0,INSERT_VALUES);CHKERRQ(ierr);
644     }
645   }
646   ierr = MatAssemblyBegin(Pscalar,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
647   ierr = MatAssemblyEnd(Pscalar,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
648   ierr = MatCreateMAIJ(Pscalar,ndof,&P);CHKERRQ(ierr);
649   ierr = MatDestroy(&Pscalar);CHKERRQ(ierr);
650   ctx->permutation = P;
651   PetscFunctionReturn(0);
652 }
653 
654 PetscErrorCode PCTelescopeSetUp_dmda_scatters(PC pc,PC_Telescope sred,PC_Telescope_DMDACtx *ctx)
655 {
656   PetscErrorCode ierr;
657   Vec            xred,yred,xtmp,x,xp;
658   VecScatter     scatter;
659   IS             isin;
660   Mat            B;
661   PetscInt       m,bs,st,ed;
662   MPI_Comm       comm;
663 
664   PetscFunctionBegin;
665   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
666   ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr);
667   ierr = MatCreateVecs(B,&x,NULL);CHKERRQ(ierr);
668   ierr = MatGetBlockSize(B,&bs);CHKERRQ(ierr);
669   ierr = VecDuplicate(x,&xp);CHKERRQ(ierr);
670   m = 0;
671   xred = NULL;
672   yred = NULL;
673   if (PCTelescope_isActiveRank(sred)) {
674     ierr = DMCreateGlobalVector(ctx->dmrepart,&xred);CHKERRQ(ierr);
675     ierr = VecDuplicate(xred,&yred);CHKERRQ(ierr);
676     ierr = VecGetOwnershipRange(xred,&st,&ed);CHKERRQ(ierr);
677     ierr = ISCreateStride(comm,ed-st,st,1,&isin);CHKERRQ(ierr);
678     ierr = VecGetLocalSize(xred,&m);CHKERRQ(ierr);
679   } else {
680     ierr = VecGetOwnershipRange(x,&st,&ed);CHKERRQ(ierr);
681     ierr = ISCreateStride(comm,0,st,1,&isin);CHKERRQ(ierr);
682   }
683   ierr = ISSetBlockSize(isin,bs);CHKERRQ(ierr);
684   ierr = VecCreate(comm,&xtmp);CHKERRQ(ierr);
685   ierr = VecSetSizes(xtmp,m,PETSC_DECIDE);CHKERRQ(ierr);
686   ierr = VecSetBlockSize(xtmp,bs);CHKERRQ(ierr);
687   ierr = VecSetType(xtmp,((PetscObject)x)->type_name);CHKERRQ(ierr);
688   ierr = VecScatterCreate(x,isin,xtmp,NULL,&scatter);CHKERRQ(ierr);
689   sred->xred    = xred;
690   sred->yred    = yred;
691   sred->isin    = isin;
692   sred->scatter = scatter;
693   sred->xtmp    = xtmp;
694 
695   ctx->xp       = xp;
696   ierr = VecDestroy(&x);CHKERRQ(ierr);
697   PetscFunctionReturn(0);
698 }
699 
700 PetscErrorCode PCTelescopeSetUp_dmda(PC pc,PC_Telescope sred)
701 {
702   PC_Telescope_DMDACtx *ctx;
703   PetscInt             dim;
704   DM                   dm;
705   MPI_Comm             comm;
706   PetscErrorCode       ierr;
707 
708   PetscFunctionBegin;
709   ierr = PetscInfo(pc,"PCTelescope: setup (DMDA)\n");CHKERRQ(ierr);
710   PetscMalloc1(1,&ctx);
711   PetscMemzero(ctx,sizeof(PC_Telescope_DMDACtx));
712   sred->dm_ctx = (void*)ctx;
713 
714   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
715   ierr = PCGetDM(pc,&dm);CHKERRQ(ierr);
716   ierr = DMDAGetInfo(dm,&dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
717 
718   PCTelescopeSetUp_dmda_repart(pc,sred,ctx);
719   PCTelescopeSetUp_dmda_repart_coors(pc,sred,ctx);
720   switch (dim) {
721   case 1: SETERRQ(comm,PETSC_ERR_SUP,"Telescope: DMDA (1D) repartitioning not provided");
722     break;
723   case 2: ierr = PCTelescopeSetUp_dmda_permutation_2d(pc,sred,ctx);CHKERRQ(ierr);
724     break;
725   case 3: ierr = PCTelescopeSetUp_dmda_permutation_3d(pc,sred,ctx);CHKERRQ(ierr);
726     break;
727   }
728   ierr = PCTelescopeSetUp_dmda_scatters(pc,sred,ctx);CHKERRQ(ierr);
729   PetscFunctionReturn(0);
730 }
731 
732 PetscErrorCode PCTelescopeMatCreate_dmda_dmactivefalse(PC pc,PC_Telescope sred,MatReuse reuse,Mat *A)
733 {
734   PetscErrorCode       ierr;
735   PC_Telescope_DMDACtx *ctx;
736   MPI_Comm             comm,subcomm;
737   Mat                  Bperm,Bred,B,P;
738   PetscInt             nr,nc;
739   IS                   isrow,iscol;
740   Mat                  Blocal,*_Blocal;
741 
742   PetscFunctionBegin;
743   ierr = PetscInfo(pc,"PCTelescope: updating the redundant preconditioned operator (DMDA)\n");CHKERRQ(ierr);
744   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
745   subcomm = PetscSubcommChild(sred->psubcomm);
746   ctx = (PC_Telescope_DMDACtx*)sred->dm_ctx;
747 
748   ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr);
749   ierr = MatGetSize(B,&nr,&nc);CHKERRQ(ierr);
750 
751   P = ctx->permutation;
752   ierr = MatPtAP(B,P,MAT_INITIAL_MATRIX,1.1,&Bperm);CHKERRQ(ierr);
753 
754   /* Get submatrices */
755   isrow = sred->isin;
756   ierr = ISCreateStride(comm,nc,0,1,&iscol);CHKERRQ(ierr);
757 
758   ierr = MatCreateSubMatrices(Bperm,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&_Blocal);CHKERRQ(ierr);
759   Blocal = *_Blocal;
760   Bred = NULL;
761   if (PCTelescope_isActiveRank(sred)) {
762     PetscInt mm;
763 
764     if (reuse != MAT_INITIAL_MATRIX) {Bred = *A;}
765     ierr = MatGetSize(Blocal,&mm,NULL);CHKERRQ(ierr);
766     /* ierr = MatCreateMPIMatConcatenateSeqMat(subcomm,Blocal,PETSC_DECIDE,reuse,&Bred);CHKERRQ(ierr); */
767     ierr = MatCreateMPIMatConcatenateSeqMat(subcomm,Blocal,mm,reuse,&Bred);CHKERRQ(ierr);
768   }
769   *A = Bred;
770 
771   ierr = ISDestroy(&iscol);CHKERRQ(ierr);
772   ierr = MatDestroy(&Bperm);CHKERRQ(ierr);
773   ierr = MatDestroyMatrices(1,&_Blocal);CHKERRQ(ierr);
774   PetscFunctionReturn(0);
775 }
776 
777 PetscErrorCode PCTelescopeMatCreate_dmda(PC pc,PC_Telescope sred,MatReuse reuse,Mat *A)
778 {
779   PetscErrorCode ierr;
780   DM             dm;
781   PetscErrorCode (*dmksp_func)(KSP,Mat,Mat,void*);
782   void           *dmksp_ctx;
783 
784   PetscFunctionBegin;
785   ierr = PCGetDM(pc,&dm);CHKERRQ(ierr);
786   ierr = DMKSPGetComputeOperators(dm,&dmksp_func,&dmksp_ctx);CHKERRQ(ierr);
787   /* We assume that dmksp_func = NULL, is equivalent to dmActive = PETSC_FALSE */
788   if (dmksp_func && !sred->ignore_kspcomputeoperators) {
789     DM  dmrepart;
790     Mat Ak;
791 
792     *A = NULL;
793     if (PCTelescope_isActiveRank(sred)) {
794       ierr = KSPGetDM(sred->ksp,&dmrepart);CHKERRQ(ierr);
795       if (reuse == MAT_INITIAL_MATRIX) {
796         ierr = DMCreateMatrix(dmrepart,&Ak);CHKERRQ(ierr);
797         *A = Ak;
798       } else if (reuse == MAT_REUSE_MATRIX) {
799         Ak = *A;
800       }
801       /*
802        There is no need to explicitly assemble the operator now,
803        the sub-KSP will call the method provided to KSPSetComputeOperators() during KSPSetUp()
804       */
805     }
806   } else {
807     ierr = PCTelescopeMatCreate_dmda_dmactivefalse(pc,sred,reuse,A);CHKERRQ(ierr);
808   }
809   PetscFunctionReturn(0);
810 }
811 
812 PetscErrorCode PCTelescopeSubNullSpaceCreate_dmda_Telescope(PC pc,PC_Telescope sred,MatNullSpace nullspace,MatNullSpace *sub_nullspace)
813 {
814   PetscErrorCode       ierr;
815   PetscBool            has_const;
816   PetscInt             i,k,n = 0;
817   const Vec            *vecs;
818   Vec                  *sub_vecs = NULL;
819   MPI_Comm             subcomm;
820   PC_Telescope_DMDACtx *ctx;
821 
822   PetscFunctionBegin;
823   ctx = (PC_Telescope_DMDACtx*)sred->dm_ctx;
824   subcomm = PetscSubcommChild(sred->psubcomm);
825   ierr = MatNullSpaceGetVecs(nullspace,&has_const,&n,&vecs);CHKERRQ(ierr);
826 
827   if (PCTelescope_isActiveRank(sred)) {
828     /* create new vectors */
829     if (n) {
830       ierr = VecDuplicateVecs(sred->xred,n,&sub_vecs);CHKERRQ(ierr);
831     }
832   }
833 
834   /* copy entries */
835   for (k=0; k<n; k++) {
836     const PetscScalar *x_array;
837     PetscScalar       *LA_sub_vec;
838     PetscInt          st,ed;
839 
840     /* permute vector into ordering associated with re-partitioned dmda */
841     ierr = MatMultTranspose(ctx->permutation,vecs[k],ctx->xp);CHKERRQ(ierr);
842 
843     /* pull in vector x->xtmp */
844     ierr = VecScatterBegin(sred->scatter,ctx->xp,sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
845     ierr = VecScatterEnd(sred->scatter,ctx->xp,sred->xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
846 
847     /* copy vector entries into xred */
848     ierr = VecGetArrayRead(sred->xtmp,&x_array);CHKERRQ(ierr);
849     if (sub_vecs) {
850       if (sub_vecs[k]) {
851         ierr = VecGetOwnershipRange(sub_vecs[k],&st,&ed);CHKERRQ(ierr);
852         ierr = VecGetArray(sub_vecs[k],&LA_sub_vec);CHKERRQ(ierr);
853         for (i=0; i<ed-st; i++) {
854           LA_sub_vec[i] = x_array[i];
855         }
856         ierr = VecRestoreArray(sub_vecs[k],&LA_sub_vec);CHKERRQ(ierr);
857       }
858     }
859     ierr = VecRestoreArrayRead(sred->xtmp,&x_array);CHKERRQ(ierr);
860   }
861 
862   if (PCTelescope_isActiveRank(sred)) {
863     /* create new (near) nullspace for redundant object */
864     ierr = MatNullSpaceCreate(subcomm,has_const,n,sub_vecs,sub_nullspace);CHKERRQ(ierr);
865     ierr = VecDestroyVecs(n,&sub_vecs);CHKERRQ(ierr);
866     if (nullspace->remove) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Propagation of custom remove callbacks not supported when propagating (near) nullspaces with PCTelescope");
867     if (nullspace->rmctx) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Propagation of custom remove callback context not supported when propagating (near) nullspaces with PCTelescope");
868   }
869   PetscFunctionReturn(0);
870 }
871 
872 PetscErrorCode PCTelescopeMatNullSpaceCreate_dmda(PC pc,PC_Telescope sred,Mat sub_mat)
873 {
874   PetscErrorCode ierr;
875   Mat            B;
876 
877   PetscFunctionBegin;
878   ierr = PCGetOperators(pc,NULL,&B);CHKERRQ(ierr);
879   {
880     MatNullSpace nullspace,sub_nullspace;
881     ierr = MatGetNullSpace(B,&nullspace);CHKERRQ(ierr);
882     if (nullspace) {
883       ierr = PetscInfo(pc,"PCTelescope: generating nullspace (DMDA)\n");CHKERRQ(ierr);
884       ierr = PCTelescopeSubNullSpaceCreate_dmda_Telescope(pc,sred,nullspace,&sub_nullspace);CHKERRQ(ierr);
885       if (PCTelescope_isActiveRank(sred)) {
886         ierr = MatSetNullSpace(sub_mat,sub_nullspace);CHKERRQ(ierr);
887         ierr = MatNullSpaceDestroy(&sub_nullspace);CHKERRQ(ierr);
888       }
889     }
890   }
891   {
892     MatNullSpace nearnullspace,sub_nearnullspace;
893     ierr = MatGetNearNullSpace(B,&nearnullspace);CHKERRQ(ierr);
894     if (nearnullspace) {
895       ierr = PetscInfo(pc,"PCTelescope: generating near nullspace (DMDA)\n");CHKERRQ(ierr);
896       ierr = PCTelescopeSubNullSpaceCreate_dmda_Telescope(pc,sred,nearnullspace,&sub_nearnullspace);CHKERRQ(ierr);
897       if (PCTelescope_isActiveRank(sred)) {
898         ierr = MatSetNearNullSpace(sub_mat,sub_nearnullspace);CHKERRQ(ierr);
899         ierr = MatNullSpaceDestroy(&sub_nearnullspace);CHKERRQ(ierr);
900       }
901     }
902   }
903   PetscFunctionReturn(0);
904 }
905 
906 PetscErrorCode PCApply_Telescope_dmda(PC pc,Vec x,Vec y)
907 {
908   PC_Telescope         sred = (PC_Telescope)pc->data;
909   PetscErrorCode       ierr;
910   Mat                  perm;
911   Vec                  xtmp,xp,xred,yred;
912   PetscInt             i,st,ed;
913   VecScatter           scatter;
914   PetscScalar          *array;
915   const PetscScalar    *x_array;
916   PC_Telescope_DMDACtx *ctx;
917 
918   ctx = (PC_Telescope_DMDACtx*)sred->dm_ctx;
919   xtmp    = sred->xtmp;
920   scatter = sred->scatter;
921   xred    = sred->xred;
922   yred    = sred->yred;
923   perm    = ctx->permutation;
924   xp      = ctx->xp;
925 
926   PetscFunctionBegin;
927   ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr);
928 
929   /* permute vector into ordering associated with re-partitioned dmda */
930   ierr = MatMultTranspose(perm,x,xp);CHKERRQ(ierr);
931 
932   /* pull in vector x->xtmp */
933   ierr = VecScatterBegin(scatter,xp,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
934   ierr = VecScatterEnd(scatter,xp,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
935 
936   /* copy vector entires into xred */
937   ierr = VecGetArrayRead(xtmp,&x_array);CHKERRQ(ierr);
938   if (xred) {
939     PetscScalar *LA_xred;
940     ierr = VecGetOwnershipRange(xred,&st,&ed);CHKERRQ(ierr);
941 
942     ierr = VecGetArray(xred,&LA_xred);CHKERRQ(ierr);
943     for (i=0; i<ed-st; i++) {
944       LA_xred[i] = x_array[i];
945     }
946     ierr = VecRestoreArray(xred,&LA_xred);CHKERRQ(ierr);
947   }
948   ierr = VecRestoreArrayRead(xtmp,&x_array);CHKERRQ(ierr);
949 
950   /* solve */
951   if (PCTelescope_isActiveRank(sred)) {
952     ierr = KSPSolve(sred->ksp,xred,yred);CHKERRQ(ierr);
953     ierr = KSPCheckSolve(sred->ksp,pc,yred);CHKERRQ(ierr);
954   }
955 
956   /* return vector */
957   ierr = VecGetArray(xtmp,&array);CHKERRQ(ierr);
958   if (yred) {
959     const PetscScalar *LA_yred;
960     ierr = VecGetOwnershipRange(yred,&st,&ed);CHKERRQ(ierr);
961     ierr = VecGetArrayRead(yred,&LA_yred);CHKERRQ(ierr);
962     for (i=0; i<ed-st; i++) {
963       array[i] = LA_yred[i];
964     }
965     ierr = VecRestoreArrayRead(yred,&LA_yred);CHKERRQ(ierr);
966   }
967   ierr = VecRestoreArray(xtmp,&array);CHKERRQ(ierr);
968   ierr = VecScatterBegin(scatter,xtmp,xp,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
969   ierr = VecScatterEnd(scatter,xtmp,xp,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
970   ierr = MatMult(perm,xp,y);CHKERRQ(ierr);
971   PetscFunctionReturn(0);
972 }
973 
974 PetscErrorCode PCApplyRichardson_Telescope_dmda(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool zeroguess,PetscInt *outits,PCRichardsonConvergedReason *reason)
975 {
976   PC_Telescope         sred = (PC_Telescope)pc->data;
977   PetscErrorCode       ierr;
978   Mat                  perm;
979   Vec                  xtmp,xp,yred;
980   PetscInt             i,st,ed;
981   VecScatter           scatter;
982   const PetscScalar    *x_array;
983   PetscBool            default_init_guess_value = PETSC_FALSE;
984   PC_Telescope_DMDACtx *ctx;
985 
986   PetscFunctionBegin;
987   ctx = (PC_Telescope_DMDACtx*)sred->dm_ctx;
988   xtmp    = sred->xtmp;
989   scatter = sred->scatter;
990   yred    = sred->yred;
991   perm    = ctx->permutation;
992   xp      = ctx->xp;
993 
994   if (its > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCApplyRichardson_Telescope_dmda only supports max_it = 1");
995   *reason = (PCRichardsonConvergedReason)0;
996 
997   if (!zeroguess) {
998     ierr = PetscInfo(pc,"PCTelescopeDMDA: Scattering y for non-zero-initial guess\n");CHKERRQ(ierr);
999     /* permute vector into ordering associated with re-partitioned dmda */
1000     ierr = MatMultTranspose(perm,y,xp);CHKERRQ(ierr);
1001 
1002     /* pull in vector x->xtmp */
1003     ierr = VecScatterBegin(scatter,xp,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1004     ierr = VecScatterEnd(scatter,xp,xtmp,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1005 
1006     /* copy vector entires into xred */
1007     ierr = VecGetArrayRead(xtmp,&x_array);CHKERRQ(ierr);
1008     if (yred) {
1009       PetscScalar *LA_yred;
1010       ierr = VecGetOwnershipRange(yred,&st,&ed);CHKERRQ(ierr);
1011       ierr = VecGetArray(yred,&LA_yred);CHKERRQ(ierr);
1012       for (i=0; i<ed-st; i++) {
1013         LA_yred[i] = x_array[i];
1014       }
1015       ierr = VecRestoreArray(yred,&LA_yred);CHKERRQ(ierr);
1016     }
1017     ierr = VecRestoreArrayRead(xtmp,&x_array);CHKERRQ(ierr);
1018   }
1019 
1020   if (PCTelescope_isActiveRank(sred)) {
1021     ierr = KSPGetInitialGuessNonzero(sred->ksp,&default_init_guess_value);CHKERRQ(ierr);
1022     if (!zeroguess) ierr = KSPSetInitialGuessNonzero(sred->ksp,PETSC_TRUE);CHKERRQ(ierr);
1023   }
1024 
1025   ierr = PCApply_Telescope_dmda(pc,x,y);CHKERRQ(ierr);
1026 
1027   if (PCTelescope_isActiveRank(sred)) {
1028     ierr = KSPSetInitialGuessNonzero(sred->ksp,default_init_guess_value);CHKERRQ(ierr);
1029   }
1030 
1031   if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
1032   *outits = 1;
1033   PetscFunctionReturn(0);
1034 }
1035 
1036 PetscErrorCode PCReset_Telescope_dmda(PC pc)
1037 {
1038   PetscErrorCode       ierr;
1039   PC_Telescope         sred = (PC_Telescope)pc->data;
1040   PC_Telescope_DMDACtx *ctx;
1041 
1042   PetscFunctionBegin;
1043   ctx = (PC_Telescope_DMDACtx*)sred->dm_ctx;
1044   ierr = VecDestroy(&ctx->xp);CHKERRQ(ierr);
1045   ierr = MatDestroy(&ctx->permutation);CHKERRQ(ierr);
1046   ierr = DMDestroy(&ctx->dmrepart);CHKERRQ(ierr);
1047   ierr = PetscFree(ctx->range_i_re);CHKERRQ(ierr);
1048   ierr = PetscFree(ctx->range_j_re);CHKERRQ(ierr);
1049   ierr = PetscFree(ctx->range_k_re);CHKERRQ(ierr);
1050   ierr = PetscFree(ctx->start_i_re);CHKERRQ(ierr);
1051   ierr = PetscFree(ctx->start_j_re);CHKERRQ(ierr);
1052   ierr = PetscFree(ctx->start_k_re);CHKERRQ(ierr);
1053   PetscFunctionReturn(0);
1054 }
1055 
1056 PetscErrorCode DMView_DA_Short_3d(DM dm,PetscViewer v)
1057 {
1058   PetscInt       M,N,P,m,n,p,ndof,nsw;
1059   MPI_Comm       comm;
1060   PetscMPIInt    size;
1061   const char*    prefix;
1062   PetscErrorCode ierr;
1063 
1064   PetscFunctionBegin;
1065   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
1066   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
1067   ierr = DMGetOptionsPrefix(dm,&prefix);CHKERRQ(ierr);
1068   ierr = DMDAGetInfo(dm,NULL,&M,&N,&P,&m,&n,&p,&ndof,&nsw,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
1069   if (prefix) {ierr = PetscViewerASCIIPrintf(v,"DMDA Object:    (%s)    %d MPI processes\n",prefix,size);CHKERRQ(ierr);}
1070   else {ierr = PetscViewerASCIIPrintf(v,"DMDA Object:    %d MPI processes\n",size);CHKERRQ(ierr);}
1071   ierr = PetscViewerASCIIPrintf(v,"  M %D N %D P %D m %D n %D p %D dof %D overlap %D\n",M,N,P,m,n,p,ndof,nsw);CHKERRQ(ierr);
1072   PetscFunctionReturn(0);
1073 }
1074 
1075 PetscErrorCode DMView_DA_Short_2d(DM dm,PetscViewer v)
1076 {
1077   PetscInt       M,N,m,n,ndof,nsw;
1078   MPI_Comm       comm;
1079   PetscMPIInt    size;
1080   const char*    prefix;
1081   PetscErrorCode ierr;
1082 
1083   PetscFunctionBegin;
1084   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
1085   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
1086   ierr = DMGetOptionsPrefix(dm,&prefix);CHKERRQ(ierr);
1087   ierr = DMDAGetInfo(dm,NULL,&M,&N,NULL,&m,&n,NULL,&ndof,&nsw,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
1088   if (prefix) {PetscViewerASCIIPrintf(v,"DMDA Object:    (%s)    %d MPI processes\n",prefix,size);CHKERRQ(ierr);}
1089   else {ierr = PetscViewerASCIIPrintf(v,"DMDA Object:    %d MPI processes\n",size);CHKERRQ(ierr);}
1090   ierr = PetscViewerASCIIPrintf(v,"  M %D N %D m %D n %D dof %D overlap %D\n",M,N,m,n,ndof,nsw);CHKERRQ(ierr);
1091   PetscFunctionReturn(0);
1092 }
1093 
1094 PetscErrorCode DMView_DA_Short(DM dm,PetscViewer v)
1095 {
1096   PetscErrorCode ierr;
1097   PetscInt       dim;
1098 
1099   PetscFunctionBegin;
1100   ierr = DMDAGetInfo(dm,&dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
1101   switch (dim) {
1102   case 2: ierr = DMView_DA_Short_2d(dm,v);CHKERRQ(ierr);
1103     break;
1104   case 3: ierr = DMView_DA_Short_3d(dm,v);CHKERRQ(ierr);
1105     break;
1106   }
1107   PetscFunctionReturn(0);
1108 }
1109 
1110