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