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 PetscErrorCode ierr; /* error checking */ 51 52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 53 Initialize program and set problem parameters 54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 55 ierr = PetscInitialize(&argc, &argv, (char*)0, help);if (ierr) return ierr; 56 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 57 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 58 59 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "ex22 DMDA tutorial example options", "DMDA");CHKERRQ(ierr); 60 CHKERRQ(PetscOptionsRangeInt("-Mx", "dimension along x-axis", "ex22.c", Mx, &Mx, NULL, 0, PETSC_MAX_INT)); 61 CHKERRQ(PetscOptionsRangeInt("-My", "dimension along y-axis", "ex22.c", My, &My, NULL, 0, PETSC_MAX_INT)); 62 CHKERRQ(PetscOptionsRangeInt("-Mz", "dimension along z-axis", "ex22.c", Mz, &Mz, NULL, 0, PETSC_MAX_INT)); 63 CHKERRQ(PetscOptionsEnum("-sliceaxis","axis along which 2D slice is extracted from : X, Y, Z","",sliceaxes,(PetscEnum)sliceaxis,(PetscEnum*)&sliceaxis,NULL)); 64 CHKERRQ(PetscOptionsRangeInt("-gp", "index along sliceaxis at which 2D slice is extracted", "ex22.c", gp, &gp, NULL, 0, PETSC_MAX_INT)); 65 ierr = PetscOptionsEnd();CHKERRQ(ierr); 66 67 /* Ensure that the requested slice is not out of bounds for the selected axis */ 68 if (sliceaxis==DM_X) { 69 PetscCheckFalse(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 PetscCheckFalse(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 PetscCheckFalse(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 ierr = DMDACreate3d(PETSC_COMM_WORLD, 80 DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 81 DMDA_STENCIL_STAR, 82 Mx, My, Mz, 83 PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 84 1, 1, 85 NULL, NULL, NULL, 86 &da3D);CHKERRQ(ierr); 87 CHKERRQ(DMSetFromOptions(da3D)); 88 CHKERRQ(DMSetUp(da3D)); 89 90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 91 Create the parent vector 92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 93 CHKERRQ(DMCreateGlobalVector(da3D, &vec_full)); 94 CHKERRQ(PetscObjectSetName((PetscObject) vec_full, "full_vector")); 95 96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97 Populate the 3D vector 98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 99 CHKERRQ(DMDAGetCorners(da3D, &ixs, &iys, &izs, &ixm, &iym, &izm)); 100 CHKERRQ(DMDAVecGetArray(da3D, vec_full, &vecdata3d)); 101 for (k=izs; k<izs+izm; k++) { 102 for (j=iys; j<iys+iym; j++) { 103 for (i=ixs; i<ixs+ixm; i++) { 104 vecdata3d[k][j][i] = ((i-Mx/2.0)*(j+Mx/2.0))+k*100; 105 } 106 } 107 } 108 CHKERRQ(DMDAVecRestoreArray(da3D, vec_full, &vecdata3d)); 109 110 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 111 Get an IS corresponding to a 2D slice 112 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 113 if (sliceaxis==DM_X) { 114 lower.i = gp; lower.j = 0; lower.k = 0; 115 upper.i = gp; upper.j = My; upper.k = Mz; 116 } else if (sliceaxis==DM_Y) { 117 lower.i = 0; lower.j = gp; lower.k = 0; 118 upper.i = Mx; upper.j = gp; upper.k = Mz; 119 } else if (sliceaxis==DM_Z) { 120 lower.i = 0; lower.j = 0; lower.k = gp; 121 upper.i = Mx; upper.j = My; upper.k = gp; 122 } 123 CHKERRQ(DMDACreatePatchIS(da3D, &lower, &upper, &patchis_3d, patchis_offproc)); 124 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "\n IS to select slice from 3D DMDA vector : \n")); 125 CHKERRQ(ISView(patchis_3d, PETSC_VIEWER_STDOUT_WORLD)); 126 127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 128 Use the obtained IS to extract the slice as a subvector 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 CHKERRQ(VecGetSubVector(vec_full, patchis_3d, &vec_extracted)); 131 132 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 133 View the extracted subvector 134 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 135 CHKERRQ(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_DENSE)); 136 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "\n Extracted slice vector, in DMDA ordering : \n")); 137 CHKERRQ(VecView(vec_extracted, PETSC_VIEWER_STDOUT_WORLD)); 138 CHKERRQ(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD)); 139 140 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 141 Query 3D DMDA layout, get the subset MPI communicator 142 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 143 if (sliceaxis==DM_X) { 144 CHKERRQ(DMDAGetInfo(da3D, NULL, NULL, NULL, NULL, NULL, &m1, &m2, NULL, NULL, NULL, NULL, NULL, NULL)); 145 CHKERRQ(DMDAGetOwnershipRanges(da3D, NULL, &l1, &l2)); 146 M1 = My; M2 = Mz; 147 CHKERRQ(DMDAGetProcessorSubset(da3D, DM_X, gp, &subset_mpi_comm)); 148 } else if (sliceaxis==DM_Y) { 149 CHKERRQ(DMDAGetInfo(da3D, NULL, NULL, NULL, NULL, &m1, NULL, &m2, NULL, NULL, NULL, NULL, NULL, NULL)); 150 CHKERRQ(DMDAGetOwnershipRanges(da3D, &l1, NULL, &l2)); 151 M1 = Mx; M2 = Mz; 152 CHKERRQ(DMDAGetProcessorSubset(da3D, DM_Y, gp, &subset_mpi_comm)); 153 } else if (sliceaxis==DM_Z) { 154 CHKERRQ(DMDAGetInfo(da3D, NULL, NULL, NULL, NULL, &m1, &m2, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 155 CHKERRQ(DMDAGetOwnershipRanges(da3D, &l1, &l2, NULL)); 156 M1 = Mx; M2 = My; 157 CHKERRQ(DMDAGetProcessorSubset(da3D, DM_Z, gp, &subset_mpi_comm)); 158 } 159 160 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 161 Create 2D DMDA object, 162 vector (that will hold the slice as a column major flattened array) & 163 index set (that will be used for scattering to the column major 164 indexed slice vector) 165 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 166 if (subset_mpi_comm != MPI_COMM_NULL) { 167 CHKERRMPI(MPI_Comm_size(subset_mpi_comm, &size)); 168 CHKERRQ(PetscSynchronizedPrintf(subset_mpi_comm, "subset MPI subcomm size is : %d, includes global rank : %d \n", size, rank)); 169 CHKERRQ(PetscSynchronizedFlush(subset_mpi_comm, PETSC_STDOUT)); 170 ierr = DMDACreate2d(subset_mpi_comm, 171 DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, 172 DMDA_STENCIL_STAR, 173 M1, M2, 174 m1, m2, 175 1, 1, 176 l1, l2, 177 &da2D);CHKERRQ(ierr); 178 CHKERRQ(DMSetFromOptions(da2D)); 179 CHKERRQ(DMSetUp(da2D)); 180 181 /* Create a 2D patch IS for the slice */ 182 lower.i = 0; lower.j = 0; 183 upper.i = M1; upper.j = M2; 184 CHKERRQ(DMDACreatePatchIS(da2D, &lower, &upper, &patchis_2d, patchis_offproc)); 185 186 /* Convert the 2D patch IS to natural indexing (column major flattened) */ 187 CHKERRQ(ISDuplicate(patchis_2d, &scatis_natural_slice)); 188 CHKERRQ(DMDAGetAO(da2D, &da2D_ao)); 189 CHKERRQ(AOPetscToApplicationIS(da2D_ao, scatis_natural_slice)); 190 CHKERRQ(ISGetIndices(scatis_natural_slice, &is_array)); 191 CHKERRQ(ISGetLocalSize(scatis_natural_slice, &in)); 192 193 /* Create an aliased IS on the 3D DMDA's communicator */ 194 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD, in, is_array, PETSC_USE_POINTER, &scatis_natural_slice_g)); 195 CHKERRQ(ISRestoreIndices(scatis_natural_slice, &is_array)); 196 197 /* Create a 2D DMDA global vector */ 198 CHKERRQ(DMCreateGlobalVector(da2D, &vec_slice)); 199 CHKERRQ(PetscObjectSetName((PetscObject) vec_slice, "slice_vector_natural")); 200 CHKERRQ(VecGetLocalSize(vec_slice ,&vn)); 201 CHKERRQ(VecGetArrayRead(vec_slice, &array)); 202 203 /* Create an aliased version of the above on the 3D DMDA's communicator */ 204 CHKERRQ(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, vn, M1*M2, array, &vec_slice_g)); 205 CHKERRQ(VecRestoreArrayRead(vec_slice, &array)); 206 } else { 207 /* Ranks not part of the subset MPI communicator provide no entries, but the routines for creating 208 the IS and Vec on the 3D DMDA's communicator still need to called, since they are collective routines */ 209 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD, 0, NULL, PETSC_USE_POINTER, &scatis_natural_slice_g)); 210 CHKERRQ(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, 0, M1*M2, NULL, &vec_slice_g)); 211 } 212 CHKERRQ(PetscBarrier(NULL)); 213 214 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 215 Create IS that maps from the extracted slice vector 216 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 217 CHKERRQ(VecGetOwnershipRange(vec_extracted, &low, &high)); 218 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD, high-low, low, 1, &scatis_extracted_slice)); 219 220 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 221 Scatter extracted subvector -> natural 2D slice vector 222 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 223 CHKERRQ(VecScatterCreate(vec_extracted, scatis_extracted_slice, vec_slice_g, scatis_natural_slice_g, &vscat)); 224 CHKERRQ(VecScatterBegin(vscat, vec_extracted, vec_slice_g, INSERT_VALUES, SCATTER_FORWARD)); 225 CHKERRQ(VecScatterEnd(vscat, vec_extracted, vec_slice_g, INSERT_VALUES, SCATTER_FORWARD)); 226 227 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 228 View the natural 2D slice vector 229 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 230 CHKERRQ(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_DENSE)); 231 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "\n Extracted slice vector, in natural ordering : \n")); 232 CHKERRQ(VecView(vec_slice_g, PETSC_VIEWER_STDOUT_WORLD)); 233 CHKERRQ(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD)); 234 235 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 236 Restore subvector 237 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 238 CHKERRQ(VecRestoreSubVector(vec_full, patchis_3d, &vec_extracted)); 239 240 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 241 Destroy data structures and exit. 242 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 243 CHKERRQ(VecDestroy(&vec_full)); 244 CHKERRQ(VecScatterDestroy(&vscat)); 245 CHKERRQ(ISDestroy(&scatis_extracted_slice)); 246 CHKERRQ(ISDestroy(&scatis_natural_slice_g)); 247 CHKERRQ(VecDestroy(&vec_slice_g)); 248 CHKERRQ(ISDestroy(&patchis_3d)); 249 CHKERRQ(DMDestroy(&da3D)); 250 251 if (subset_mpi_comm != MPI_COMM_NULL) { 252 CHKERRQ(ISDestroy(&patchis_2d)); 253 CHKERRQ(ISDestroy(&scatis_natural_slice)); 254 CHKERRQ(VecDestroy(&vec_slice)); 255 CHKERRQ(DMDestroy(&da2D)); 256 CHKERRMPI(MPI_Comm_free(&subset_mpi_comm)); 257 } 258 259 ierr = PetscFinalize(); 260 return ierr; 261 } 262 263 /*TEST 264 265 test: 266 nsize: 1 267 args: -sliceaxis X -gp 0 268 269 test: 270 suffix: 2 271 nsize: 2 272 args: -sliceaxis Y -gp 1 273 filter: grep -v "subset MPI subcomm" 274 275 test: 276 suffix: 3 277 nsize: 3 278 args: -sliceaxis Z -gp 2 279 filter: grep -v "subset MPI subcomm" 280 281 test: 282 suffix: 4 283 nsize: 4 284 args: -sliceaxis X -gp 2 285 filter: grep -v "subset MPI subcomm" 286 287 test: 288 suffix: 5 289 nsize: 4 290 args: -sliceaxis Z -gp 1 291 filter: grep -v "subset MPI subcomm" 292 293 TEST*/ 294