write.f90 Source File


This file depends on

sourcefile~~write.f90~~EfferentGraph sourcefile~write.f90 write.f90 sourcefile~interface.f90 interface.f90 sourcefile~write.f90->sourcefile~interface.f90

Files dependent on this one

sourcefile~~write.f90~~AfferentGraph sourcefile~write.f90 write.f90 sourcefile~writer_lt.in.f90 writer_lt.in.f90 sourcefile~writer_lt.in.f90->sourcefile~write.f90 sourcefile~writer_nd.f90 writer_nd.f90 sourcefile~writer_nd.f90->sourcefile~write.f90 sourcefile~writer.in.f90 writer.in.f90 sourcefile~writer.in.f90->sourcefile~write.f90

Contents

Source Code


Source Code

submodule (h5fortran) write
!! This submodule is for writing HDF5 data via child submodules
use hdf5, only: &
h5screate_f, H5S_SCALAR_F, &
h5dcreate_f, &
h5pset_chunk_f, h5pset_layout_f, h5pset_deflate_f, h5pset_shuffle_f, h5pset_fletcher32_f, h5pcreate_f, h5pclose_f, &
H5P_DATASET_CREATE_F, &
h5gopen_f, &
H5Lcreate_soft_f

use H5LT, only: h5ltmake_dataset_string_f, h5ltpath_valid_f

implicit none (type, external)

contains


module procedure hdf_create

logical :: exists
integer :: ierr
integer(HID_T) :: pid, space_id, ds_id

if(.not.self%is_open) error stop 'h5fortran:write: file handle is not open: ' // self%filename

call h5ltpath_valid_f(self%lid, dname, .true., exists, ierr)
!! h5lexists_f can false error with groups--just use h5ltpath_valid

!! stricter than self%exists() since we're creating and/or writing variable
if (ierr /= 0) error stop 'ERROR:h5fortran:create: variable path invalid: ' // dname // ' in ' // self%filename

if(self%debug) print *,'h5fortran:TRACE:create:exists: ' // dname, exists

if(exists) then
  if (.not.present(istart)) call hdf_shape_check(self, dname, dims)
  !! FIXME: read and write slice shape not checked; but should check in future versions
  !> open dataset
  call h5dopen_f(self%lid, dname, ds_id, ierr)
  if (ierr /= 0) error stop 'ERROR:h5fortran:create: could not open ' // dname // ' in ' // self%filename

  if(present(did)) did = ds_id
  if(present(sid)) then
    call h5dget_space_f(ds_id, sid, ierr)
    if(ierr /= 0) error stop 'h5fortran:create could not get dataset ' // dname // ' in ' // self%filename
  end if
  return
endif

if(self%debug) print *, 'h5fortran:TRACE1: ' // dname

!> Only new datasets go past this point

call self%write_group(dname, ierr)
if (ierr /= 0) error stop 'ERROR:h5fortran:create: could not create group for ' // dname // ' in ' // self%filename


!> create properties
pid = 0 !< sentinel

if(size(dims) >= 2) then
  if(self%debug) print *, 'h5fortran:TRACE:create: deflate: ' // dname
  call set_deflate(self, dims, pid, ierr, chunk_size)
  if (ierr/=0) error stop 'ERROR:h5fortran:create: problem setting deflate on ' // dname
endif

if(present(compact)) then
! print *, "TRACE1: COMPACT", compact
!! don't set COMPACT after CHUNKED, will fail. And it's either or anyway.
if(compact .and. pid == 0 .and. product(dims) * 8 < 60000)  then
!! 64000 byte limit, here we assumed 8 bytes / element
  call h5pcreate_f(H5P_DATASET_CREATE_F, pid, ierr)
  if (check(ierr, self%filename)) error stop "h5fortran:h5pcreate: " // dname

  call h5pset_layout_f(pid, H5D_COMPACT_F, ierr)
  if (check(ierr, self%filename)) error stop "h5fortran:h5pset_layout: " // dname
endif
endif

!> create dataspace
if(size(dims) == 0) then
  call h5screate_f(H5S_SCALAR_F, space_id, ierr)
else
  call h5screate_simple_f(size(dims), dims, space_id, ierr)
endif
if (check(ierr, self%filename, dname)) error stop "h5fortran:h5screate: " // dname

!> create dataset
if(pid == 0) then
  call h5dcreate_f(self%lid, dname, dtype, space_id, ds_id, ierr)
else
  call h5dcreate_f(self%lid, dname, dtype, space_id, ds_id, ierr, pid)
  if (check(ierr, self%filename, dname)) error stop "h5fortran:h5dcreate: " // dname
  call h5pclose_f(pid, ierr)
  if (check(ierr, self%filename, dname)) error stop "h5fortran:h5pclose: " // dname
endif
if (check(ierr, self%filename, dname)) error stop "h5fortran:h5dcreate: " // dname

if(.not.(present(did) .and. present(sid))) then
  if(self%debug) print *, 'h5fortran:TRACE:create: closing dataset ', dname
  call hdf_wrapup(ds_id, space_id, ierr)
