无法从 R 调用 Fortran 子例程

Question

我试图从 R 调用一个简单的 Fortran 子例程，但出了点问题。我设法编译了 Fortran 代码（我准确地说我只是这门语言的初学者）但是当我在 R 中调用子例程时失败发生了。

Fortran 代码是一个简单的子例程，它计算函数在（定义的）范围内每个点的值。然后将结果存储在大小为 2x100 000 的矩阵（fortran 中的数组）中。一行存储函数 f(x) 的值，另一行存储相应的变量 (x)。

Subroutine algo(n, tho, c, phi, ydata, eps, results)

    ! n : number of observations
    ! tho : number of steps in the range
    ! phi and c are given parameters 
    ! ydata : the vector of data
    ! eps : the iteration step
    ! results : array which stores the results

    IMPLICIT NONE

    INTEGER                                 :: n, i, j, tho
    DOUBLE PRECISION                        :: c, phi, sigma2, eps, ll 
    DOUBLE PRECISION, DIMENSION(1:n)        :: ydata 
    DOUBLE PRECISION, DIMENSION(1:tho)      :: vecteps
    DOUBLE PRECISION, DIMENSION(1:2,1:tho)  :: results
    DOUBLE PRECISION, PARAMETER             :: pi=acos(-1.d0)

    ! vecteps is the vector of all the x 
    vecteps(1)=0

    do i=1,tho
        vecteps(i)=vecteps(i)+eps
    end do

    do j=1,tho
        sigma2=vecteps(j)

        ll=-( (-(n-1)/2)*log(2*pi)-((n-1)/2)*log(sigma2)-
        sum((ydata(2:n)-c-phi*ydata(1:n-1))**2)/(2*sigma2) )

        results(1,j)=ll
        results(2,j)=vecteps(j)
    end do

end subroutine algo

然后是R的电话

y=rnorm(200,0,1)
y[1]=4/(1-0.6)
for(i in 2:length(y)){
  y[i]=4+0.6*y[i-1]+rnorm(1,0,1)
}

results=matrix(0, nrow=2, ncol=100000)

dyn.load("algo.dll")
    .Fortran('algo',n=as.integer(200), tho=as.integer(100000), c=as.double(4), phi=as.double(0.6), ydata=as.double(y), eps=as.double(0.0001), results=as.double(results) )

从 R 给出的数组 "results" 的 1000 个观察值中，绝大多数的数字都是相同的，只有少数 NaN ：

$results
   [1]  8.921136e+05  1.000000e-04  8.921136e+05  1.000000e-04  8.921136e+05  1.000000e-04
   [7]  8.921136e+05  1.000000e-04  8.921136e+05  1.000000e-04  8.921136e+05  1.000000e-04
  [13]           NaN           NaN  8.921136e+05  1.000000e-04  8.921136e+05  1.000000e-04

我希望他给我一个输出 [1] f(x1) 0.0001 f(x2) 0.0002 ...

对我来说，问题可能来自 3 个不同的问题

Fortran 代码中的数学错误
错误使用 matrix/arrays
R ( .Fortran(..., n=as... ) ) 的子程序调用

但是我没能解决。任何帮助将不胜感激。谢谢！

Answer 1

您可以使用 inline 包 -- 这是它的第一个示例：

R> x <- as.numeric(1:10)
R> n <- as.integer(10)
R> code <- "
+       integer i
+       do 1 i=1, n(1)
+     1 x(i) = x(i)**3
+ "
R> cubefn <- cfunction(signature(n="integer", x="numeric"), code, 
+                      convention=".Fortran")
R> print(cubefn)
An object of class 'CFunc'
function (n, x) 
.Primitive(".Fortran")(<pointer: 0x7fc8a8db7660>, n = as.integer(n), 
    x = as.double(x))
<environment: 0x8193e98>
code:
  1: 
  2:        SUBROUTINE file2e6b876b50a29 ( n, x )
  3:       INTEGER n(*)
  4:       DOUBLE PRECISION x(*)
  5: 
  6:       integer i
  7:       do 1 i=1, n(1)
  8:     1 x(i) = x(i)**3
  9: 
 10:       RETURN
 11:       END
 12: 
R> cubefn(n, x)$x
 [1]    1    8   27   64  125  216  343  512  729 1000
R>

我会用它来确保您的代码构建和运行，一旦完成，建议创建一个包。

无法从 R 调用 Fortran 子例程

Fail in calling a Fortran subroutine from R

fortran

r

call

subroutine