我有一個Fortran子程序,我從別人的Fortran程序轉換而來。我想通過R的.Fortran
函數來調用它。當我從Fortan程序調用子程序時,該子程序立即工作,但是當我嘗試從R中調用它時,沒有任何反應(事實上,當我輸入時,R仍然運行這個子程序)。Fortran子程序卡住了(但僅在通過R調用時)
這裏是Fortran程序(也包含子程序):
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
(對不起,我從我原來的程序轉換到子程序留下的評論)。
該計劃被稱爲midpss1_prog.f
和調用子程序midpss1.f
我編譯和執行調用程序如下:
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>
這是工作就好了!
當我調用子程序,我執行以下操作:
在命令行,我把這種:
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
然後,我進入R端和做到這一點:
> 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))
而這最後一次撥打.Fortran
永遠不會返回任何東西。它只是卡住...
任何幫助搞清楚這裏發生了什麼將不勝感激,謝謝。
能不能請你看在m儘可能減少你的例子,同時仍然保持這個問題?這樣做不僅有助於你自己解決問題的努力,而且對我們來說閱讀起來也不那麼重要。 [如果沒有其他問題,請刪除那些大量註釋掉的代碼。]請參見[mcve]。 – francescalus
在前面的問題中,答案建議你使用'implicit none'。你在這裏沒有,並且一些虛擬變量(至少)是隱式聲明的。我們怎樣才能確定「真實(4)」與默認真實相同? – francescalus
你可能會對http://www.inside-r.org/node/54298感興趣... –