xref: /petsc/src/vec/is/sf/tests/ex9.c (revision 327415f76d85372a4417cf1aaa14db707d4d6c04)
197929ea7SJunchao Zhang 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\
297929ea7SJunchao Zhang   2) how to transfer vectors from a subcommunicator to vectors on another subcommunicator. The two subcommunicators are not\n\
38a53a0a4SSajid Ali   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\
48a53a0a4SSajid Ali   To run any example with VECCUDA vectors, add -vectype cuda to the argument list\n\n";
597929ea7SJunchao Zhang 
697929ea7SJunchao Zhang #include <petscvec.h>
797929ea7SJunchao Zhang int main(int argc,char **argv)
897929ea7SJunchao Zhang {
997929ea7SJunchao Zhang   PetscMPIInt    nproc,grank,mycolor;
1097929ea7SJunchao Zhang   PetscInt       i,n,N = 20,low,high;
1197929ea7SJunchao Zhang   MPI_Comm       subcomm;
128a53a0a4SSajid Ali   Vec            x  = PETSC_NULL; /* global vectors on PETSC_COMM_WORLD */
138a53a0a4SSajid Ali   Vec            yg = PETSC_NULL; /* global vectors on PETSC_COMM_WORLD */
1497929ea7SJunchao Zhang   VecScatter     vscat;
1597929ea7SJunchao Zhang   IS             ix,iy;
168a53a0a4SSajid Ali   PetscBool      iscuda = PETSC_FALSE;      /* Option to use VECCUDA vectors */
178a53a0a4SSajid Ali   PetscBool      optionflag, compareflag;
188a53a0a4SSajid Ali   char           vectypename[PETSC_MAX_PATH_LEN];
1997929ea7SJunchao Zhang   PetscBool      world2sub  = PETSC_FALSE;  /* Copy a vector from WORLD to a subcomm? */
2097929ea7SJunchao Zhang   PetscBool      sub2sub    = PETSC_FALSE;  /* Copy a vector from a subcomm to another subcomm? */
2197929ea7SJunchao Zhang   PetscBool      world2subs = PETSC_FALSE;  /* Copy a vector from WORLD to multiple subcomms? */
2297929ea7SJunchao Zhang 
23*327415f7SBarry Smith   PetscFunctionBeginUser;
249566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc,&argv,(char*)0,help));
259566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&nproc));
269566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&grank));
2797929ea7SJunchao Zhang 
2808401ef6SPierre Jolivet   PetscCheck(nproc >= 2,PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"This test must have at least two processes to run");
2997929ea7SJunchao Zhang 
309566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,0,"-world2sub",&world2sub,NULL));
319566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,0,"-sub2sub",&sub2sub,NULL));
329566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL,0,"-world2subs",&world2subs,NULL));
339566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetString(NULL,NULL,"-vectype",vectypename,sizeof(vectypename),&optionflag));
348a53a0a4SSajid Ali   if (optionflag) {
359566063dSJacob Faibussowitsch     PetscCall(PetscStrncmp(vectypename, "cuda", (size_t)4, &compareflag));
368a53a0a4SSajid Ali     if (compareflag) iscuda = PETSC_TRUE;
378a53a0a4SSajid Ali   }
3897929ea7SJunchao Zhang 
3997929ea7SJunchao Zhang   /* Split PETSC_COMM_WORLD into three subcomms. Each process can only see the subcomm it belongs to */
4097929ea7SJunchao Zhang   mycolor = grank % 3;
419566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD,mycolor,grank,&subcomm));
4297929ea7SJunchao Zhang 
4397929ea7SJunchao Zhang   /*===========================================================================
4497929ea7SJunchao Zhang    *  Transfer a vector x defined on PETSC_COMM_WORLD to a vector y defined on
4597929ea7SJunchao Zhang    *  a subcommunicator of PETSC_COMM_WORLD and vice versa.
4697929ea7SJunchao Zhang    *===========================================================================*/
4797929ea7SJunchao Zhang   if (world2sub) {
489566063dSJacob Faibussowitsch     PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
499566063dSJacob Faibussowitsch     PetscCall(VecSetSizes(x, PETSC_DECIDE, N));
508a53a0a4SSajid Ali     if (iscuda) {
519566063dSJacob Faibussowitsch       PetscCall(VecSetType(x, VECCUDA));
528a53a0a4SSajid Ali     } else {
539566063dSJacob Faibussowitsch       PetscCall(VecSetType(x, VECSTANDARD));
548a53a0a4SSajid Ali     }
559566063dSJacob Faibussowitsch     PetscCall(VecSetUp(x));
569566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)x,"x_commworld")); /* Give a name to view x clearly */
5797929ea7SJunchao Zhang 
5897929ea7SJunchao Zhang     /* Initialize x to [-0.0, -1.0, -2.0, ..., -19.0] */
599566063dSJacob Faibussowitsch     PetscCall(VecGetOwnershipRange(x,&low,&high));
6097929ea7SJunchao Zhang     for (i=low; i<high; i++) {
6197929ea7SJunchao Zhang       PetscScalar val = -i;
629566063dSJacob Faibussowitsch       PetscCall(VecSetValue(x,i,val,INSERT_VALUES));
6397929ea7SJunchao Zhang     }
649566063dSJacob Faibussowitsch     PetscCall(VecAssemblyBegin(x));
659566063dSJacob Faibussowitsch     PetscCall(VecAssemblyEnd(x));
6697929ea7SJunchao Zhang 
6797929ea7SJunchao Zhang     /* Transfer x to a vector y only defined on subcomm0 and vice versa */
6897929ea7SJunchao Zhang     if (mycolor == 0) { /* subcomm0 contains ranks 0, 3, 6, ... in PETSC_COMM_WORLD */
6997929ea7SJunchao Zhang       Vec         y;
7097929ea7SJunchao Zhang       PetscScalar *yvalue;
719566063dSJacob Faibussowitsch        PetscCall(VecCreate(subcomm, &y));
729566063dSJacob Faibussowitsch       PetscCall(VecSetSizes(y, PETSC_DECIDE, N));
738a53a0a4SSajid Ali       if (iscuda) {
749566063dSJacob Faibussowitsch         PetscCall(VecSetType(y, VECCUDA));
758a53a0a4SSajid Ali       } else {
769566063dSJacob Faibussowitsch         PetscCall(VecSetType(y, VECSTANDARD));
778a53a0a4SSajid Ali       }
789566063dSJacob Faibussowitsch       PetscCall(VecSetUp(y));
799566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)y,"y_subcomm_0")); /* Give a name to view y clearly */
809566063dSJacob Faibussowitsch       PetscCall(VecGetLocalSize(y,&n));
818a53a0a4SSajid Ali       if (iscuda) {
828a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
839566063dSJacob Faibussowitsch           PetscCall(VecCUDAGetArray(y,&yvalue));
848a53a0a4SSajid Ali         #endif
858a53a0a4SSajid Ali       } else {
869566063dSJacob Faibussowitsch         PetscCall(VecGetArray(y,&yvalue));
878a53a0a4SSajid Ali       }
8897929ea7SJunchao Zhang       /* Create yg on PETSC_COMM_WORLD and alias yg with y. They share the memory pointed by yvalue.
8997929ea7SJunchao Zhang         Note this is a collective call. All processes have to call it and supply consistent N.
9097929ea7SJunchao Zhang       */
918a53a0a4SSajid Ali       if (iscuda) {
928a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
939566063dSJacob Faibussowitsch           PetscCall(VecCreateMPICUDAWithArray(PETSC_COMM_WORLD,1,n,N,yvalue,&yg));
948a53a0a4SSajid Ali         #endif
958a53a0a4SSajid Ali       } else {
969566063dSJacob Faibussowitsch         PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,yvalue,&yg));
978a53a0a4SSajid Ali       }
9897929ea7SJunchao Zhang 
9997929ea7SJunchao Zhang       /* Create an identity map that makes yg[i] = x[i], i=0..N-1 */
1009566063dSJacob Faibussowitsch       PetscCall(VecGetOwnershipRange(yg,&low,&high)); /* low, high are global indices */
1019566063dSJacob Faibussowitsch       PetscCall(ISCreateStride(PETSC_COMM_SELF,high-low,low,1,&ix));
1029566063dSJacob Faibussowitsch       PetscCall(ISDuplicate(ix,&iy));
10397929ea7SJunchao Zhang 
10497929ea7SJunchao Zhang       /* Union of ix's on subcomm0 covers the full range of [0,N) */
1059566063dSJacob Faibussowitsch       PetscCall(VecScatterCreate(x,ix,yg,iy,&vscat));
1069566063dSJacob Faibussowitsch       PetscCall(VecScatterBegin(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
1079566063dSJacob Faibussowitsch       PetscCall(VecScatterEnd(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
10897929ea7SJunchao Zhang 
10997929ea7SJunchao Zhang       /* Once yg got the data from x, we return yvalue to y so that we can use y in other operations.
11097929ea7SJunchao Zhang         VecGetArray must be paired with VecRestoreArray.
11197929ea7SJunchao Zhang       */
1128a53a0a4SSajid Ali       if (iscuda) {
1138a53a0a4SSajid Ali          #if defined(PETSC_HAVE_CUDA)
1149566063dSJacob Faibussowitsch            PetscCall(VecCUDARestoreArray(y,&yvalue));
1158a53a0a4SSajid Ali          #endif
1168a53a0a4SSajid Ali       } else {
1179566063dSJacob Faibussowitsch         PetscCall(VecRestoreArray(y,&yvalue));
1188a53a0a4SSajid Ali       }
11997929ea7SJunchao Zhang 
12097929ea7SJunchao Zhang       /* Libraries on subcomm0 can safely use y now, for example, view and scale it */
1219566063dSJacob Faibussowitsch       PetscCall(VecView(y,PETSC_VIEWER_STDOUT_(subcomm)));
1229566063dSJacob Faibussowitsch       PetscCall(VecScale(y,2.0));
12397929ea7SJunchao Zhang 
12497929ea7SJunchao Zhang       /* Send the new y back to x */
1259566063dSJacob Faibussowitsch       PetscCall(VecGetArray(y,&yvalue)); /* If VecScale is done on GPU, Petsc will prepare a valid yvalue for access */
12697929ea7SJunchao Zhang       /* Supply new yvalue to yg without memory copying */
1279566063dSJacob Faibussowitsch       PetscCall(VecPlaceArray(yg,yvalue));
1289566063dSJacob Faibussowitsch       PetscCall(VecScatterBegin(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE));
1299566063dSJacob Faibussowitsch       PetscCall(VecScatterEnd(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE));
1309566063dSJacob Faibussowitsch       PetscCall(VecResetArray(yg));
1318a53a0a4SSajid Ali       if (iscuda) {
1328a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
1339566063dSJacob Faibussowitsch           PetscCall(VecCUDARestoreArray(y,&yvalue));
1348a53a0a4SSajid Ali         #endif
1358a53a0a4SSajid Ali       } else {
1369566063dSJacob Faibussowitsch         PetscCall(VecRestoreArray(y,&yvalue));
1378a53a0a4SSajid Ali       }
13897929ea7SJunchao Zhang 
1399566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&y));
14097929ea7SJunchao Zhang     } else {
14197929ea7SJunchao Zhang       /* Ranks outside of subcomm0 do not supply values to yg */
1428a53a0a4SSajid Ali       if (iscuda) {
1438a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
1449566063dSJacob Faibussowitsch           PetscCall(VecCreateMPICUDAWithArray(PETSC_COMM_WORLD,1,0/*n*/,N,NULL,&yg));
1458a53a0a4SSajid Ali         #endif
1468a53a0a4SSajid Ali       } else {
1479566063dSJacob Faibussowitsch         PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,0/*n*/,N,NULL,&yg));
1488a53a0a4SSajid Ali       }
14997929ea7SJunchao Zhang 
15097929ea7SJunchao Zhang       /* Ranks in subcomm0 already specified the full range of the identity map. The remaining
15197929ea7SJunchao Zhang         ranks just need to create empty ISes to cheat VecScatterCreate.
15297929ea7SJunchao Zhang       */
1539566063dSJacob Faibussowitsch       PetscCall(ISCreateStride(PETSC_COMM_SELF,0,0,1,&ix));
1549566063dSJacob Faibussowitsch       PetscCall(ISDuplicate(ix,&iy));
15597929ea7SJunchao Zhang 
1569566063dSJacob Faibussowitsch       PetscCall(VecScatterCreate(x,ix,yg,iy,&vscat));
1579566063dSJacob Faibussowitsch       PetscCall(VecScatterBegin(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
1589566063dSJacob Faibussowitsch       PetscCall(VecScatterEnd(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
15997929ea7SJunchao Zhang 
16097929ea7SJunchao Zhang       /* Send the new y back to x. Ranks outside of subcomm0 actually have nothing to send.
16197929ea7SJunchao Zhang         But they have to call VecScatterBegin/End since these routines are collective.
16297929ea7SJunchao Zhang       */
1639566063dSJacob Faibussowitsch       PetscCall(VecScatterBegin(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE));
1649566063dSJacob Faibussowitsch       PetscCall(VecScatterEnd(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE));
16597929ea7SJunchao Zhang     }
16697929ea7SJunchao Zhang 
1679566063dSJacob Faibussowitsch     PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));
1689566063dSJacob Faibussowitsch     PetscCall(ISDestroy(&ix));
1699566063dSJacob Faibussowitsch     PetscCall(ISDestroy(&iy));
1709566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&x));
1719566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&yg));
1729566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&vscat));
17397929ea7SJunchao Zhang   } /* world2sub */
17497929ea7SJunchao Zhang 
17597929ea7SJunchao Zhang   /*===========================================================================
17697929ea7SJunchao Zhang    *  Transfer a vector x defined on subcomm0 to a vector y defined on
17797929ea7SJunchao Zhang    *  subcomm1. The two subcomms are not overlapping and their union is
17897929ea7SJunchao Zhang    *  not necessarily equal to PETSC_COMM_WORLD.
17997929ea7SJunchao Zhang    *===========================================================================*/
18097929ea7SJunchao Zhang   if (sub2sub) {
18197929ea7SJunchao Zhang     if (mycolor == 0) {
18297929ea7SJunchao Zhang       /* Intentionally declare N as a local variable so that processes in subcomm1 do not know its value */
18397929ea7SJunchao Zhang       PetscInt          n,N = 22;
18497929ea7SJunchao Zhang       Vec               x,xg,yg;
18597929ea7SJunchao Zhang       IS                ix,iy;
18697929ea7SJunchao Zhang       VecScatter        vscat;
18797929ea7SJunchao Zhang       const PetscScalar *xvalue;
18897929ea7SJunchao Zhang       MPI_Comm          intercomm,parentcomm;
18997929ea7SJunchao Zhang       PetscMPIInt       lrank;
19097929ea7SJunchao Zhang 
1919566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Comm_rank(subcomm,&lrank));
1928a53a0a4SSajid Ali       /* x is on subcomm */
1939566063dSJacob Faibussowitsch       PetscCall(VecCreate(subcomm, &x));
1949566063dSJacob Faibussowitsch       PetscCall(VecSetSizes(x, PETSC_DECIDE, N));
1958a53a0a4SSajid Ali       if (iscuda) {
1969566063dSJacob Faibussowitsch         PetscCall(VecSetType(x, VECCUDA));
1978a53a0a4SSajid Ali       } else {
1989566063dSJacob Faibussowitsch         PetscCall(VecSetType(x, VECSTANDARD));
1998a53a0a4SSajid Ali       }
2009566063dSJacob Faibussowitsch       PetscCall(VecSetUp(x));
2019566063dSJacob Faibussowitsch       PetscCall(VecGetOwnershipRange(x,&low,&high));
20297929ea7SJunchao Zhang 
20397929ea7SJunchao Zhang       /* initialize x = [0.0, 1.0, 2.0, ..., 21.0] */
20497929ea7SJunchao Zhang       for (i=low; i<high; i++) {
20597929ea7SJunchao Zhang         PetscScalar val = i;
2069566063dSJacob Faibussowitsch         PetscCall(VecSetValue(x,i,val,INSERT_VALUES));
20797929ea7SJunchao Zhang       }
2089566063dSJacob Faibussowitsch       PetscCall(VecAssemblyBegin(x));
2099566063dSJacob Faibussowitsch       PetscCall(VecAssemblyEnd(x));
21097929ea7SJunchao Zhang 
2119566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Intercomm_create(subcomm,0,PETSC_COMM_WORLD/*peer_comm*/,1,100/*tag*/,&intercomm));
21297929ea7SJunchao Zhang 
21397929ea7SJunchao Zhang       /* Tell rank 0 of subcomm1 the global size of x */
2149566063dSJacob Faibussowitsch       if (!lrank) PetscCallMPI(MPI_Send(&N,1,MPIU_INT,0/*receiver's rank in remote comm, i.e., subcomm1*/,200/*tag*/,intercomm));
21597929ea7SJunchao Zhang 
21697929ea7SJunchao Zhang       /* Create an intracomm Petsc can work on. Ranks in subcomm0 are ordered before ranks in subcomm1 in parentcomm.
21797929ea7SJunchao Zhang         But this order actually does not matter, since what we care is vector y, which is defined on subcomm1.
21897929ea7SJunchao Zhang       */
2199566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Intercomm_merge(intercomm,0/*low*/,&parentcomm));
22097929ea7SJunchao Zhang 
22197929ea7SJunchao Zhang       /* Create a vector xg on parentcomm, which shares memory with x */
2229566063dSJacob Faibussowitsch       PetscCall(VecGetLocalSize(x,&n));
2238a53a0a4SSajid Ali       if (iscuda) {
2248a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
2259566063dSJacob Faibussowitsch           PetscCall(VecCUDAGetArrayRead(x,&xvalue));
2269566063dSJacob Faibussowitsch           PetscCall(VecCreateMPICUDAWithArray(parentcomm,1,n,N,xvalue,&xg));
2278a53a0a4SSajid Ali         #endif
2288a53a0a4SSajid Ali       } else {
2299566063dSJacob Faibussowitsch         PetscCall(VecGetArrayRead(x,&xvalue));
2309566063dSJacob Faibussowitsch         PetscCall(VecCreateMPIWithArray(parentcomm,1,n,N,xvalue,&xg));
2318a53a0a4SSajid Ali       }
23297929ea7SJunchao Zhang 
23397929ea7SJunchao Zhang       /* Ranks in subcomm 0 have nothing on yg, so they simply have n=0, array=NULL */
2348a53a0a4SSajid Ali       if (iscuda) {
2358a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
2369566063dSJacob Faibussowitsch           PetscCall(VecCreateMPICUDAWithArray(parentcomm,1,0/*n*/,N,NULL/*array*/,&yg));
2378a53a0a4SSajid Ali         #endif
2388a53a0a4SSajid Ali       } else {
2399566063dSJacob Faibussowitsch         PetscCall(VecCreateMPIWithArray(parentcomm,1,0/*n*/,N,NULL/*array*/,&yg));
2408a53a0a4SSajid Ali       }
24197929ea7SJunchao Zhang 
24297929ea7SJunchao Zhang       /* Create the vecscatter, which does identity map by setting yg[i] = xg[i], i=0..N-1. */
2439566063dSJacob Faibussowitsch       PetscCall(VecGetOwnershipRange(xg,&low,&high)); /* low, high are global indices of xg */
2449566063dSJacob Faibussowitsch       PetscCall(ISCreateStride(PETSC_COMM_SELF,high-low,low,1,&ix));
2459566063dSJacob Faibussowitsch       PetscCall(ISDuplicate(ix,&iy));
2469566063dSJacob Faibussowitsch       PetscCall(VecScatterCreate(xg,ix,yg,iy,&vscat));
24797929ea7SJunchao Zhang 
24897929ea7SJunchao Zhang       /* Scatter values from xg to yg */
2499566063dSJacob Faibussowitsch       PetscCall(VecScatterBegin(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD));
2509566063dSJacob Faibussowitsch       PetscCall(VecScatterEnd(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD));
25197929ea7SJunchao Zhang 
25297929ea7SJunchao Zhang       /* After the VecScatter is done, xg is idle so we can safely return xvalue to x */
2538a53a0a4SSajid Ali       if (iscuda) {
2548a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
2559566063dSJacob Faibussowitsch           PetscCall(VecCUDARestoreArrayRead(x,&xvalue));
2568a53a0a4SSajid Ali         #endif
2578a53a0a4SSajid Ali       } else {
2589566063dSJacob Faibussowitsch         PetscCall(VecRestoreArrayRead(x,&xvalue));
2598a53a0a4SSajid Ali       }
2609566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&x));
2619566063dSJacob Faibussowitsch       PetscCall(ISDestroy(&ix));
2629566063dSJacob Faibussowitsch       PetscCall(ISDestroy(&iy));
2639566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&xg));
2649566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&yg));
2659566063dSJacob Faibussowitsch       PetscCall(VecScatterDestroy(&vscat));
2669566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Comm_free(&intercomm));
2679566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Comm_free(&parentcomm));
26897929ea7SJunchao Zhang     } else if (mycolor == 1) { /* subcomm 1, containing ranks 1, 4, 7, ... in PETSC_COMM_WORLD */
26997929ea7SJunchao Zhang       PetscInt    n,N;
27097929ea7SJunchao Zhang       Vec         y,xg,yg;
27197929ea7SJunchao Zhang       IS          ix,iy;
27297929ea7SJunchao Zhang       VecScatter  vscat;
27397929ea7SJunchao Zhang       PetscScalar *yvalue;
27497929ea7SJunchao Zhang       MPI_Comm    intercomm,parentcomm;
27597929ea7SJunchao Zhang       PetscMPIInt lrank;
27697929ea7SJunchao Zhang 
2779566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Comm_rank(subcomm,&lrank));
2789566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Intercomm_create(subcomm,0,PETSC_COMM_WORLD/*peer_comm*/,0/*remote_leader*/,100/*tag*/,&intercomm));
27997929ea7SJunchao Zhang 
28097929ea7SJunchao Zhang       /* Two rank-0 are talking */
2819566063dSJacob Faibussowitsch       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));
28297929ea7SJunchao Zhang       /* Rank 0 of subcomm1 bcasts N to its members */
2839566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Bcast(&N,1,MPIU_INT,0/*local root*/,subcomm));
28497929ea7SJunchao Zhang 
28597929ea7SJunchao Zhang       /* Create a intracomm Petsc can work on */
2869566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Intercomm_merge(intercomm,1/*high*/,&parentcomm));
28797929ea7SJunchao Zhang 
28897929ea7SJunchao Zhang       /* Ranks in subcomm1 have nothing on xg, so they simply have n=0, array=NULL.*/
2898a53a0a4SSajid Ali       if (iscuda) {
2908a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
2919566063dSJacob Faibussowitsch           PetscCall(VecCreateMPICUDAWithArray(parentcomm,1/*bs*/,0/*n*/,N,NULL/*array*/,&xg));
2928a53a0a4SSajid Ali         #endif
2938a53a0a4SSajid Ali       } else {
2949566063dSJacob Faibussowitsch         PetscCall(VecCreateMPIWithArray(parentcomm,1/*bs*/,0/*n*/,N,NULL/*array*/,&xg));
2958a53a0a4SSajid Ali       }
29697929ea7SJunchao Zhang 
2979566063dSJacob Faibussowitsch       PetscCall(VecCreate(subcomm, &y));
2989566063dSJacob Faibussowitsch       PetscCall(VecSetSizes(y, PETSC_DECIDE, N));
2998a53a0a4SSajid Ali       if (iscuda) {
3009566063dSJacob Faibussowitsch         PetscCall(VecSetType(y, VECCUDA));
3018a53a0a4SSajid Ali       } else {
3029566063dSJacob Faibussowitsch         PetscCall(VecSetType(y, VECSTANDARD));
3038a53a0a4SSajid Ali       }
3049566063dSJacob Faibussowitsch       PetscCall(VecSetUp(y));
3058a53a0a4SSajid Ali 
3069566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)y,"y_subcomm_1")); /* Give a name to view y clearly */
3079566063dSJacob Faibussowitsch       PetscCall(VecGetLocalSize(y,&n));
3088a53a0a4SSajid Ali       if (iscuda) {
3098a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
3109566063dSJacob Faibussowitsch           PetscCall(VecCUDAGetArray(y,&yvalue));
3118a53a0a4SSajid Ali         #endif
3128a53a0a4SSajid Ali       } else {
3139566063dSJacob Faibussowitsch         PetscCall(VecGetArray(y,&yvalue));
3148a53a0a4SSajid Ali       }
31597929ea7SJunchao Zhang       /* Create a vector yg on parentcomm, which shares memory with y. xg and yg must be
31697929ea7SJunchao Zhang         created in the same order in subcomm0/1. For example, we can not reverse the order of
31797929ea7SJunchao Zhang         creating xg and yg in subcomm1.
31897929ea7SJunchao Zhang       */
3198a53a0a4SSajid Ali       if (iscuda) {
3208a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
3219566063dSJacob Faibussowitsch           PetscCall(VecCreateMPICUDAWithArray(parentcomm,1/*bs*/,n,N,yvalue,&yg));
3228a53a0a4SSajid Ali         #endif
3238a53a0a4SSajid Ali       } else {
3249566063dSJacob Faibussowitsch         PetscCall(VecCreateMPIWithArray(parentcomm,1/*bs*/,n,N,yvalue,&yg));
3258a53a0a4SSajid Ali       }
32697929ea7SJunchao Zhang 
32797929ea7SJunchao Zhang       /* Ranks in subcomm0 already specified the full range of the identity map.
32897929ea7SJunchao Zhang         ranks in subcomm1 just need to create empty ISes to cheat VecScatterCreate.
32997929ea7SJunchao Zhang       */
3309566063dSJacob Faibussowitsch       PetscCall(ISCreateStride(PETSC_COMM_SELF,0,0,1,&ix));
3319566063dSJacob Faibussowitsch       PetscCall(ISDuplicate(ix,&iy));
3329566063dSJacob Faibussowitsch       PetscCall(VecScatterCreate(xg,ix,yg,iy,&vscat));
33397929ea7SJunchao Zhang 
33497929ea7SJunchao Zhang       /* Scatter values from xg to yg */
3359566063dSJacob Faibussowitsch       PetscCall(VecScatterBegin(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD));
3369566063dSJacob Faibussowitsch       PetscCall(VecScatterEnd(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD));
33797929ea7SJunchao Zhang 
33897929ea7SJunchao Zhang       /* After the VecScatter is done, values in yg are available. y is our interest, so we return yvalue to y */
3398a53a0a4SSajid Ali       if (iscuda) {
3408a53a0a4SSajid Ali         #if defined(PETSC_HAVE_CUDA)
3419566063dSJacob Faibussowitsch           PetscCall(VecCUDARestoreArray(y,&yvalue));
3428a53a0a4SSajid Ali         #endif
3438a53a0a4SSajid Ali       } else {
3449566063dSJacob Faibussowitsch         PetscCall(VecRestoreArray(y,&yvalue));
3458a53a0a4SSajid Ali       }
34697929ea7SJunchao Zhang 
34797929ea7SJunchao Zhang       /* Libraries on subcomm1 can safely use y now, for example, view it */
3489566063dSJacob Faibussowitsch       PetscCall(VecView(y,PETSC_VIEWER_STDOUT_(subcomm)));
34997929ea7SJunchao Zhang 
3509566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&y));
3519566063dSJacob Faibussowitsch       PetscCall(ISDestroy(&ix));
3529566063dSJacob Faibussowitsch       PetscCall(ISDestroy(&iy));
3539566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&xg));
3549566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&yg));
3559566063dSJacob Faibussowitsch       PetscCall(VecScatterDestroy(&vscat));
3569566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Comm_free(&intercomm));
3579566063dSJacob Faibussowitsch       PetscCallMPI(MPI_Comm_free(&parentcomm));
35897929ea7SJunchao Zhang     } else if (mycolor == 2) { /* subcomm2 */
35997929ea7SJunchao Zhang       /* Processes in subcomm2 do not participate in the VecScatter. They can freely do unrelated things on subcomm2 */
36097929ea7SJunchao Zhang     }
36197929ea7SJunchao Zhang   } /* sub2sub */
36297929ea7SJunchao Zhang 
36397929ea7SJunchao Zhang   /*===========================================================================
36497929ea7SJunchao Zhang    *  Copy a vector x defined on PETSC_COMM_WORLD to vectors y defined on
36597929ea7SJunchao Zhang    *  every subcommunicator of PETSC_COMM_WORLD. We could use multiple transfers
36697929ea7SJunchao Zhang    *  as we did in case 1, but that is not efficient. Instead, we use one vecscatter
36797929ea7SJunchao Zhang    *  to achieve that.
36897929ea7SJunchao Zhang    *===========================================================================*/
36997929ea7SJunchao Zhang   if (world2subs) {
3708a53a0a4SSajid Ali     Vec         y;
37197929ea7SJunchao Zhang     PetscInt    n,N=15,xstart,ystart,low,high;
37297929ea7SJunchao Zhang     PetscScalar *yvalue;
37397929ea7SJunchao Zhang 
37497929ea7SJunchao Zhang     /* Initialize x to [0, 1, 2, 3, ..., N-1] */
3759566063dSJacob Faibussowitsch     PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
3769566063dSJacob Faibussowitsch     PetscCall(VecSetSizes(x, PETSC_DECIDE, N));
3778a53a0a4SSajid Ali     if (iscuda) {
3789566063dSJacob Faibussowitsch       PetscCall(VecSetType(x, VECCUDA));
3798a53a0a4SSajid Ali     } else {
3809566063dSJacob Faibussowitsch       PetscCall(VecSetType(x, VECSTANDARD));
3818a53a0a4SSajid Ali     }
3829566063dSJacob Faibussowitsch     PetscCall(VecSetUp(x));
3839566063dSJacob Faibussowitsch     PetscCall(VecGetOwnershipRange(x,&low,&high));
3849566063dSJacob Faibussowitsch     for (i=low; i<high; i++) PetscCall(VecSetValue(x,i,(PetscScalar)i,INSERT_VALUES));
3859566063dSJacob Faibussowitsch     PetscCall(VecAssemblyBegin(x));
3869566063dSJacob Faibussowitsch     PetscCall(VecAssemblyEnd(x));
38797929ea7SJunchao Zhang 
38897929ea7SJunchao Zhang     /* Every subcomm has a y as long as x */
3899566063dSJacob Faibussowitsch     PetscCall(VecCreate(subcomm, &y));
3909566063dSJacob Faibussowitsch     PetscCall(VecSetSizes(y, PETSC_DECIDE, N));
3918a53a0a4SSajid Ali     if (iscuda) {
3929566063dSJacob Faibussowitsch       PetscCall(VecSetType(y, VECCUDA));
3938a53a0a4SSajid Ali     } else {
3949566063dSJacob Faibussowitsch       PetscCall(VecSetType(y, VECSTANDARD));
3958a53a0a4SSajid Ali     }
3969566063dSJacob Faibussowitsch     PetscCall(VecSetUp(y));
3979566063dSJacob Faibussowitsch     PetscCall(VecGetLocalSize(y,&n));
39897929ea7SJunchao Zhang 
39997929ea7SJunchao Zhang     /* Create a global vector yg on PETSC_COMM_WORLD using y's memory. yg's global size = N*(number of subcommunicators).
40097929ea7SJunchao Zhang        Eeach rank in subcomms contributes a piece to construct the global yg. Keep in mind that pieces from a subcomm are not
40197929ea7SJunchao Zhang        necessarily consecutive in yg. That depends on how PETSC_COMM_WORLD is split. In our case, subcomm0 is made of rank
40297929ea7SJunchao Zhang        0, 3, 6 etc from PETSC_COMM_WORLD. So subcomm0's pieces are interleaved with pieces from other subcomms in yg.
40397929ea7SJunchao Zhang     */
4048a53a0a4SSajid Ali     if (iscuda) {
4058a53a0a4SSajid Ali       #if defined(PETSC_HAVE_CUDA)
4069566063dSJacob Faibussowitsch         PetscCall(VecCUDAGetArray(y,&yvalue));
4079566063dSJacob Faibussowitsch         PetscCall(VecCreateMPICUDAWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,yvalue,&yg));
4088a53a0a4SSajid Ali       #endif
4098a53a0a4SSajid Ali     } else {
4109566063dSJacob Faibussowitsch       PetscCall(VecGetArray(y,&yvalue));
4119566063dSJacob Faibussowitsch       PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,yvalue,&yg));
4128a53a0a4SSajid Ali     }
4139566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetName((PetscObject)yg,"yg_on_subcomms")); /* Give a name to view yg clearly */
41497929ea7SJunchao Zhang 
41597929ea7SJunchao Zhang     /* The following two lines are key. From xstart, we know where to pull entries from x. Note that we get xstart from y,
41697929ea7SJunchao Zhang        since first entry of y on this rank is from x[xstart]. From ystart, we know where ot put entries to yg.
41797929ea7SJunchao Zhang      */
4189566063dSJacob Faibussowitsch     PetscCall(VecGetOwnershipRange(y,&xstart,NULL));
4199566063dSJacob Faibussowitsch     PetscCall(VecGetOwnershipRange(yg,&ystart,NULL));
42097929ea7SJunchao Zhang 
4219566063dSJacob Faibussowitsch     PetscCall(ISCreateStride(PETSC_COMM_SELF,n,xstart,1,&ix));
4229566063dSJacob Faibussowitsch     PetscCall(ISCreateStride(PETSC_COMM_SELF,n,ystart,1,&iy));
4239566063dSJacob Faibussowitsch     PetscCall(VecScatterCreate(x,ix,yg,iy,&vscat));
4249566063dSJacob Faibussowitsch     PetscCall(VecScatterBegin(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
4259566063dSJacob Faibussowitsch     PetscCall(VecScatterEnd(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD));
42697929ea7SJunchao Zhang 
42797929ea7SJunchao Zhang     /* View yg on PETSC_COMM_WORLD before destroying it. We shall see the interleaving effect in output. */
4289566063dSJacob Faibussowitsch     PetscCall(VecView(yg,PETSC_VIEWER_STDOUT_WORLD));
4299566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&yg));
43097929ea7SJunchao Zhang 
43197929ea7SJunchao Zhang     /* Restory yvalue so that processes in subcomm can use y from now on. */
4328a53a0a4SSajid Ali     if (iscuda) {
4338a53a0a4SSajid Ali       #if defined(PETSC_HAVE_CUDA)
4349566063dSJacob Faibussowitsch         PetscCall(VecCUDARestoreArray(y,&yvalue));
4358a53a0a4SSajid Ali       #endif
4368a53a0a4SSajid Ali     } else {
4379566063dSJacob Faibussowitsch       PetscCall(VecRestoreArray(y,&yvalue));
4388a53a0a4SSajid Ali     }
4399566063dSJacob Faibussowitsch     PetscCall(VecScale(y,3.0));
44097929ea7SJunchao Zhang 
4419566063dSJacob Faibussowitsch     PetscCall(ISDestroy(&ix)); /* One can also destroy ix, iy immediately after VecScatterCreate() */
4429566063dSJacob Faibussowitsch     PetscCall(ISDestroy(&iy));
4439566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&x));
4449566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&y));
4459566063dSJacob Faibussowitsch     PetscCall(VecScatterDestroy(&vscat));
44697929ea7SJunchao Zhang   } /* world2subs */
44797929ea7SJunchao Zhang 
4489566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_free(&subcomm));
4499566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
450b122ec5aSJacob Faibussowitsch   return 0;
45197929ea7SJunchao Zhang }
45297929ea7SJunchao Zhang 
45397929ea7SJunchao Zhang /*TEST
45497929ea7SJunchao Zhang 
45597929ea7SJunchao Zhang    build:
456dfd57a17SPierre Jolivet      requires: !defined(PETSC_HAVE_MPIUNI)
457c6c723bbSStefano Zampini 
45897929ea7SJunchao Zhang    testset:
45997929ea7SJunchao Zhang      nsize: 7
46097929ea7SJunchao Zhang 
46197929ea7SJunchao Zhang      test:
46297929ea7SJunchao Zhang        suffix: 1
46397929ea7SJunchao Zhang        args: -world2sub
46497929ea7SJunchao Zhang 
46597929ea7SJunchao Zhang      test:
46697929ea7SJunchao Zhang        suffix: 2
46797929ea7SJunchao Zhang        args: -sub2sub
468f424265bSStefano Zampini        # deadlocks with NECMPI and INTELMPI (20210400300)
469f424265bSStefano Zampini        requires: !defined(PETSC_HAVE_NECMPI) !defined(PETSC_HAVE_I_MPI_NUMVERSION)
47097929ea7SJunchao Zhang 
47197929ea7SJunchao Zhang      test:
47297929ea7SJunchao Zhang        suffix: 3
47397929ea7SJunchao Zhang        args: -world2subs
47497929ea7SJunchao Zhang 
47597929ea7SJunchao Zhang      test:
47697929ea7SJunchao Zhang        suffix: 4
4778a53a0a4SSajid Ali        args: -world2sub -vectype cuda
4788a53a0a4SSajid Ali        requires: cuda
4798a53a0a4SSajid Ali 
4808a53a0a4SSajid Ali      test:
4818a53a0a4SSajid Ali        suffix: 5
4828a53a0a4SSajid Ali        args: -sub2sub -vectype cuda
4838a53a0a4SSajid Ali        requires: cuda
4848a53a0a4SSajid Ali 
4858a53a0a4SSajid Ali      test:
4868a53a0a4SSajid Ali       suffix: 6
4878a53a0a4SSajid Ali       args: -world2subs -vectype cuda
4888a53a0a4SSajid Ali       requires: cuda
4898a53a0a4SSajid Ali 
4908a53a0a4SSajid Ali      test:
4918a53a0a4SSajid Ali        suffix: 7
49297929ea7SJunchao Zhang        args: -world2sub -sf_type neighbor
49397929ea7SJunchao Zhang        output_file: output/ex9_1.out
494c6c723bbSStefano Zampini        # 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)
495c6c723bbSStefano Zampini        # segfaults with NECMPI
496c6c723bbSStefano Zampini        requires: defined(PETSC_HAVE_MPI_NEIGHBORHOOD_COLLECTIVES) !defined(PETSC_HAVE_OMPI_MAJOR_VERSION) !defined(PETSC_HAVE_NECMPI)
49797929ea7SJunchao Zhang TEST*/
498