#include <stdio.h>
#include <petscvec.h>

int main(int argc, char **argv){
  PetscErrorCode ierr;
  PetscInt n = 5, loc_tgt_size;
  PetscMPIInt rank, size;
  PetscReal vnorm;
  Vec src_vec, tgt_vec;
  PetscRandom rctx;
  IS is;
  VecScatter sctx;

  ierr = PetscInitialize(&argc, &argv, NULL, NULL); CHKERRQ(ierr);

  MPI_Comm_size(PETSC_COMM_WORLD, &size);
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

  // Create parallel source vector with random entries
  ierr = VecCreate(PETSC_COMM_WORLD, &src_vec); CHKERRQ(ierr);
  ierr = VecSetType(src_vec, VECSTANDARD); CHKERRQ(ierr);
  ierr = VecSetSizes(src_vec, PETSC_DECIDE, n); CHKERRQ(ierr);
  ierr = PetscRandomCreate(PETSC_COMM_WORLD, &rctx); CHKERRQ(ierr);
  ierr = VecSetRandom(src_vec, rctx); CHKERRQ(ierr);
  ierr = PetscRandomDestroy(&rctx); CHKERRQ(ierr);

  // Create sequential target vector w. zero length on all but one process, zero entries.
  loc_tgt_size = 0;
  if (rank == (size-1)){
    loc_tgt_size = n;
  }
  ierr = VecCreateSeq(PETSC_COMM_SELF, loc_tgt_size, &tgt_vec); CHKERRQ(ierr);
  ierr = VecZeroEntries(tgt_vec); CHKERRQ(ierr);

  // Scatter source vector to target vector on one process
  ierr = ISCreateStride(PETSC_COMM_SELF, loc_tgt_size, 0, 1, &is); CHKERRQ(ierr);
  ierr = VecScatterCreate(src_vec, is, tgt_vec, is, &sctx); CHKERRQ(ierr);
  ierr = VecScatterBegin(sctx, src_vec, tgt_vec, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
  ierr = VecScatterEnd(sctx, src_vec, tgt_vec, INSERT_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);

  // Get and print norm of target vector on each process
  ierr = VecNorm(tgt_vec, NORM_2, &vnorm); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_SELF, "[%i] |v| = %e\n", rank, vnorm); CHKERRQ(ierr);

  ierr = VecDestroy(&src_vec); CHKERRQ(ierr);
  ierr = VecDestroy(&tgt_vec); CHKERRQ(ierr);
  ierr = ISDestroy(&is); CHKERRQ(ierr);
  ierr = VecScatterDestroy(&sctx); CHKERRQ(ierr);

  ierr = PetscFinalize(); CHKERRQ(ierr);

  return(ierr);
}
