R help - Sep 2009 - call Fortran from R

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

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


 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

-- 
-------------------------------------------------------
Giacomo Santini PhD
Dipartimento di Biologia Evoluzionistica "Leo Pardi"
Universita' degli Studi di Firenze
Via Romana 17
I-50125 Firenze
Italy

Tel: +39 055 2288288 (DBE) - +39 0574 447727 (CESPRO)
Fax: +39 055 2288289
www.dbe.unifi.it/santini

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
>