1 static char help[]= " Test VecScatter with x, y on different communicators\n\n"; 2 3 #include <petscvec.h> 4 5 int main(int argc,char **argv) 6 { 7 PetscInt i,n=5,rstart; 8 PetscScalar *val; 9 const PetscScalar *dat; 10 Vec x,y1,y2; 11 MPI_Comm newcomm; 12 PetscMPIInt nproc,rank; 13 IS ix; 14 VecScatter vscat1,vscat2; 15 16 PetscFunctionBegin; 17 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 18 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&nproc)); 19 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 20 21 PetscCheckFalse(nproc != 2,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This test must run with exactly two MPI ranks"); 22 23 /* Create MPI vectors x and y, which are on the same comm (i.e., MPI_IDENT) */ 24 CHKERRQ(VecCreateMPI(PETSC_COMM_WORLD,n,PETSC_DECIDE,&x)); 25 CHKERRQ(VecDuplicate(x,&y1)); 26 CHKERRQ(VecGetOwnershipRange(x,&rstart,NULL)); 27 28 /* Set x's value locally. x would be {0., 1., 2., ..., 9.} */ 29 CHKERRQ(VecGetArray(x,&val)); 30 for (i=0; i<n; i++) val[i] = rstart + i; 31 CHKERRQ(VecRestoreArray(x,&val)); 32 33 /* Create index set ix = {0, 1, 2, ..., 9} */ 34 CHKERRQ(ISCreateStride(PETSC_COMM_WORLD,n,rstart,1,&ix)); 35 36 /* Create newcomm that reverses processes in x's comm, and then create y2 on it*/ 37 CHKERRMPI(MPI_Comm_split(PETSC_COMM_WORLD,0/*color*/,-rank/*key*/,&newcomm)); 38 CHKERRQ(VecCreateMPI(newcomm,n,PETSC_DECIDE,&y2)); 39 40 /* It looks vscat1/2 are the same, but actually not. y2 is on a different communicator than x */ 41 CHKERRQ(VecScatterCreate(x,ix,y1,ix,&vscat1)); 42 CHKERRQ(VecScatterCreate(x,ix,y2,ix,&vscat2)); 43 44 CHKERRQ(VecScatterBegin(vscat1,x,y1,INSERT_VALUES,SCATTER_FORWARD)); 45 CHKERRQ(VecScatterBegin(vscat2,x,y2,INSERT_VALUES,SCATTER_FORWARD)); 46 CHKERRQ(VecScatterEnd (vscat1,x,y1,INSERT_VALUES,SCATTER_FORWARD)); 47 CHKERRQ(VecScatterEnd (vscat2,x,y2,INSERT_VALUES,SCATTER_FORWARD)); 48 49 /* View on rank 0 of x's comm, which is PETSC_COMM_WORLD */ 50 if (rank == 0) { 51 /* Print the part of x on rank 0, which is 0 1 2 3 4 */ 52 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"\nOn rank 0 of PETSC_COMM_WORLD, x = ")); 53 CHKERRQ(VecGetArrayRead(x,&dat)); 54 for (i=0; i<n; i++) CHKERRQ(PetscPrintf(PETSC_COMM_SELF," %.0g",(double)PetscRealPart(dat[i]))); 55 CHKERRQ(VecRestoreArrayRead(x,&dat)); 56 57 /* Print the part of y1 on rank 0, which is 0 1 2 3 4 */ 58 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"\nOn rank 0 of PETSC_COMM_WORLD, y1 = ")); 59 CHKERRQ(VecGetArrayRead(y1,&dat)); 60 for (i=0; i<n; i++) CHKERRQ(PetscPrintf(PETSC_COMM_SELF," %.0g",(double)PetscRealPart(dat[i]))); 61 CHKERRQ(VecRestoreArrayRead(y1,&dat)); 62 63 /* Print the part of y2 on rank 0, which is 5 6 7 8 9 since y2 swapped the processes of PETSC_COMM_WORLD */ 64 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"\nOn rank 0 of PETSC_COMM_WORLD, y2 = ")); 65 CHKERRQ(VecGetArrayRead(y2,&dat)); 66 for (i=0; i<n; i++) CHKERRQ(PetscPrintf(PETSC_COMM_SELF," %.0g",(double)PetscRealPart(dat[i]))); 67 CHKERRQ(VecRestoreArrayRead(y2,&dat)); 68 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"\n")); 69 } 70 71 CHKERRQ(ISDestroy(&ix)); 72 CHKERRQ(VecDestroy(&x)); 73 CHKERRQ(VecDestroy(&y1)); 74 CHKERRQ(VecDestroy(&y2)); 75 CHKERRQ(VecScatterDestroy(&vscat1)); 76 CHKERRQ(VecScatterDestroy(&vscat2)); 77 CHKERRMPI(MPI_Comm_free(&newcomm)); 78 CHKERRQ(PetscFinalize()); 79 return 0; 80 } 81 82 /*TEST 83 84 test: 85 nsize: 2 86 requires: double 87 TEST*/ 88