Hi,
I have a code that reads in an hdf5 file.
The code to read the file works with GNU and Intel with -O3.
With PGI 17.4 or 17.7 (I have not tested other versions), the code only works with -O1 or lower, but with -O2 and -O3 it either segfaults or says it is trying to reallocated an allocated array.
I used valgrind to pinpoint the routine.
Since the code works with -O1, the problem is not in the PGI compilation of the HDF5 library itself.
The code is given below.
First, the definition of the modules is given, then the code that reads the file.
c#######################################################################
module rp1d_def
c
c-----------------------------------------------------------------------
c ****** Define a structure to hold a REAL 1D pointer.
c-----------------------------------------------------------------------
c
use number_types
c
implicit none
c
type :: rp1d
real(r_typ), dimension(:), pointer :: f
end type
c
end module
c#######################################################################
module sds_def
c
c-----------------------------------------------------------------------
c ****** Definition of the SDS data structure.
c-----------------------------------------------------------------------
c
use number_types
use rp1d_def
c
implicit none
c
integer, parameter, private :: mxdim=3
c
type :: sds
integer :: ndim
integer, dimension(mxdim) :: dims
logical :: scale
logical :: hdf32
type(rp1d), dimension(mxdim) :: scales
real(r_typ), dimension(:,:,:), pointer :: f
end type
c
end module
c#######################################################################
The routine itself is:
c#######################################################################
subroutine rdh5 (fname,s,ierr)
c
c-----------------------------------------------------------------------
c
c ****** Read a 1D, 2D, or 3D scientific data set from an HDF5 file.
c ****** The HDF5 file is currently assumed to contain only one
c ****** dataset (1D,2d,3D), with or without scales, in group "/",
c ****** and has no other data members.
c
c-----------------------------------------------------------------------
c
c ****** This routine allocates the required memory and returns
c ****** pointers to the data and scale arrays.
c
c-----------------------------------------------------------------------
c
c ****** Input arguments:
c
c FNAME : [character(*)]
c HDF5 data file name to read from.
c
c ****** Output arguments:
c
c S : [structure of type SDS]
c A structure that holds the field, its
c dimensions, and the scales, with the
c components described below.
c
c IERR : [integer]
c IERR=0 is returned if the data set was read
c successfully. Otherwise, IERR is set to a
c nonzero value.
c
c ****** Components of structure S:
c
c NDIM : [integer]
c Number of dimensions found in the data set.
c
c DIMS : [integer, dimension(3)]
c Number of points in the data set dimensions.
c For a 1D data set, DIMS(2)=DIMS(3)=1.
c For a 2D data set, DIMS(3)=1.
c
c SCALE : [logical]
c Flag to indicate the presence of scales (axes)
c in the data set. SCALE=.false. means that scales
c were not found; SCALE=.true. means that scales
c were found.
c
c HDF32 : [logical]
c Flag to indicate the precision of the data set
c read in. HDF32=.true. means that the data is
c 32-bit; HDF32=.false. means that the data is
c 64-bit.
c
c SCALES : [structure of type RP1D, dimension(3)]
c This array holds the pointers to the scales
c when SCALE=.true., and is undefined otherwise.
c
c F : [real, pointer to a rank-3 array]
c This array holds the data set values.
c
c ****** The storage for the arrays pointed to by F, and the
c ****** scales (if present) in structure SCALES, is allocated by
c ****** this routine.
c
c-----------------------------------------------------------------------
c
use number_types
use sds_def
use hdf5
use h5ds
c
c-----------------------------------------------------------------------
c
implicit none
c
c-----------------------------------------------------------------------
c
type(sds) :: s
character(*) :: fname
c
c-----------------------------------------------------------------------
c
integer :: ierr
c
c-----------------------------------------------------------------------
c
integer :: i,obj_type,n_members
c
integer(HID_T) :: file_id ! File identifier
integer(HID_T) :: dset_id ! Dataset identifier
integer(HID_T) :: dspace_id ! Dataspace identifier
integer(HID_T) :: dim_id ! Dimension identifiers
integer(HID_T) :: datatype_id ! Datatype identifiers
c
integer(SIZE_T) :: prec
c
integer(HSIZE_T),dimension(3) :: s_dims
integer(HSIZE_T),dimension(1) :: s_dims_i
integer(HSIZE_T),dimension(1) :: totpts
c
real(KIND_REAL_4), dimension(:,:,:), allocatable :: f4
real(KIND_REAL_4), dimension(:), allocatable :: f4dim
real(KIND_REAL_8), dimension(:,:,:), allocatable :: f8
real(KIND_REAL_8), dimension(:), allocatable :: f8dim
c
character(512) :: obj_name
character(4), parameter :: cname='RDH5'
c
logical :: is_scale,is_attached
c
c-----------------------------------------------------------------------
c
c
c ****** Initialize dimension count and arrays.
c
s%ndim=0
s%dims(:)=1
s_dims(:)=1
c
c ****** Initialize hdf5 interface.
c
call h5open_f (ierr)
c
c ****** Open hdf5 file.
c
call h5Fopen_f (trim(fname),H5F_ACC_RDONLY_F,file_id,ierr)
c
c ****** Get information about the hdf5 file.
c
call h5Gn_members_f (file_id,"/",n_members,ierr)
c
c ****** Make sure there is (at maximum) one 3D dataset with scales.
c
if (n_members.eq.0.or.n_members.gt.4) then
write (*,*)
write (*,*) '### ERROR in ',cname,':'
write (*,*) '### Input file contains too few/many datasets.'
write (*,*) 'File name: ',trim(fname)
return
endif
c
c ****** Assume the Dataset is in index 0 and get its name.
c
call h5Gget_obj_info_idx_f (file_id,"/",0,obj_name,obj_type,ierr)
c
c ****** Open Dataset.
c
call h5Dopen_f (file_id,trim(obj_name),dset_id,ierr)
c
c ****** Make sure the Dataset is not a scale.
c
call h5DSis_scale_f(dset_id,is_scale,ierr)
if (is_scale) then
write (*,*)
write (*,*) '### ERROR in ',cname,':'
write (*,*) '### Input file Dataset at index 0 is a scale.'
write (*,*) 'File name: ',trim(fname)
return
endif
c
c ****** Get dimensions (need s_dims array for format requirements).
c
call h5Dget_space_f (dset_id,dspace_id,ierr)
call h5Sget_simple_extent_ndims_f (dspace_id,s%ndim,ierr)
call h5Sget_simple_extent_dims_f (dspace_id,s_dims,totpts,ierr)
s%dims(:)=s_dims(:)
c
c ****** Get the floating-point precision of the data and set flag.
c
call h5Dget_type_f (dset_id,datatype_id,ierr)
call h5Tget_precision_f (datatype_id,prec,ierr)
c
if (prec.eq.32) then
s%hdf32=.true.
elseif (prec.eq.64) then
s%hdf32=.false.
end if
c
c ****** Allocate the memory for the Dataset array in s.
c
allocate (s%f(s%dims(1),s%dims(2),s%dims(3)))
c
c ****** Need to read the file in its own datatype, and then convert
c ****** to datatype of s%f.
c
if (s%hdf32) then
allocate (f4(s%dims(1),s%dims(2),s%dims(3)))
call h5Dread_f (dset_id,datatype_id,f4,s_dims,ierr)
s%f(:,:,:)=f4(:,:,:)
deallocate (f4)
else
allocate (f8(s%dims(1),s%dims(2),s%dims(3)))
call h5Dread_f (dset_id,datatype_id,f8,s_dims,ierr)
s%f(:,:,:)=f8(:,:,:)
deallocate (f8)
end if
c
if (ierr.ne.0) then
write (*,*)
write (*,*) '### ERROR in RDH5:'
write (*,*) '### Error while reading the dataset.'
write (*,*) 'File name: ',trim(fname)
write (*,*) '[Error return (from H5DREAD_F) = ',ierr,']'
ierr=4
return
end if
c
c ****** Close the hdf5 type descriptor.
c
call h5Tclose_f (datatype_id,ierr)
c
c ****** Check if there might be scales present, if so, read them.
c
if (n_members.gt.1) then
c
c ***** First check that the number of scale datasets match the # dim.
c
if (n_members-1.ne.s%ndim) then
write (*,*)
write (*,*) '### ERROR in RDH5:'
write (*,*) '### # scales does not match # dims.'
write (*,*) 'File name: ',trim(fname)
return
end if
c
s%scale=.true.
c
c ****** Loop through scales, make sure each is a scale, and read them.
c
do i=1,s%ndim
c
c ****** Get the name of scale dataset.
c
call h5Gget_obj_info_idx_f (file_id,"/",i,
& obj_name,obj_type,ierr)
c
c ****** Open scale dataset.
c
call h5Dopen_f (file_id,trim(obj_name),dim_id,ierr)
c
c ****** Make sure the scale is a scale.
c
call h5DSis_scale_f (dim_id,is_scale,ierr)
if (.not.is_scale) then
write (*,*)
write (*,*) '### ERROR in RDH5:'
write (*,*) '### Scale is not a scale.'
write (*,*) 'File name: ',trim(fname)
return
end if
c
c ****** Make sure the scale is attached to the dim we think it is.
c [This is buggy, using 0-idx or 1-idx on some systems.
c until worked out, ignore this error check.]
c
c call h5DSis_attached_f (dim_id,dset_id,i,is_attached,ierr)
c if (.not.is_attached) then
c write (*,*)
c write (*,*) '### ERROR in RDH5:'
c write (*,*) '### Scale is not attached to dataset.'
c write (*,*) 'File name: ',trim(fname)
c return
c end if
c
c ****** Get dimension of scale.
c
s_dims_i(1)=s_dims(i)
c
c ****** Allocate scale.
c
allocate (s%scales(i)%f(s_dims_i(1)))
c
c ****** Get the floating-point precision of the scale.
c
call h5Dget_type_f (dim_id,datatype_id,ierr)
call h5Tget_precision_f (datatype_id,prec,ierr)
c
c ****** Read in the scale data.
c
if (s%hdf32) then
allocate (f4dim(s_dims_i(1)))
call h5Dread_f (dim_id,datatype_id,f4dim,s_dims_i,ierr)
s%scales(i)%f(:)=f4dim(:)
deallocate (f4dim)
else
allocate (f8dim(s_dims_i(1)))
call h5Dread_f (dim_id,datatype_id,f8dim,s_dims_i,ierr)
s%scales(i)%f(:)=f8dim(:)
deallocate (f8dim)
end if
c
c ****** Close the scale dataset.
c
call h5Dclose_f (dim_id,ierr)
c
enddo
c
c ****** Allocate dummy scales (of length 1) for empty dimensions.
c
do i=s%ndim+1,3
allocate (s%scales(i)%f(1))
enddo
else
c
c ****** If scales are not present, allocate dummy
c ****** scales (of length 1) so that the pointers to the scales
c ****** are valid.
c
s%scale=.false.
c
allocate (s%scales(1)%f(1))
allocate (s%scales(2)%f(1))
allocate (s%scales(3)%f(1))
end if
c
c ****** Close the dataset.
c
call h5Dclose_f (dset_id,ierr)
c
c ****** Close the file.
c
call h5Fclose_f (file_id,ierr)
c
c ****** Close FORTRAN interface.
c
call h5close_f (ierr)
c
return
end subroutine
c#######################################################################