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