[R] call Fortran from R

[Duncan Murdoch]

On 11/09/2009 4:39 AM, Giacomo Santini wrote:
> Dear R users,
> 
> 
> I have to call fortran program from within R (R 2.8.1 on ubuntu 8.10 
> machine).
> Suppose I have a fortran code like this (this is only a toy model, my 
> working model is far more complex, but input/output is similar)
> 
> 
>       DOUBLE PRECISION FUNCTION model(times, alfa, beta)
>       DOUBLE PRECISION alfa, beta, times
>       model=beta*sin(times)+alfa*cos(times)
>       END FUNCTION
> 
> which is saved as model.f.
> 
> I wrote a wapper like this (saved as wrapper.f)
> 
>        SUBROUTINE model_wrapper(times, alfa, beta, answer)
>        DOUBLE PRECISION times, alfa, beta, answer
>        EXTERNAL model
>        answer = model(times, alfa, beta)
>        END SUBROUTINE

It's been a long time since I worked in Fortran, but don't you need to 
declare that model returns a double precision value in the wrapper?  By 
default in ancient Fortran model would be an integer value.  Since you 
compile the two functions separately, the wrapper code won't see the 
declaration in model.f.

> 
> Then I compiled all this stuff
> 
> g77 -fno-second-underscore -c -fPIC model.f
> g77 -fno-second-underscore -c -fPIC wrapper.f
> g77 -fno-second-underscore -shared -o model.so model.o wrapper.o

I would use R CMD SHLIB model.f wrapper.f

to be sure the options match what R is expecting.

Duncan Murdoch

> 
> 
>  From within R
> 
> dyn.load('model.so')
> model <- function(times, alfa, beta) {
>       returned_data = .Fortran('model_wrapper', times=as.double(times), 
> alfa=as.double(alfa),
>       beta=as.double(beta), result=double(1))
>       return(returned_data) }
> 
> # example run
> test_1<-model(1.0,0.2,0.3)
> 
> which gives
> 
> test_1$times
> [1] 1.0
> $alfa
> [1] 0.2
> $beta
> [1] 0.3
> $result
> [1] 147456887
> 
> where $result is clearly wrong.
> 
> I suppose I made some mistake with the handling of  data types, but I am 
> not able to figure out  where.
> 
> Can someboby help me?
> 
> 
> Giacomo
>