Calling a Fortran subroutine in R

Question

EDITED

I am trying to integrate R and Fortran, calling a Fortran subroutine in R. I'm new to Fortran, but I'm not new to R or to programming.

I wrote a subroutine in Fortran to calculate the sinc function, and I want to call it in R. This is the subroutine:

SUBROUTINE SINC(X,Y)
    IMPLICIT NONE
    DOUBLE PRECISION , INTENT(IN) :: X
    DOUBLE PRECISION , PARAMETER :: PI = 4.D0*ATAN(1.0)
    DOUBLE PRECISION , INTENT(OUT) :: Y
    
    IF(X == 0) THEN
        Y = 1
    ELSE
        Y = (SIN(PI*X))/(PI*X)
    END IF
END SUBROUTINE SINC

I wrote a Fortran code to test it, and it runs nicely.

PROGRAM DEMO
DOUBLE PRECISION A, B
PRINT *, "WHAT IS THE NUMBER YOU WANT TO SINC?"
READ *, A
CALL SINC(A, B)
PRINT *, "THE SINC IS", B
END PROGRAM DEMO

SUBROUTINE SINC(X,Y)
    IMPLICIT NONE
    DOUBLE PRECISION , INTENT(IN) :: X
    DOUBLE PRECISION , PARAMETER :: PI = 4.D0*ATAN(1.0)
    DOUBLE PRECISION , INTENT(OUT) :: Y
    
    IF(X == 0) THEN
        PRINT *, "X IS NULL"
        Y = 1
    ELSE
        PRINT *, "X IS NON NULL"
        Y = (SIN(PI*X))/(PI*X)
    END IF
END SUBROUTINE SINC

I compiled the subroutine (not the program) to a .dll file using the following command

R CMD SHLIB sinc.f95

If I defined my variables with REAL instead of DOUBLE PRECISION, when I ran it passing the parameter for Y, it basically returned me the parameter for Y I passed and I don't know why.

dyn.load("D:\\Fortran\\sinc.dll")
.Fortran("sinc", X=as.double(0), Y=as.double(1))

$X
[1] 0

$Y
[1] 1

The subroutine works perfectly fine for values that are non-integers, such as 1.1 or 2.5. However, when I try computing it for integer values, such as 1,2,3,4... the subroutine returns me -2.7827534378485793E-008, whereas in R i get 3.898043e-17. I know both are zeroes, but the precision is different in that case.

What am I doing wrong? How can I fix this to call my subroutine in R? I'm using gfortran as my compiler.

do you need `DOUBLE PRECISION` rather than `REAL`? Looking at https://cran.r-project.org/doc/manuals/r-release/R-exts.html#Interface-functions-_002eC-and-_002eFortran ... — Ben Bolker, Jun 15 '22 at 03:25
Thanks! Oddly it worked out! But I don't know why changing it to ```DOUBLE PRECISION``` would make my code work. — YetAnotherUsr, Jun 15 '22 at 03:53
You know that it works using `double precision` on the Fortran side instead of `real`, but `real` means something very different from `double precision` and one wouldn't expect things to work using `integer` or `character` instead, just as there's no reason for `real` to be right. Given that, can you explain what explanation you require for why you see the behaviour with `real`? — francescalus, Jun 15 '22 at 08:47
@francescalus thanks, I was reading more about Fortran declarations, I had misinterpreted the meaning of ```REAL```. I believe it is because ```REAL``` refers to a single precision floating point. Still, I don't know how to make my code return the same precision I have in R (2.7E-8 in FORTRAN as compared to 3.8E-17 in R). — YetAnotherUsr, Jun 15 '22 at 10:22
Your `PI = 4.D0*ATAN(1.0)` simply isn't as precise. You should be using `PI = 4.D0*ATAN(1.0D0)` or similar. The first has `atan` as single precision instead of the double precision of the second. — francescalus, Jun 15 '22 at 10:55
Between the two linked questions you should hopefully get a clear understanding of why the imprecision happens. If it isn't so clear, do let us know and we can try to clarify further — francescalus, Jun 15 '22 at 11:09
By the way, if you've got your answer from the comments you're welcome (encouraged) to write up the answer to your own question and post it. — Ben Bolker, Jun 15 '22 at 13:35

Calling a Fortran subroutine in R

0 Answers0