[petsc-users] VecGetArrayReadF90 with gfortran

[Samuel Lanthaler]

Hi all,

I am having trouble using the function VecGetArrayReadF90 under gfortran. I
have created a minimal example and put it in the attachment. Basically, it
appears that the input pointer is empty upon return from
VecGetArrayReadF90. The code runs fine, when compiled with ifort, so I
don't know what's going on.

I am attaching the code and makefile that I used, as well as the command
line output that I can see. Though the errors are just to do with the fact
that the output pointer is empty, so the command line output is not all
that helpful...

Thanks in advance for your help!
Sam

 setting values
 assembling vector
 copying vector to array
 N =           10
 m =           10
 getting vecX size
 getting ownership range
 ownership range:            0          10
 copy values to global array
 vecX_pt:
 SIZE =            0
At line 127 of file allreduce.f90
Fortran runtime error: Array bound mismatch for dimension 1 of array 
'copy_buffer' (10/0)

Error termination. Backtrace:
#0  0x402621 in petscvec2array_
        at /home/lanthale/Progs/test_PetscVec2Array/allreduce.f90:127
#1  0x403550 in allreduce
        at /home/lanthale/Progs/test_PetscVec2Array/allreduce.f90:58
#2  0x4035d8 in main
        at /home/lanthale/Progs/test_PetscVec2Array/allreduce.f90:4

makefile
Description: Binary data

PROGRAM allreduce
#include "slepc/finclude/slepc.h"
  ! ------
  USE slepcsys
  USE slepceps
  IMPLICIT NONE
  
  ! --- pure PETSc
  PetscErrorCode :: ierr
  PetscScalar, ALLOCATABLE :: vals(:)
  PetscInt, ALLOCATABLE :: idxn(:)
  PetscInt :: m,i
  ! --- SLEPc
  Vec :: vecX                ! eigenvector (real/imaginary parts)
  !
  COMPLEX, ALLOCATABLE :: vecXcopy(:)
  INTEGER :: N

  ! initialize SLEPc
  CALL SlepcInitialize(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 SlepcFinalize(ierr); 


END PROGRAM allreduce


  !> Copy Petsc vector to fortran array
  SUBROUTINE PetscVec2Array(vecX,buffer,N)
#include "slepc/finclude/slepc.h"
    ! ------
!    USE petscsys
!    USE petscvec
    ! ------
    USE slepcsys
    USE slepceps
    IMPLICIT NONE
    Vec, INTENT(in) :: vecX
    COMPLEX, INTENT(inout) :: buffer(0:N-1)
    INTEGER, INTENT(in) :: N
    !
    INTEGER :: istart, iend, Np, i,j
    INTEGER :: me2,nprocs
    COMPLEX, ALLOCATABLE :: copy_buffer(:), copy_buffer2(:)
    !
    PetscScalar, POINTER :: vecX_pt(:) => NULL()
    PetscInt :: petsc_istart, petsc_iend, 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_DOUBLE_COMPLEX,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