program vecscatter_test

  implicit none

#include "finclude/petscsys.h"
#include "finclude/petscis.h"
#include "finclude/petscvec.h"

  MPI_Comm    :: comm
  PetscMPIInt :: myrank, mysize
  PetscErrorCode :: ierr
  PetscInt    :: ii,offset
  PetscViewer :: viewer

  PetscInt,pointer :: tmp_int_array(:),tmp_int_array1(:),tmp_int_array2(:)
  PetscInt     :: ncells_from,ncells_to
  IS           :: is_to, is_from
  Vec :: vec_from, vec_to
  VecScatter :: vscat

  call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

  comm = PETSC_COMM_WORLD
  call MPI_Comm_rank(comm,myrank,ierr)
  call MPI_Comm_size(comm,mysize,ierr)

  if(myrank==0) then
     ncells_from = 1
     ncells_to   = 2
  endif
  if(myrank==1) then
     ncells_from = 2
     ncells_to   = 1
  endif

  offset = 0
  call MPI_Exscan(ncells_from,offset,1,MPIU_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
  allocate(tmp_int_array(ncells_from))
  do ii = 1,ncells_from
     tmp_int_array(ii) = ii - 1 + offset
  enddo
  call ISCreateGeneral(PETSC_COMM_WORLD,ncells_from,tmp_int_array, &
       PETSC_COPY_VALUES,is_from,ierr)
  deallocate(tmp_int_array)

  offset = 0
  call MPI_Exscan(ncells_to,offset,1, &
                  MPIU_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
  allocate(tmp_int_array(ncells_to))
  do ii = 1,ncells_to
     tmp_int_array(ii) = ii - 1+offset
  enddo
  call ISCreateGeneral(PETSC_COMM_WORLD,ncells_to,tmp_int_array, &
       PETSC_COPY_VALUES,is_to,ierr)
  deallocate(tmp_int_array)

  call ISView(is_from,PETSC_VIEWER_STDOUT_WORLD,ierr)
  call ISView(is_to,PETSC_VIEWER_STDOUT_WORLD,ierr)


  call VecCreateMPI(PETSC_COMM_WORLD,ncells_from,PETSC_DECIDE,vec_from,ierr)
  call VecCreateMPI(PETSC_COMM_WORLD,ncells_to  ,PETSC_DECIDE,vec_to  ,ierr)

  call VecView(vec_from,PETSC_VIEWER_STDOUT_WORLD,ierr)
  call VecView(vec_to,PETSC_VIEWER_STDOUT_WORLD,ierr)

  call VecScatterCreate(vec_from,is_from,vec_to,is_to,vscat,ierr)
  call PetscViewerASCIIOpen(PETSC_COMM_WORLD, 'vscat.out', viewer, ierr)
  !call VecScatterView(vscat,viewer,ierr)
  call PetscViewerDestroy(viewer,ierr)

  call ISDestroy(is_from,ierr)
  call ISDestroy(is_to,ierr)
  call VecDestroy(vec_from,ierr)
  call VecDestroy(vec_to,ierr)

end program vecscatter_test