fortran array reshape without temporary array

659 Views Asked by At

I'm passing some data that's read in by a C function to a Fortran procedure which will do all the number crunching. The array is naturally viewed as having shape (2, nn) in Fortran. Since C does not have multidimensional arrays, the array is allocated and used in C with length 2 * nn.

I'd like to reshape the array in Fortran so that I can use the convenient indexing in the numberical code. The Fortran routines would look something like this:

subroutine fortran_glue_code(n, array) bind(c)
    use iso_c_binding, only: c_int
    integer(c_int), intent(in) :: n, array(2 * n)
    real(double_precision) :: result

    call do_the_number_crunching(result, reshape(array, [2, n]))
    { do stuff with the result... }
end subroutine

subroutine do_the_number_crunching(result, array)    ! not C-interoperable
    real(double_precision), intent(out) :: result
    integer, intent(in) :: array(:,:)

    if (size(array, 1) /= 2) then
        print *, "Wrong size!"
        call exit(1)
    end if

    { crunch crunch crunch... }
end subroutine

What I'm not happy with is that reshape unnecessarily creates a temporary array, and with the data I'll be using, a very large one at that.

The procedure in question treats array as read-only, so I would think that, rather than make an array copy, the compiler could simply make a new Fortran array header that refers to the same location in memory for the array contents, just with different dimensions. Is there any way to avoid the array copy with reshape?

2

There are 2 best solutions below

1
On BEST ANSWER

I do not get your problem, why don't you just do

subroutine fortran_glue_code(n, array) bind(c)
    use iso_c_binding, only: c_int
    integer(c_int), intent(in) :: n, array(2, n)
    real(double_precision) :: result

    call do_the_number_crunching(result, array)
    { do stuff with the result... }
end subroutine

?

Alternatively, you could also use sequence association

subroutine fortran_glue_code(n, array) bind(c)
    use iso_c_binding, only: c_int
    integer(c_int), intent(in) :: n, array(2 * n)
    real(double_precision) :: result

    call do_the_number_crunching(result, array, n)
    { do stuff with the result... }
end subroutine

subroutine do_the_number_crunching(result, array, n)    ! not C-interoperable
    real(double_precision), intent(out) :: result
    integer, intent(in) :: array(2,n)
    integer, intent(in) :: n

    if (size(array, 1) /= 2) then
        print *, "Wrong size!"
        call exit(1)
    end if

    { crunch crunch crunch... }
end subroutine

but that is an unnecessary complication.

0
On

If you used int ** in the C code, then type (c_ptr) in the Fortran, you could use intrinsic subroutine c_f_pointer to associate the C-pointer with a Fortran array with pointer attribute, designating the dimensions of the Fortran array with the third argument of c_f_pointer. Guessing, it seems to me that the compiler would not need to perform a copy, but would instead create the Fortran pointer/structure, using the C-pointer address to locate the data.