#include #include #include #include /*I "petscksp.h" I*/ #include #include #include "../src/ksp/pc/impls/telescope/telescope.h" static PetscBool cited = PETSC_FALSE; static const char citation[] = "@inproceedings{MaySananRuppKnepleySmith2016,\n" " title = {Extreme-Scale Multigrid Components within PETSc},\n" " author = {Dave A. May and Patrick Sanan and Karl Rupp and Matthew G. Knepley and Barry F. Smith},\n" " booktitle = {Proceedings of the Platform for Advanced Scientific Computing Conference},\n" " series = {PASC '16},\n" " isbn = {978-1-4503-4126-4},\n" " location = {Lausanne, Switzerland},\n" " pages = {5:1--5:12},\n" " articleno = {5},\n" " numpages = {12},\n" " url = {https://doi.acm.org/10.1145/2929908.2929913},\n" " doi = {10.1145/2929908.2929913},\n" " acmid = {2929913},\n" " publisher = {ACM},\n" " address = {New York, NY, USA},\n" " keywords = {GPU, HPC, agglomeration, coarse-level solver, multigrid, parallel computing, preconditioning},\n" " year = {2016}\n" "}\n"; 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) { PetscInt pi, pj, pk, n; PetscFunctionBegin; *rank_re = -1; if (_pi) *_pi = -1; if (_pj) *_pj = -1; if (_pk) *_pk = -1; pi = pj = pk = -1; if (_pi) { for (n = 0; n < Mp; n++) { if ((i >= start_i[n]) && (i < start_i[n] + span_i[n])) { pi = n; break; } } PetscCheck(pi != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmda-ijk] pi cannot be determined : range %" PetscInt_FMT ", val %" PetscInt_FMT, Mp, i); PetscCall(PetscMPIIntCast(pi, _pi)); } if (_pj) { for (n = 0; n < Np; n++) { if ((j >= start_j[n]) && (j < start_j[n] + span_j[n])) { pj = n; break; } } PetscCheck(pj != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmda-ijk] pj cannot be determined : range %" PetscInt_FMT ", val %" PetscInt_FMT, Np, j); PetscCall(PetscMPIIntCast(pj, _pj)); } if (_pk) { for (n = 0; n < Pp; n++) { if ((k >= start_k[n]) && (k < start_k[n] + span_k[n])) { pk = n; break; } } PetscCheck(pk != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmda-ijk] pk cannot be determined : range %" PetscInt_FMT ", val %" PetscInt_FMT, Pp, k); PetscCall(PetscMPIIntCast(pk, _pk)); } switch (dim) { case 1: PetscCall(PetscMPIIntCast(pi, rank_re)); break; case 2: PetscCall(PetscMPIIntCast(pi + pj * Mp, rank_re)); break; case 3: PetscCall(PetscMPIIntCast(pi + pj * Mp + pk * (Mp * Np), rank_re)); break; } PetscFunctionReturn(PETSC_SUCCESS); } 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) { PetscInt i, j, k, start_IJK = 0; PetscInt rank_ijk; PetscFunctionBegin; switch (dim) { case 1: for (i = 0; i < Mp_re; i++) { rank_ijk = i; if (rank_ijk < rank_re) start_IJK += range_i_re[i]; } break; case 2: for (j = 0; j < Np_re; j++) { for (i = 0; i < Mp_re; i++) { rank_ijk = i + j * Mp_re; if (rank_ijk < rank_re) start_IJK += range_i_re[i] * range_j_re[j]; } } break; case 3: for (k = 0; k < Pp_re; k++) { for (j = 0; j < Np_re; j++) { for (i = 0; i < Mp_re; i++) { rank_ijk = i + j * Mp_re + k * Mp_re * Np_re; if (rank_ijk < rank_re) start_IJK += range_i_re[i] * range_j_re[j] * range_k_re[k]; } } } break; } *s0 = start_IJK; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCTelescopeSetUp_dmda_repart_coors2d(PC_Telescope sred, DM dm, DM subdm) { DM cdm; Vec coor, coor_natural, perm_coors; PetscInt i, j, si, sj, ni, nj, M, N, Ml, Nl, c, nidx; PetscInt *fine_indices; IS is_fine, is_local; VecScatter sctx; PetscFunctionBegin; PetscCall(DMGetCoordinates(dm, &coor)); if (!coor) PetscFunctionReturn(PETSC_SUCCESS); if (PCTelescope_isActiveRank(sred)) PetscCall(DMDASetUniformCoordinates(subdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); /* Get the coordinate vector from the distributed array */ PetscCall(DMGetCoordinateDM(dm, &cdm)); PetscCall(DMDACreateNaturalVector(cdm, &coor_natural)); PetscCall(DMDAGlobalToNaturalBegin(cdm, coor, INSERT_VALUES, coor_natural)); PetscCall(DMDAGlobalToNaturalEnd(cdm, coor, INSERT_VALUES, coor_natural)); /* get indices of the guys I want to grab */ PetscCall(DMDAGetInfo(dm, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); if (PCTelescope_isActiveRank(sred)) { PetscCall(DMDAGetCorners(subdm, &si, &sj, NULL, &ni, &nj, NULL)); Ml = ni; Nl = nj; } else { si = sj = 0; ni = nj = 0; Ml = Nl = 0; } PetscCall(PetscMalloc1(Ml * Nl * 2, &fine_indices)); c = 0; if (PCTelescope_isActiveRank(sred)) { for (j = sj; j < sj + nj; j++) { for (i = si; i < si + ni; i++) { nidx = (i) + (j)*M; fine_indices[c] = 2 * nidx; fine_indices[c + 1] = 2 * nidx + 1; c = c + 2; } } 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); } /* generate scatter */ PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), Ml * Nl * 2, fine_indices, PETSC_USE_POINTER, &is_fine)); PetscCall(ISCreateStride(PETSC_COMM_SELF, Ml * Nl * 2, 0, 1, &is_local)); /* scatter */ PetscCall(VecCreate(PETSC_COMM_SELF, &perm_coors)); PetscCall(VecSetSizes(perm_coors, PETSC_DECIDE, Ml * Nl * 2)); PetscCall(VecSetType(perm_coors, VECSEQ)); PetscCall(VecScatterCreate(coor_natural, is_fine, perm_coors, is_local, &sctx)); PetscCall(VecScatterBegin(sctx, coor_natural, perm_coors, INSERT_VALUES, SCATTER_FORWARD)); PetscCall(VecScatterEnd(sctx, coor_natural, perm_coors, INSERT_VALUES, SCATTER_FORWARD)); /* access */ if (PCTelescope_isActiveRank(sred)) { Vec _coors; const PetscScalar *LA_perm; PetscScalar *LA_coors; PetscCall(DMGetCoordinates(subdm, &_coors)); PetscCall(VecGetArrayRead(perm_coors, &LA_perm)); PetscCall(VecGetArray(_coors, &LA_coors)); for (i = 0; i < Ml * Nl * 2; i++) LA_coors[i] = LA_perm[i]; PetscCall(VecRestoreArray(_coors, &LA_coors)); PetscCall(VecRestoreArrayRead(perm_coors, &LA_perm)); } /* update local coords */ if (PCTelescope_isActiveRank(sred)) { DM _dmc; Vec _coors, _coors_local; PetscCall(DMGetCoordinateDM(subdm, &_dmc)); PetscCall(DMGetCoordinates(subdm, &_coors)); PetscCall(DMGetCoordinatesLocal(subdm, &_coors_local)); PetscCall(DMGlobalToLocalBegin(_dmc, _coors, INSERT_VALUES, _coors_local)); PetscCall(DMGlobalToLocalEnd(_dmc, _coors, INSERT_VALUES, _coors_local)); } PetscCall(VecScatterDestroy(&sctx)); PetscCall(ISDestroy(&is_fine)); PetscCall(PetscFree(fine_indices)); PetscCall(ISDestroy(&is_local)); PetscCall(VecDestroy(&perm_coors)); PetscCall(VecDestroy(&coor_natural)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCTelescopeSetUp_dmda_repart_coors3d(PC_Telescope sred, DM dm, DM subdm) { DM cdm; Vec coor, coor_natural, perm_coors; PetscInt i, j, k, si, sj, sk, ni, nj, nk, M, N, P, Ml, Nl, Pl, c, nidx; PetscInt *fine_indices; IS is_fine, is_local; VecScatter sctx; PetscFunctionBegin; PetscCall(DMGetCoordinates(dm, &coor)); if (!coor) PetscFunctionReturn(PETSC_SUCCESS); if (PCTelescope_isActiveRank(sred)) PetscCall(DMDASetUniformCoordinates(subdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); /* Get the coordinate vector from the distributed array */ PetscCall(DMGetCoordinateDM(dm, &cdm)); PetscCall(DMDACreateNaturalVector(cdm, &coor_natural)); PetscCall(DMDAGlobalToNaturalBegin(cdm, coor, INSERT_VALUES, coor_natural)); PetscCall(DMDAGlobalToNaturalEnd(cdm, coor, INSERT_VALUES, coor_natural)); /* get indices of the guys I want to grab */ PetscCall(DMDAGetInfo(dm, NULL, &M, &N, &P, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); if (PCTelescope_isActiveRank(sred)) { PetscCall(DMDAGetCorners(subdm, &si, &sj, &sk, &ni, &nj, &nk)); Ml = ni; Nl = nj; Pl = nk; } else { si = sj = sk = 0; ni = nj = nk = 0; Ml = Nl = Pl = 0; } PetscCall(PetscMalloc1(Ml * Nl * Pl * 3, &fine_indices)); c = 0; if (PCTelescope_isActiveRank(sred)) { for (k = sk; k < sk + nk; k++) { for (j = sj; j < sj + nj; j++) { for (i = si; i < si + ni; i++) { nidx = (i) + (j)*M + (k)*M * N; fine_indices[c] = 3 * nidx; fine_indices[c + 1] = 3 * nidx + 1; fine_indices[c + 2] = 3 * nidx + 2; c = c + 3; } } } } /* generate scatter */ PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), Ml * Nl * Pl * 3, fine_indices, PETSC_USE_POINTER, &is_fine)); PetscCall(ISCreateStride(PETSC_COMM_SELF, Ml * Nl * Pl * 3, 0, 1, &is_local)); /* scatter */ PetscCall(VecCreate(PETSC_COMM_SELF, &perm_coors)); PetscCall(VecSetSizes(perm_coors, PETSC_DECIDE, Ml * Nl * Pl * 3)); PetscCall(VecSetType(perm_coors, VECSEQ)); PetscCall(VecScatterCreate(coor_natural, is_fine, perm_coors, is_local, &sctx)); PetscCall(VecScatterBegin(sctx, coor_natural, perm_coors, INSERT_VALUES, SCATTER_FORWARD)); PetscCall(VecScatterEnd(sctx, coor_natural, perm_coors, INSERT_VALUES, SCATTER_FORWARD)); /* access */ if (PCTelescope_isActiveRank(sred)) { Vec _coors; const PetscScalar *LA_perm; PetscScalar *LA_coors; PetscCall(DMGetCoordinates(subdm, &_coors)); PetscCall(VecGetArrayRead(perm_coors, &LA_perm)); PetscCall(VecGetArray(_coors, &LA_coors)); for (i = 0; i < Ml * Nl * Pl * 3; i++) LA_coors[i] = LA_perm[i]; PetscCall(VecRestoreArray(_coors, &LA_coors)); PetscCall(VecRestoreArrayRead(perm_coors, &LA_perm)); } /* update local coords */ if (PCTelescope_isActiveRank(sred)) { DM _dmc; Vec _coors, _coors_local; PetscCall(DMGetCoordinateDM(subdm, &_dmc)); PetscCall(DMGetCoordinates(subdm, &_coors)); PetscCall(DMGetCoordinatesLocal(subdm, &_coors_local)); PetscCall(DMGlobalToLocalBegin(_dmc, _coors, INSERT_VALUES, _coors_local)); PetscCall(DMGlobalToLocalEnd(_dmc, _coors, INSERT_VALUES, _coors_local)); } PetscCall(VecScatterDestroy(&sctx)); PetscCall(ISDestroy(&is_fine)); PetscCall(PetscFree(fine_indices)); PetscCall(ISDestroy(&is_local)); PetscCall(VecDestroy(&perm_coors)); PetscCall(VecDestroy(&coor_natural)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCTelescopeSetUp_dmda_repart_coors(PC pc, PC_Telescope sred, PC_Telescope_DMDACtx *ctx) { PetscInt dim; DM dm, subdm; PetscSubcomm psubcomm; MPI_Comm comm; Vec coor; PetscFunctionBegin; PetscCall(PCGetDM(pc, &dm)); PetscCall(DMGetCoordinates(dm, &coor)); if (!coor) PetscFunctionReturn(PETSC_SUCCESS); psubcomm = sred->psubcomm; comm = PetscSubcommParent(psubcomm); subdm = ctx->dmrepart; PetscCall(PetscInfo(pc, "PCTelescope: setting up the coordinates (DMDA)\n")); PetscCall(DMDAGetInfo(dm, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); switch (dim) { case 1: SETERRQ(comm, PETSC_ERR_SUP, "Telescope: DMDA (1D) repartitioning not provided"); case 2: PetscCall(PCTelescopeSetUp_dmda_repart_coors2d(sred, dm, subdm)); break; case 3: PetscCall(PCTelescopeSetUp_dmda_repart_coors3d(sred, dm, subdm)); break; } PetscFunctionReturn(PETSC_SUCCESS); } /* setup repartitioned dm */ static PetscErrorCode PCTelescopeSetUp_dmda_repart(PC pc, PC_Telescope sred, PC_Telescope_DMDACtx *ctx) { DM dm; PetscInt dim, nx, ny, nz, ndof, nsw, sum, k; DMBoundaryType bx, by, bz; DMDAStencilType stencil; const PetscInt *_range_i_re; const PetscInt *_range_j_re; const PetscInt *_range_k_re; DMDAInterpolationType itype; PetscInt refine_x, refine_y, refine_z; MPI_Comm comm, subcomm; const char *prefix; PetscMPIInt ni; PetscFunctionBegin; comm = PetscSubcommParent(sred->psubcomm); subcomm = PetscSubcommChild(sred->psubcomm); PetscCall(PCGetDM(pc, &dm)); PetscCall(DMDAGetInfo(dm, &dim, &nx, &ny, &nz, NULL, NULL, NULL, &ndof, &nsw, &bx, &by, &bz, &stencil)); PetscCall(DMDAGetInterpolationType(dm, &itype)); PetscCall(DMDAGetRefinementFactor(dm, &refine_x, &refine_y, &refine_z)); ctx->dmrepart = NULL; _range_i_re = _range_j_re = _range_k_re = NULL; /* Create DMDA on the child communicator */ if (PCTelescope_isActiveRank(sred)) { switch (dim) { case 1: PetscCall(PetscInfo(pc, "PCTelescope: setting up the DMDA on comm subset (1D)\n")); /* PetscCall(DMDACreate1d(subcomm,bx,nx,ndof,nsw,NULL,&ctx->dmrepart)); */ ny = nz = 1; by = bz = DM_BOUNDARY_NONE; break; case 2: PetscCall(PetscInfo(pc, "PCTelescope: setting up the DMDA on comm subset (2D)\n")); /* PetscCall(DMDACreate2d(subcomm,bx,by,stencil,nx,ny, PETSC_DECIDE,PETSC_DECIDE, ndof,nsw, NULL,NULL,&ctx->dmrepart)); */ nz = 1; bz = DM_BOUNDARY_NONE; break; case 3: PetscCall(PetscInfo(pc, "PCTelescope: setting up the DMDA on comm subset (3D)\n")); /* PetscCall(DMDACreate3d(subcomm,bx,by,bz,stencil,nx,ny,nz, PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, ndof,nsw, NULL,NULL,NULL,&ctx->dmrepart)); */ break; } /* The API DMDACreate1d(), DMDACreate2d(), DMDACreate3d() does not allow us to set/append a unique option prefix for the DM, thus I prefer to expose the contents of these API's here. This allows users to control the partitioning of the subDM. */ PetscCall(DMDACreate(subcomm, &ctx->dmrepart)); /* Set unique option prefix name */ PetscCall(KSPGetOptionsPrefix(sred->ksp, &prefix)); PetscCall(DMSetOptionsPrefix(ctx->dmrepart, prefix)); PetscCall(DMAppendOptionsPrefix(ctx->dmrepart, "repart_")); /* standard setup from DMDACreate{1,2,3}d() */ PetscCall(DMSetDimension(ctx->dmrepart, dim)); PetscCall(DMDASetSizes(ctx->dmrepart, nx, ny, nz)); PetscCall(DMDASetNumProcs(ctx->dmrepart, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE)); PetscCall(DMDASetBoundaryType(ctx->dmrepart, bx, by, bz)); PetscCall(DMDASetDof(ctx->dmrepart, ndof)); PetscCall(DMDASetStencilType(ctx->dmrepart, stencil)); PetscCall(DMDASetStencilWidth(ctx->dmrepart, nsw)); PetscCall(DMDASetOwnershipRanges(ctx->dmrepart, NULL, NULL, NULL)); PetscCall(DMSetFromOptions(ctx->dmrepart)); PetscCall(DMSetUp(ctx->dmrepart)); /* Set refinement factors and interpolation type from the partent */ PetscCall(DMDASetRefinementFactor(ctx->dmrepart, refine_x, refine_y, refine_z)); PetscCall(DMDASetInterpolationType(ctx->dmrepart, itype)); PetscCall(DMDAGetInfo(ctx->dmrepart, NULL, NULL, NULL, NULL, &ctx->Mp_re, &ctx->Np_re, &ctx->Pp_re, NULL, NULL, NULL, NULL, NULL, NULL)); PetscCall(DMDAGetOwnershipRanges(ctx->dmrepart, &_range_i_re, &_range_j_re, &_range_k_re)); ctx->dmrepart->ops->creatematrix = dm->ops->creatematrix; ctx->dmrepart->ops->createdomaindecomposition = dm->ops->createdomaindecomposition; } /* generate ranges for repartitioned dm */ /* note - assume rank 0 always participates */ /* TODO: use a single MPI call */ PetscCallMPI(MPI_Bcast(&ctx->Mp_re, 1, MPIU_INT, 0, comm)); PetscCallMPI(MPI_Bcast(&ctx->Np_re, 1, MPIU_INT, 0, comm)); PetscCallMPI(MPI_Bcast(&ctx->Pp_re, 1, MPIU_INT, 0, comm)); PetscCall(PetscCalloc3(ctx->Mp_re, &ctx->range_i_re, ctx->Np_re, &ctx->range_j_re, ctx->Pp_re, &ctx->range_k_re)); if (_range_i_re) PetscCall(PetscArraycpy(ctx->range_i_re, _range_i_re, ctx->Mp_re)); if (_range_j_re) PetscCall(PetscArraycpy(ctx->range_j_re, _range_j_re, ctx->Np_re)); if (_range_k_re) PetscCall(PetscArraycpy(ctx->range_k_re, _range_k_re, ctx->Pp_re)); /* TODO: use a single MPI call */ PetscCall(PetscMPIIntCast(ctx->Mp_re, &ni)); PetscCallMPI(MPI_Bcast(ctx->range_i_re, ni, MPIU_INT, 0, comm)); PetscCall(PetscMPIIntCast(ctx->Np_re, &ni)); PetscCallMPI(MPI_Bcast(ctx->range_j_re, ni, MPIU_INT, 0, comm)); PetscCall(PetscMPIIntCast(ctx->Pp_re, &ni)); PetscCallMPI(MPI_Bcast(ctx->range_k_re, ni, MPIU_INT, 0, comm)); PetscCall(PetscMalloc3(ctx->Mp_re, &ctx->start_i_re, ctx->Np_re, &ctx->start_j_re, ctx->Pp_re, &ctx->start_k_re)); sum = 0; for (k = 0; k < ctx->Mp_re; k++) { ctx->start_i_re[k] = sum; sum += ctx->range_i_re[k]; } sum = 0; for (k = 0; k < ctx->Np_re; k++) { ctx->start_j_re[k] = sum; sum += ctx->range_j_re[k]; } sum = 0; for (k = 0; k < ctx->Pp_re; k++) { ctx->start_k_re[k] = sum; sum += ctx->range_k_re[k]; } /* attach repartitioned dm to child ksp */ { PetscErrorCode (*dmksp_func)(KSP, Mat, Mat, void *); void *dmksp_ctx; PetscCall(DMKSPGetComputeOperators(dm, &dmksp_func, &dmksp_ctx)); /* attach dm to ksp on sub communicator */ if (PCTelescope_isActiveRank(sred)) { PetscCall(KSPSetDM(sred->ksp, ctx->dmrepart)); if (!dmksp_func || sred->ignore_kspcomputeoperators) { PetscCall(KSPSetDMActive(sred->ksp, KSP_DMACTIVE_ALL, PETSC_FALSE)); } else { /* sub ksp inherits dmksp_func and context provided by user */ PetscCall(KSPSetComputeOperators(sred->ksp, dmksp_func, dmksp_ctx)); PetscCall(KSPSetDMActive(sred->ksp, KSP_DMACTIVE_ALL, PETSC_TRUE)); } } } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCTelescopeSetUp_dmda_permutation_3d(PC pc, PC_Telescope sred, PC_Telescope_DMDACtx *ctx) { DM dm; MPI_Comm comm; Mat Pscalar, P; PetscInt ndof; PetscInt i, j, k, location, startI[3], endI[3], lenI[3], nx, ny, nz; PetscInt sr, er, Mr; Vec V; PetscFunctionBegin; PetscCall(PetscInfo(pc, "PCTelescope: setting up the permutation matrix (DMDA-3D)\n")); PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); PetscCall(PCGetDM(pc, &dm)); PetscCall(DMDAGetInfo(dm, NULL, &nx, &ny, &nz, NULL, NULL, NULL, &ndof, NULL, NULL, NULL, NULL, NULL)); PetscCall(DMGetGlobalVector(dm, &V)); PetscCall(VecGetSize(V, &Mr)); PetscCall(VecGetOwnershipRange(V, &sr, &er)); PetscCall(DMRestoreGlobalVector(dm, &V)); sr = sr / ndof; er = er / ndof; Mr = Mr / ndof; PetscCall(MatCreate(comm, &Pscalar)); PetscCall(MatSetSizes(Pscalar, er - sr, er - sr, Mr, Mr)); PetscCall(MatSetType(Pscalar, MATAIJ)); PetscCall(MatSeqAIJSetPreallocation(Pscalar, 1, NULL)); PetscCall(MatMPIAIJSetPreallocation(Pscalar, 1, NULL, 1, NULL)); PetscCall(DMDAGetCorners(dm, NULL, NULL, NULL, &lenI[0], &lenI[1], &lenI[2])); PetscCall(DMDAGetCorners(dm, &startI[0], &startI[1], &startI[2], &endI[0], &endI[1], &endI[2])); endI[0] += startI[0]; endI[1] += startI[1]; endI[2] += startI[2]; for (k = startI[2]; k < endI[2]; k++) { for (j = startI[1]; j < endI[1]; j++) { for (i = startI[0]; i < endI[0]; i++) { PetscMPIInt rank_ijk_re, rank_reI[3]; PetscInt s0_re; PetscInt ii, jj, kk, local_ijk_re, mapped_ijk; PetscInt lenI_re[3]; location = (i - startI[0]) + (j - startI[1]) * lenI[0] + (k - startI[2]) * lenI[0] * lenI[1]; 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)); 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)); ii = i - ctx->start_i_re[rank_reI[0]]; PetscCheck(ii >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmdarepart-perm3d] index error ii"); jj = j - ctx->start_j_re[rank_reI[1]]; PetscCheck(jj >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmdarepart-perm3d] index error jj"); kk = k - ctx->start_k_re[rank_reI[2]]; PetscCheck(kk >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmdarepart-perm3d] index error kk"); lenI_re[0] = ctx->range_i_re[rank_reI[0]]; lenI_re[1] = ctx->range_j_re[rank_reI[1]]; lenI_re[2] = ctx->range_k_re[rank_reI[2]]; local_ijk_re = ii + jj * lenI_re[0] + kk * lenI_re[0] * lenI_re[1]; mapped_ijk = s0_re + local_ijk_re; PetscCall(MatSetValue(Pscalar, sr + location, mapped_ijk, 1.0, INSERT_VALUES)); } } } PetscCall(MatAssemblyBegin(Pscalar, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(Pscalar, MAT_FINAL_ASSEMBLY)); PetscCall(MatCreateMAIJ(Pscalar, ndof, &P)); PetscCall(MatDestroy(&Pscalar)); ctx->permutation = P; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCTelescopeSetUp_dmda_permutation_2d(PC pc, PC_Telescope sred, PC_Telescope_DMDACtx *ctx) { DM dm; MPI_Comm comm; Mat Pscalar, P; PetscInt ndof; PetscInt i, j, location, startI[2], endI[2], lenI[2], nx, ny, nz; PetscInt sr, er, Mr; Vec V; PetscFunctionBegin; PetscCall(PetscInfo(pc, "PCTelescope: setting up the permutation matrix (DMDA-2D)\n")); PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); PetscCall(PCGetDM(pc, &dm)); PetscCall(DMDAGetInfo(dm, NULL, &nx, &ny, &nz, NULL, NULL, NULL, &ndof, NULL, NULL, NULL, NULL, NULL)); PetscCall(DMGetGlobalVector(dm, &V)); PetscCall(VecGetSize(V, &Mr)); PetscCall(VecGetOwnershipRange(V, &sr, &er)); PetscCall(DMRestoreGlobalVector(dm, &V)); sr = sr / ndof; er = er / ndof; Mr = Mr / ndof; PetscCall(MatCreate(comm, &Pscalar)); PetscCall(MatSetSizes(Pscalar, er - sr, er - sr, Mr, Mr)); PetscCall(MatSetType(Pscalar, MATAIJ)); PetscCall(MatSeqAIJSetPreallocation(Pscalar, 1, NULL)); PetscCall(MatMPIAIJSetPreallocation(Pscalar, 1, NULL, 1, NULL)); PetscCall(DMDAGetCorners(dm, NULL, NULL, NULL, &lenI[0], &lenI[1], NULL)); PetscCall(DMDAGetCorners(dm, &startI[0], &startI[1], NULL, &endI[0], &endI[1], NULL)); endI[0] += startI[0]; endI[1] += startI[1]; for (j = startI[1]; j < endI[1]; j++) { for (i = startI[0]; i < endI[0]; i++) { PetscMPIInt rank_ijk_re, rank_reI[3]; PetscInt s0_re; PetscInt ii, jj, local_ijk_re, mapped_ijk; PetscInt lenI_re[3]; location = (i - startI[0]) + (j - startI[1]) * lenI[0]; 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)); 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)); ii = i - ctx->start_i_re[rank_reI[0]]; PetscCheck(ii >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmdarepart-perm2d] index error ii"); jj = j - ctx->start_j_re[rank_reI[1]]; PetscCheck(jj >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmdarepart-perm2d] index error jj"); lenI_re[0] = ctx->range_i_re[rank_reI[0]]; lenI_re[1] = ctx->range_j_re[rank_reI[1]]; local_ijk_re = ii + jj * lenI_re[0]; mapped_ijk = s0_re + local_ijk_re; PetscCall(MatSetValue(Pscalar, sr + location, mapped_ijk, 1.0, INSERT_VALUES)); } } PetscCall(MatAssemblyBegin(Pscalar, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(Pscalar, MAT_FINAL_ASSEMBLY)); PetscCall(MatCreateMAIJ(Pscalar, ndof, &P)); PetscCall(MatDestroy(&Pscalar)); ctx->permutation = P; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCTelescopeSetUp_dmda_scatters(PC pc, PC_Telescope sred, PC_Telescope_DMDACtx *ctx) { Vec xred, yred, xtmp, x, xp; VecScatter scatter; IS isin; Mat B; PetscInt m, bs, st, ed; MPI_Comm comm; PetscFunctionBegin; PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); PetscCall(PCGetOperators(pc, NULL, &B)); PetscCall(MatCreateVecs(B, &x, NULL)); PetscCall(MatGetBlockSize(B, &bs)); PetscCall(VecDuplicate(x, &xp)); m = 0; xred = NULL; yred = NULL; if (PCTelescope_isActiveRank(sred)) { PetscCall(DMCreateGlobalVector(ctx->dmrepart, &xred)); PetscCall(VecDuplicate(xred, &yred)); PetscCall(VecGetOwnershipRange(xred, &st, &ed)); PetscCall(ISCreateStride(comm, ed - st, st, 1, &isin)); PetscCall(VecGetLocalSize(xred, &m)); } else { PetscCall(VecGetOwnershipRange(x, &st, &ed)); PetscCall(ISCreateStride(comm, 0, st, 1, &isin)); } PetscCall(ISSetBlockSize(isin, bs)); PetscCall(VecCreate(comm, &xtmp)); PetscCall(VecSetSizes(xtmp, m, PETSC_DECIDE)); PetscCall(VecSetBlockSize(xtmp, bs)); PetscCall(VecSetType(xtmp, ((PetscObject)x)->type_name)); PetscCall(VecScatterCreate(x, isin, xtmp, NULL, &scatter)); sred->xred = xred; sred->yred = yred; sred->isin = isin; sred->scatter = scatter; sred->xtmp = xtmp; ctx->xp = xp; PetscCall(VecDestroy(&x)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode PCTelescopeSetUp_dmda(PC pc, PC_Telescope sred) { PC_Telescope_DMDACtx *ctx; PetscInt dim; DM dm; MPI_Comm comm; PetscFunctionBegin; PetscCall(PetscInfo(pc, "PCTelescope: setup (DMDA)\n")); PetscCall(PetscNew(&ctx)); sred->dm_ctx = (void *)ctx; PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); PetscCall(PCGetDM(pc, &dm)); PetscCall(DMDAGetInfo(dm, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); PetscCall(PCTelescopeSetUp_dmda_repart(pc, sred, ctx)); PetscCall(PCTelescopeSetUp_dmda_repart_coors(pc, sred, ctx)); switch (dim) { case 1: SETERRQ(comm, PETSC_ERR_SUP, "Telescope: DMDA (1D) repartitioning not provided"); case 2: PetscCall(PCTelescopeSetUp_dmda_permutation_2d(pc, sred, ctx)); break; case 3: PetscCall(PCTelescopeSetUp_dmda_permutation_3d(pc, sred, ctx)); break; } PetscCall(PCTelescopeSetUp_dmda_scatters(pc, sred, ctx)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCTelescopeMatCreate_dmda_dmactivefalse(PC pc, PC_Telescope sred, MatReuse reuse, Mat *A) { PC_Telescope_DMDACtx *ctx; MPI_Comm comm, subcomm; Mat Bperm, Bred, B, P; PetscInt nr, nc; IS isrow, iscol; Mat Blocal, *_Blocal; PetscFunctionBegin; PetscCall(PetscInfo(pc, "PCTelescope: updating the redundant preconditioned operator (DMDA)\n")); PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); subcomm = PetscSubcommChild(sred->psubcomm); ctx = (PC_Telescope_DMDACtx *)sred->dm_ctx; PetscCall(PCGetOperators(pc, NULL, &B)); PetscCall(MatGetSize(B, &nr, &nc)); P = ctx->permutation; PetscCall(MatPtAP(B, P, MAT_INITIAL_MATRIX, 1.1, &Bperm)); /* Get submatrices */ isrow = sred->isin; PetscCall(ISCreateStride(comm, nc, 0, 1, &iscol)); PetscCall(MatCreateSubMatrices(Bperm, 1, &isrow, &iscol, MAT_INITIAL_MATRIX, &_Blocal)); Blocal = *_Blocal; Bred = NULL; if (PCTelescope_isActiveRank(sred)) { PetscInt mm; if (reuse != MAT_INITIAL_MATRIX) Bred = *A; PetscCall(MatGetSize(Blocal, &mm, NULL)); /* PetscCall(MatCreateMPIMatConcatenateSeqMat(subcomm,Blocal,PETSC_DECIDE,reuse,&Bred)); */ PetscCall(MatCreateMPIMatConcatenateSeqMat(subcomm, Blocal, mm, reuse, &Bred)); } *A = Bred; PetscCall(ISDestroy(&iscol)); PetscCall(MatDestroy(&Bperm)); PetscCall(MatDestroyMatrices(1, &_Blocal)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode PCTelescopeMatCreate_dmda(PC pc, PC_Telescope sred, MatReuse reuse, Mat *A) { DM dm; PetscErrorCode (*dmksp_func)(KSP, Mat, Mat, void *); void *dmksp_ctx; PetscFunctionBegin; PetscCall(PCGetDM(pc, &dm)); PetscCall(DMKSPGetComputeOperators(dm, &dmksp_func, &dmksp_ctx)); /* We assume that dmksp_func = NULL, is equivalent to dmActive = PETSC_FALSE */ if (dmksp_func && !sred->ignore_kspcomputeoperators) { DM dmrepart; Mat Ak; *A = NULL; if (PCTelescope_isActiveRank(sred)) { PetscCall(KSPGetDM(sred->ksp, &dmrepart)); if (reuse == MAT_INITIAL_MATRIX) { PetscCall(DMCreateMatrix(dmrepart, &Ak)); *A = Ak; } else if (reuse == MAT_REUSE_MATRIX) { Ak = *A; } /* There is no need to explicitly assemble the operator now, the sub-KSP will call the method provided to KSPSetComputeOperators() during KSPSetUp() */ } } else { PetscCall(PCTelescopeMatCreate_dmda_dmactivefalse(pc, sred, reuse, A)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCTelescopeSubNullSpaceCreate_dmda_Telescope(PC pc, PC_Telescope sred, MatNullSpace nullspace, MatNullSpace *sub_nullspace) { PetscBool has_const; PetscInt i, k, n = 0; const Vec *vecs; Vec *sub_vecs = NULL; MPI_Comm subcomm; PC_Telescope_DMDACtx *ctx; PetscFunctionBegin; ctx = (PC_Telescope_DMDACtx *)sred->dm_ctx; subcomm = PetscSubcommChild(sred->psubcomm); PetscCall(MatNullSpaceGetVecs(nullspace, &has_const, &n, &vecs)); if (PCTelescope_isActiveRank(sred)) { /* create new vectors */ if (n) PetscCall(VecDuplicateVecs(sred->xred, n, &sub_vecs)); } /* copy entries */ for (k = 0; k < n; k++) { const PetscScalar *x_array; PetscScalar *LA_sub_vec; PetscInt st, ed; /* permute vector into ordering associated with re-partitioned dmda */ PetscCall(MatMultTranspose(ctx->permutation, vecs[k], ctx->xp)); /* pull in vector x->xtmp */ PetscCall(VecScatterBegin(sred->scatter, ctx->xp, sred->xtmp, INSERT_VALUES, SCATTER_FORWARD)); PetscCall(VecScatterEnd(sred->scatter, ctx->xp, sred->xtmp, INSERT_VALUES, SCATTER_FORWARD)); /* copy vector entries into xred */ PetscCall(VecGetArrayRead(sred->xtmp, &x_array)); if (sub_vecs) { if (sub_vecs[k]) { PetscCall(VecGetOwnershipRange(sub_vecs[k], &st, &ed)); PetscCall(VecGetArray(sub_vecs[k], &LA_sub_vec)); for (i = 0; i < ed - st; i++) LA_sub_vec[i] = x_array[i]; PetscCall(VecRestoreArray(sub_vecs[k], &LA_sub_vec)); } } PetscCall(VecRestoreArrayRead(sred->xtmp, &x_array)); } if (PCTelescope_isActiveRank(sred)) { /* create new (near) nullspace for redundant object */ PetscCall(MatNullSpaceCreate(subcomm, has_const, n, sub_vecs, sub_nullspace)); PetscCall(VecDestroyVecs(n, &sub_vecs)); PetscCheck(!nullspace->remove, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Propagation of custom remove callbacks not supported when propagating (near) nullspaces with PCTelescope"); PetscCheck(!nullspace->rmctx, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Propagation of custom remove callback context not supported when propagating (near) nullspaces with PCTelescope"); } PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode PCTelescopeMatNullSpaceCreate_dmda(PC pc, PC_Telescope sred, Mat sub_mat) { Mat B; PetscFunctionBegin; PetscCall(PCGetOperators(pc, NULL, &B)); { MatNullSpace nullspace, sub_nullspace; PetscCall(MatGetNullSpace(B, &nullspace)); if (nullspace) { PetscCall(PetscInfo(pc, "PCTelescope: generating nullspace (DMDA)\n")); PetscCall(PCTelescopeSubNullSpaceCreate_dmda_Telescope(pc, sred, nullspace, &sub_nullspace)); if (PCTelescope_isActiveRank(sred)) { PetscCall(MatSetNullSpace(sub_mat, sub_nullspace)); PetscCall(MatNullSpaceDestroy(&sub_nullspace)); } } } { MatNullSpace nearnullspace, sub_nearnullspace; PetscCall(MatGetNearNullSpace(B, &nearnullspace)); if (nearnullspace) { PetscCall(PetscInfo(pc, "PCTelescope: generating near nullspace (DMDA)\n")); PetscCall(PCTelescopeSubNullSpaceCreate_dmda_Telescope(pc, sred, nearnullspace, &sub_nearnullspace)); if (PCTelescope_isActiveRank(sred)) { PetscCall(MatSetNearNullSpace(sub_mat, sub_nearnullspace)); PetscCall(MatNullSpaceDestroy(&sub_nearnullspace)); } } } PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode PCApply_Telescope_dmda(PC pc, Vec x, Vec y) { PC_Telescope sred = (PC_Telescope)pc->data; Mat perm; Vec xtmp, xp, xred, yred; PetscInt i, st, ed; VecScatter scatter; PetscScalar *array; const PetscScalar *x_array; PC_Telescope_DMDACtx *ctx; ctx = (PC_Telescope_DMDACtx *)sred->dm_ctx; xtmp = sred->xtmp; scatter = sred->scatter; xred = sred->xred; yred = sred->yred; perm = ctx->permutation; xp = ctx->xp; PetscFunctionBegin; PetscCall(PetscCitationsRegister(citation, &cited)); /* permute vector into ordering associated with re-partitioned dmda */ PetscCall(MatMultTranspose(perm, x, xp)); /* pull in vector x->xtmp */ PetscCall(VecScatterBegin(scatter, xp, xtmp, INSERT_VALUES, SCATTER_FORWARD)); PetscCall(VecScatterEnd(scatter, xp, xtmp, INSERT_VALUES, SCATTER_FORWARD)); /* copy vector entries into xred */ PetscCall(VecGetArrayRead(xtmp, &x_array)); if (xred) { PetscScalar *LA_xred; PetscCall(VecGetOwnershipRange(xred, &st, &ed)); PetscCall(VecGetArray(xred, &LA_xred)); for (i = 0; i < ed - st; i++) LA_xred[i] = x_array[i]; PetscCall(VecRestoreArray(xred, &LA_xred)); } PetscCall(VecRestoreArrayRead(xtmp, &x_array)); /* solve */ if (PCTelescope_isActiveRank(sred)) { PetscCall(KSPSolve(sred->ksp, xred, yred)); PetscCall(KSPCheckSolve(sred->ksp, pc, yred)); } /* return vector */ PetscCall(VecGetArray(xtmp, &array)); if (yred) { const PetscScalar *LA_yred; PetscCall(VecGetOwnershipRange(yred, &st, &ed)); PetscCall(VecGetArrayRead(yred, &LA_yred)); for (i = 0; i < ed - st; i++) array[i] = LA_yred[i]; PetscCall(VecRestoreArrayRead(yred, &LA_yred)); } PetscCall(VecRestoreArray(xtmp, &array)); PetscCall(VecScatterBegin(scatter, xtmp, xp, INSERT_VALUES, SCATTER_REVERSE)); PetscCall(VecScatterEnd(scatter, xtmp, xp, INSERT_VALUES, SCATTER_REVERSE)); PetscCall(MatMult(perm, xp, y)); PetscFunctionReturn(PETSC_SUCCESS); } 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) { PC_Telescope sred = (PC_Telescope)pc->data; Mat perm; Vec xtmp, xp, yred; PetscInt i, st, ed; VecScatter scatter; const PetscScalar *x_array; PetscBool default_init_guess_value = PETSC_FALSE; PC_Telescope_DMDACtx *ctx; PetscFunctionBegin; ctx = (PC_Telescope_DMDACtx *)sred->dm_ctx; xtmp = sred->xtmp; scatter = sred->scatter; yred = sred->yred; perm = ctx->permutation; xp = ctx->xp; PetscCheck(its <= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCApplyRichardson_Telescope_dmda only supports max_it = 1"); *reason = (PCRichardsonConvergedReason)0; if (!zeroguess) { PetscCall(PetscInfo(pc, "PCTelescopeDMDA: Scattering y for non-zero-initial guess\n")); /* permute vector into ordering associated with re-partitioned dmda */ PetscCall(MatMultTranspose(perm, y, xp)); /* pull in vector x->xtmp */ PetscCall(VecScatterBegin(scatter, xp, xtmp, INSERT_VALUES, SCATTER_FORWARD)); PetscCall(VecScatterEnd(scatter, xp, xtmp, INSERT_VALUES, SCATTER_FORWARD)); /* copy vector entries into xred */ PetscCall(VecGetArrayRead(xtmp, &x_array)); if (yred) { PetscScalar *LA_yred; PetscCall(VecGetOwnershipRange(yred, &st, &ed)); PetscCall(VecGetArray(yred, &LA_yred)); for (i = 0; i < ed - st; i++) LA_yred[i] = x_array[i]; PetscCall(VecRestoreArray(yred, &LA_yred)); } PetscCall(VecRestoreArrayRead(xtmp, &x_array)); } if (PCTelescope_isActiveRank(sred)) { PetscCall(KSPGetInitialGuessNonzero(sred->ksp, &default_init_guess_value)); if (!zeroguess) PetscCall(KSPSetInitialGuessNonzero(sred->ksp, PETSC_TRUE)); } PetscCall(PCApply_Telescope_dmda(pc, x, y)); if (PCTelescope_isActiveRank(sred)) PetscCall(KSPSetInitialGuessNonzero(sred->ksp, default_init_guess_value)); if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS; *outits = 1; PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode PCReset_Telescope_dmda(PC pc) { PC_Telescope sred = (PC_Telescope)pc->data; PC_Telescope_DMDACtx *ctx; PetscFunctionBegin; ctx = (PC_Telescope_DMDACtx *)sred->dm_ctx; PetscCall(VecDestroy(&ctx->xp)); PetscCall(MatDestroy(&ctx->permutation)); PetscCall(DMDestroy(&ctx->dmrepart)); PetscCall(PetscFree3(ctx->range_i_re, ctx->range_j_re, ctx->range_k_re)); PetscCall(PetscFree3(ctx->start_i_re, ctx->start_j_re, ctx->start_k_re)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMView_DA_Short_3d(DM dm, PetscViewer v) { PetscInt M, N, P, m, n, p, ndof, nsw; MPI_Comm comm; PetscMPIInt size; const char *prefix; PetscFunctionBegin; PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); PetscCallMPI(MPI_Comm_size(comm, &size)); PetscCall(DMGetOptionsPrefix(dm, &prefix)); PetscCall(DMDAGetInfo(dm, NULL, &M, &N, &P, &m, &n, &p, &ndof, &nsw, NULL, NULL, NULL, NULL)); if (prefix) PetscCall(PetscViewerASCIIPrintf(v, "DMDA Object: (%s) %d MPI processes\n", prefix, size)); else PetscCall(PetscViewerASCIIPrintf(v, "DMDA Object: %d MPI processes\n", size)); 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)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMView_DA_Short_2d(DM dm, PetscViewer v) { PetscInt M, N, m, n, ndof, nsw; MPI_Comm comm; PetscMPIInt size; const char *prefix; PetscFunctionBegin; PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); PetscCallMPI(MPI_Comm_size(comm, &size)); PetscCall(DMGetOptionsPrefix(dm, &prefix)); PetscCall(DMDAGetInfo(dm, NULL, &M, &N, NULL, &m, &n, NULL, &ndof, &nsw, NULL, NULL, NULL, NULL)); if (prefix) PetscCall(PetscViewerASCIIPrintf(v, "DMDA Object: (%s) %d MPI processes\n", prefix, size)); else PetscCall(PetscViewerASCIIPrintf(v, "DMDA Object: %d MPI processes\n", size)); 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)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMView_DA_Short(DM dm, PetscViewer v) { PetscInt dim; PetscFunctionBegin; PetscCall(DMDAGetInfo(dm, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); switch (dim) { case 2: PetscCall(DMView_DA_Short_2d(dm, v)); break; case 3: PetscCall(DMView_DA_Short_3d(dm, v)); break; } PetscFunctionReturn(PETSC_SUCCESS); }