// 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();
}