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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55: PetscInitialize(&argc, &argv, (char *)0, help);
56: MPI_Comm_size(PETSC_COMM_WORLD, &size);
57: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
59: PetscOptionsBegin(PETSC_COMM_WORLD, "", "ex22 DMDA tutorial example options", "DMDA");
60: PetscOptionsRangeInt("-Mx", "dimension along x-axis", "ex22.c", Mx, &Mx, NULL, 0, PETSC_MAX_INT);
61: PetscOptionsRangeInt("-My", "dimension along y-axis", "ex22.c", My, &My, NULL, 0, PETSC_MAX_INT);
62: PetscOptionsRangeInt("-Mz", "dimension along z-axis", "ex22.c", Mz, &Mz, NULL, 0, PETSC_MAX_INT);
63: PetscOptionsEnum("-sliceaxis", "axis along which 2D slice is extracted from : X, Y, Z", "", sliceaxes, (PetscEnum)sliceaxis, (PetscEnum *)&sliceaxis, NULL);
64: 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) {
70: } else if (sliceaxis == DM_Y) {
72: } else if (sliceaxis == DM_Z) {
74: }
76: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: Create 3D DMDA object.
78: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
79: 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: DMSetFromOptions(da3D);
81: DMSetUp(da3D);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create the parent vector
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
86: DMCreateGlobalVector(da3D, &vec_full);
87: PetscObjectSetName((PetscObject)vec_full, "full_vector");
89: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: Populate the 3D vector
91: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: DMDAGetCorners(da3D, &ixs, &iys, &izs, &ixm, &iym, &izm);
93: 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: 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: DMDACreatePatchIS(da3D, &lower, &upper, &patchis_3d, patchis_offproc);
127: PetscPrintf(PETSC_COMM_WORLD, "\n IS to select slice from 3D DMDA vector : \n");
128: ISView(patchis_3d, PETSC_VIEWER_STDOUT_WORLD);
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Use the obtained IS to extract the slice as a subvector
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: VecGetSubVector(vec_full, patchis_3d, &vec_extracted);
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: View the extracted subvector
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DENSE);
139: PetscPrintf(PETSC_COMM_WORLD, "\n Extracted slice vector, in DMDA ordering : \n");
140: VecView(vec_extracted, PETSC_VIEWER_STDOUT_WORLD);
141: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Query 3D DMDA layout, get the subset MPI communicator
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
146: if (sliceaxis == DM_X) {
147: DMDAGetInfo(da3D, NULL, NULL, NULL, NULL, NULL, &m1, &m2, NULL, NULL, NULL, NULL, NULL, NULL);
148: DMDAGetOwnershipRanges(da3D, NULL, &l1, &l2);
149: M1 = My;
150: M2 = Mz;
151: DMDAGetProcessorSubset(da3D, DM_X, gp, &subset_mpi_comm);
152: } else if (sliceaxis == DM_Y) {
153: DMDAGetInfo(da3D, NULL, NULL, NULL, NULL, &m1, NULL, &m2, NULL, NULL, NULL, NULL, NULL, NULL);
154: DMDAGetOwnershipRanges(da3D, &l1, NULL, &l2);
155: M1 = Mx;
156: M2 = Mz;
157: DMDAGetProcessorSubset(da3D, DM_Y, gp, &subset_mpi_comm);
158: } else if (sliceaxis == DM_Z) {
159: DMDAGetInfo(da3D, NULL, NULL, NULL, NULL, &m1, &m2, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
160: DMDAGetOwnershipRanges(da3D, &l1, &l2, NULL);
161: M1 = Mx;
162: M2 = My;
163: 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: MPI_Comm_size(subset_mpi_comm, &size);
174: PetscSynchronizedPrintf(subset_mpi_comm, "subset MPI subcomm size is : %d, includes global rank : %d \n", size, rank);
175: PetscSynchronizedFlush(subset_mpi_comm, PETSC_STDOUT);
176: DMDACreate2d(subset_mpi_comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M1, M2, m1, m2, 1, 1, l1, l2, &da2D);
177: DMSetFromOptions(da2D);
178: 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: DMDACreatePatchIS(da2D, &lower, &upper, &patchis_2d, patchis_offproc);
187: /* Convert the 2D patch IS to natural indexing (column major flattened) */
188: ISDuplicate(patchis_2d, &scatis_natural_slice);
189: DMDAGetAO(da2D, &da2D_ao);
190: AOPetscToApplicationIS(da2D_ao, scatis_natural_slice);
191: ISGetIndices(scatis_natural_slice, &is_array);
192: ISGetLocalSize(scatis_natural_slice, &in);
194: /* Create an aliased IS on the 3D DMDA's communicator */
195: ISCreateGeneral(PETSC_COMM_WORLD, in, is_array, PETSC_USE_POINTER, &scatis_natural_slice_g);
196: ISRestoreIndices(scatis_natural_slice, &is_array);
198: /* Create a 2D DMDA global vector */
199: DMCreateGlobalVector(da2D, &vec_slice);
200: PetscObjectSetName((PetscObject)vec_slice, "slice_vector_natural");
201: VecGetLocalSize(vec_slice, &vn);
202: VecGetArrayRead(vec_slice, &array);
204: /* Create an aliased version of the above on the 3D DMDA's communicator */
205: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, vn, M1 * M2, array, &vec_slice_g);
206: 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: ISCreateGeneral(PETSC_COMM_WORLD, 0, NULL, PETSC_USE_POINTER, &scatis_natural_slice_g);
211: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, 0, M1 * M2, NULL, &vec_slice_g);
212: }
213: PetscBarrier(NULL);
215: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216: Create IS that maps from the extracted slice vector
217: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
218: VecGetOwnershipRange(vec_extracted, &low, &high);
219: ISCreateStride(PETSC_COMM_WORLD, high - low, low, 1, &scatis_extracted_slice);
221: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222: Scatter extracted subvector -> natural 2D slice vector
223: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
224: VecScatterCreate(vec_extracted, scatis_extracted_slice, vec_slice_g, scatis_natural_slice_g, &vscat);
225: VecScatterBegin(vscat, vec_extracted, vec_slice_g, INSERT_VALUES, SCATTER_FORWARD);
226: VecScatterEnd(vscat, vec_extracted, vec_slice_g, INSERT_VALUES, SCATTER_FORWARD);
228: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229: View the natural 2D slice vector
230: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
231: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DENSE);
232: PetscPrintf(PETSC_COMM_WORLD, "\n Extracted slice vector, in natural ordering : \n");
233: VecView(vec_slice_g, PETSC_VIEWER_STDOUT_WORLD);
234: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
236: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237: Restore subvector
238: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
239: VecRestoreSubVector(vec_full, patchis_3d, &vec_extracted);
241: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
242: Destroy data structures and exit.
243: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
244: VecDestroy(&vec_full);
245: VecScatterDestroy(&vscat);
246: ISDestroy(&scatis_extracted_slice);
247: ISDestroy(&scatis_natural_slice_g);
248: VecDestroy(&vec_slice_g);
249: ISDestroy(&patchis_3d);
250: DMDestroy(&da3D);
252: if (subset_mpi_comm != MPI_COMM_NULL) {
253: ISDestroy(&patchis_2d);
254: ISDestroy(&scatis_natural_slice);
255: VecDestroy(&vec_slice);
256: DMDestroy(&da2D);
257: MPI_Comm_free(&subset_mpi_comm);
258: }
260: 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*/