Actual source code: telescope_dmda.c


  2: #include <petsc/private/matimpl.h>
  3: #include <petsc/private/pcimpl.h>
  4: #include <petsc/private/dmimpl.h>
  5: #include <petscksp.h>
  6: #include <petscdm.h>
  7: #include <petscdmda.h>

  9: #include "../src/ksp/pc/impls/telescope/telescope.h"

 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";

 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;

 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:   }

 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:   }

 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:   }

 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: }

 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;

 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: }

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;

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));

141:   PetscCall(DMDAGlobalToNaturalBegin(cdm, coor, INSERT_VALUES, coor_natural));
142:   PetscCall(DMDAGlobalToNaturalEnd(cdm, coor, INSERT_VALUES, coor_natural));

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:   }

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:   }

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));

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));

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;

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:   }

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: }

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;

224:   PetscFunctionBegin;
225:   PetscCall(DMGetCoordinates(dm, &coor));
226:   if (!coor) PetscFunctionReturn(PETSC_SUCCESS);

228:   if (PCTelescope_isActiveRank(sred)) PetscCall(DMDASetUniformCoordinates(subdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));

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));

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));

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:   }

250:   PetscCall(PetscMalloc1(Ml * Nl * Pl * 3, &fine_indices));

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:   }

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));

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));

279:   /* access */
280:   if (PCTelescope_isActiveRank(sred)) {
281:     Vec                _coors;
282:     const PetscScalar *LA_perm;
283:     PetscScalar       *LA_coors;

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:   }

293:   /* update local coords */
294:   if (PCTelescope_isActiveRank(sred)) {
295:     DM  _dmc;
296:     Vec _coors, _coors_local;

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:   }

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: }

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;

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;

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: }

345: /* setup repartitioned dm */
346: 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;

360:   PetscFunctionBegin;
361:   comm    = PetscSubcommParent(sred->psubcomm);
362:   subcomm = PetscSubcommChild(sred->psubcomm);
363:   PetscCall(PCGetDM(pc, &dm));

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));

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));

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));

421:     ctx->dmrepart->ops->creatematrix              = dm->ops->creatematrix;
422:     ctx->dmrepart->ops->createdomaindecomposition = dm->ops->createdomaindecomposition;
423:   }

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));

432:   PetscCall(PetscCalloc3(ctx->Mp_re, &ctx->range_i_re, ctx->Np_re, &ctx->range_j_re, ctx->Pp_re, &ctx->range_k_re));

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));

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));

443:   PetscCall(PetscMalloc3(ctx->Mp_re, &ctx->start_i_re, ctx->Np_re, &ctx->start_j_re, ctx->Pp_re, &ctx->start_k_re));

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:   }

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:   }

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:   }

463:   /* attach repartitioned dm to child ksp */
464:   {
465:     PetscErrorCode (*dmksp_func)(KSP, Mat, Mat, void *);
466:     void *dmksp_ctx;

468:     PetscCall(DMKSPGetComputeOperators(dm, &dmksp_func, &dmksp_ctx));

470:     /* attach dm to ksp on sub communicator */
471:     if (PCTelescope_isActiveRank(sred)) {
472:       PetscCall(KSPSetDM(sred->ksp, ctx->dmrepart));

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: }

486: 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;

496:   PetscFunctionBegin;
497:   PetscCall(PetscInfo(pc, "PCTelescope: setting up the permutation matrix (DMDA-3D)\n"));
498:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));

500:   PetscCall(PCGetDM(pc, &dm));
501:   PetscCall(DMDAGetInfo(dm, NULL, &nx, &ny, &nz, NULL, NULL, NULL, &ndof, NULL, NULL, NULL, NULL, NULL));

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;

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));

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];

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];

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: }

557: 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;

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;

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));

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];

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];

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));

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));

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");

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: }

623: 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;

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;

663:   ctx->xp = xp;
664:   PetscCall(VecDestroy(&x));
665:   PetscFunctionReturn(PETSC_SUCCESS);
666: }

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;

