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