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