r - 包装fortran程序以在R中使用

Question

我正在使用 R，但需要做很多数字运算，我想在 fortran 中做。我对 R 比较陌生，对 fortran 还是个新手……我已经有一个可以工作的 R 程序，我想对其进行优化。我创建了一个解决 ODE 系统的 Fortran 程序，我将其保存在一个子程序中。我还使用了一个名为 aux.f90 的模块来存储参数和一个函数，该函数创建一个输入方程的信号。这按预期工作，数据保存在 .txt 文件中。

我现在想做的是创建一个 R 前端，它向 Fortran 程序扔参数，例如模拟的长度或解决方案中使用的步骤数。然后 Fortran 完成繁重的工作，将结果保存在文件中，我可以使用 R 来可视化文件中的数据。请参阅下面的 Fortran 代码:

! The auxiliary module contains all parameters
module aux
implicit none

integer,parameter :: n = 2**8       ! number of steps 
real(kind=4) :: jout = 0.5          ! for normal metabolism
real(kind=4) :: alp = 4.0           ! factor for growth
real(kind=4) :: bet = 1.0           ! benefit value
real(kind=4) :: etay = 0.1          ! cost value y
real(kind=4) :: etaz = 0.10         ! cost value z
real(kind=4) :: h  = 3.0            ! Hill coefficient
real(kind=4) :: Kx = 0.07           ! saturation coefficient
real(kind=4) :: t1 = 0.0, t2 = 30.0 ! start and end point of the simulation

contains                            ! function f(t) to create a step signal

real(kind=4) function f(t)
    implicit none
    real(kind=4), intent(in)  :: t  ! taking time value
    real(kind=4)              :: tt 
    !real(kind=4), intent(out) :: s  ! giving out the signal 
    real(kind=4) :: period = 5      ! defining the period length

    tt = MODULO(t,period)

        ! Signal function
        if (tt > 0.5*period) then 
            f = 1
        else
            f = 0
        endif
end function

end module aux

! The program solving the ODE system and giving the output as a fileprogram ffl
program ffl
use aux         ! Use module aux
implicit none
integer                     :: m,j          ! iteration variable
real(kind=4), dimension(6)  :: b =(/0.0, 0.2, 0.4, 0.6, 0.8, 1.0/) ! expression
real(kind=4)                :: dt                                  ! time resolution
real(kind=4), dimension(n)  :: t                                   ! time vector
real(kind=4), dimension(4)  :: x_new, x_aux, y_new, y_aux, q0      ! vectors 

! computing the time vector 
dt=(t2-t1)/real(n) ! calculating time resolution
t(1) = t1          ! setting first time value to t1 = 0

do m = 1,n         ! filling the time vector
    t(m) = t(m-1)+dt
end do
open(unit = 10,file = 'ffl.txt', status = 'unknown')
do j=1,6
    k = b(j)
    ! initial conditions
    q0(1) = 0   ! x
    q0(2) = k   ! y
    q0(3) = 0   ! z
    q0(4) = 0   ! w

    !open(unit = 10,file = 'ffl.txt', status = 'unknown')
    x_new = q0                ! set initial conditions
    write(10,*)t(1),x_new(1),x_new(2),x_new(3),x_new(4) ! saving data

    do m = 2,n
        ! Solving with a second order method 
        call derivate(t(m-1),x_new,y_new) ! call the derivate routine
        x_aux = x_new + dt*y_new

        call derivate(t(m),x_aux,y_aux)
        x_new = x_new + 0.5*dt*(y_new+y_aux)
        write(10,*)t(m),x_new(1),x_new(2),x_new(3),x_new(4) ! saving data

    end do
 end do
 close(10)
end program ffl

! The subroutine derivate gives the system of ODE's to be solved

subroutine derivate(time,y,z)
use aux              ! Use module aux
implicit none
real(kind=4), intent(in)               :: time ! input: time vector  
real(kind=4), dimension(4), intent(in) :: y    ! input: initial conditions vector
real(kind=4), dimension(4), intent(out):: z    ! output: results vector

