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