PROGRAM allreduce
#include <petsc/finclude/petsc.h>
  USE petsc
  IMPLICIT NONE
  
  ! --- pure PETSc
  PetscErrorCode :: ierr
  PetscScalar, ALLOCATABLE :: vals(:)
  PetscInt, ALLOCATABLE :: idxn(:)
  PetscInt :: m,i
  Vec :: vecX                ! eigenvector (real/imaginary parts)
  !
  REAL*8, ALLOCATABLE :: vecXcopy(:)
  INTEGER :: N

  ! initialize PETSc
  CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)

  ! Set up a new matrix
  m = 10
  ALLOCATE(vals(10),idxn(10))

  !
  CALL VecCreate(PETSC_COMM_WORLD,vecX,ierr); 
  CALL VecSetType(vecX,VECMPI,ierr); 
  CALL VecSetSizes(vecX,PETSC_DECIDE,m,ierr); 

  ! set values of vector
  DO i=1,m
     vals(i) = i
     idxn(i) = i-1
  END DO
  !
  PRINT*,'setting values'
  CALL VecSetValues(vecX,m,idxn,vals,INSERT_VALUES,ierr); 

  ! assemble vector
  PRINT*,'assembling vector'
  CALL VecAssemblyBegin(vecX,ierr); 
  CALL VecAssemblyEnd(vecX,ierr); 

  ! 
  PRINT*,'copying vector to array'
  N = m
  PRINT*,'N = ',N
  PRINT*,'m = ',m
  ALLOCATE(vecXcopy(N))
  CALL PetscVec2Array(vecX,vecXcopy,N)

  !
  CALL VecDestroy(vecX,ierr); 
  CALL PetscFinalize(ierr); 


END PROGRAM allreduce


  !> Copy Petsc vector to fortran array
  SUBROUTINE PetscVec2Array(vecX,buffer,N)
    ! ------
#include <petsc/finclude/petsc.h>
    USE petsc
    ! ------
    IMPLICIT NONE
    Vec, INTENT(in) :: vecX
    REAL*8, INTENT(inout) :: buffer(0:N-1)
    INTEGER, INTENT(in) :: N
    !
    INTEGER :: istart, iend, Np, i,j
    INTEGER :: me2,nprocs
    REAL*8, ALLOCATABLE :: copy_buffer(:)
    !
    PetscScalar, POINTER :: vecX_pt(:) => NULL()
    PetscInt :: petsc_N
    PetscErrorCode :: ierr

    !
    CALL MPI_COMM_RANK(MPI_COMM_WORLD,me2,ierr);
    CALL MPI_COMM_SIZE(MPI_COMM_WORLD,Nprocs,ierr);

    ! initailize buffer to 0.
    buffer = 0.
    ALLOCATE(copy_buffer(0:N-1))
   

    ! check length of petsc vector is N
    PRINT*,'getting vecX size'
    CALL VecGetSize(vecX,petsc_N,ierr); CHKERRQ(ierr)
    Np = petsc_N
    IF(N.NE.Np) THEN
       PRINT*,'ERROR: In PetscVec2Array:'
       PRINT*,'petsc_N = ',petsc_N, Np
       PRINT*,'N       = ',N
       CALL ABORT
    END IF

    ! vecX_pt will contain the element in the local range [istart,iend-1]
    PRINT*,'getting ownership range'
    CALL VecGetOwnershipRange(vecX,istart,iend,ierr); CHKERRQ(ierr)
    PRINT*,'ownership range: ',istart,iend

    ! initialize the pointer
    CALL VecGetArrayReadF90(vecX,vecX_pt,ierr); CHKERRQ(ierr)
    ! copy values to global array
    PRINT*,'copy values to global array'

    DO i=1,Nprocs
       !
       CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
       !
       IF(me2.EQ.i-1) THEN
          PRINT*,'vecX_pt:'
          PRINT*,'SIZE = ',SIZE(VecX_pt)
          DO j=1,SIZE(vecX_pt)
             PRINT*,'  ',vecX_pt(j)
          END DO
       END IF
    END DO


    copy_buffer(istart:iend-1) = vecX_pt(:)
    ! free pointer
    CALL VecRestoreArrayReadF90(vecX,vecX_pt,ierr); CHKERRQ(ierr)

    CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

    DO i=1,Nprocs
       !
       CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
       !
       IF(me2.EQ.i-1) THEN
          PRINT*,'copy_buffer = '
          DO j=istart,iend-1
             PRINT*,'  ',copy_buffer(j)
          END DO
       END IF
    END DO

    CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

    ! combine all components from all processes
    PRINT*,'calling allreduce'
    CALL MPI_ALLREDUCE(copy_buffer,buffer,N,MPI_REAL,MPI_SUM,MPI_COMM_WORLD,ierr)
    CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

    DO i=1,Nprocs
       !
       CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
       !
       IF(me2.EQ.i-1) THEN
          PRINT*,'buffer = '
          DO j=0,N-1
             PRINT*,'  ',buffer(j)
          END DO
       END IF
    END DO

    
    ! free local memory
    DEALLOCATE(copy_buffer)

  END SUBROUTINE PetscVec2Array