z(1) = f(time)-y(1)                                                             !dx/dt 
z(2) = k+(1-k)*((1+Kx)*(y(1)**h))/((y(1)**h)+Kx)-y(2)                           !dy/dt
z(3) = ((1+Kx)*(y(1)**h)/((y(1)**h)+Kx)*((1+Kx)*(y(2)**h))/((y(2)**h)+Kx)-y(3)) !dz/dt
z(4) = f(time)*y(3)-etay*(k+(1-k)*((1+Kx)*(y(1)**h))/((y(1)**h)+Kx)) &          !dw/dt
       -etaz*(((1+Kx)*(y(1)**h))/((y(1)**h)+Kx)*((1+Kx)*(y(2)**h))/((y(2)**h)+Kx)) 

结束子程序派生

我已阅读“编写 R 扩展”文档，但没有发现它很有帮助...

现在问题:由于 R 需要一个 Fortran 子例程，我想在 fortran 中创建一个包装子例程，它利用我现有的文件，然后我可以从 R 调用它。但是，我找不到创建这个的方法首先是包装子程序。甚至可以在子程序中调用实际程序吗？我在网上找不到任何有用的东西。

Answer 1

程序应该作为可执行文件链接，因此您不能从子程序调用它 - 或者您调用可执行文件(在 gfortran 中使用 SYSTEM)，但您可以这样做 directly from R .

从 R 调用 Fortran 的简单方法是 .Fortran R 函数，调用 Fortran subroutine (不是 function ，也不是 program )。

基本步骤是:

编译 Fortran DLL，导出您需要的子程序(当然它们可能是其他子程序或函数的包装器)

将 DLL 放入系统路径中的目录

从 R 中，使用 dyn.load 加载 DLL

使用 .Fortran 调用您的子程序.

如果你使用 gfortran，你可以只安装 Rtools，它有你需要的一切。如果你想使用另一个编译器，你可能会遇到一些麻烦，尤其是名称。

从您的评论到 user2188538 的回答，我看到您已经知道所有这些步骤，但要非常小心符号名称。来自 .Fortran帮助:使用 .Fortran 时要注意编译的 Fortran 9x 代码:如果使用的 Fortran 9x 编译器与配置 R 时使用的 Fortran 77 编译器不同，它可能不起作用，尤其是在子例程名称不是小写或包含下划线的情况下。也可以使用 .C 并自己进行任何必要的符号名称转换。

另外，我怀疑您的包装子例程不应该驻留在模块中，否则您可能会在名称方面遇到额外的麻烦。但这只是对包装函数的限制，它必须从 R 中可见。

您可以在 DLL 中检查导出的名称(将 objdump -x your.so 发送到文件并查找导出的符号)。并在 R 中使用 is.loaded("your.symbol") 进行检查，加载DLL后。请注意，gfortran 通常会在名称后附加一个额外的下划线，而在调用 .Fortran 时不需要它。来自 R。如上所述，您可以使用 .C相反(但是，记住 Fortran 参数是通过引用传递的)。

为了检查你是否理解整个过程，我建议你用一个简单的例子来测试它，比如一个独特的子程序 mysub(x,y,z)只有 z=x+y .当这个运行时，您可以详细说明它以调用更复杂的例程。

编辑
当您将数组参数从 R 传递给 Fortran 时，不应使用假定形状或延迟形状数组，而只能使用假定大小的数组，即在 Fortran 77 中传递的通常数组。这是因为 R 只知道如何传递指向原始数据的指针，而假设形状和延迟形状需要更多信息，而 R 不知道执行此操作的数据结构。

例如，您可以这样做:

subroutine mysub(n, a)
    real :: a(n, n)
    ...
end subroutine

但这几乎肯定会失败:

subroutine mysub(a)
    real :: a(:, :)
    ...
end subroutine

此外，您不能将函数参数从 R 传递到 Fortran，因为这需要一个特殊的数据结构用于回调(在底层，R 是一个 Scheme 方言，它使用 S 表达式)。您可以通过 .C 在 C 中执行此操作或 .Call (参见 .Call 和 R Internals 的帮助)。