Actual source code: ex9.c
1: static char help[] = "This example shows 1) how to transfer vectors from a parent communicator to vectors on a child communicator and vice versa;\n\
2: 2) how to transfer vectors from a subcommunicator to vectors on another subcommunicator. The two subcommunicators are not\n\
3: required to cover all processes in PETSC_COMM_WORLD; 3) how to copy a vector from a parent communicator to vectors on its child communicators.\n\
4: To run any example with VECCUDA vectors, add -vectype cuda to the argument list\n\n";
6: #include <petscvec.h>
7: int main(int argc, char **argv)
8: {
9: PetscMPIInt nproc, grank, mycolor;
10: PetscInt i, n, N = 20, low, high;
11: MPI_Comm subcomm;
12: Vec x = NULL; /* global vectors on PETSC_COMM_WORLD */
13: Vec yg = NULL; /* global vectors on PETSC_COMM_WORLD */
14: VecScatter vscat;
15: IS ix, iy;
16: PetscBool iscuda = PETSC_FALSE; /* Option to use VECCUDA vectors */
17: PetscBool optionflag, compareflag;
18: char vectypename[PETSC_MAX_PATH_LEN];
19: PetscBool world2sub = PETSC_FALSE; /* Copy a vector from WORLD to a subcomm? */
20: PetscBool sub2sub = PETSC_FALSE; /* Copy a vector from a subcomm to another subcomm? */
21: PetscBool world2subs = PETSC_FALSE; /* Copy a vector from WORLD to multiple subcomms? */
23: PetscFunctionBeginUser;
24: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
25: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &nproc));
26: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &grank));
28: PetscCheck(nproc >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "This test must have at least two processes to run");
30: PetscCall(PetscOptionsGetBool(NULL, 0, "-world2sub", &world2sub, NULL));
31: PetscCall(PetscOptionsGetBool(NULL, 0, "-sub2sub", &sub2sub, NULL));
32: PetscCall(PetscOptionsGetBool(NULL, 0, "-world2subs", &world2subs, NULL));
33: PetscCall(PetscOptionsGetString(NULL, NULL, "-vectype", vectypename, sizeof(vectypename), &optionflag));
34: if (optionflag) {
35: PetscCall(PetscStrncmp(vectypename, "cuda", (size_t)4, &compareflag));
36: if (compareflag) iscuda = PETSC_TRUE;
37: }
39: /* Split PETSC_COMM_WORLD into three subcomms. Each process can only see the subcomm it belongs to */
40: mycolor = grank % 3;
41: PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, mycolor, grank, &subcomm));
43: /*===========================================================================
44: * Transfer a vector x defined on PETSC_COMM_WORLD to a vector y defined on
45: * a subcommunicator of PETSC_COMM_WORLD and vice versa.
46: *===========================================================================*/
47: if (world2sub) {
48: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
49: PetscCall(VecSetSizes(x, PETSC_DECIDE, N));
50: if (iscuda) PetscCall(VecSetType(x, VECCUDA));
51: else PetscCall(VecSetType(x, VECSTANDARD));
52: PetscCall(VecSetUp(x));
53: PetscCall(PetscObjectSetName((PetscObject)x, "x_commworld")); /* Give a name to view x clearly */
55: /* Initialize x to [-0.0, -1.0, -2.0, ..., -19.0] */
56: PetscCall(VecGetOwnershipRange(x, &low, &high));
57: for (i = low; i < high; i++) {
58: PetscScalar val = -i;
59: PetscCall(VecSetValue(x, i, val, INSERT_VALUES));
60: }
61: PetscCall(VecAssemblyBegin(x));
62: PetscCall(VecAssemblyEnd(x));
64: /* Transfer x to a vector y only defined on subcomm0 and vice versa */
65: if (mycolor == 0) { /* subcomm0 contains ranks 0, 3, 6, ... in PETSC_COMM_WORLD */
66: Vec y;
67: PetscScalar *yvalue;
68: PetscCall(VecCreate(subcomm, &y));
69: PetscCall(VecSetSizes(y, PETSC_DECIDE, N));
70: if (iscuda) PetscCall(VecSetType(y, VECCUDA));
71: else PetscCall(VecSetType(y, VECSTANDARD));
72: PetscCall(VecSetUp(y));
73: PetscCall(PetscObjectSetName((PetscObject)y, "y_subcomm_0")); /* Give a name to view y clearly */
74: PetscCall(VecGetLocalSize(y, &n));
75: if (iscuda) PetscCall(VecCUDAGetArray(y, &yvalue));
76: else PetscCall(VecGetArray(y, &yvalue));
77: /* Create yg on PETSC_COMM_WORLD and alias yg with y. They share the memory pointed by yvalue.
78: Note this is a collective call. All processes have to call it and supply consistent N.
79: */
80: PetscCall(VecCreateMPIWithArrayAndMemType(PETSC_COMM_WORLD, iscuda ? PETSC_MEMTYPE_CUDA : PETSC_MEMTYPE_HOST, 1, n, N, yvalue, &yg));
82: /* Create an identity map that makes yg[i] = x[i], i=0..N-1 */
83: PetscCall(VecGetOwnershipRange(yg, &low, &high)); /* low, high are global indices */
84: PetscCall(ISCreateStride(PETSC_COMM_SELF, high - low, low, 1, &ix));
85: PetscCall(ISDuplicate(ix, &iy));
87: /* Union of ix's on subcomm0 covers the full range of [0,N) */
88: PetscCall(VecScatterCreate(x, ix, yg, iy, &vscat));
89: PetscCall(VecScatterBegin(vscat, x, yg, INSERT_VALUES, SCATTER_FORWARD));
90: PetscCall(VecScatterEnd(vscat, x, yg, INSERT_VALUES, SCATTER_FORWARD));
92: /* Once yg got the data from x, we return yvalue to y so that we can use y in other operations.
93: VecGetArray must be paired with VecRestoreArray.
94: */
95: if (iscuda) PetscCall(VecCUDARestoreArray(y, &yvalue));
96: else PetscCall(VecRestoreArray(y, &yvalue));
98: /* Libraries on subcomm0 can safely use y now, for example, view and scale it */
99: PetscCall(VecView(y, PETSC_VIEWER_STDOUT_(subcomm)));
100: PetscCall(VecScale(y, 2.0));
102: /* Send the new y back to x */
103: PetscCall(VecGetArray(y, &yvalue)); /* If VecScale is done on GPU, PETSc will prepare a valid yvalue for access */
104: /* Supply new yvalue to yg without memory copying */
105: PetscCall(VecPlaceArray(yg, yvalue));
106: PetscCall(VecScatterBegin(vscat, yg, x, INSERT_VALUES, SCATTER_REVERSE));
107: PetscCall(VecScatterEnd(vscat, yg, x, INSERT_VALUES, SCATTER_REVERSE));
108: PetscCall(VecResetArray(yg));
109: PetscCall(VecRestoreArray(y, &yvalue));
110: PetscCall(VecDestroy(&y));
111: } else {
112: /* Ranks outside of subcomm0 do not supply values to yg */
113: PetscCall(VecCreateMPIWithArrayAndMemType(PETSC_COMM_WORLD, iscuda ? PETSC_MEMTYPE_CUDA : PETSC_MEMTYPE_HOST, 1, 0 /*n*/, N, NULL, &yg));
115: /* Ranks in subcomm0 already specified the full range of the identity map. The remaining
116: ranks just need to create empty ISes to cheat VecScatterCreate.
117: */
118: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &ix));
119: PetscCall(ISDuplicate(ix, &iy));
121: PetscCall(VecScatterCreate(x, ix, yg, iy, &vscat));
122: PetscCall(VecScatterBegin(vscat, x, yg, INSERT_VALUES, SCATTER_FORWARD));
123: PetscCall(VecScatterEnd(vscat, x, yg, INSERT_VALUES, SCATTER_FORWARD));
125: /* Send the new y back to x. Ranks outside of subcomm0 actually have nothing to send.
126: But they have to call VecScatterBegin/End since these routines are collective.
127: */
128: PetscCall(VecScatterBegin(vscat, yg, x, INSERT_VALUES, SCATTER_REVERSE));
129: PetscCall(VecScatterEnd(vscat, yg, x, INSERT_VALUES, SCATTER_REVERSE));
130: }
132: PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
133: PetscCall(ISDestroy(&ix));
134: PetscCall(ISDestroy(&iy));
135: PetscCall(VecDestroy(&x));
136: PetscCall(VecDestroy(&yg));
137: PetscCall(VecScatterDestroy(&vscat));
138: } /* world2sub */
140: /*===========================================================================
141: * Transfer a vector x defined on subcomm0 to a vector y defined on
142: * subcomm1. The two subcomms are not overlapping and their union is
143: * not necessarily equal to PETSC_COMM_WORLD.
144: *===========================================================================*/
145: if (sub2sub) {
146: if (mycolor == 0) {
147: /* Intentionally declare N as a local variable so that processes in subcomm1 do not know its value */
148: PetscInt n, N = 22;
149: Vec x, xg, yg;
150: IS ix, iy;
151: VecScatter vscat;
152: const PetscScalar *xvalue;
153: MPI_Comm intercomm, parentcomm;
154: PetscMPIInt lrank;
156: PetscCallMPI(MPI_Comm_rank(subcomm, &lrank));
157: /* x is on subcomm */
158: PetscCall(VecCreate(subcomm, &x));
159: PetscCall(VecSetSizes(x, PETSC_DECIDE, N));
160: if (iscuda) PetscCall(VecSetType(x, VECCUDA));
161: else PetscCall(VecSetType(x, VECSTANDARD));
162: PetscCall(VecSetUp(x));
163: PetscCall(VecGetOwnershipRange(x, &low, &high));
165: /* initialize x = [0.0, 1.0, 2.0, ..., 21.0] */
166: for (i = low; i < high; i++) {
167: PetscScalar val = i;
168: PetscCall(VecSetValue(x, i, val, INSERT_VALUES));
169: }
170: PetscCall(VecAssemblyBegin(x));
171: PetscCall(VecAssemblyEnd(x));
173: PetscCallMPI(MPI_Intercomm_create(subcomm, 0, PETSC_COMM_WORLD /*peer_comm*/, 1, 100 /*tag*/, &intercomm));
175: /* Tell rank 0 of subcomm1 the global size of x */
176: if (!lrank) PetscCallMPI(MPI_Send(&N, 1, MPIU_INT, 0 /*receiver's rank in remote comm, i.e., subcomm1*/, 200 /*tag*/, intercomm));
178: /* Create an intracomm PETSc can work on. Ranks in subcomm0 are ordered before ranks in subcomm1 in parentcomm.
179: But this order actually does not matter, since what we care is vector y, which is defined on subcomm1.
180: */
181: PetscCallMPI(MPI_Intercomm_merge(intercomm, 0 /*low*/, &parentcomm));
183: /* Create a vector xg on parentcomm, which shares memory with x */
184: PetscCall(VecGetLocalSize(x, &n));
185: if (iscuda) PetscCall(VecCUDAGetArrayRead(x, &xvalue));
186: else PetscCall(VecGetArrayRead(x, &xvalue));
187: PetscCall(VecCreateMPIWithArrayAndMemType(parentcomm, iscuda ? PETSC_MEMTYPE_CUDA : PETSC_MEMTYPE_HOST, 1, n, N, xvalue, &xg));
189: /* Ranks in subcomm 0 have nothing on yg, so they simply have n=0, array=NULL */
190: PetscCall(VecCreateMPIWithArrayAndMemType(parentcomm, iscuda ? PETSC_MEMTYPE_CUDA : PETSC_MEMTYPE_HOST, 1, 0 /*n*/, N, NULL /*array*/, &yg));
192: /* Create the vecscatter, which does identity map by setting yg[i] = xg[i], i=0..N-1. */
193: PetscCall(VecGetOwnershipRange(xg, &low, &high)); /* low, high are global indices of xg */
194: PetscCall(ISCreateStride(PETSC_COMM_SELF, high - low, low, 1, &ix));
195: PetscCall(ISDuplicate(ix, &iy));
196: PetscCall(VecScatterCreate(xg, ix, yg, iy, &vscat));
198: /* Scatter values from xg to yg */
199: PetscCall(VecScatterBegin(vscat, xg, yg, INSERT_VALUES, SCATTER_FORWARD));
200: PetscCall(VecScatterEnd(vscat, xg, yg, INSERT_VALUES, SCATTER_FORWARD));
202: /* After the VecScatter is done, xg is idle so we can safely return xvalue to x */
203: if (iscuda) PetscCall(VecCUDARestoreArrayRead(x, &xvalue));
204: else PetscCall(VecRestoreArrayRead(x, &xvalue));
205: PetscCall(VecDestroy(&x));
206: PetscCall(ISDestroy(&ix));
207: PetscCall(ISDestroy(&iy));
208: PetscCall(VecDestroy(&xg));
209: PetscCall(VecDestroy(&yg));
210: PetscCall(VecScatterDestroy(&vscat));
211: PetscCallMPI(MPI_Comm_free(&intercomm));
212: PetscCallMPI(MPI_Comm_free(&parentcomm));
213: } else if (mycolor == 1) { /* subcomm 1, containing ranks 1, 4, 7, ... in PETSC_COMM_WORLD */
214: PetscInt n, N;
215: Vec y, xg, yg;
216: IS ix, iy;
217: VecScatter vscat;
218: PetscScalar *yvalue;
219: MPI_Comm intercomm, parentcomm;
220: PetscMPIInt lrank;
222: PetscCallMPI(MPI_Comm_rank(subcomm, &lrank));
223: PetscCallMPI(MPI_Intercomm_create(subcomm, 0, PETSC_COMM_WORLD /*peer_comm*/, 0 /*remote_leader*/, 100 /*tag*/, &intercomm));
225: /* Two rank-0 are talking */
226: if (!lrank) PetscCallMPI(MPI_Recv(&N, 1, MPIU_INT, 0 /*sender's rank in remote comm, i.e. subcomm0*/, 200 /*tag*/, intercomm, MPI_STATUS_IGNORE));
227: /* Rank 0 of subcomm1 bcasts N to its members */
228: PetscCallMPI(MPI_Bcast(&N, 1, MPIU_INT, 0 /*local root*/, subcomm));
230: /* Create a intracomm PETSc can work on */
231: PetscCallMPI(MPI_Intercomm_merge(intercomm, 1 /*high*/, &parentcomm));
233: /* Ranks in subcomm1 have nothing on xg, so they simply have n=0, array=NULL.*/
234: PetscCall(VecCreateMPIWithArrayAndMemType(parentcomm, iscuda ? PETSC_MEMTYPE_CUDA : PETSC_MEMTYPE_HOST, 1 /*bs*/, 0 /*n*/, N, NULL /*array*/, &xg));
236: PetscCall(VecCreate(subcomm, &y));
237: PetscCall(VecSetSizes(y, PETSC_DECIDE, N));
238: if (iscuda) PetscCall(VecSetType(y, VECCUDA));
239: else PetscCall(VecSetType(y, VECSTANDARD));
240: PetscCall(VecSetUp(y));
242: PetscCall(PetscObjectSetName((PetscObject)y, "y_subcomm_1")); /* Give a name to view y clearly */
243: PetscCall(VecGetLocalSize(y, &n));
244: if (iscuda) PetscCall(VecCUDAGetArray(y, &yvalue));
245: else PetscCall(VecGetArray(y, &yvalue));
246: /* Create a vector yg on parentcomm, which shares memory with y. xg and yg must be
247: created in the same order in subcomm0/1. For example, we can not reverse the order of
248: creating xg and yg in subcomm1.
249: */
250: PetscCall(VecCreateMPIWithArrayAndMemType(parentcomm, iscuda ? PETSC_MEMTYPE_CUDA : PETSC_MEMTYPE_HOST, 1 /*bs*/, n, N, yvalue, &yg));
252: /* Ranks in subcomm0 already specified the full range of the identity map.
253: ranks in subcomm1 just need to create empty ISes to cheat VecScatterCreate.
254: */
255: PetscCall(ISCreateStride(PETSC_COMM_SELF, 0, 0, 1, &ix));
256: PetscCall(ISDuplicate(ix, &iy));
257: PetscCall(VecScatterCreate(xg, ix, yg, iy, &vscat));
259: /* Scatter values from xg to yg */
260: PetscCall(VecScatterBegin(vscat, xg, yg, INSERT_VALUES, SCATTER_FORWARD));
261: PetscCall(VecScatterEnd(vscat, xg, yg, INSERT_VALUES, SCATTER_FORWARD));
263: /* After the VecScatter is done, values in yg are available. y is our interest, so we return yvalue to y */
264: if (iscuda) PetscCall(VecCUDARestoreArray(y, &yvalue));
265: else PetscCall(VecRestoreArray(y, &yvalue));
267: /* Libraries on subcomm1 can safely use y now, for example, view it */
268: PetscCall(VecView(y, PETSC_VIEWER_STDOUT_(subcomm)));
270: PetscCall(VecDestroy(&y));
271: PetscCall(ISDestroy(&ix));
272: PetscCall(ISDestroy(&iy));
273: PetscCall(VecDestroy(&xg));
274: PetscCall(VecDestroy(&yg));
275: PetscCall(VecScatterDestroy(&vscat));
276: PetscCallMPI(MPI_Comm_free(&intercomm));
277: PetscCallMPI(MPI_Comm_free(&parentcomm));
278: } else if (mycolor == 2) { /* subcomm2 */
279: /* Processes in subcomm2 do not participate in the VecScatter. They can freely do unrelated things on subcomm2 */
280: }
281: } /* sub2sub */
283: /*===========================================================================
284: * Copy a vector x defined on PETSC_COMM_WORLD to vectors y defined on
285: * every subcommunicator of PETSC_COMM_WORLD. We could use multiple transfers
286: * as we did in case 1, but that is not efficient. Instead, we use one vecscatter
287: * to achieve that.
288: *===========================================================================*/
289: if (world2subs) {
290: Vec y;
291: PetscInt n, N = 15, xstart, ystart, low, high;
292: PetscScalar *yvalue;
294: /* Initialize x to [0, 1, 2, 3, ..., N-1] */
295: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
296: PetscCall(VecSetSizes(x, PETSC_DECIDE, N));
297: if (iscuda) PetscCall(VecSetType(x, VECCUDA));
298: else PetscCall(VecSetType(x, VECSTANDARD));
299: PetscCall(VecSetUp(x));
300: PetscCall(VecGetOwnershipRange(x, &low, &high));
301: for (i = low; i < high; i++) PetscCall(VecSetValue(x, i, (PetscScalar)i, INSERT_VALUES));
302: PetscCall(VecAssemblyBegin(x));
303: PetscCall(VecAssemblyEnd(x));
305: /* Every subcomm has a y as long as x */
306: PetscCall(VecCreate(subcomm, &y));
307: PetscCall(VecSetSizes(y, PETSC_DECIDE, N));
308: if (iscuda) PetscCall(VecSetType(y, VECCUDA));
309: else PetscCall(VecSetType(y, VECSTANDARD));
310: PetscCall(VecSetUp(y));
311: PetscCall(VecGetLocalSize(y, &n));
313: /* Create a global vector yg on PETSC_COMM_WORLD using y's memory. yg's global size = N*(number of subcommunicators).
314: Eeach rank in subcomms contributes a piece to construct the global yg. Keep in mind that pieces from a subcomm are not
315: necessarily consecutive in yg. That depends on how PETSC_COMM_WORLD is split. In our case, subcomm0 is made of rank
316: 0, 3, 6 etc from PETSC_COMM_WORLD. So subcomm0's pieces are interleaved with pieces from other subcomms in yg.
317: */
318: if (iscuda) PetscCall(VecCUDAGetArray(y, &yvalue));
319: else PetscCall(VecGetArray(y, &yvalue));
320: PetscCall(VecCreateMPIWithArrayAndMemType(PETSC_COMM_WORLD, iscuda ? PETSC_MEMTYPE_CUDA : PETSC_MEMTYPE_HOST, 1, n, PETSC_DECIDE, yvalue, &yg));
321: PetscCall(PetscObjectSetName((PetscObject)yg, "yg_on_subcomms")); /* Give a name to view yg clearly */
323: /* The following two lines are key. From xstart, we know where to pull entries from x. Note that we get xstart from y,
324: since first entry of y on this rank is from x[xstart]. From ystart, we know where ot put entries to yg.
325: */
326: PetscCall(VecGetOwnershipRange(y, &xstart, NULL));
327: PetscCall(VecGetOwnershipRange(yg, &ystart, NULL));
329: PetscCall(ISCreateStride(PETSC_COMM_SELF, n, xstart, 1, &ix));
330: PetscCall(ISCreateStride(PETSC_COMM_SELF, n, ystart, 1, &iy));
331: PetscCall(VecScatterCreate(x, ix, yg, iy, &vscat));
332: PetscCall(VecScatterBegin(vscat, x, yg, INSERT_VALUES, SCATTER_FORWARD));
333: PetscCall(VecScatterEnd(vscat, x, yg, INSERT_VALUES, SCATTER_FORWARD));
335: /* View yg on PETSC_COMM_WORLD before destroying it. We shall see the interleaving effect in output. */
336: PetscCall(VecView(yg, PETSC_VIEWER_STDOUT_WORLD));
337: PetscCall(VecDestroy(&yg));
339: /* Restory yvalue so that processes in subcomm can use y from now on. */
340: if (iscuda) PetscCall(VecCUDARestoreArray(y, &yvalue));
341: else PetscCall(VecRestoreArray(y, &yvalue));
342: PetscCall(VecScale(y, 3.0));
344: PetscCall(ISDestroy(&ix)); /* One can also destroy ix, iy immediately after VecScatterCreate() */
345: PetscCall(ISDestroy(&iy));
346: PetscCall(VecDestroy(&x));
347: PetscCall(VecDestroy(&y));
348: PetscCall(VecScatterDestroy(&vscat));
349: } /* world2subs */
351: PetscCallMPI(MPI_Comm_free(&subcomm));
352: PetscCall(PetscFinalize());
353: return 0;
354: }
356: /*TEST
358: build:
359: requires: !defined(PETSC_HAVE_MPIUNI)
361: testset:
362: nsize: 7
364: test:
365: suffix: 1
366: args: -world2sub
368: test:
369: suffix: 2
370: args: -sub2sub
371: # deadlocks with NECMPI and INTELMPI (20210400300)
372: requires: !defined(PETSC_HAVE_NECMPI) !defined(PETSC_HAVE_I_MPI)
374: test:
375: suffix: 3
376: args: -world2subs
378: test:
379: suffix: 4
380: args: -world2sub -vectype cuda
381: requires: cuda
383: test:
384: suffix: 5
385: args: -sub2sub -vectype cuda
386: requires: cuda
388: test:
389: suffix: 6
390: args: -world2subs -vectype cuda
391: requires: cuda
393: testset:
394: nsize: 7
395: args: -world2sub -sf_type neighbor
396: output_file: output/ex9_1.out
397: # segfaults with NECMPI
398: requires: defined(PETSC_HAVE_MPI_NEIGHBORHOOD_COLLECTIVES) !defined(PETSC_HAVE_NECMPI)
400: test:
401: suffix: 71
403: test:
404: suffix: 72
405: requires: defined(PETSC_HAVE_MPI_PERSISTENT_NEIGHBORHOOD_COLLECTIVES)
406: args: -sf_neighbor_persistent
408: TEST*/