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