"C pointer trickery" to allow mismatched Fortran array ranks

Question

I'm writing a HDF5 wrapper subroutine that will read/write a double precision array of any shape from/to a dataset inside a HDF5 file. To achieve this, I use some C pointer trickery such that the subroutine takes in only the first element of the array as val, but it actually reads/writes the whole array using the temporary buffer buf(1:sz_buf).

So far I have the following for the read subroutine (after removing error checks to keep it concise):

SUBROUTINE hdf5_read_array_d( fname, path, name, val, dims )
  USE ISO_C_BINDING, ONLY: C_SIZE_T, C_LOC, C_F_POINTER

  ! Input arguments
  CHARACTER(LEN=*), INTENT(IN) :: fname, path, name
  REAL(KIND(1.D0)), TARGET, INTENT(OUT) :: val
  INTEGER, DIMENSION(:), INTENT(IN) :: dims

  ! Internal variables
  INTEGER(KIND=HID_T) :: h5root, h5path, h5dset
  INTEGER(KIND=HSIZE_T), DIMENSION(SIZE(dims)) :: h5dims
  REAL(KIND(1.D0)), DIMENSION(:), POINTER :: buf
  INTEGER(KIND=C_SIZE_T) :: sz_buf
  INTEGER :: dim

  ! Open the file in read-only mode
  CALL h5fopen_f( TRIM(fname), H5F_ACC_RDONLY_F, h5root, ierr )

  ! Open the pre-existing path in the file as a group
  CALL h5gopen_f( h5root, TRIM(path), h5path, ierr )

  ! Open the dataset
  CALL h5dopen_f( h5path, TRIM(name), h5dset, ierr )

  ! Convert dims to HSIZE_T
  h5dims(:) = dims(:)

  ! C pointer trickery: cast double -> void* -> double*
  sz_buf = PRODUCT(dims)
  ALLOCATE( buf( sz_buf ) )
  CALL C_F_POINTER( C_LOC(val), buf, (/sz_buf/) )

  ! Read data from dataset through buffer
  CALL h5dread_f( h5dset, H5T_NATIVE_DOUBLE, buf, h5dims, ierr )

  ! Clean up and close HDF5 file
  NULLIFY(buf)
  CALL h5dclose_f( h5dset, ierr )
  CALL h5gclose_f( h5path, ierr )
  CALL h5fclose_f( h5root, ierr )

  RETURN
END SUBROUTINE hdf5_read_array_d

Now, the question is, do I need to also put in DEALLOCATE(buf) in addition to / in place of the NULLIFY(buf)?

Any help would be appreciated.

Note: I am aware that Fortran 2018 includes assumed-rank arrays val(..) that will elegantly solve this problem. But again, it's a newer feature that might not be implemented by all compilers yet.

Edit: On C_F_POINTER(), here's a screenshot of Metcalf, Reid, and Cohen (4th Edition, not the newest one that has Fortran 2018 stuff):

user_3157076 · Answer

You can use C-style pointer trickery to do what you want, but you have some things to address in your approach:

you have a memory leak with allocate(buf)
you are (subtly) lying about the scalar nature of val
you'll horribly confuse anyone reading your code

The reason why this is horribly confusing, is because you don't need to do this trickery. That's also why I won't show you how to do it, or to address the question "do I need to deallocate as well as nullify?".

You know that you have an array val to stuff n values in, in a contiguous lump. You worry that that can't do that because you (without using an assumed-rank dummy) have to match array shape. Worry not.

  integer :: a(2,2,2,2), b(4,2,2), c(4,4)

are all arrays with 16 elements. So is

  integer :: d(16)

You can associate actual arguments a, b and c with dummy argument d. Let's see that in action:

  implicit none

  integer :: a(2,2,2,2), b(4,2,2), c(4,4)

  call set_them(a, SHAPE(a))
  call set_them(b, SHAPE(b))
  call set_them(c, SHAPE(c))

  print '(16I3)', a, b, c

contains

  subroutine set_them(d, dims)
    integer, intent(in)  :: dims(:)
    integer, intent(out) :: d(PRODUCT(dims))

    integer i
    d=[(i,i=1,SIZE(d))]
  end subroutine

end program

You can even associate array sections in this way to define portions.

You can see several other questions around here about this sequence association, in particular looking at changing shapes of arrays. This answer is more of a motivation of what to look for when tempted to do something complicated instead.

"C pointer trickery" to allow mismatched Fortran array ranks

1 Answers