2016-01-06 152 views
1

我有一個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永遠不會返回任何東西。它只是卡住...

任何幫助搞清楚這裏發生了什麼將不勝感激,謝謝。

+1

能不能請你看在m儘可能減少你的例子,同時仍然保持這個問題?這樣做不僅有助於你自己解決問題的努力,而且對我們來說閱讀起來也不那麼重要。 [如果沒有其他問題,請刪除那些大量註釋掉的代碼。]請參見[mcve]。 – francescalus

+0

在前面的問題中,答案建議你使用'implicit none'。你在這裏沒有,並且一些虛擬變量(至少)是隱式聲明的。我們怎樣才能確定「真實(4)」與默認真實相同? – francescalus

+0

你可能會對http://www.inside-r.org/node/54298感興趣... –

回答

4

R似乎是將雙精度浮點數發送到Fortran子程序,所以我們可能需要相應地聲明相應的僞參數。由於您的程序在子程序頂部沒有implicit none,因此僞參數w,xy被隱式視爲單精度,導致R和Fortran之間的參數類型不一致(因此導致掛起)。爲了解決這個問題,簡單地宣佈他們明確(這裏我們假設real(8)對應於雙精度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... 
      .... 

那麼我們可以得到預期的結果(在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 

順便說一句,這樣的問題可以通過使用implicit none(如多次建議)是可以避免的,因爲所有的變量都需要顯式聲明,例如:

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 
     ... 
+1

似乎有不同的頁面解釋瞭如何將浮點數傳遞給C/Fortran(例如https://stat.ethz.ch/R-manual/R-devel/library/base/html/Foreign.html )和一些選項似乎有必要將它們作爲單精度傳遞。 – roygvib