HDF5 for data files written with fortran

Question

The HDF5 data storage uses the C convention, i.e. if I am storing a matrix A(N,M,K) in a binary file, the fastest changing dimension of the stored data will have size N. Apparently when I use the Fortran wrapper of HDF5, HDF5 automatically transposes the matrix, to be consistent with C.

I have a data of size (256 by 128 by 256) stored in an unformatted binary file written by fortran. I am trying to convert it into h5 format by using a program given below. But the final output is giving me the dimensions of the stored matrix as (128,256,256). I have no idea what to do to make sure that the final hd5 file can be rightly visualized in the visualizing software (Paraview).

PROGRAM H5_RDWT

 USE HDF5 ! This module contains all necessary modules

 IMPLICIT NONE


 CHARACTER(LEN=6), parameter :: out_file = "out.h5"  ! File name
 CHARACTER(LEN=6), parameter :: in_file  = "in.dat" ! File name
 CHARACTER(LEN=4), parameter :: dsetname =  "vort"! Dataset name
 CHARACTER(LEN=50) :: len

 INTEGER(HID_T) :: in_file_id   ! File identifier
 INTEGER(HID_T) :: out_file_id  ! File identifier
 INTEGER(HID_T) :: dset_id      ! Dataset identifier
 INTEGER(HID_T) :: dspace_id    ! Dataspace identifier

 INTEGER :: in_file_id = 23

 INTEGER     :: nx = 256, ny=128, nz=256

 INTEGER(HSIZE_T), DIMENSION(3) :: dims             ! Dataset dimensions
 INTEGER     ::   rank = 3                          ! Dataset rank

 INTEGER     ::   error                  ! Error flag
 INTEGER     ::   i, j, k, ii, jj, kk    ! Indices

 REAL,    allocatable :: buff_r(:,:,:)   ! buffer for reading from input file

 dims(1) = nx
 dims(2) = ny
 dims(3) = nz
 allocate(buff_r(nx,ny,nz))

 ! Read the input data.
 open (in_file_id,FILE=in_file,form='unformatted',access='direct',recl=4*nx*ny*nz)  
 read (in_file_id,rec=1) buff_r


 ! Initialize FORTRAN interface of HDF5.
 CALL h5open_f(error)

 ! Create a new file.
 CALL h5fcreate_f (out_file, H5F_ACC_TRUNC_F, out_file_id, error)

 ! Create the dataspace.
 CALL h5screate_simple_f(rank, dims, dspace_id, error)


 ! Create the dataset with default properties.
 CALL h5dcreate_f(out_file_id, dsetname, H5T_NATIVE_REAL, dspace_id, &
         dset_id, error)

 ! Write the dataset.
 CALL h5dwrite_f(dset_id, H5T_NATIVE_REAL, buff_r, dims, error)

 ! End access to the dataset and release resources used by it.
 CALL h5dclose_f(dset_id, error)

 ! Terminate access to the data space.
 CALL h5sclose_f(dspace_id, error)

 ! Close the file.
 CALL h5fclose_f(out_file_id, error)

 ! Close FORTRAN interface.
 CALL h5close_f(error)

 deallocate(buff_r)

 END PROGRAM H5_RDWT

---

To illustrate what is happening, I am generating a sample data file using the following script:

  program main

  !-------- initialize variables -------------
  character(8) :: fname
  integer, parameter:: n = 32
  real*8, dimension(n,n,2*n) :: re
  integer i,j,k, recl
  Inquire( iolength =  recl ) re

  !------ fill in the array with sample data ----

  do k = 1, 2*n
     do j = 1, n
        do i = 1, n
           re(i,j,k) = 1.0
        end do
     end do
  end do

  !------ write in data in a file -----------
  write(fname, "(A)") "data.dat"
  open (10, file=fname, form='unformatted', access='direct', recl=recl)
  write(10,rec=1) re
  close(10)

  stop
  end program main

I copy pasted the program by Ian Bush and changed the values of nx, ny and nz to 32, 32 and 64 respectively. I would expect the generated h5 file to have dimensions (32,32,64). But it is coming out to be (64,32,32). Here is what is happening in my machine:

[pradeep@laptop]$gfortran generate_data.f90 
[pradeep@laptop]$./a.out 
[pradeep@laptop]$ls -l data.dat 
-rw-r--r--  1 pradeep  staff  524288 Mar 12 14:04 data.dat
[pradeep@laptop]$h5fc convert_to_h5.f90 
[pradeep@laptop]$./a.out 
[pradeep@laptop]$ls -l out.h5 
-rw-r--r--  1 pradeep  staff  526432 Mar 12 14:05 out.h5
[pradeep@laptop]$h5dump -H out.h5 
HDF5 "out.h5" {
GROUP "/" {
   DATASET "data" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 64, 32, 32 ) / ( 64, 32, 32 ) }
   }
}
}

Please confirm with me if you are seeing the same thing.

Answer 1

I have also run into trouble with viewing HDF5 files that I've written with a Fortran application. The basic issue is that Fortran and C store multidimensional arrays differently (Fortran is column-major, C is row-major), and since the Fortran HDF5 libraries are interfaces into the C HDF5 libraries, the Fortran wrapper transposes the dimensions before passing the data into the C code. Likewise, when a Fortran application reads an HDF5 file, the Fortran wrapper transposes the dimensions again.

So if you do all your writing and reading with Fortran applications, you shouldn't notice any discrepancies. If you write the file with a Fortran app and then read it with a C app (such as h5dump), the dimensions will appear transposed. That's not a bug, it's just how it works.

If you want to display the data correctly, either read the data with a Fortran application or use a C app and transpose the data first. (Or you could transpose the data before writing it in the first place.)

As already mentioned, this is all explained fairly well in section 7.3.2.5 of the documentation: http://www.hdfgroup.org/HDF5/doc/UG/UG_frame12Dataspaces.html