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