675:   PetscFunctionBegin;
676:   PetscCall(PetscInfo(pc, "PCTelescope: setup (DMDA)\n"));
677:   PetscCall(PetscNew(&ctx));
678:   sred->dm_ctx = (void *)ctx;

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));

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: }

700: 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;

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;

715:   PetscCall(PCGetOperators(pc, NULL, &B));
716:   PetscCall(MatGetSize(B, &nr, &nc));

718:   P = ctx->permutation;
719:   PetscCall(MatPtAP(B, P, MAT_INITIAL_MATRIX, 1.1, &Bperm));

721:   /* Get submatrices */
722:   isrow = sred->isin;
723:   PetscCall(ISCreateStride(comm, nc, 0, 1, &iscol));

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;

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;

738:   PetscCall(ISDestroy(&iscol));
739:   PetscCall(MatDestroy(&Bperm));
740:   PetscCall(MatDestroyMatrices(1, &_Blocal));
741:   PetscFunctionReturn(PETSC_SUCCESS);
742: }

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;

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;

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: }

778: 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;

787:   PetscFunctionBegin;
788:   ctx     = (PC_Telescope_DMDACtx *)sred->dm_ctx;
789:   subcomm = PetscSubcommChild(sred->psubcomm);
790:   PetscCall(MatNullSpaceGetVecs(nullspace, &has_const, &n, &vecs));

792:   if (PCTelescope_isActiveRank(sred)) {
793:     /* create new vectors */
794:     if (n) PetscCall(VecDuplicateVecs(sred->xred, n, &sub_vecs));
795:   }

797:   /* copy entries */
798:   for (k = 0; k < n; k++) {
799:     const PetscScalar *x_array;
800:     PetscScalar       *LA_sub_vec;
801:     PetscInt           st, ed;

803:     /* permute vector into ordering associated with re-partitioned dmda */
804:     PetscCall(MatMultTranspose(ctx->permutation, vecs[k], ctx->xp));

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));

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:   }

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: }

833: PetscErrorCode PCTelescopeMatNullSpaceCreate_dmda(PC pc, PC_Telescope sred, Mat sub_mat)
834: {
835:   Mat B;

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: }

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;

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;

885:   PetscFunctionBegin;
886:   PetscCall(PetscCitationsRegister(citation, &cited));

888:   /* permute vector into ordering associated with re-partitioned dmda */
889:   PetscCall(MatMultTranspose(perm, x, xp));

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));

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));

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));

907:   /* solve */
908:   if (PCTelescope_isActiveRank(sred)) {
909:     PetscCall(KSPSolve(sred->ksp, xred, yred));
910:     PetscCall(KSPCheckSolve(sred->ksp, pc, yred));
911:   }

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: }

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;

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;

948:   PetscCheck(its <= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCApplyRichardson_Telescope_dmda only supports max_it = 1");
949:   *reason = (PCRichardsonConvergedReason)0;

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));

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));

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:   }

972:   if (PCTelescope_isActiveRank(sred)) {
973:     PetscCall(KSPGetInitialGuessNonzero(sred->ksp, &default_init_guess_value));
974:     if (!zeroguess) PetscCall(KSPSetInitialGuessNonzero(sred->ksp, PETSC_TRUE));
975:   }

977:   PetscCall(PCApply_Telescope_dmda(pc, x, y));

979:   if (PCTelescope_isActiveRank(sred)) PetscCall(KSPSetInitialGuessNonzero(sred->ksp, default_init_guess_value));

981:   if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
982:   *outits = 1;
983:   PetscFunctionReturn(PETSC_SUCCESS);
984: }

986: PetscErrorCode PCReset_Telescope_dmda(PC pc)
987: {
988:   PC_Telescope          sred = (PC_Telescope)pc->data;
989:   PC_Telescope_DMDACtx *ctx;

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: }

1001: 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;

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: }

1019: 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;

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: }

1037: PetscErrorCode DMView_DA_Short(DM dm, PetscViewer v)
1038: {
1039:   PetscInt dim;

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: }