1 static char help[] ="Extract a 2D slice in natural ordering from a 3D vector, Command line options : \n\ 2 Mx/My/Mz - set the dimensions of the parent vector \n\ 3 sliceaxis - string describing the axis along which the sice will be selected : X, Y, Z \n\ 4 gp - global grid point number along the sliceaxis direction where the slice will be extracted from the parent vector \n"; 5 6 /* 7 This example shows to extract a 2D slice in natural ordering 8 from a 3D DMDA vector (first by extracting the slice and then 9 by converting it to natural ordering) 10 */ 11 12 #include <petscdmda.h> 13 14 const char *const sliceaxes[] = {"X","Y","Z","sliceaxis","DM_",NULL}; 15 16 int main(int argc,char **argv) 17 { 18 DM da3D; /* 3D DMDA object */ 19 DM da2D; /* 2D DMDA object */ 20 Vec vec_full; /* Parent vector */ 21 Vec vec_extracted; /* Extracted slice vector (in DMDA ordering) */ 22 Vec vec_slice; /* vector in natural ordering */ 23 Vec vec_slice_g; /* aliased vector in natural ordering */ 24 IS patchis_3d; /* IS to select slice and extract subvector */ 25 IS patchis_2d; /* Patch IS for 2D vector, will be converted to application ordering */ 26 IS scatis_extracted_slice; /* PETSc indexed IS for extracted slice */ 27 IS scatis_natural_slice; /* natural/application ordered IS for slice*/ 28 IS scatis_natural_slice_g; /* aliased natural/application ordered IS for slice */ 29 VecScatter vscat; /* scatter slice in DMDA ordering <-> slice in column major ordering */ 30 AO da2D_ao; /* AO associated with 2D DMDA */ 31 MPI_Comm subset_mpi_comm=MPI_COMM_NULL; /* MPI communicator where the slice lives */ 32 PetscScalar ***vecdata3d; /* Pointer to access 3d parent vector */ 33 const PetscScalar *array; /* pointer to create aliased Vec */ 34 PetscInt Mx=4,My=4,Mz=4; /* Dimensions for 3D DMDA */ 35 const PetscInt *l1,*l2; /* 3D DMDA layout */ 36 PetscInt M1=-1,M2=-1; /* Dimensions for 2D DMDA */ 37 PetscInt m1=-1,m2=-1; /* Layouts for 2D DMDA */ 38 PetscInt gp=2; /* grid point along sliceaxis to pick the slice */ 39 DMDirection sliceaxis=DM_X; /* Select axis along which the slice will be extracted */ 40 PetscInt i,j,k; /* Iteration indices */ 41 PetscInt ixs,iys,izs; /* Corner indices for 3D vector */ 42 PetscInt ixm,iym,izm; /* Widths of parent vector */ 43 PetscInt low, high; /* ownership range indices */ 44 PetscInt in; /* local size index for IS*/ 45 PetscInt vn; /* local size index */ 46 const PetscInt *is_array; /* pointer to create aliased IS */ 47 MatStencil lower, upper; /* Stencils to select slice for Vec */ 48 PetscBool patchis_offproc = PETSC_FALSE; /* flag to DMDACreatePatchIS indicating that off-proc values are to be ignored */ 49 PetscMPIInt rank,size; /* MPI rank and size */ 50 51 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 52 Initialize program and set problem parameters 53 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 54 PetscFunctionBeginUser; 55 PetscCall(PetscInitialize(&argc, &argv, (char*)0, help)); 56 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 57 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 58 59 PetscOptionsBegin(PETSC_COMM_WORLD, "", "ex22 DMDA tutorial example options", "DMDA"); 60 PetscCall(PetscOptionsRangeInt("-Mx", "dimension along x-axis", "ex22.c", Mx, &Mx, NULL, 0, PETSC_MAX_INT)); 61 PetscCall(PetscOptionsRangeInt("-My", "dimension along y-axis", "ex22.c", My, &My, NULL, 0, PETSC_MAX_INT)); 62 PetscCall(PetscOptionsRangeInt("-Mz", "dimension along z-axis", "ex22.c", Mz, &Mz, NULL, 0, PETSC_MAX_INT)); 63 PetscCall(PetscOptionsEnum("-sliceaxis","axis along which 2D slice is extracted from : X, Y, Z","",sliceaxes,(PetscEnum)sliceaxis,(PetscEnum*)&sliceaxis,NULL)); 64 PetscCall(PetscOptionsRangeInt("-gp", "index along sliceaxis at which 2D slice is extracted", "ex22.c", gp, &gp, NULL, 0, PETSC_MAX_INT)); 65 PetscOptionsEnd(); 66 67 /* Ensure that the requested slice is not out of bounds for the selected axis */ 68 if (sliceaxis==DM_X) { 69 PetscCheck(gp<=Mx,PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "grid point along sliceaxis is larger than largest index!"); 70 } else if (sliceaxis==DM_Y) { 71 PetscCheck(gp<=My,PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "grid point along sliceaxis is larger than largest index!"); 72 } else if (sliceaxis==DM_Z) { 73 PetscCheck(gp<=Mz,PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "grid point along sliceaxis is larger than largest index!"); 74 } 75 76 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 77 Create 3D DMDA object. 78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 79 PetscCall(DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,Mx, My, Mz, 80 PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,1, 1,NULL, NULL, NULL,&da3D)); 81 PetscCall(DMSetFromOptions(da3D)); 82 PetscCall(DMSetUp(da3D)); 83 84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 85 Create the parent vector 86 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 87 PetscCall(DMCreateGlobalVector(da3D, &vec_full)); 88 PetscCall(PetscObjectSetName((PetscObject) vec_full, "full_vector")); 89 90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 91 Populate the 3D vector 92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 93 PetscCall(DMDAGetCorners(da3D, &ixs, &iys, &izs, &ixm, &iym, &izm)); 94 PetscCall(DMDAVecGetArray(da3D, vec_full, &vecdata3d)); 95 for (k=izs; k<izs+izm; k++) { 96 for (j=iys; j<iys+iym; j++) { 97 for (i=ixs; i<ixs+ixm; i++) { 98 vecdata3d[k][j][i] = ((i-Mx/2.0)*(j+Mx/2.0))+k*100; 99 } 100 } 101 } 102 PetscCall(DMDAVecRestoreArray(da3D, vec_full, &vecdata3d)); 103 104 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 105 Get an IS corresponding to a 2D slice 106 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 107 if (sliceaxis==DM_X) { 108 lower.i = gp; lower.j = 0; lower.k = 0; 109 upper.i = gp; upper.j = My; upper.k = Mz; 110 } else if (sliceaxis==DM_Y) { 111 lower.i = 0; lower.j = gp; lower.k = 0; 112 upper.i = Mx; upper.j = gp; upper.k = Mz; 113 } else if (sliceaxis==DM_Z) { 114 lower.i = 0; lower.j = 0; lower.k = gp; 115 upper.i = Mx; upper.j = My; upper.k = gp; 116 } 117 PetscCall(DMDACreatePatchIS(da3D, &lower, &upper, &patchis_3d, patchis_offproc)); 118 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n IS to select slice from 3D DMDA vector : \n")); 119 PetscCall(ISView(patchis_3d, PETSC_VIEWER_STDOUT_WORLD)); 120 121 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 122 Use the obtained IS to extract the slice as a subvector 123 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 124 PetscCall(VecGetSubVector(vec_full, patchis_3d, &vec_extracted)); 125 126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127 View the extracted subvector 128 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 129 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_DENSE)); 130 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n Extracted slice vector, in DMDA ordering : \n")); 131 PetscCall(VecView(vec_extracted, PETSC_VIEWER_STDOUT_WORLD)); 132 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD)); 133 134 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 135 Query 3D DMDA layout, get the subset MPI communicator 136 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 137 if (sliceaxis==DM_X) { 138 PetscCall(DMDAGetInfo(da3D, NULL, NULL, NULL, NULL, NULL, &m1, &m2, NULL, NULL, NULL, NULL, NULL, NULL)); 139 PetscCall(DMDAGetOwnershipRanges(da3D, NULL, &l1, &l2)); 140 M1 = My; M2 = Mz; 141 PetscCall(DMDAGetProcessorSubset(da3D, DM_X, gp, &subset_mpi_comm)); 142 } else if (sliceaxis==DM_Y) { 143 PetscCall(DMDAGetInfo(da3D, NULL, NULL, NULL, NULL, &m1, NULL, &m2, NULL, NULL, NULL, NULL, NULL, NULL)); 144 PetscCall(DMDAGetOwnershipRanges(da3D, &l1, NULL, &l2)); 145 M1 = Mx; M2 = Mz; 146 PetscCall(DMDAGetProcessorSubset(da3D, DM_Y, gp, &subset_mpi_comm)); 147 } else if (sliceaxis==DM_Z) { 148 PetscCall(DMDAGetInfo(da3D, NULL, NULL, NULL, NULL, &m1, &m2, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 149 PetscCall(DMDAGetOwnershipRanges(da3D, &l1, &l2, NULL)); 150 M1 = Mx; M2 = My; 151 PetscCall(DMDAGetProcessorSubset(da3D, DM_Z, gp, &subset_mpi_comm)); 152 } 153 154 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 155 Create 2D DMDA object, 156 vector (that will hold the slice as a column major flattened array) & 157 index set (that will be used for scattering to the column major 158 indexed slice vector) 159 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 160 if (subset_mpi_comm != MPI_COMM_NULL) { 161 PetscCallMPI(MPI_Comm_size(subset_mpi_comm, &size)); 162 PetscCall(PetscSynchronizedPrintf(subset_mpi_comm, "subset MPI subcomm size is : %d, includes global rank : %d \n", size, rank)); 163 PetscCall(PetscSynchronizedFlush(subset_mpi_comm, PETSC_STDOUT)); 164 PetscCall(DMDACreate2d(subset_mpi_comm,DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M1, M2,m1, m2,1, 1,l1, l2,&da2D)); 165 PetscCall(DMSetFromOptions(da2D)); 166 PetscCall(DMSetUp(da2D)); 167 168 /* Create a 2D patch IS for the slice */ 169 lower.i = 0; lower.j = 0; 170 upper.i = M1; upper.j = M2; 171 PetscCall(DMDACreatePatchIS(da2D, &lower, &upper, &patchis_2d, patchis_offproc)); 172 173 /* Convert the 2D patch IS to natural indexing (column major flattened) */ 174 PetscCall(ISDuplicate(patchis_2d, &scatis_natural_slice)); 175 PetscCall(DMDAGetAO(da2D, &da2D_ao)); 176 PetscCall(AOPetscToApplicationIS(da2D_ao, scatis_natural_slice)); 177 PetscCall(ISGetIndices(scatis_natural_slice, &is_array)); 178 PetscCall(ISGetLocalSize(scatis_natural_slice, &in)); 179 180 /* Create an aliased IS on the 3D DMDA's communicator */ 181 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, in, is_array, PETSC_USE_POINTER, &scatis_natural_slice_g)); 182 PetscCall(ISRestoreIndices(scatis_natural_slice, &is_array)); 183 184 /* Create a 2D DMDA global vector */ 185 PetscCall(DMCreateGlobalVector(da2D, &vec_slice)); 186 PetscCall(PetscObjectSetName((PetscObject) vec_slice, "slice_vector_natural")); 187 PetscCall(VecGetLocalSize(vec_slice ,&vn)); 188 PetscCall(VecGetArrayRead(vec_slice, &array)); 189 190 /* Create an aliased version of the above on the 3D DMDA's communicator */ 191 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, vn, M1*M2, array, &vec_slice_g)); 192 PetscCall(VecRestoreArrayRead(vec_slice, &array)); 193 } else { 194 /* Ranks not part of the subset MPI communicator provide no entries, but the routines for creating 195 the IS and Vec on the 3D DMDA's communicator still need to called, since they are collective routines */ 196 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 0, NULL, PETSC_USE_POINTER, &scatis_natural_slice_g)); 197 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, 0, M1*M2, NULL, &vec_slice_g)); 198 } 199 PetscCall(PetscBarrier(NULL)); 200 201 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 202 Create IS that maps from the extracted slice vector 203 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 204 PetscCall(VecGetOwnershipRange(vec_extracted, &low, &high)); 205 PetscCall(ISCreateStride(PETSC_COMM_WORLD, high-low, low, 1, &scatis_extracted_slice)); 206 207 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 208 Scatter extracted subvector -> natural 2D slice vector 209 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 210 PetscCall(VecScatterCreate(vec_extracted, scatis_extracted_slice, vec_slice_g, scatis_natural_slice_g, &vscat)); 211 PetscCall(VecScatterBegin(vscat, vec_extracted, vec_slice_g, INSERT_VALUES, SCATTER_FORWARD)); 212 PetscCall(VecScatterEnd(vscat, vec_extracted, vec_slice_g, INSERT_VALUES, SCATTER_FORWARD)); 213 214 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 215 View the natural 2D slice vector 216 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 217 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_DENSE)); 218 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n Extracted slice vector, in natural ordering : \n")); 219 PetscCall(VecView(vec_slice_g, PETSC_VIEWER_STDOUT_WORLD)); 220 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD)); 221 222 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 223 Restore subvector 224 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 225 PetscCall(VecRestoreSubVector(vec_full, patchis_3d, &vec_extracted)); 226 227 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 228 Destroy data structures and exit. 229 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 230 PetscCall(VecDestroy(&vec_full)); 231 PetscCall(VecScatterDestroy(&vscat)); 232 PetscCall(ISDestroy(&scatis_extracted_slice)); 233 PetscCall(ISDestroy(&scatis_natural_slice_g)); 234 PetscCall(VecDestroy(&vec_slice_g)); 235 PetscCall(ISDestroy(&patchis_3d)); 236 PetscCall(DMDestroy(&da3D)); 237 238 if (subset_mpi_comm != MPI_COMM_NULL) { 239 PetscCall(ISDestroy(&patchis_2d)); 240 PetscCall(ISDestroy(&scatis_natural_slice)); 241 PetscCall(VecDestroy(&vec_slice)); 242 PetscCall(DMDestroy(&da2D)); 243 PetscCallMPI(MPI_Comm_free(&subset_mpi_comm)); 244 } 245 246 PetscCall(PetscFinalize()); 247 return 0; 248 } 249 250 /*TEST 251 252 test: 253 nsize: 1 254 args: -sliceaxis X -gp 0 255 256 test: 257 suffix: 2 258 nsize: 2 259 args: -sliceaxis Y -gp 1 260 filter: grep -v "subset MPI subcomm" 261 262 test: 263 suffix: 3 264 nsize: 3 265 args: -sliceaxis Z -gp 2 266 filter: grep -v "subset MPI subcomm" 267 268 test: 269 suffix: 4 270 nsize: 4 271 args: -sliceaxis X -gp 2 272 filter: grep -v "subset MPI subcomm" 273 274 test: 275 suffix: 5 276 nsize: 4 277 args: -sliceaxis Z -gp 1 278 filter: grep -v "subset MPI subcomm" 279 280 TEST*/ 281