endif
if (check(ierr, self%filename, dname)) error stop "h5fortran:wrapup: " // dname

if(present(sid)) sid = space_id
if(present(did)) did = ds_id

end procedure hdf_create


module procedure create_softlink
!! HDF5 soft link -- to variables in same file
!! target need not exist (dangling link)
!! linking to external files requires an external link (different function required)

integer :: ierr

call H5Lcreate_soft_f(target, self%lid, link, ierr)
if (check(ierr, self%filename)) return

end procedure create_softlink


subroutine set_deflate(self, dims, pid, ierr, chunk_size)
class(hdf5_file), intent(inout) :: self
integer(HSIZE_T), intent(in) :: dims(:)
integer(HID_T), intent(out) :: pid
integer, intent(out) :: ierr
integer, intent(in), optional :: chunk_size(:)

integer(HSIZE_T) :: cs(size(dims))

pid = 0
if (self%comp_lvl < 1 .or. self%comp_lvl > 9) return

if (present(chunk_size)) then
  cs = chunk_size
  where (cs > dims) cs = dims
  if(self%debug) print *,'TRACE: user request chunk_size ',cs
else
  !! guess chunk size, keeping in mind 1 Megabyte recommended maximum chunk size
  call guess_chunk_size(dims, cs)
endif

if(any(cs < 1)) return

if(self%debug) print *,'DEBUG:set_deflate: dims: ',dims,'chunk size: ', cs

call h5pcreate_f(H5P_DATASET_CREATE_F, pid, ierr)
if (check(ierr, self%filename)) return

call h5pset_chunk_f(pid, size(dims), cs, ierr)
if (check(ierr, self%filename)) return

call h5pset_shuffle_f(pid, ierr)
if (check(ierr, self%filename)) return

call h5pset_fletcher32_f(pid, ierr)
if (check(ierr, self%filename)) return

call h5pset_deflate_f(pid, self%comp_lvl, ierr)
if (check(ierr, self%filename)) return

if(self%debug) print *,'TRACE:set_deflate done'

end subroutine set_deflate


subroutine guess_chunk_size(dims, chunk_size)
!! based on https://github.com/h5py/h5py/blob/master/h5py/_hl/filters.py
!! refer to https://support.hdfgroup.org/HDF5/Tutor/layout.html
integer(HSIZE_T), intent(in) :: dims(:)
integer(HSIZE_T), intent(out) :: chunk_size(:)

integer(hsize_t), parameter :: &
CHUNK_BASE = 16000, &    !< Multiplier by which chunks are adjusted
CHUNK_MIN = 8000, &      !< lower limit: 8 kbyte
CHUNK_MAX = 1000000, &   !< upper limit: 1 Mbyte
TYPESIZE = 8             !< bytes, assume real64 for simplicity

integer(hsize_t) :: dset_size, target_size, chunk_bytes, i, j, ndims

if (product(dims) * TYPESIZE < CHUNK_MIN) then
  chunk_size = 0
  return
endif

ndims = size(chunk_size)
chunk_size = dims

dset_size = product(chunk_size) * TYPESIZE
target_size = int(CHUNK_BASE * (2**log10(real(dset_size) / 1e6)), hsize_t)
if (target_size > CHUNK_MAX) target_size = CHUNK_MAX

! print *,'target_size [bytes]: ',target_size

i = 0
do
  !! Repeatedly loop over the axes, dividing them by 2.
  !! Stop when:
  !!   1a. We're smaller than the target chunk size, OR
  !!   1b. We're within 50% of the target chunk size, AND
  !!    2. The chunk is smaller than the maximum chunk size

  chunk_bytes = product(chunk_size) * TYPESIZE

  if ((chunk_bytes < target_size .or. 2*(abs(chunk_bytes-target_size) / target_size) < 1) .and. &
     chunk_bytes < CHUNK_MAX) exit

  if (product(chunk_size) == 1) exit
  !! Element size larger than CHUNK_MAX
  j = int(modulo(i, ndims), hsize_t) + 1
  if (j < 1 .or. j > ndims) error stop 'h5fortran: auto index bounds error'
  chunk_size(j) = ceiling(real(chunk_size(j)) / 2.0)
  i = i+1
end do

end subroutine guess_chunk_size


module procedure hdf_open_group

integer :: ier

if(.not.self%is_open) error stop 'h5fortran:open_group: file handle is not open: ' // self%filename

call h5gopen_f(self%lid, gname, self%gid, ier)

if (present(ierr)) ierr = ier
if (check(ier, self%filename, gname)) then
  if (present(ierr)) return
  error stop
endif

self%glid = self%lid
self%lid  = self%gid

end procedure hdf_open_group


module procedure hdf_close_group

integer :: ier

if(.not.self%is_open) error stop 'h5fortran:close_group: file handle is not open: ' // self%filename

call h5gclose_f(self%gid, ier)

if (present(ierr)) ierr = ier
if (check(ier, self%filename)) then
  if (present(ierr)) return
  error stop
endif

self%lid = self%glid

end procedure hdf_close_group


end submodule write