Fortran子程序卡住了（但仅在通过R调用时） [英] Fortran subroutine gets stuck (but only when calling via R)
本文介绍了Fortran子程序卡住了（但仅在通过R调用时）的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我有一个Fortran子程序，我从别人的Fortran程序转换而来。我想通过R的 .Fortran 函数调用它。当我从Fortan程序中调用子程序时，子程序立即工作，但是当我尝试从R中调用它时，什么事都没有发生（事实上，当我输入时，R仍然运行这个子程序）。
以下是Fortran程序（也包含子程序）：
PROGRAM blep
整数a
实数（4）b，c，d
b = 0.9
c = 0.1
d = 0.99
a = 0
call midpss（b，c，d，a）
4格式（'计算样本量为'，i6）
print 4，a
end
子程序midpss（w，x，y，numbr）
c这是从fosgate_original_working.f中获得的，然后被转换为
real（8）probA，probB，part1，part2，part3，part4
real（8）totprA，totprB，factt，resp
整数numbr
c字符resp
1格式（'输入比例'，$）
'格式（'输入错误限制'，$）
3格式（'输入置信水平'，$）
4格式（'计算样本大小为'，i6）
5格式'精确的mid-P with'，f7.5，'2-tail概率'）
6格式（'抱歉，无法用数学方法解决这个问题。'）
7格式（'Reported sample size is '）
8格式（'输入q退出'，$）
9格式（'分配的实际限制'，f5.3，' - '，f5.3）
print *，'精确采样'
print *，'使用Mid-P方法'
print *，'Geoff Fosgate DVM PhD'
print *，'兽医学院'
print *，'Texas A& M University'
print *
10 prop1 = w
range = x
conlev = y
c换算比例小于0.5算法
if（prop1 .lt。 0.5）则
prop = 1 - prop1
nprop = 1
else
prop = prop1
nprop = 0
如果
结束slimit =最大（（道具范围），0.0001）
supper = min（（道具+范围），0.9999）
c对于p = 0和p = 1
alpha =（1 - conlev）
if（alpha .gt。1.0）转到10
if（alpha .lt。0.0）转到10
if（prop .gt。1.0）转到10
if（prop .lt。0.0）转到10
numbr =（1 /（1-prop）） - 1
c定义和初始化变量
c注意基于Fortran的变量名称77规则
c起始样本大小基于估计的比例
c所得到的样本量必须足够大以获得该比例
100 numbr = numbr + 1
numx =（numbr * prop）+ 0.001
c这是如果（numx .eq。numbr）达到100时导致比例
的二项成功的数目
if（numx .lt。 1）转到100
totprA = slimit ** numbr
totprB = supper ** numbr
do 130 loop1 = numx，（numbr - 1）
c必须初始化循环内的变量
factt = 1.0
probA = 0.0
probB = 0.0
part1 = 0.0
part2 = 0.0
part3 = 0.0
part4 = 0.0
c开始循环计算二项概率的阶乘成分
c请注意，由于取消
而不需要进行完整的阶乘计算
do 110 loop2 =（loop1 + 1），numbr
factt = factt *（loop2 ）/（numbr - （loop2 - 1））
110继续
c计算这个特定数目成功的概率
c总概率是一个运行总数
c请注意，实际变量必须高精度和构成
c的多个字节，因为阶乘组件可能非常大
c和指数化组件可能非常小
c程序w如果任何组件被识别为零或无穷大，则会失败
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
如果
c则结束这是错误处理。不打印，设置NUMBR = -1
c ************************************* ****************************
if（probA .eq。0.0）then
c print 6
c print 7
c print *
c转到150
numbr = -1
如果
结束if（probB .eq。0.0）then
c print 6
c print 7
c print *
c转到150
numbr = -1
如果
c为************** ************************************************** *
130继续
140如果（（totprA +（1 - totprB））.gt。alpha）转到100
c转到开头，如果没有$ b $，则将样本大小增加1 bc达到了指定的信心水平
c IE如果输入比例小于0.5
c（我不认为这是必要的 - 它只是打印结果）
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我们是否需要这部分？
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
准确样本量
使用Mid-P方法
Geoff Fosgate DVM博士
兽医学院
德克萨斯农工大学
计算d样本大小为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-编译/ 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 时不会返回任何内容。它只是卡住了......
任何帮助了解这里发生的事情将不胜感激，谢谢。
隐式无，所以伪参数
w ， x 和 y 被隐式视为单精度，导致R和Fortran之间的参数类型不一致（导致挂起）。为了解决这个问题，只需明确声明它们（这里我们假设 real（8）对应于R中的双精度）： 子程序midpss（w，x，y，numbr）
real（8）:: w，x，y！< --- insert this line
！ double precision :: w，x，y！<---或这行（但不是两个）
!!无需修改剩余部分...
....
然后我们获得预期的结果（在Linux x86_64上）：
> .fortran（midpss，w = as.numeric（0.9），x = as.numeric（0.1），y = as.numeric（0.99），numbr = as.integer（0））
精确的抽样比例
使用Mid-P方法
Geoff Fosgate DVM博士
兽医学院
德克萨斯农工大学
$ w
[1] 0.9
$ x
[1] 0.1
$ y
[1] 0.99
$ numbr
[1] 80
顺便说一下，这种问题可以通过使用 implicit none （重复建议），因为所有变量都需要显式声明，例如：
<$ c （8）:: w，x，y
real（8）:: prop，prop1，范围（$）$ c>子程序midpss（w，x，y，numbr）
隐式无
real ，conlev，slimit，supper，alpha
integer :: loop1，loop2，numx，nprop
...
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.
解决方案
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
...
这篇关于Fortran子程序卡住了（但仅在通过R调用时）的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
