在Julia中调用Fortran代码：进阶篇（动态数组）

在基础篇中，我们解决了下面的四个问题。

• Fortran function的调用
• Fortran subroutine的调用
• 定长数组
• 全局变量

在本篇中，我们将继续探索更有难度的情形。

例六：Fortran动态数组

在前面的例子中，我们处理了定长数组的情形。那么，如果遇到通过allocate方式动态分配的数组，应当如何处理呢？

Fortran代码：

module foo
    implicit none
    integer, allocatable :: arr(:)

contains
    subroutine bar(n)
        integer :: i, n
        if(allocated(arr)) deallocate(arr)
        allocate(arr(n))
        do i = 1,n
            arr(i) = (-1)**i
        end do
    end
end module foo

可以先尝试刚才的方法。

这里我们使用了unsafe_wrap(Type, Ptr, Size)函数来将一个Fortran的指针转换为Julia中的数组。

ccall((:__foo_MOD_bar, "f6"), Cvoid, (Ref{Int32},), 10)
arr = cglobal((:__foo_MOD_arr, "f6"), Int32)
println(unsafe_wrap(Array{Int32}, arr, (10,)))

但我们得到的结果是这样的（每次运行结果可能不一致）：

Int32[-739398176, 32666, -1, -1, 4, 0, 0, 257, 4, 0]

和我们的预期不一致。这是为什么呢？可以参考Introduction to the gfortran array descriptor
这篇文章对Fortran中数组描述符的介绍。我们利用cglobal拿到的指针，并不是指向我们希望得到的数据，而是指向了数组描述符中的地址元。也就是说，我们需要多走一步，首先利用指针得到地址，再把地址当作指针，进一步拿到数据。

ccall((:__foo_MOD_bar, "f6"), Cvoid, (Ref{Int32},), 10)
arr_addr_ptr = cglobal((:__foo_MOD_arr, "f6"), UInt64)
arr_addr = convert(Ptr{Int32}, unsafe_load(arr_addr_ptr))
println(unsafe_wrap(Array{Int32}, arr_addr, (10,)))

注意这里arr_addr_ptr的类型为UInt64。（使用Int64也能得到相同的结果，但可能存在风险）

输出的结果为：

Int32[-1, 1, -1, 1, -1, 1, -1, 1, -1, 1]

与我们的预期相符。

例七：Fortran动态数组（二维）

接下来，看看二维的情形。

Fortran代码：

module foo
    implicit none
    integer, allocatable :: arr(:, :)

contains
    subroutine bar(n)
        integer :: i, j, n
        if(allocated(arr)) deallocate(arr)
        allocate(arr(n, n))
        do i = 1,n
            do j = 1,n
                arr(i, j) = i * j
            end do
        end do
    end
end module foo

Julia代码：

ccall((:__foo_MOD_bar, "f7"), Cvoid, (Ref{Int32},), 9)
arr_addr_ptr = cglobal((:__foo_MOD_arr, "f7"), UInt64)
arr_addr = convert(Ptr{Int32}, unsafe_load(arr_addr_ptr))
println(unsafe_wrap(Array{Int32,2}, arr_addr, (9, 9)))

输出：

Int32[1 2 3 4 5 6 7 8 9; 2 4 6 8 10 12 14 16 18; 3 6 9 12 15 18 21 24 27; 4 8 12 16 20 24 28 32 36; 5 10 15 20 25 30 35 40 45; 6 12 18 24 30 36 42 48 54; 7 14 21 28 35 42 49 56 63; 8 16 24 32 40 48 56 64 72; 9 18 27 36 45 54 63 72 81]

可以看到，我们顺利地输出了九九乘法表。

---

我是吴子寒，欢迎关注我的GitHub账号@lucifer1004，以及我的个人站点CP-Wiki。