1 static char help[] = "Scatters from a sequential vector to a parallel vector. \n\ 2 uses block index sets\n\n"; 3 4 #include <petscvec.h> 5 6 int main(int argc, char **argv) 7 { 8 PetscInt bs = 1, n = 5, i, low; 9 PetscInt ix0[3] = {5, 7, 9}, iy0[3] = {1, 2, 4}, ix1[3] = {2, 3, 4}, iy1[3] = {0, 1, 3}; 10 PetscMPIInt size, rank; 11 PetscScalar *array; 12 Vec x, y; 13 IS isx, isy; 14 VecScatter ctx; 15 PetscViewer sviewer; 16 17 PetscFunctionBeginUser; 18 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 19 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 20 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 21 22 PetscCheck(size >= 2, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Must run more than one processor"); 23 24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL)); 25 n = bs * n; 26 27 /* Create vector x over shared memory */ 28 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 29 PetscCall(VecSetSizes(x, n, PETSC_DECIDE)); 30 PetscCall(VecSetFromOptions(x)); 31 32 PetscCall(VecGetOwnershipRange(x, &low, NULL)); 33 PetscCall(VecGetArray(x, &array)); 34 for (i = 0; i < n; i++) array[i] = (PetscScalar)(i + low); 35 PetscCall(VecRestoreArray(x, &array)); 36 37 /* Create a sequential vector y */ 38 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &y)); 39 PetscCall(VecGetArray(y, &array)); 40 for (i = 0; i < n; i++) array[i] = (PetscScalar)(i + 100 * rank); 41 PetscCall(VecRestoreArray(y, &array)); 42 43 /* Create two index sets */ 44 if (rank == 0) { 45 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, 3, ix0, PETSC_COPY_VALUES, &isx)); 46 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, 3, iy0, PETSC_COPY_VALUES, &isy)); 47 } else { 48 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, 3, ix1, PETSC_COPY_VALUES, &isx)); 49 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, 3, iy1, PETSC_COPY_VALUES, &isy)); 50 } 51 52 if (rank == 10) { 53 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n[%d] isx:\n", rank)); 54 PetscCall(ISView(isx, PETSC_VIEWER_STDOUT_SELF)); 55 } 56 57 PetscCall(VecScatterCreate(y, isy, x, isx, &ctx)); 58 PetscCall(VecScatterSetFromOptions(ctx)); 59 60 /* Test forward vecscatter */ 61 PetscCall(VecSet(x, 0.0)); 62 PetscCall(VecScatterBegin(ctx, y, x, ADD_VALUES, SCATTER_FORWARD)); 63 PetscCall(VecScatterEnd(ctx, y, x, ADD_VALUES, SCATTER_FORWARD)); 64 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 65 66 /* Test reverse vecscatter */ 67 PetscCall(VecScale(x, -1.0)); 68 PetscCall(VecScatterBegin(ctx, x, y, ADD_VALUES, SCATTER_REVERSE)); 69 PetscCall(VecScatterEnd(ctx, x, y, ADD_VALUES, SCATTER_REVERSE)); 70 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer)); 71 if (rank == 1) PetscCall(VecView(y, sviewer)); 72 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer)); 73 74 /* Free spaces */ 75 PetscCall(VecScatterDestroy(&ctx)); 76 PetscCall(ISDestroy(&isx)); 77 PetscCall(ISDestroy(&isy)); 78 PetscCall(VecDestroy(&x)); 79 PetscCall(VecDestroy(&y)); 80 PetscCall(PetscFinalize()); 81 return 0; 82 } 83 84 /*TEST 85 86 test: 87 nsize: 3 88 89 TEST*/ 90