Reputation: 8140
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:
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
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
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
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