chw21
chw21

Reputation: 8140

Best way to pass multi-dimensional array as 1D array

I have a subroutine that populates a rank 1 array of arbitrary length.

subroutine populate_array_1D(array)
    implicit none
    integer, dimension(:), intent(out) :: array
    array(:) = 1
end subroutine populate_array_1D

I now want to modify this routine using overloading to accept higher-dimensioned arrays.

interface populate_array
    module procedure populate_array_1D, populate_array_2D, populate_array_3D
end interface populate_array

What I want to do is to refer the higher-dimensioned subroutines back to the 1D subroutine.

My ideas so far:

  1. Slice the array, and pass the slices one-by-one:

    subroutine populate_array_2D(array)
        implicit none
        integer, dimension(:,:), intent(out) :: array
        integer :: i
        do i = lbound(array, 2), ubound(array, 2)
            call populate_array_1D(array(:, i))
        end do
    end subroutine populate_array_2D
    
  2. Create a new, automatic 1D array, and then use reshape:

    subroutine populate_array_3D(array)
        implicit none
        integer, dimension(:,:,:), intent(out) :: array
        integer :: tmp_array(size(array))
        call populate_array_1D(tmp_array)
        array = reshape(tmp_array, shape(array))
    end subroutine populate_array_3D
    

Assume that, like in this example, the result of either version will be okay, I'm only worried about the performance. (I should probably add that I can't just make an elemental subroutine because it can't be pure.)

Or is there an even better version, where I can pass a multidimensional array into a subroutine pretending it's a 1D array?

(Ultimately I want to replace the call to RANDOM_NUMBER with a PRNG that is as deterministic as possible across different compilers and/or machines, assuming they use the same seed.)

Upvotes: 2

Views: 797

Answers (1)

knia
knia

Reputation: 749

As francescalus pointed out, elemental procedures can be impure (as of Fortran 2008).

Alternatively, you could use Fortran's intrinsic C interoperability to gain low-level access to the memory location of the 2D array. The following subroutine builds a 1D pointer to the 2D array:

subroutine populate_array_2D(array)
    use iso_c_binding, only: C_F_POINTER, C_LOC

    integer, dimension(:,:), target, contiguous, intent(out) :: array

    integer, dimension(:), pointer :: ptr_1d

    call C_F_POINTER (C_LOC(array), ptr_1d, [size(array)])
    call populate_array_1D(ptr_1d)
end subroutine populate_array_2D

Edit:

The call to c_f_pointer can be replaced by the equivalent bound-remapping list in the pointer: ptr_1d(1:size(array)) => array (though compiling with ifort18 gives an erroneous warning message that this is not valid Fortran 2008).

Upvotes: 1

Related Questions