Reputation: 1355
Fortran subroutine gets stuck (but only when calling via R)
I have a Fortran subroutine that I converted from someone else's Fortran program. I would like to call this via R's .Fortran function. The subroutine works (instantly) when I call it from a Fortan program, but when I try to call it from R, nothing happens (In fact, R is still running this subroutine as I am typing this).
Here is the Fortran program (also containing the subroutine):
PROGRAM blep
integer a
real(4) b, c, d
b = 0.9
c = 0.1
d = 0.99
a = 0
call midpss(b, c, d, a)
4 format ('Calculated sample size is ',i6)
print 4, a
end
subroutine midpss(w, x, y, numbr)
c THIS IS TAKEN FROM "fosgate_original_working.f" AND THEN CONVERTED
real(8) probA,probB,part1,part2,part3,part4
real(8) totprA,totprB,factt, resp
integer numbr
c character resp
1 format ('Enter proportion ',$)
2 format ('Enter error limit ',$)
3 format ('Enter confidence level ',$)
4 format ('Calculated sample size is ',i6)
5 format ('Exact mid-P with ',f7.5,' 2-tail probability')
6 format ('Sorry, unable to mathmatically solve this problem.')
7 format ('Reported sample size is not accuarate.')
8 format ('Enter q to quit ',$)
9 format ('Actual limits for distribution ',f5.3,' - ',f5.3)
print *, 'Exact sampleroportions'
print *, 'Using Mid-P methods'
print *, 'Geoff Fosgate DVM PhD'
print *, 'College of Veterinary Medicine'
print *, 'Texas A&M University'
print *
10 prop1 = w
range = x
conlev = y
c Convert proportions less than 0.5 for algorithm
if (prop1 .lt. 0.5) then
prop = 1 - prop1
nprop = 1
else
prop = prop1
nprop = 0
end if
slimit = max ((prop - range) , 0.0001)
supper = min ((prop + range) , 0.9999)
c Probabilities cannot be calculated for p=0 and p=1
alpha = (1 - conlev)
if (alpha .gt. 1.0) go to 10
if (alpha .lt. 0.0) go to 10
if (prop .gt. 1.0) go to 10
if (prop .lt. 0.0) go to 10
numbr = (1 / (1 - prop)) - 1
c Define and initialize variables
c Note names of variables based on Fortran 77 rules
c Starting sample size is based on estimated proportion
c Resulting sample size must be large enough to obtain this proportion
100 numbr = numbr + 1
numx = (numbr * prop) + 0.001
c This is the number of binomial "successes" resulting in the proportion
if (numx .eq. numbr) go to 100
if (numx .lt. 1) go to 100
totprA = slimit**numbr
totprB = supper**numbr
do 130 loop1 = numx, (numbr - 1)
c Must initialize variables within loop
factt = 1.0
probA = 0.0
probB = 0.0
part1 = 0.0
part2 = 0.0
part3 = 0.0
part4 = 0.0
c Start loop to calculate factorial component of binomial probability
c Note that complete factorial calculations not necessary due to cancellations
do 110 loop2 = (loop1 + 1) , numbr
factt = factt * (loop2) / (numbr - (loop2 - 1))
110 continue
c Calculate probability for this particular number of successes
c Total probability is a running total
c Note that real variables must have high precision and be comprised
c of multiple bytes because factorial component can be very large
c and exponentiated component can be very small
c Program will fail if any component is recognized as zero or infinity
part1 = slimit**loop1
part2 = (1.0-slimit)**(numbr-loop1)
part3 = supper**loop1
part4 = (1.0-supper)**(numbr-loop1)
if (part1 .eq. 0.0) part1 = 1.0D-307
if (part2 .eq. 0.0) part2 = 1.0D-307
if (part3 .eq. 0.0) part3 = 1.0D-307
if (part4 .eq. 0.0) part4 = 1.0D-307
if (factt .gt. 1.0D308) factt = 1.0D308
probA = part1 * part2 * factt
probB = part3 * part4 * factt
if (loop1 .eq. numx) then
totprA = totprA + (0.5 * probA)
totprB = totprB + (0.5 * probB)
else
totprA = totprA + probA
totprB = totprB + probB
end if
c THIS IS ERROR HANDLING. INSTEAD OF PRINTING, SET NUMBR = -1
c *****************************************************************
if (probA .eq. 0.0) then
c print 6
c print 7
c print *
c go to 150
numbr = -1
end if
if (probB .eq. 0.0) then
c print 6
c print 7
c print *
c go to 150
numbr = -1
end if
c *****************************************************************
130 continue
140 if ((totprA + (1 - totprB)) .gt. alpha) go to 100
c go to beginning and increase sample size by 1 if have not
c reached specified level of confidence
c I.E. IF INPUT PROPORTION IS LESS THAN 0.5
c (I DONT THINK THIS IS NECESSARY -- IT JUST PRINTS THE RESULTS)
c150 if (nprop .eq. 1) then
c print 4,numbr
c print 9, (1-supper),(1-slimit)
c else
c print 4,numbr
c print 9, slimit,supper
c end if
c DO WE NEED THIS PART????
c *****************************************************************
c if (totprA+(1-totprB) .lt. alpha) print 5,(totprA+(1-totprB))
c print *
c print 8
c result = resp
c print *
c if (resp .ne. 'q') go to 10
c print *
c print *
998 return
999 end
(Sorry for the comments left over from my converting the original program to a subroutine).
The program is called midpss1_prog.f and the subroutine is called midpss1.f
I compile and call the program by doing the following:
C:\Users\panterasBox>gfortran midpss1_prog.f
C:\Users\panterasBox>a.exe
Exact sampleroportions
Using Mid-P methods
Geoff Fosgate DVM PhD
College of Veterinary Medicine
Texas A&M University
Calculated sample size is 80
C:\Users\panterasBox>
This is working just fine!
When I call the subroutine, I do the following:
In the command line, I call this:
C:\Users\panterasBox>R CMD SHLIB midpss1.f
gfortran -m64 -O2 -mtune=core2 -c midpss1.f -o midpss1.o
gcc -m64 -shared -s -static-libgcc -o midpss1.dll tmp.def midpss1.o -Ld:/RCompil
e/r-compiling/local/local320/lib/x64 -Ld:/RCompile/r-compiling/local/local320/li
b -lgfortran -LC:/Users/panterasBox/Documents/R/R-3.2.2/bin/x64 -lR
Then, I go into the R terminal and do this:
> setwd("C:/Users/panterasBox")
> dyn.load("midpss1.dll")
> is.loaded("midpss")
[1] TRUE
> .Fortran("midpss", w=as.numeric(0.9), x=as.numeric(0.1), y=as.numeric(0.90), numbr=as.integer(0))
And this last call to .Fortran never returns anything. It is just stuck...
Any help figuring out what is going on here would be greatly appreciated, thank you.
Upvotes: 1
Views: 226
Answers (1)
Reputation: 7395
R seems to be sending floating-point numbers to Fortran subroutines as double-precision, so we probably need to declare the corresponding dummy arguments accordingly. Because your program has no implicit none at the top of the subroutine, the dummy arguments w, x, and y are regarded implicitly as single-precision, making argument types inconsistent between R and Fortran (so resulting in hung-up). To fix this, simply declare them explicitly (here we assume that real(8) corresponds to double precision in R):
subroutine midpss(w, x, y, numbr)
real(8) :: w, x, y !<--- insert this line
!! double precision :: w, x, y !<--- or this line (but not both)
!! No need to modify the remaining part...
....
then we obtain the expected result (on Linux x86_64):
> .Fortran("midpss", w=as.numeric(0.9), x=as.numeric(0.1), y=as.numeric(0.99), numbr=as.integer(0))
Exact sampleroportions
Using Mid-P methods
Geoff Fosgate DVM PhD
College of Veterinary Medicine
Texas A&M University
$w
[1] 0.9
$x
[1] 0.1
$y
[1] 0.99
$numbr
[1] 80
Btw, this kind of problem may be avoided by using implicit none (as suggested repeatedly), because all the variables need to be declared explicitly, for example:
subroutine midpss (w, x, y, numbr)
implicit none
real(8) :: w, x, y
real(8) :: prop, prop1, range, conlev, slimit, supper, alpha
integer :: loop1, loop2, numx, nprop
...
Upvotes: 4