从 R 调用 Fortran 子例程失败答案

【问题标题】：Fail in calling a Fortran subroutine from R从 R 调用 Fortran 子例程失败
【发布时间】：2017-03-14 21:26:51
【问题描述】：

我试图从 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 给出的对数组“结果”的 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 代码中的数学错误
矩阵/数组的错误使用
来自 R 的子例程调用 (.Fortran(..., n=as...))

但我无法解决它。任何帮助将不胜感激。谢谢！

【问题讨论】：

我刚刚编辑了最初的帖子以使结果更加清晰。
我没有计算机来检查这个，但在第一个循环中vecteps(i)=vecteps(i)+eps 看起来非常可疑：到目前为止，只有vecteps(1) 被定义。也许你的意思是vecteps=0（或vecteps(:)=0），尽管之后的循环对我来说都没有多大意义？
通过这个循环，我尝试生成一个类似于由这个 R 代码创建的向量的向量：“vecteps=seq(from=0.0001, to=10, by=0.0001)”.. . 但你是对的，它似乎在这里失败了。
This question's answers 提供了如何在 Fortran 中创建/构造这样一个数组的示例。

标签： r fortran call subroutine

【解决方案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>

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

【讨论】：