[R] call Fortran from R

Giacomo Santini giacomo.santini at unifi.it
Fri Sep 11 11:19:35 CEST 2009


Thanks! now it works !

Giacomo
 

Duncan Murdoch wrote:
> 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
>>
>
>

-- 
-------------------------------------------------------
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




More information about the R-help mailing list