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