xref: /petsc/src/vec/is/sf/tests/ex9.c (revision 40badf4fbc550ac1f60bd080eaff6de6d55b946d)
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";
5 
6 #include <petscvec.h>
7 int main(int argc,char **argv)
8 {
9   PetscErrorCode ierr;
10   PetscMPIInt    nproc,grank,mycolor;
11   PetscInt       i,n,N = 20,low,high;
12   MPI_Comm       subcomm;
13   Vec            x  = PETSC_NULL; /* global vectors on PETSC_COMM_WORLD */
14   Vec            yg = PETSC_NULL; /* global vectors on PETSC_COMM_WORLD */
15   VecScatter     vscat;
16   IS             ix,iy;
17   PetscBool      iscuda = PETSC_FALSE;      /* Option to use VECCUDA vectors */
18   PetscBool      optionflag, compareflag;
19   char           vectypename[PETSC_MAX_PATH_LEN];
20   PetscBool      world2sub  = PETSC_FALSE;  /* Copy a vector from WORLD to a subcomm? */
21   PetscBool      sub2sub    = PETSC_FALSE;  /* Copy a vector from a subcomm to another subcomm? */
22   PetscBool      world2subs = PETSC_FALSE;  /* Copy a vector from WORLD to multiple subcomms? */
23 
24   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
25   CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&nproc));
26   CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&grank));
27 
28   PetscCheckFalse(nproc < 2,PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"This test must have at least two processes to run");
29 
30   CHKERRQ(PetscOptionsGetBool(NULL,0,"-world2sub",&world2sub,NULL));
31   CHKERRQ(PetscOptionsGetBool(NULL,0,"-sub2sub",&sub2sub,NULL));
32   CHKERRQ(PetscOptionsGetBool(NULL,0,"-world2subs",&world2subs,NULL));
33   CHKERRQ(PetscOptionsGetString(NULL,NULL,"-vectype",vectypename,sizeof(vectypename),&optionflag));
34   if (optionflag) {
35     CHKERRQ(PetscStrncmp(vectypename, "cuda", (size_t)4, &compareflag));
36     if (compareflag) iscuda = PETSC_TRUE;
37   }
38 
39   /* Split PETSC_COMM_WORLD into three subcomms. Each process can only see the subcomm it belongs to */
40   mycolor = grank % 3;
41   CHKERRMPI(MPI_Comm_split(PETSC_COMM_WORLD,mycolor,grank,&subcomm));
42 
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     CHKERRQ(VecCreate(PETSC_COMM_WORLD, &x));
49     CHKERRQ(VecSetSizes(x, PETSC_DECIDE, N));
50     if (iscuda) {
51       CHKERRQ(VecSetType(x, VECCUDA));
52     } else {
53       CHKERRQ(VecSetType(x, VECSTANDARD));
54     }
55     CHKERRQ(VecSetUp(x));
56     CHKERRQ(PetscObjectSetName((PetscObject)x,"x_commworld")); /* Give a name to view x clearly */
57 
58     /* Initialize x to [-0.0, -1.0, -2.0, ..., -19.0] */
59     CHKERRQ(VecGetOwnershipRange(x,&low,&high));
60     for (i=low; i<high; i++) {
61       PetscScalar val = -i;
62       CHKERRQ(VecSetValue(x,i,val,INSERT_VALUES));
63     }
64     CHKERRQ(VecAssemblyBegin(x));
65     CHKERRQ(VecAssemblyEnd(x));
66 
67     /* Transfer x to a vector y only defined on subcomm0 and vice versa */
68     if (mycolor == 0) { /* subcomm0 contains ranks 0, 3, 6, ... in PETSC_COMM_WORLD */
69       Vec         y;
70       PetscScalar *yvalue;
71        CHKERRQ(VecCreate(subcomm, &y));
72       CHKERRQ(VecSetSizes(y, PETSC_DECIDE, N));
73       if (iscuda) {
74         CHKERRQ(VecSetType(y, VECCUDA));
75       } else {
76         CHKERRQ(VecSetType(y, VECSTANDARD));
77       }
78       CHKERRQ(VecSetUp(y));
79       CHKERRQ(PetscObjectSetName((PetscObject)y,"y_subcomm_0")); /* Give a name to view y clearly */
80       CHKERRQ(VecGetLocalSize(y,&n));
81       if (iscuda) {
82         #if defined(PETSC_HAVE_CUDA)
83           CHKERRQ(VecCUDAGetArray(y,&yvalue));
84         #endif
85       } else {
86         CHKERRQ(VecGetArray(y,&yvalue));
87       }
88       /* Create yg on PETSC_COMM_WORLD and alias yg with y. They share the memory pointed by yvalue.
89         Note this is a collective call. All processes have to call it and supply consistent N.
90       */
91       if (iscuda) {
92         #if defined(PETSC_HAVE_CUDA)
93           CHKERRQ(VecCreateMPICUDAWithArray(PETSC_COMM_WORLD,1,n,N,yvalue,&yg));
94         #endif
95       } else {
96         CHKERRQ(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,yvalue,&yg));
97       }
98 
99       /* Create an identity map that makes yg[i] = x[i], i=0..N-1 */
100       CHKERRQ(VecGetOwnershipRange(yg,&low,&high)); /* low, high are global indices */
101       CHKERRQ(ISCreateStride(PETSC_COMM_SELF,high-low,low,1,&ix));
102       CHKERRQ(ISDuplicate(ix,&iy));
103 
104       /* Union of ix's on subcomm0 covers the full range of [0,N) */
105       CHKERRQ(VecScatterCreate(x,ix,yg,iy,&vscat));
106       CHKERRQ(VecScatterBegin(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
107       CHKERRQ(VecScatterEnd(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
108 
109       /* Once yg got the data from x, we return yvalue to y so that we can use y in other operations.
110         VecGetArray must be paired with VecRestoreArray.
111       */
112       if (iscuda) {
113          #if defined(PETSC_HAVE_CUDA)
114            CHKERRQ(VecCUDARestoreArray(y,&yvalue));
115          #endif
116       } else {
117         CHKERRQ(VecRestoreArray(y,&yvalue));
118       }
119 
120       /* Libraries on subcomm0 can safely use y now, for example, view and scale it */
121       CHKERRQ(VecView(y,PETSC_VIEWER_STDOUT_(subcomm)));
122       CHKERRQ(VecScale(y,2.0));
123 
124       /* Send the new y back to x */
125       CHKERRQ(VecGetArray(y,&yvalue)); /* If VecScale is done on GPU, Petsc will prepare a valid yvalue for access */
126       /* Supply new yvalue to yg without memory copying */
127       CHKERRQ(VecPlaceArray(yg,yvalue));
128       CHKERRQ(VecScatterBegin(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE));
129       CHKERRQ(VecScatterEnd(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE));
130       CHKERRQ(VecResetArray(yg));
131       if (iscuda) {
132         #if defined(PETSC_HAVE_CUDA)
133           CHKERRQ(VecCUDARestoreArray(y,&yvalue));
134         #endif
135       } else {
136         CHKERRQ(VecRestoreArray(y,&yvalue));
137       }
138 
139       CHKERRQ(VecDestroy(&y));
140     } else {
141       /* Ranks outside of subcomm0 do not supply values to yg */
142       if (iscuda) {
143         #if defined(PETSC_HAVE_CUDA)
144           CHKERRQ(VecCreateMPICUDAWithArray(PETSC_COMM_WORLD,1,0/*n*/,N,NULL,&yg));
145         #endif
146       } else {
147         CHKERRQ(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,0/*n*/,N,NULL,&yg));
148       }
149 
150       /* Ranks in subcomm0 already specified the full range of the identity map. The remaining
151         ranks just need to create empty ISes to cheat VecScatterCreate.
152       */
153       CHKERRQ(ISCreateStride(PETSC_COMM_SELF,0,0,1,&ix));
154       CHKERRQ(ISDuplicate(ix,&iy));
155 
156       CHKERRQ(VecScatterCreate(x,ix,yg,iy,&vscat));
157       CHKERRQ(VecScatterBegin(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
158       CHKERRQ(VecScatterEnd(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
159 
160       /* Send the new y back to x. Ranks outside of subcomm0 actually have nothing to send.
161         But they have to call VecScatterBegin/End since these routines are collective.
162       */
163       CHKERRQ(VecScatterBegin(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE));
164       CHKERRQ(VecScatterEnd(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE));
165     }
166 
167     CHKERRQ(VecView(x,PETSC_VIEWER_STDOUT_WORLD));
168     CHKERRQ(ISDestroy(&ix));
169     CHKERRQ(ISDestroy(&iy));
170     CHKERRQ(VecDestroy(&x));
171     CHKERRQ(VecDestroy(&yg));
172     CHKERRQ(VecScatterDestroy(&vscat));
173   } /* world2sub */
174 
175   /*===========================================================================
176    *  Transfer a vector x defined on subcomm0 to a vector y defined on
177    *  subcomm1. The two subcomms are not overlapping and their union is
178    *  not necessarily equal to PETSC_COMM_WORLD.
179    *===========================================================================*/
180   if (sub2sub) {
181     if (mycolor == 0) {
182       /* Intentionally declare N as a local variable so that processes in subcomm1 do not know its value */
183       PetscInt          n,N = 22;
184       Vec               x,xg,yg;
185       IS                ix,iy;
186       VecScatter        vscat;
187       const PetscScalar *xvalue;
188       MPI_Comm          intercomm,parentcomm;
189       PetscMPIInt       lrank;
190 
191       CHKERRMPI(MPI_Comm_rank(subcomm,&lrank));
192       /* x is on subcomm */
193       CHKERRQ(VecCreate(subcomm, &x));
194       CHKERRQ(VecSetSizes(x, PETSC_DECIDE, N));
195       if (iscuda) {
196         CHKERRQ(VecSetType(x, VECCUDA));
197       } else {
198         CHKERRQ(VecSetType(x, VECSTANDARD));
199       }
200       CHKERRQ(VecSetUp(x));
201       CHKERRQ(VecGetOwnershipRange(x,&low,&high));
202 
203       /* initialize x = [0.0, 1.0, 2.0, ..., 21.0] */
204       for (i=low; i<high; i++) {
205         PetscScalar val = i;
206         CHKERRQ(VecSetValue(x,i,val,INSERT_VALUES));
207       }
208       CHKERRQ(VecAssemblyBegin(x));
209       CHKERRQ(VecAssemblyEnd(x));
210 
211       CHKERRMPI(MPI_Intercomm_create(subcomm,0,PETSC_COMM_WORLD/*peer_comm*/,1,100/*tag*/,&intercomm));
212 
213       /* Tell rank 0 of subcomm1 the global size of x */
214       if (!lrank) CHKERRMPI(MPI_Send(&N,1,MPIU_INT,0/*receiver's rank in remote comm, i.e., subcomm1*/,200/*tag*/,intercomm));
215 
216       /* Create an intracomm Petsc can work on. Ranks in subcomm0 are ordered before ranks in subcomm1 in parentcomm.
217         But this order actually does not matter, since what we care is vector y, which is defined on subcomm1.
218       */
219       CHKERRMPI(MPI_Intercomm_merge(intercomm,0/*low*/,&parentcomm));
220 
221       /* Create a vector xg on parentcomm, which shares memory with x */
222       CHKERRQ(VecGetLocalSize(x,&n));
223       if (iscuda) {
224         #if defined(PETSC_HAVE_CUDA)
225           CHKERRQ(VecCUDAGetArrayRead(x,&xvalue));
226           CHKERRQ(VecCreateMPICUDAWithArray(parentcomm,1,n,N,xvalue,&xg));
227         #endif
228       } else {
229         CHKERRQ(VecGetArrayRead(x,&xvalue));
230         CHKERRQ(VecCreateMPIWithArray(parentcomm,1,n,N,xvalue,&xg));
231       }
232 
233       /* Ranks in subcomm 0 have nothing on yg, so they simply have n=0, array=NULL */
234       if (iscuda) {
235         #if defined(PETSC_HAVE_CUDA)
236           CHKERRQ(VecCreateMPICUDAWithArray(parentcomm,1,0/*n*/,N,NULL/*array*/,&yg));
237         #endif
238       } else {
239         CHKERRQ(VecCreateMPIWithArray(parentcomm,1,0/*n*/,N,NULL/*array*/,&yg));
240       }
241 
242       /* Create the vecscatter, which does identity map by setting yg[i] = xg[i], i=0..N-1. */
243       CHKERRQ(VecGetOwnershipRange(xg,&low,&high)); /* low, high are global indices of xg */
244       CHKERRQ(ISCreateStride(PETSC_COMM_SELF,high-low,low,1,&ix));
245       CHKERRQ(ISDuplicate(ix,&iy));
246       CHKERRQ(VecScatterCreate(xg,ix,yg,iy,&vscat));
247 
248       /* Scatter values from xg to yg */
249       CHKERRQ(VecScatterBegin(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD));
250       CHKERRQ(VecScatterEnd(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD));
251 
252       /* After the VecScatter is done, xg is idle so we can safely return xvalue to x */
253       if (iscuda) {
254         #if defined(PETSC_HAVE_CUDA)
255           CHKERRQ(VecCUDARestoreArrayRead(x,&xvalue));
256         #endif
257       } else {
258         CHKERRQ(VecRestoreArrayRead(x,&xvalue));
259       }
260       CHKERRQ(VecDestroy(&x));
261       CHKERRQ(ISDestroy(&ix));
262       CHKERRQ(ISDestroy(&iy));
263       CHKERRQ(VecDestroy(&xg));
264       CHKERRQ(VecDestroy(&yg));
265       CHKERRQ(VecScatterDestroy(&vscat));
266       CHKERRMPI(MPI_Comm_free(&intercomm));
267       CHKERRMPI(MPI_Comm_free(&parentcomm));
268     } else if (mycolor == 1) { /* subcomm 1, containing ranks 1, 4, 7, ... in PETSC_COMM_WORLD */
269       PetscInt    n,N;
270       Vec         y,xg,yg;
271       IS          ix,iy;
272       VecScatter  vscat;
273       PetscScalar *yvalue;
274       MPI_Comm    intercomm,parentcomm;
275       PetscMPIInt lrank;
276 
277       CHKERRMPI(MPI_Comm_rank(subcomm,&lrank));
278       CHKERRMPI(MPI_Intercomm_create(subcomm,0,PETSC_COMM_WORLD/*peer_comm*/,0/*remote_leader*/,100/*tag*/,&intercomm));
279 
280       /* Two rank-0 are talking */
281       if (!lrank) CHKERRMPI(MPI_Recv(&N,1,MPIU_INT,0/*sender's rank in remote comm, i.e. subcomm0*/,200/*tag*/,intercomm,MPI_STATUS_IGNORE));
282       /* Rank 0 of subcomm1 bcasts N to its members */
283       CHKERRMPI(MPI_Bcast(&N,1,MPIU_INT,0/*local root*/,subcomm));
284 
285       /* Create a intracomm Petsc can work on */
286       CHKERRMPI(MPI_Intercomm_merge(intercomm,1/*high*/,&parentcomm));
287 
288       /* Ranks in subcomm1 have nothing on xg, so they simply have n=0, array=NULL.*/
289       if (iscuda) {
290         #if defined(PETSC_HAVE_CUDA)
291           CHKERRQ(VecCreateMPICUDAWithArray(parentcomm,1/*bs*/,0/*n*/,N,NULL/*array*/,&xg));
292         #endif
293       } else {
294         CHKERRQ(VecCreateMPIWithArray(parentcomm,1/*bs*/,0/*n*/,N,NULL/*array*/,&xg));
295       }
296 
297       CHKERRQ(VecCreate(subcomm, &y));
298       CHKERRQ(VecSetSizes(y, PETSC_DECIDE, N));
299       if (iscuda) {
300         CHKERRQ(VecSetType(y, VECCUDA));
301       } else {
302         CHKERRQ(VecSetType(y, VECSTANDARD));
303       }
304       CHKERRQ(VecSetUp(y));
305 
306       CHKERRQ(PetscObjectSetName((PetscObject)y,"y_subcomm_1")); /* Give a name to view y clearly */
307       CHKERRQ(VecGetLocalSize(y,&n));
308       if (iscuda) {
309         #if defined(PETSC_HAVE_CUDA)
310           CHKERRQ(VecCUDAGetArray(y,&yvalue));
311         #endif
312       } else {
313         CHKERRQ(VecGetArray(y,&yvalue));
314       }
315       /* Create a vector yg on parentcomm, which shares memory with y. xg and yg must be
316         created in the same order in subcomm0/1. For example, we can not reverse the order of
317         creating xg and yg in subcomm1.
318       */
319       if (iscuda) {
320         #if defined(PETSC_HAVE_CUDA)
321           CHKERRQ(VecCreateMPICUDAWithArray(parentcomm,1/*bs*/,n,N,yvalue,&yg));
322         #endif
323       } else {
324         CHKERRQ(VecCreateMPIWithArray(parentcomm,1/*bs*/,n,N,yvalue,&yg));
325       }
326 
327       /* Ranks in subcomm0 already specified the full range of the identity map.
328         ranks in subcomm1 just need to create empty ISes to cheat VecScatterCreate.
329       */
330       CHKERRQ(ISCreateStride(PETSC_COMM_SELF,0,0,1,&ix));
331       CHKERRQ(ISDuplicate(ix,&iy));
332       CHKERRQ(VecScatterCreate(xg,ix,yg,iy,&vscat));
333 
334       /* Scatter values from xg to yg */
335       CHKERRQ(VecScatterBegin(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD));
336       CHKERRQ(VecScatterEnd(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD));
337 
338       /* After the VecScatter is done, values in yg are available. y is our interest, so we return yvalue to y */
339       if (iscuda) {
340         #if defined(PETSC_HAVE_CUDA)
341           CHKERRQ(VecCUDARestoreArray(y,&yvalue));
342         #endif
343       } else {
344         CHKERRQ(VecRestoreArray(y,&yvalue));
345       }
346 
347       /* Libraries on subcomm1 can safely use y now, for example, view it */
348       CHKERRQ(VecView(y,PETSC_VIEWER_STDOUT_(subcomm)));
349 
350       CHKERRQ(VecDestroy(&y));
351       CHKERRQ(ISDestroy(&ix));
352       CHKERRQ(ISDestroy(&iy));
353       CHKERRQ(VecDestroy(&xg));
354       CHKERRQ(VecDestroy(&yg));
355       CHKERRQ(VecScatterDestroy(&vscat));
356       CHKERRMPI(MPI_Comm_free(&intercomm));
357       CHKERRMPI(MPI_Comm_free(&parentcomm));
358     } else if (mycolor == 2) { /* subcomm2 */
359       /* Processes in subcomm2 do not participate in the VecScatter. They can freely do unrelated things on subcomm2 */
360     }
361   } /* sub2sub */
362 
363   /*===========================================================================
364    *  Copy a vector x defined on PETSC_COMM_WORLD to vectors y defined on
365    *  every subcommunicator of PETSC_COMM_WORLD. We could use multiple transfers
366    *  as we did in case 1, but that is not efficient. Instead, we use one vecscatter
367    *  to achieve that.
368    *===========================================================================*/
369   if (world2subs) {
370     Vec         y;
371     PetscInt    n,N=15,xstart,ystart,low,high;
372     PetscScalar *yvalue;
373 
374     /* Initialize x to [0, 1, 2, 3, ..., N-1] */
375     CHKERRQ(VecCreate(PETSC_COMM_WORLD, &x));
376     CHKERRQ(VecSetSizes(x, PETSC_DECIDE, N));
377     if (iscuda) {
378       CHKERRQ(VecSetType(x, VECCUDA));
379     } else {
380       CHKERRQ(VecSetType(x, VECSTANDARD));
381     }
382     CHKERRQ(VecSetUp(x));
383     CHKERRQ(VecGetOwnershipRange(x,&low,&high));
384     for (i=low; i<high; i++) CHKERRQ(VecSetValue(x,i,(PetscScalar)i,INSERT_VALUES));
385     CHKERRQ(VecAssemblyBegin(x));
386     CHKERRQ(VecAssemblyEnd(x));
387 
388     /* Every subcomm has a y as long as x */
389     CHKERRQ(VecCreate(subcomm, &y));
390     CHKERRQ(VecSetSizes(y, PETSC_DECIDE, N));
391     if (iscuda) {
392       CHKERRQ(VecSetType(y, VECCUDA));
393     } else {
394       CHKERRQ(VecSetType(y, VECSTANDARD));
395     }
396     CHKERRQ(VecSetUp(y));
397     CHKERRQ(VecGetLocalSize(y,&n));
398 
399     /* Create a global vector yg on PETSC_COMM_WORLD using y's memory. yg's global size = N*(number of subcommunicators).
400        Eeach rank in subcomms contributes a piece to construct the global yg. Keep in mind that pieces from a subcomm are not
401        necessarily consecutive in yg. That depends on how PETSC_COMM_WORLD is split. In our case, subcomm0 is made of rank
402        0, 3, 6 etc from PETSC_COMM_WORLD. So subcomm0's pieces are interleaved with pieces from other subcomms in yg.
403     */
404     if (iscuda) {
405       #if defined(PETSC_HAVE_CUDA)
406         CHKERRQ(VecCUDAGetArray(y,&yvalue));
407         CHKERRQ(VecCreateMPICUDAWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,yvalue,&yg));
408       #endif
409     } else {
410       CHKERRQ(VecGetArray(y,&yvalue));
411       CHKERRQ(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,yvalue,&yg));
412     }
413     CHKERRQ(PetscObjectSetName((PetscObject)yg,"yg_on_subcomms")); /* Give a name to view yg clearly */
414 
415     /* The following two lines are key. From xstart, we know where to pull entries from x. Note that we get xstart from y,
416        since first entry of y on this rank is from x[xstart]. From ystart, we know where ot put entries to yg.
417      */
418     CHKERRQ(VecGetOwnershipRange(y,&xstart,NULL));
419     CHKERRQ(VecGetOwnershipRange(yg,&ystart,NULL));
420 
421     CHKERRQ(ISCreateStride(PETSC_COMM_SELF,n,xstart,1,&ix));
422     CHKERRQ(ISCreateStride(PETSC_COMM_SELF,n,ystart,1,&iy));
423     CHKERRQ(VecScatterCreate(x,ix,yg,iy,&vscat));
424     CHKERRQ(VecScatterBegin(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
425     CHKERRQ(VecScatterEnd(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
426 
427     /* View yg on PETSC_COMM_WORLD before destroying it. We shall see the interleaving effect in output. */
428     CHKERRQ(VecView(yg,PETSC_VIEWER_STDOUT_WORLD));
429     CHKERRQ(VecDestroy(&yg));
430 
431     /* Restory yvalue so that processes in subcomm can use y from now on. */
432     if (iscuda) {
433       #if defined(PETSC_HAVE_CUDA)
434         CHKERRQ(VecCUDARestoreArray(y,&yvalue));
435       #endif
436     } else {
437       CHKERRQ(VecRestoreArray(y,&yvalue));
438     }
439     CHKERRQ(VecScale(y,3.0));
440 
441     CHKERRQ(ISDestroy(&ix)); /* One can also destroy ix, iy immediately after VecScatterCreate() */
442     CHKERRQ(ISDestroy(&iy));
443     CHKERRQ(VecDestroy(&x));
444     CHKERRQ(VecDestroy(&y));
445     CHKERRQ(VecScatterDestroy(&vscat));
446   } /* world2subs */
447 
448   CHKERRMPI(MPI_Comm_free(&subcomm));
449   ierr = PetscFinalize();
450   return ierr;
451 }
452 
453 /*TEST
454 
455    build:
456      requires: !defined(PETSC_HAVE_MPIUNI)
457 
458    testset:
459      nsize: 7
460 
461      test:
462        suffix: 1
463        args: -world2sub
464 
465      test:
466        suffix: 2
467        args: -sub2sub
468        # deadlocks with NECMPI and INTELMPI (20210400300)
469        requires: !defined(PETSC_HAVE_NECMPI) !defined(PETSC_HAVE_I_MPI_NUMVERSION)
470 
471      test:
472        suffix: 3
473        args: -world2subs
474 
475      test:
476        suffix: 4
477        args: -world2sub -vectype cuda
478        requires: cuda
479 
480      test:
481        suffix: 5
482        args: -sub2sub -vectype cuda
483        requires: cuda
484 
485      test:
486       suffix: 6
487       args: -world2subs -vectype cuda
488       requires: cuda
489 
490      test:
491        suffix: 7
492        args: -world2sub -sf_type neighbor
493        output_file: output/ex9_1.out
494        # OpenMPI has a bug wrt MPI_Neighbor_alltoallv etc (https://github.com/open-mpi/ompi/pull/6782). Once the patch is in, we can remove !define(PETSC_HAVE_OMPI_MAJOR_VERSION)
495        # segfaults with NECMPI
496        requires: defined(PETSC_HAVE_MPI_NEIGHBORHOOD_COLLECTIVES) !defined(PETSC_HAVE_OMPI_MAJOR_VERSION) !defined(PETSC_HAVE_NECMPI)
497 TEST*/
498