converting a fortran program to a subroutine to be called from R

Question

I previously posted a question about a question here about trying to reproduce sample sizes from a statistics paper. The author has provided me the code which was used in the paper which is a Fortran program. I would like to call this code using the R package my team at work uses. To do this, I need to convert this Fortran program to a subroutine. This is my first experience with Fortran and I seem to have broken the program somehow and get wrong numbers that do not change when I change the input. I should also note that I had to make a few small changes from the original code that the author provided to get the program to work.

The "original" program that produces correct results when compiled and then ran:

    PROGRAM Exact_Mid_P_binomial_sample_size
          real(8) probA,probB,part1,part2,part3,part4
          real(8) totprA,totprB,factt, resp
!       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      print *
        print 1
          read *, prop1
        print 2
          read *,range
        print 3
          read *,conlev
        print *
!           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)
!           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
!           Define and initialize variables
!             Note names of variables based on Fortran 77 rules
!             Starting sample size is based on estimated proportion
!             Resulting sample size must be large enough to obtain this proportion
100     numbr = numbr + 1
        numx = (numbr * prop) + 0.001
!             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)
!               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
!                Start loop to calculate factorial component of binomial probability
!                Note that complete factorial calculations not necessary due to cancellations
            do 110 loop2 = (loop1 + 1) , numbr
               factt = factt * (loop2) / (numbr - (loop2 - 1))
110         continue
!                Calculate probability for this particular number of successes
!                Total probability is a running total
!                Note that real variables must have high precision and be comprised
!                of multiple bytes because factorial component can be very large
!                and exponentiated component can be very small
!                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
            if (probA .eq. 0.0) then 
                    print 6
                    print 7
                    print *
                    go to 150
                end if
            if (probB .eq. 0.0) then
                    print 6
                    print 7
                    print *
                    go to 150
                end if
130       continue
140     if ((totprA + (1 - totprB)) .gt. alpha) go to 100
!             go to beginning and increase sample size by 1 if have not
!             reached specified level of confidence
150         if (nprop .eq. 1) then
                print 4,numbr
                print 9, (1-supper),(1-slimit)
            else
                print 4,numbr
                print 9, slimit,supper
            end if
            if (totprA+(1-totprB) .lt. alpha) print 5,(totprA+(1-totprB))
            print *
            print 8
            result = resp
!           print *
!       if (resp .ne. 'q') go to 10
      print *
      print *
999   end

The program above turned into a subroutine:

       subroutine midpss(prop1, range, conlev, numbr)

       integer numbr, nprop
       real(8) prop1, range, conlev, prop, slimit, supper
       real(8) probA,probB,part1,part2,part3,part4,factt
       real(8) totprA,totprB, resp
c         character  resp

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)

           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

130      continue
140      if ((totprA + (1 - totprB)) .gt. alpha) go to 100


       return
       end

To compile the original program, I run this in the command line:

gfortran midpSS_original.f

To compile the subroutine, I run this in the command line:

R CMD SHLIB midpSS_subroutine.f

and then from the R console, I run this:

> dyn.load("midpSS_subroutine.dll")
> is.loaded("midpss")
[1] TRUE
> .Fortran("midpss", prop1=as.numeric(0.9), range=as.numeric(0.1), conlev=as.numeric(0.90), numbr=as.integer(0)) # numbr should be 29
$prop1
[1] 0.9

$range
[1] 0.1

$conlev
[1] 0.9

$numbr
[1] 2091

> .Fortran("midpss", prop1=as.numeric(0.9), range=as.numeric(0.1), conlev=as.numeric(0.95), numbr=as.integer(0)) # numbr should be 47
$prop1
[1] 0.9

$range
[1] 0.1

$conlev
[1] 0.95

$numbr
[1] 2091

In the first call to the subroutine, the numbr should be 29. In the second one, numbr should be 47. The results are correct when compiling the "original" fortran program and running the same parameters. I cannot figure out what is going on here. Any help would be greatly appreciated.

EhsanSafari · Accepted Answer

The following block is checking the input data:

    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

This block makes the code go to the entry position to read the input data again. You might get wrong results because your input data into the subroutine might be out of range. Please check this first!

Another missing block is

if (probA .eq. 0.0) then 
  print 6
  print 7
  print *
 go to 150
end if

if (probB .eq. 0.0) then
 print 6
 print 7
 print *
 go to 150
end if

and the label 150 as well. Are you sure that omitting these two blocks have not changed the accuracy and functionality of the original code?

Converting a program into a subroutine is very simple. The main steps are as follows:

Convert the syntax program into subroutine
Convert the syntax end program into end subroutine
Put a parenthesis ‘()’ after the name of your subroutine
Put all the input and output of the program inside that parenthesis ‘()’

It is also a good practice to write implicit none after the program or subroutine

I convert the program into a subroutine of my own. Please note that the result = resp at the end of the subroutine might be wrong because as far as I know the command ‘result’ is for a function not a subroutine.

In the subroutine you have to specify the output. Since I don’t know which variable is the output, I have not included that variable inside the parenthesis ‘()’ after the subroutine name. please write it there. My subroutine is like the following:

!    PROGRAM Exact_Mid_P_binomial_sample_size
    subroutine midpss (prop1, range, conlev, output)
    implicit none
          real(8) probA,probB,part1,part2,part3,part4
          real(8) totprA,totprB,factt, resp
          integer numbr, nprop,loop1,loop2
          real(8) prop1,prop,range,conlev,slimit,supper,alpha,numx,output

!       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      print *
!        print 1
!          read *, prop1
!        print 2
!         read *, range
!       print 3
!         read *, conlev
!       print *
!           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)
!           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
!           Define and initialize variables
!             Note names of variables based on Fortran 77 rules
!             Starting sample size is based on estimated proportion
!             Resulting sample size must be large enough to obtain this proportion
100     numbr = numbr + 1
        numx = (numbr * prop) + 0.001
!             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)
!               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
!                Start loop to calculate factorial component of binomial probability
!                Note that complete factorial calculations not necessary due to cancellations
            do 110 loop2 = (loop1 + 1) , numbr
               factt = factt * (loop2) / (numbr - (loop2 - 1))
110         continue
!                Calculate probability for this particular number of successes
!                Total probability is a running total
!                Note that real variables must have high precision and be comprised
!                of multiple bytes because factorial component can be very large
!                and exponentiated component can be very small
!                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
            if (probA .eq. 0.0) then 
                    print 6
                    print 7
                    print *
                    go to 150
                end if
            if (probB .eq. 0.0) then
                    print 6
                    print 7
                    print *
                    go to 150
                end if
130       continue
140     if ((totprA + (1 - totprB)) .gt. alpha) go to 100
!             go to beginning and increase sample size by 1 if have not
!             reached specified level of confidence
150         if (nprop .eq. 1) then
                print 4,numbr
                print 9, (1-supper),(1-slimit)
            else
                print 4,numbr
                print 9, slimit,supper
            end if
            if (totprA+(1-totprB) .lt. alpha) print 5,(totprA+(1-totprB))
            print *
            print 8
            !result = resp
!           print *
!       if (resp .ne. 'q') go to 10
!      output=***     !write the output variable instead of ***
      print *
      print *
999   end subroutine midpss

converting a fortran program to a subroutine to be called from R

Answers (1)

Related Questions