I tried to call a simple Fortran subroutine from R but something get wrong. I managed to compile the Fortran code (I precise that I'm only a beginner in this language) but the fail happens when I call the subroutine in R.
The Fortran code is a simple subroutine which compute the value of a fonction at each point of an (defined) range. Then the results are stored in a matrix (array in fortran) of size 2x100 000. One row stores the value of the fonction f(x) and the other row stores the corresponding variable (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
Then the call from 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) )
From the 1000 observations of the array "results" that R gives, the numbers are for the vast majority all the same with few 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
I would like that he gives me an output with [1] f(x1) 0.0001 f(x2) 0.0002 ...
To me, the problem may come from 3 different issues
but I fail to solve it. Any help would be greatly appreciated. Thanks!
You could use the inline package -- here is its first example:
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>
I would use this to make sure your code builds and runs, and once that is accomplished, suggest to create a package.