[petsc-users] How to VecScatter from global to local vector, and then, VecGather back ?

[Franck Houssen]

How to VecScatter from global to local vector, and then, VecGather back ? 

This is a very simple use case: I need to split a global vector in local 
(possibly overlapping) pieces, then I need to modify each local piece (x2), and 
finally I need to assemble (+=) back local parts into a global vector. 
Read the doc and went through examples... But still can't make this work: can I 
get some help on this ? 

Note: running petsc-3.7.6 on debian with gcc-6.3 

Thanks, 

Franck 

~> head -n 12 vecScatterGather.cpp 
// How to VecScatter from global to local vector, and then, VecGather back ? 
// 
// global vector: 3x1 2 overlapping local vector: 2x1 global vector: 3x1 
// 
// x2 
// |1 -> |2 
// |1 scatter |1 |2 gather |2 
// |1 -> -> |4 
// |1 |1 -> |2 |2 
// |1 |2 
// 
// ~> g++ -o vecScatterGather.exe vecScatterGather.cpp -lpetsc -lm; mpirun -n 2 
vecScatterGather.exe 



// How to VecScatter from global to local vector, and then, VecGather back ?
//
//  global vector: 3x1           2 overlapping local vector: 2x1                    global vector: 3x1
//
//                                            x2
//                                      |1    ->    |2
//           |1        scatter          |1          |2                gather               |2
//           |1          ->                                             ->                 |4
//           |1                         |1    ->    |2                                     |2
//                                      |1          |2
//
// ~> g++ -o vecScatterGather.exe vecScatterGather.cpp -lpetsc -lm; mpirun -n 2 vecScatterGather.exe

#include "petsc.h"

int main(int argc,char **argv) {
  PetscInitialize(&argc, &argv, NULL, NULL);
  int size = 0; MPI_Comm_size(MPI_COMM_WORLD, &size); if (size != 2) return 1;
  int rank = 0; MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  PetscInt globSize = 3;
  Vec globVec; VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, globSize, &globVec); VecSet(globVec, 1.);
  VecView(globVec, PETSC_VIEWER_STDOUT_WORLD); PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);

  PetscInt locSize = 2;
  Vec locVec; VecCreateSeq(PETSC_COMM_SELF, locSize, &locVec); VecSet(globVec, -1.);

  PetscInt locIdx[2] = {0, 0};
  if (rank == 0) {locIdx[0] = 0; locIdx[1] = 1;} // 1st (local) domain is {0, 1} (global index)
  else           {locIdx[0] = 1; locIdx[1] = 2;} // 2nd (local) domain is {1, 2} (global index)
  IS is; ISCreateGeneral(PETSC_COMM_WORLD, locSize, locIdx, PETSC_COPY_VALUES, &is);
  ISView(is, PETSC_VIEWER_STDOUT_WORLD);

  VecScatter scatCtx; VecScatterCreate(globVec, is, locVec, NULL, &scatCtx);
  VecScatterBegin(scatCtx, globVec, locVec, INSERT_VALUES, SCATTER_FORWARD);
  VecScatterEnd  (scatCtx, globVec, locVec, INSERT_VALUES, SCATTER_FORWARD);

  if (rank == 0) {VecView(locVec, PETSC_VIEWER_STDOUT_SELF); PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);}
  MPI_Barrier(MPI_COMM_WORLD);
  if (rank == 1) {VecView(locVec, PETSC_VIEWER_STDOUT_SELF); PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);}
  MPI_Barrier(MPI_COMM_WORLD);

  VecScale(locVec, 2.);

  VecScatterBegin(scatCtx, locVec, globVec, ADD_VALUES, SCATTER_REVERSE);
  VecScatterEnd  (scatCtx, locVec, globVec, ADD_VALUES, SCATTER_REVERSE);
  VecView(globVec, PETSC_VIEWER_STDOUT_WORLD); PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);

  PetscFinalize();
}

Vec Object: 2 MPI processes
  type: mpi
Process [0]
1.
1.
Process [1]
1.
IS Object: 2 MPI processes
  type: general
[0] Number of indices in set 2
[0] 0 0
[0] 1 1
[1] Number of indices in set 2
[1] 0 1
[1] 1 2
Vec Object: 1 MPI processes
  type: seq
-1.
-1.
Vec Object: 1 MPI processes
  type: seq
-1.
-1.
Vec Object: 2 MPI processes
  type: mpi
Process [0]
-3.
-5.
Process [1]
-3.

vecScatterGather.log.expected
Description: Binary data