Code works with -O1, segfaults with -O2 or -O3

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#######################################################################

Hello,

Please send the entire failing program to trs@pgroup.com,
along with any files needed to recreate the segfault failure at runtime.

Also indicate how hdf5 was built, so we can duplicate.

thanks,
dave

Hi,

This problem can be avoided by using “-Mlarge_arrays” when compiling the code (using 17.9). I am not sure why since the arrays I was using were less than 2GB.

  • Ron

-Mlarge_arrays may be necessary when the total amount if array space
is greater than 2gb , not just individual arrays.

dave

HI

As I was just gathering together a reproducible example of this issue, I tried it with -O1 and -O3 and it seems to work now.

I am now using PGI 18.10.

Was this issue fixed since 17.9 and I missed the bug report?

Hi Ron,

Did you send something to Dave and have a TPR# I can check?

If not then I’d still need a reproducer to determine what the original problem was as well as when and how it got fixed. Or at the very least, I can determine when it started working.

It’s possible another user sent in a similar issue that we fixed, one of our engineers ran upon it independently, or possibly it was fixed when a separate problem was addressed.

-Mat

Hi Matt,

I have a very large suspicion that the issue in this post is due to the bug you found in the RDH5 routine.
Maybe -O1 was putting things on the heap and -O3 on the stack?
I do not think it is a coincidence that it is this exact routine that had the bug in it.

If you think this is a reasonable explanation, you can tell whoever has the TPR to close it.

Thanks!

  • Ron