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