
static char help[] = "Attempt at mpi-vec scatter, with ordering reversed*\n\n";

#include <petscis.h>
#include <petscvec.h>
#include <petscviewer.h>

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  PetscMPIInt    rank;
  Vec            vec1,vec2;
  PetscInt       N = 10,i;
  PetscInt       start1,end1,start2,end2;
  PetscScalar    v;


  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  /*-----------------------------------------------
    Create and get ownership ranges of vectors
    -----------------------------------------------*/ 
  ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,N,&vec1);
  ierr = VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,N,&vec2);

  ierr = VecGetOwnershipRange(vec1,&start1,&end1);CHKERRQ(ierr);
  ierr = VecGetOwnershipRange(vec2,&start2,&end2);CHKERRQ(ierr);


  /*-----------------------------------------------
    Fill data in vectors
    -----------------------------------------------*/ 
  for(i=start1;i<end1;i++){
	v = (PetscScalar) i;
	ierr = VecSetValues(vec1,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
  	}
  v = 5;
  ierr = VecSet(vec2,v);CHKERRQ(ierr);


  /*-----------------------------------------------
    View vectors before scatter
    -----------------------------------------------*/ 
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Viewing vec1, before scatter\n");CHKERRQ(ierr);
  ierr = VecView(vec1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);
  ierr = PetscBarrier(NULL);CHKERRQ(ierr);

  ierr = PetscPrintf(PETSC_COMM_WORLD,"Viewing vec2, before scatter\n");CHKERRQ(ierr);
  ierr = VecView(vec2,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);
  ierr = PetscBarrier(NULL);CHKERRQ(ierr);


  /*-----------------------------------------------
    Create two IS which will be used for vecscatter
    -----------------------------------------------*/
  IS             is1,is2;
  
  PetscInt       *indices1;
  PetscInt       *indices2;
  
  PetscInt idx_size1 = end1 - start1;
  PetscInt idx_size2 = end2 - start2;
  
  ierr = PetscMalloc1(idx_size1,&indices1);CHKERRQ(ierr);
  ierr = PetscMalloc1(idx_size2,&indices2);CHKERRQ(ierr);
  
  for(i=0;i<idx_size1;i++){indices1[i] = rank*(PetscInt)(N/2) + i;}
  for(i=0;i<idx_size2;i++){indices2[i] = (N-1) - (rank*(PetscInt)(N/2) + i);}
  
  ierr = ISCreateGeneral(PETSC_COMM_SELF,idx_size1,indices1,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
  ierr = ISCreateGeneral(PETSC_COMM_SELF,idx_size2,indices2,PETSC_COPY_VALUES,&is2);CHKERRQ(ierr);
  
  ierr = PetscFree(indices1);CHKERRQ(ierr);
  ierr = PetscFree(indices2);CHKERRQ(ierr);

  /*Optionally, view the index sets via stdout*/
  //ierr = ISView(is1,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
  //ierr = ISView(is2,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);


  /*-----------------------------------------------
    Scatter, define and execute*
    -----------------------------------------------*/ 
  VecScatter newctx;
  ierr = VecScatterCreate(vec1,is1,vec2,is2,&newctx);CHKERRQ(ierr);
  ierr = VecScatterBegin(newctx,vec1,vec2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd(newctx,vec1,vec2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
 

  /*-----------------------------------------------
    View vectors after scatter
    -----------------------------------------------*/ 
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Viewing vec1, after scatter\n");CHKERRQ(ierr);
  ierr = VecView(vec1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);
  ierr = PetscBarrier(NULL);CHKERRQ(ierr);

  ierr = PetscPrintf(PETSC_COMM_WORLD,"Viewing vec2, after scatter\n");CHKERRQ(ierr);
  ierr = VecView(vec2,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);
  ierr = PetscBarrier(NULL);CHKERRQ(ierr);


  /*Free space*/
  ierr = VecDestroy(&vec1);CHKERRQ(ierr);
  ierr = VecDestroy(&vec2);CHKERRQ(ierr);
  ierr = ISDestroy(&is1);CHKERRQ(ierr);
  ierr = ISDestroy(&is2);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}
