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 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&nproc);CHKERRMPI(ierr); 26 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&grank);CHKERRMPI(ierr); 27 28 if (nproc < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ,"This test must have at least two processes to run"); 29 30 ierr = PetscOptionsGetBool(NULL,0,"-world2sub",&world2sub,NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsGetBool(NULL,0,"-sub2sub",&sub2sub,NULL);CHKERRQ(ierr); 32 ierr = PetscOptionsGetBool(NULL,0,"-world2subs",&world2subs,NULL);CHKERRQ(ierr); 33 ierr = PetscOptionsGetString(NULL,NULL,"-vectype",vectypename,sizeof(vectypename),&optionflag);CHKERRQ(ierr); 34 if (optionflag) { 35 ierr = PetscStrncmp(vectypename, "cuda", (size_t)4, &compareflag);CHKERRQ(ierr); 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 ierr = MPI_Comm_split(PETSC_COMM_WORLD,mycolor,grank,&subcomm);CHKERRMPI(ierr); 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 ierr = VecCreate(PETSC_COMM_WORLD, &x);CHKERRQ(ierr); 49 ierr = VecSetSizes(x, PETSC_DECIDE, N);CHKERRQ(ierr); 50 if (iscuda) { 51 ierr = VecSetType(x, VECCUDA);CHKERRQ(ierr); 52 } else { 53 ierr = VecSetType(x, VECSTANDARD);CHKERRQ(ierr); 54 } 55 ierr = VecSetUp(x);CHKERRQ(ierr); 56 ierr = PetscObjectSetName((PetscObject)x,"x_commworld");CHKERRQ(ierr); /* Give a name to view x clearly */ 57 58 /* Initialize x to [-0.0, -1.0, -2.0, ..., -19.0] */ 59 ierr = VecGetOwnershipRange(x,&low,&high);CHKERRQ(ierr); 60 for (i=low; i<high; i++) { 61 PetscScalar val = -i; 62 ierr = VecSetValue(x,i,val,INSERT_VALUES);CHKERRQ(ierr); 63 } 64 ierr = VecAssemblyBegin(x);CHKERRQ(ierr); 65 ierr = VecAssemblyEnd(x);CHKERRQ(ierr); 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 ierr = VecCreate(subcomm, &y);CHKERRQ(ierr); 72 ierr = VecSetSizes(y, PETSC_DECIDE, N);CHKERRQ(ierr); 73 if (iscuda) { 74 ierr = VecSetType(y, VECCUDA);CHKERRQ(ierr); 75 } else { 76 ierr = VecSetType(y, VECSTANDARD);CHKERRQ(ierr); 77 } 78 ierr = VecSetUp(y);CHKERRQ(ierr); 79 ierr = PetscObjectSetName((PetscObject)y,"y_subcomm_0");CHKERRQ(ierr); /* Give a name to view y clearly */ 80 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 81 if (iscuda) { 82 #if defined(PETSC_HAVE_CUDA) 83 ierr = VecCUDAGetArray(y,&yvalue);CHKERRQ(ierr); 84 #endif 85 } else { 86 ierr = VecGetArray(y,&yvalue);CHKERRQ(ierr); 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 ierr = VecCreateMPICUDAWithArray(PETSC_COMM_WORLD,1,n,N,yvalue,&yg);CHKERRQ(ierr); 94 #endif 95 } else { 96 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,yvalue,&yg);CHKERRQ(ierr); 97 } 98 99 /* Create an identity map that makes yg[i] = x[i], i=0..N-1 */ 100 ierr = VecGetOwnershipRange(yg,&low,&high);CHKERRQ(ierr); /* low, high are global indices */ 101 ierr = ISCreateStride(PETSC_COMM_SELF,high-low,low,1,&ix);CHKERRQ(ierr); 102 ierr = ISDuplicate(ix,&iy);CHKERRQ(ierr); 103 104 /* Union of ix's on subcomm0 covers the full range of [0,N) */ 105 ierr = VecScatterCreate(x,ix,yg,iy,&vscat);CHKERRQ(ierr); 106 ierr = VecScatterBegin(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 107 ierr = VecScatterEnd(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 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 ierr = VecCUDARestoreArray(y,&yvalue);CHKERRQ(ierr); 115 #endif 116 } else { 117 ierr = VecRestoreArray(y,&yvalue);CHKERRQ(ierr); 118 } 119 120 /* Libraries on subcomm0 can safely use y now, for example, view and scale it */ 121 ierr = VecView(y,PETSC_VIEWER_STDOUT_(subcomm));CHKERRQ(ierr); 122 ierr = VecScale(y,2.0);CHKERRQ(ierr); 123 124 /* Send the new y back to x */ 125 ierr = VecGetArray(y,&yvalue);CHKERRQ(ierr); /* If VecScale is done on GPU, Petsc will prepare a valid yvalue for access */ 126 /* Supply new yvalue to yg without memory copying */ 127 ierr = VecPlaceArray(yg,yvalue);CHKERRQ(ierr); 128 ierr = VecScatterBegin(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 129 ierr = VecScatterEnd(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 130 ierr = VecResetArray(yg);CHKERRQ(ierr); 131 if (iscuda) { 132 #if defined(PETSC_HAVE_CUDA) 133 ierr = VecCUDARestoreArray(y,&yvalue);CHKERRQ(ierr); 134 #endif 135 } else { 136 ierr = VecRestoreArray(y,&yvalue);CHKERRQ(ierr); 137 } 138 139 ierr = VecDestroy(&y);CHKERRQ(ierr); 140 } else { 141 /* Ranks outside of subcomm0 do not supply values to yg */ 142 if (iscuda) { 143 #if defined(PETSC_HAVE_CUDA) 144 ierr = VecCreateMPICUDAWithArray(PETSC_COMM_WORLD,1,0/*n*/,N,NULL,&yg);CHKERRQ(ierr); 145 #endif 146 } else { 147 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,0/*n*/,N,NULL,&yg);CHKERRQ(ierr); 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 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&ix);CHKERRQ(ierr); 154 ierr = ISDuplicate(ix,&iy);CHKERRQ(ierr); 155 156 ierr = VecScatterCreate(x,ix,yg,iy,&vscat);CHKERRQ(ierr); 157 ierr = VecScatterBegin(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 158 ierr = VecScatterEnd(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 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 ierr = VecScatterBegin(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 164 ierr = VecScatterEnd(vscat,yg,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 165 } 166 167 ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 168 ierr = ISDestroy(&ix);CHKERRQ(ierr); 169 ierr = ISDestroy(&iy);CHKERRQ(ierr); 170 ierr = VecDestroy(&x);CHKERRQ(ierr); 171 ierr = VecDestroy(&yg);CHKERRQ(ierr); 172 ierr = VecScatterDestroy(&vscat);CHKERRQ(ierr); 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 ierr = MPI_Comm_rank(subcomm,&lrank);CHKERRMPI(ierr); 192 /* x is on subcomm */ 193 ierr = VecCreate(subcomm, &x);CHKERRQ(ierr); 194 ierr = VecSetSizes(x, PETSC_DECIDE, N);CHKERRQ(ierr); 195 if (iscuda) { 196 ierr = VecSetType(x, VECCUDA);CHKERRQ(ierr); 197 } else { 198 ierr = VecSetType(x, VECSTANDARD);CHKERRQ(ierr); 199 } 200 ierr = VecSetUp(x);CHKERRQ(ierr); 201 ierr = VecGetOwnershipRange(x,&low,&high);CHKERRQ(ierr); 202 203 /* initialize x = [0.0, 1.0, 2.0, ..., 21.0] */ 204 for (i=low; i<high; i++) { 205 PetscScalar val = i; 206 ierr = VecSetValue(x,i,val,INSERT_VALUES);CHKERRQ(ierr); 207 } 208 ierr = VecAssemblyBegin(x);CHKERRQ(ierr); 209 ierr = VecAssemblyEnd(x);CHKERRQ(ierr); 210 211 ierr = MPI_Intercomm_create(subcomm,0,PETSC_COMM_WORLD/*peer_comm*/,1,100/*tag*/,&intercomm);CHKERRMPI(ierr); 212 213 /* Tell rank 0 of subcomm1 the global size of x */ 214 if (!lrank) {ierr = MPI_Send(&N,1,MPIU_INT,0/*receiver's rank in remote comm, i.e., subcomm1*/,200/*tag*/,intercomm);CHKERRMPI(ierr);} 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 ierr = MPI_Intercomm_merge(intercomm,0/*low*/,&parentcomm);CHKERRMPI(ierr); 220 221 /* Create a vector xg on parentcomm, which shares memory with x */ 222 ierr = VecGetLocalSize(x,&n);CHKERRQ(ierr); 223 if (iscuda) { 224 #if defined(PETSC_HAVE_CUDA) 225 ierr = VecCUDAGetArrayRead(x,&xvalue);CHKERRQ(ierr); 226 ierr = VecCreateMPICUDAWithArray(parentcomm,1,n,N,xvalue,&xg);CHKERRQ(ierr); 227 #endif 228 } else { 229 ierr = VecGetArrayRead(x,&xvalue);CHKERRQ(ierr); 230 ierr = VecCreateMPIWithArray(parentcomm,1,n,N,xvalue,&xg);CHKERRQ(ierr); 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 ierr = VecCreateMPICUDAWithArray(parentcomm,1,0/*n*/,N,NULL/*array*/,&yg);CHKERRQ(ierr); 237 #endif 238 } else { 239 ierr = VecCreateMPIWithArray(parentcomm,1,0/*n*/,N,NULL/*array*/,&yg);CHKERRQ(ierr); 240 } 241 242 /* Create the vecscatter, which does identity map by setting yg[i] = xg[i], i=0..N-1. */ 243 ierr = VecGetOwnershipRange(xg,&low,&high);CHKERRQ(ierr); /* low, high are global indices of xg */ 244 ierr = ISCreateStride(PETSC_COMM_SELF,high-low,low,1,&ix);CHKERRQ(ierr); 245 ierr = ISDuplicate(ix,&iy);CHKERRQ(ierr); 246 ierr = VecScatterCreate(xg,ix,yg,iy,&vscat);CHKERRQ(ierr); 247 248 /* Scatter values from xg to yg */ 249 ierr = VecScatterBegin(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 250 ierr = VecScatterEnd(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 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 ierr = VecCUDARestoreArrayRead(x,&xvalue);CHKERRQ(ierr); 256 #endif 257 } else { 258 ierr = VecRestoreArrayRead(x,&xvalue);CHKERRQ(ierr); 259 } 260 ierr = VecDestroy(&x);CHKERRQ(ierr); 261 ierr = ISDestroy(&ix);CHKERRQ(ierr); 262 ierr = ISDestroy(&iy);CHKERRQ(ierr); 263 ierr = VecDestroy(&xg);CHKERRQ(ierr); 264 ierr = VecDestroy(&yg);CHKERRQ(ierr); 265 ierr = VecScatterDestroy(&vscat);CHKERRQ(ierr); 266 ierr = MPI_Comm_free(&intercomm);CHKERRMPI(ierr); 267 ierr = MPI_Comm_free(&parentcomm);CHKERRMPI(ierr); 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 ierr = MPI_Comm_rank(subcomm,&lrank);CHKERRMPI(ierr); 278 ierr = MPI_Intercomm_create(subcomm,0,PETSC_COMM_WORLD/*peer_comm*/,0/*remote_leader*/,100/*tag*/,&intercomm);CHKERRMPI(ierr); 279 280 /* Two rank-0 are talking */ 281 if (!lrank) {ierr = MPI_Recv(&N,1,MPIU_INT,0/*sender's rank in remote comm, i.e. subcomm0*/,200/*tag*/,intercomm,MPI_STATUS_IGNORE);CHKERRMPI(ierr);} 282 /* Rank 0 of subcomm1 bcasts N to its members */ 283 ierr = MPI_Bcast(&N,1,MPIU_INT,0/*local root*/,subcomm);CHKERRMPI(ierr); 284 285 /* Create a intracomm Petsc can work on */ 286 ierr = MPI_Intercomm_merge(intercomm,1/*high*/,&parentcomm);CHKERRMPI(ierr); 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 ierr = VecCreateMPICUDAWithArray(parentcomm,1/*bs*/,0/*n*/,N,NULL/*array*/,&xg);CHKERRQ(ierr); 292 #endif 293 } else { 294 ierr = VecCreateMPIWithArray(parentcomm,1/*bs*/,0/*n*/,N,NULL/*array*/,&xg);CHKERRQ(ierr); 295 } 296 297 ierr = VecCreate(subcomm, &y);CHKERRQ(ierr); 298 ierr = VecSetSizes(y, PETSC_DECIDE, N);CHKERRQ(ierr); 299 if (iscuda) { 300 ierr = VecSetType(y, VECCUDA);CHKERRQ(ierr); 301 } else { 302 ierr = VecSetType(y, VECSTANDARD);CHKERRQ(ierr); 303 } 304 ierr = VecSetUp(y);CHKERRQ(ierr); 305 306 ierr = PetscObjectSetName((PetscObject)y,"y_subcomm_1");CHKERRQ(ierr); /* Give a name to view y clearly */ 307 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 308 if (iscuda) { 309 #if defined(PETSC_HAVE_CUDA) 310 ierr = VecCUDAGetArray(y,&yvalue);CHKERRQ(ierr); 311 #endif 312 } else { 313 ierr = VecGetArray(y,&yvalue);CHKERRQ(ierr); 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 ierr = VecCreateMPICUDAWithArray(parentcomm,1/*bs*/,n,N,yvalue,&yg);CHKERRQ(ierr); 322 #endif 323 } else { 324 ierr = VecCreateMPIWithArray(parentcomm,1/*bs*/,n,N,yvalue,&yg);CHKERRQ(ierr); 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 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&ix);CHKERRQ(ierr); 331 ierr = ISDuplicate(ix,&iy);CHKERRQ(ierr); 332 ierr = VecScatterCreate(xg,ix,yg,iy,&vscat);CHKERRQ(ierr); 333 334 /* Scatter values from xg to yg */ 335 ierr = VecScatterBegin(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 336 ierr = VecScatterEnd(vscat,xg,yg,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 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 ierr = VecCUDARestoreArray(y,&yvalue);CHKERRQ(ierr); 342 #endif 343 } else { 344 ierr = VecRestoreArray(y,&yvalue);CHKERRQ(ierr); 345 } 346 347 /* Libraries on subcomm1 can safely use y now, for example, view it */ 348 ierr = VecView(y,PETSC_VIEWER_STDOUT_(subcomm));CHKERRQ(ierr); 349 350 ierr = VecDestroy(&y);CHKERRQ(ierr); 351 ierr = ISDestroy(&ix);CHKERRQ(ierr); 352 ierr = ISDestroy(&iy);CHKERRQ(ierr); 353 ierr = VecDestroy(&xg);CHKERRQ(ierr); 354 ierr = VecDestroy(&yg);CHKERRQ(ierr); 355 ierr = VecScatterDestroy(&vscat);CHKERRQ(ierr); 356 ierr = MPI_Comm_free(&intercomm);CHKERRMPI(ierr); 357 ierr = MPI_Comm_free(&parentcomm);CHKERRMPI(ierr); 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 ierr = VecCreate(PETSC_COMM_WORLD, &x);CHKERRQ(ierr); 376 ierr = VecSetSizes(x, PETSC_DECIDE, N);CHKERRQ(ierr); 377 if (iscuda) { 378 ierr = VecSetType(x, VECCUDA);CHKERRQ(ierr); 379 } else { 380 ierr = VecSetType(x, VECSTANDARD);CHKERRQ(ierr); 381 } 382 ierr = VecSetUp(x);CHKERRQ(ierr); 383 ierr = VecGetOwnershipRange(x,&low,&high);CHKERRQ(ierr); 384 for (i=low; i<high; i++) {ierr = VecSetValue(x,i,(PetscScalar)i,INSERT_VALUES);CHKERRQ(ierr);} 385 ierr = VecAssemblyBegin(x);CHKERRQ(ierr); 386 ierr = VecAssemblyEnd(x);CHKERRQ(ierr); 387 388 /* Every subcomm has a y as long as x */ 389 ierr = VecCreate(subcomm, &y);CHKERRQ(ierr); 390 ierr = VecSetSizes(y, PETSC_DECIDE, N);CHKERRQ(ierr); 391 if (iscuda) { 392 ierr = VecSetType(y, VECCUDA);CHKERRQ(ierr); 393 } else { 394 ierr = VecSetType(y, VECSTANDARD);CHKERRQ(ierr); 395 } 396 ierr = VecSetUp(y);CHKERRQ(ierr); 397 ierr = VecGetLocalSize(y,&n);CHKERRQ(ierr); 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 ierr = VecCUDAGetArray(y,&yvalue);CHKERRQ(ierr); 407 ierr = VecCreateMPICUDAWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,yvalue,&yg);CHKERRQ(ierr); 408 #endif 409 } else { 410 ierr = VecGetArray(y,&yvalue);CHKERRQ(ierr); 411 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,yvalue,&yg);CHKERRQ(ierr); 412 } 413 ierr = PetscObjectSetName((PetscObject)yg,"yg_on_subcomms");CHKERRQ(ierr); /* 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 ierr = VecGetOwnershipRange(y,&xstart,NULL);CHKERRQ(ierr); 419 ierr = VecGetOwnershipRange(yg,&ystart,NULL);CHKERRQ(ierr); 420 421 ierr = ISCreateStride(PETSC_COMM_SELF,n,xstart,1,&ix);CHKERRQ(ierr); 422 ierr = ISCreateStride(PETSC_COMM_SELF,n,ystart,1,&iy);CHKERRQ(ierr); 423 ierr = VecScatterCreate(x,ix,yg,iy,&vscat);CHKERRQ(ierr); 424 ierr = VecScatterBegin(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 425 ierr = VecScatterEnd(vscat,x,yg,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 426 427 /* View yg on PETSC_COMM_WORLD before destroying it. We shall see the interleaving effect in output. */ 428 ierr = VecView(yg,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 429 ierr = VecDestroy(&yg);CHKERRQ(ierr); 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 ierr = VecCUDARestoreArray(y,&yvalue);CHKERRQ(ierr); 435 #endif 436 } else { 437 ierr = VecRestoreArray(y,&yvalue);CHKERRQ(ierr); 438 } 439 ierr = VecScale(y,3.0);CHKERRQ(ierr); 440 441 ierr = ISDestroy(&ix);CHKERRQ(ierr); /* One can also destroy ix, iy immediately after VecScatterCreate() */ 442 ierr = ISDestroy(&iy);CHKERRQ(ierr); 443 ierr = VecDestroy(&x);CHKERRQ(ierr); 444 ierr = VecDestroy(&y);CHKERRQ(ierr); 445 ierr = VecScatterDestroy(&vscat);CHKERRQ(ierr); 446 } /* world2subs */ 447 448 ierr = MPI_Comm_free(&subcomm);CHKERRMPI(ierr); 449 ierr = PetscFinalize(); 450 return ierr; 451 } 452 453 /*TEST 454 455 build: 456 requires: !define(PETSC_HAVE_MPIUNI) 457 testset: 458 nsize: 7 459 460 test: 461 suffix: 1 462 args: -world2sub 463 464 test: 465 suffix: 2 466 args: -sub2sub 467 468 test: 469 suffix: 3 470 args: -world2subs 471 472 test: 473 suffix: 4 474 args: -world2sub -vectype cuda 475 requires: cuda 476 477 test: 478 suffix: 5 479 args: -sub2sub -vectype cuda 480 requires: cuda 481 482 test: 483 suffix: 6 484 args: -world2subs -vectype cuda 485 requires: cuda 486 487 test: 488 suffix: 7 489 args: -world2sub -sf_type neighbor 490 output_file: output/ex9_1.out 491 # 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) 492 requires: define(PETSC_HAVE_MPI_NEIGHBORHOOD_COLLECTIVES) !define(PETSC_HAVE_OMPI_MAJOR_VERSION) 493 TEST*/ 494 495