I'm attaching the modified code - that appears to work fine for me. My changes
are:
1. use .F90 suffix - so the compiler does the f90 freeform/preprocessing
automatically
2. remove slepc references and use petsc only
3. compile with a petsc formatted makefile
4. replace complex with real*8 - thats compatible with PetscScalar
So perhaps you can modify your version of test code - to be closer to
this one [or vice versa] to see what is triggering the difference in
behavior.
Satish
On Fri, 23 Mar 2018, Samuel Lanthaler wrote:
> 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
>
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
CFLAGS =
FFLAGS =
CPPFLAGS =
FPPFLAGS =
include ${PETSC_DIR}/lib/petsc/conf/variables
include ${PETSC_DIR}/lib/petsc/conf/rules
ex1f: ex1f.o chkopts
-${FLINKER} -o ex1f ex1f.o ${PETSC_KSP_LIB}
${RM} ex1f.o
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 = 10
1.0000000000000000
2.0000000000000000
3.0000000000000000
4.0000000000000000
5.0000000000000000
6.0000000000000000
7.0000000000000000
8.0000000000000000
9.0000000000000000
10.000000000000000
copy_buffer =
1.0000000000000000
2.0000000000000000
3.0000000000000000
4.0000000000000000
5.0000000000000000
6.0000000000000000
7.0000000000000000
8.0000000000000000
9.0000000000000000
10.000000000000000
calling allreduce
buffer =
1.0000000000000000
2.0000000000000000
3.0000000000000000
4.0000000000000000
5.0000000000000000
0.0000000000000000
0.0000000000000000
0.0000000000000000
0.0000000000000000
0.0000000000000000
