Recursion – Fortran finally gets it!

Fortran gave recursion a wide berth until Fortran 90. More than 25 years after ALGOL first allowed recursion, Fortran was finally allowing. Modern Fortran is one of the few languages where the function must be explicitly defined recursive. The prefix recursive must be used before the keyword function.

recursive function fibonacci (n) result (fib_result)
   integer, intent (in) :: n
   integer :: fib_result
   if (n <= 2) then
      fib_result = 1
   else
      fib_result = fibonacci (n-1) + fibonacci (n-2)
   end if
end function fibonacci

In Fortran 2018, the specs changed again. Here procedures without an explicit recursive attribute behave as if recursive is specified. F2018 provides a new attribute non_recursive, which can be used to mark a procedure that may not be called recursively.

Ditching a tricky arithmetic if

A three-way branch known as the arithmetic-if, can be tricky to remove sometimes. Consider the code below:

    do 6 i = 1,3
       if (dif - c(i)) 7,7,6
  6 continue
    i = 4
  7 ffm = b(i)*exp(a(i)*dif)

The arithmetic-if works by branching based on evaluating the expression: the result is either less-than zero, zero, or greater-than zero. Here if the value of dif-c(i) is ≤0 then the program branches to label 7, otherwise if >0 it branches to label 6, which is the next iteration of the loop. Basically the loop exits when dif is less than or equal to c(i), and the value of i is used in the calculation of ffm. If the condition is never met, then the loop ends, and the value of i is set to 4, and ffm is calculated. 

The first thing to do when re-engineering this, is to update the do-loop to F90+ standards, and reconfigure how the if statement works, essentially transforming it to a regular logical if, that uses goto to branch to label 7 if the difference is ≤ 0.

    do i = 1,3
       if (dif - c(i) .le. 0) then
          go to 7
       end if
    end do
    i = 4
  7 ffm = b(i)*exp(a(i)*dif)

The problem here is removing the goto and still being able to skip the statement i=4. The trick is to incorporate i=4 into the loop. For example:

do i = 1,4
   if (i == 4) then
      exit
   else if (dif <= c(i)) then
      exit
   end if
end do
ffm = b(i)*exp(a(i)*dif)

The loop runs through 4 times, when the index i becomes 4, the loop exits, and ffm is calculated. Otherwise for i is 1…3, the loop only exits if dif is less than or equal to c(i). This is only one approach, there are others of course.

2D strings in Fortran

So elsewhere we have talked about how in Fortran, strings are different to character arrays, so what about 2D strings? Strings that hold multiple words, such as in a dictionary. Here we are effectively looking at an array of strings. In early versions of Fortran, they were created in the following manner:

     CHARACTER*20 DICT(12)
     DATA DICT/'bantha','dewback','ronto','womprat','woodoo','worrt',
   M           'sarlacc','rancor','kraytdragon','tauntaun','wampa','galoomp'/

Here DICT was an array of 20 strings, each of which could be 12 characters in length. This works fine, except the DATA statement for specifying initial values is now somewhat obsolete. Now a 2D string can be created in the following manner in F95:

character(len=20),dimension(12) :: dict

This create a character array named dict in which there are 12 strings each which can contain a maximum of 20 characters. So values can be assigned to dict in the following manner:

dict = [character(20) :: 'bantha','dewback','ronto','womprat', &
                         'woodoo','worrt','sarlacc','rancor', &
                         'kraytdragon','tauntaun','wampa','galoomp']

Note the use of a type specification in the constructor – character(20). If this is omitted, then each string must be of the same size. For example, if the above constructor were to look like:

dict = ['bantha','dewback','ronto','womprat', &
        'woodoo','worrt','sarlacc','rancor', &
        'kraytdragon','tauntaun','wampa','galoomp']

Then when compiled, the following error would appear:

dict = ['bantha','dewback','ronto','womprat', &
                1
Error: Different CHARACTER lengths (6/7) in array constructor at (1)

However the error would not occur if the constructor looked like this (each string is padded with blanks to the size of the maximum string in the group:

dict = ['bantha     ','dewback    ','ronto      ','womprat    ', &
        'woodoo     ','worrt      ','sarlacc    ','rancor     ', &
        'kraytdragon','tauntaun   ','wampa      ','galoomp    ']

Defensive programming: Validating input in C and Fortran

As mentioned previously, defensive programming often involves validating user input. In many languages this involves checking what the user has entered to make sure it is valid within the context of what is being stored. If a program expects an integer, will an integer be entered? How will invalid data affect the programs behaviour? To look at two differing perspectives of validation input, let’s consider C and Fortran.

In C, the scanf() function does return a value which indicates something – the total number of inputs scanned successfully. So if a scanf() reads in one variable, and returns the value 1, then all is good. The problem then lies in creating a piece of code that actually works. Consider a simple loop which reads in an integer into the variable a. If the user enters 3, the program functions properly, however if the user enters 3.7, the loop executes infinitely. Why? Because of C’s wonderful issues with its input stream.

while (scanf("%d",&a) != 1)
   printf("invalid input\n");

The better way in C is to create a function to read the input as a string, and parse the string to determine if it is a valid integer, returning the integer. Below is such a function getInt(), which basically checks each element of the string to make sure it is a digit.

int getInt(){
   int valid, len, i;
   char s[100];
   valid = 0;
   while (valid == 0){
      fgets(s, sizeof(s), stdin);
      len = strlen(s);
      while (len > 0 && isspace(s[len-1]))
         len = len - 1;
      if (len > 0){
         valid = 1;
         for (i=0; i<len; i=i+1)
            if (!isdigit(s[i])){
               valid = 0;
               break;
            }
      }
      if (valid == 1)
         return atoi(s);
      else
         printf("invalid input\n");
   }
}

Of course you could create a library of such functions for integers, and floats etc. which is fun, but should you have to? Other languages have better built-in functionality to deal with invalid input. Fortran has a functionality in the read() function which returns a validity value – iostat. Basically if iostat returns a value of 0, then read() was executed flawlessly and all variables have received their valid input values. Below is a much simpler version of getInt() in Fortran, without the need to parse a string!

integer function getInt()
   integer :: n, istat
   istat = 1
   chklp: do
      read(*,*,iostat=istat) n
      if (istat == 0) then
         exit chklp
      else
         write(*,*) "invalid input"
      end if
   end do chklp
   getInt = n
end function getInt

Fortran Reengineering: entombing legacy subprograms

A couple of years ago, I termed the reengineering concept of entombing. This basically takes legacy subprograms and encapsulates the code within a modern framework. The interface block can be very useful in this process, so we will term this interface entombing. Consider the following example of a F77 subroutine. It calculates the average thrust of a rocket engine, given the weight of the rocket, the thrust (F) of the engines, and the velocity (V) of the exhaust gases. Below is the pseudocode for the F77 function SUBACL.

• Compute the mass lost after t seconds: m_loss = t (F/−V)
• Compute the net mass remaining after t seconds: m_{t (remaining}) = (W/gravity) − m_loss
• Compute the weight at W_t: W_t = m_{t (remaining}) gravity
• Compute the unbalanced force that causes acceleration: F_t=F−W_t
• Compute the acceleration: a = F_t / m_t

Here is Fortran 77 function (subacl.for).

      REAL FUNCTION SUBACL(RCKTWT,THRUST,VELGAS,TIME)
      REAL RCKTWT,THRUST,VELGAS,TIME,FTUNB,MT,UBWT
      REAL FUNFOR, FUNACL
      REAL GRAVTY
      PARAMETER(GRAVTY=32.0)
      FUNFOR(THRUST,VELGAS) = THRUST / (-VELGAS)
      FUNACL(FTUNB,MT) = FTUNB/MT
      MT = RCKTWT / GRAVTY - FUNFOR(THRUST,VELGAS) * TIME
      UBWT = MT * GRAVTY
      FTUNB = THRUST - UBWT
      SUBACL = FUNACL(FTUNB,MT)
      RETURN
      END

Now we may not want to change anything here, for whatever reason. Note that the variables are all explicitly declared, and the function retains a number of legacy structures. Its functionality could be accessed by a Fortran 95 program which specifies it using an interface block. This basically adds a modern layer of interface between the new program, and the legacy function.

program entomb1
   implicit none

   interface
      real function subacl(rcktwt,thrust,velgas,time)
         real, intent(in) :: rcktwt, thrust, velgas, time
      end function subacl
   end interface

   real :: weight, thrust, velgas, time

   write(*,*) 'Input data: '
   write(*,*) 'Weight of rocket (lbs): '
   read(*,*) weight
   write(*,*) 'Thrust of the engine (lbs): '
   read(*,*) thrust
   write(*,*) 'Velocity of exhaust gas (ft/sec): '
   read(*,*) velgas
   write(*,*) 'Duration thrust (sec): '
   read(*,*) time
   write(*,10) weight, thrust, velgas
10 format(/' ',4x,'Output:' &
          /' ',8x,'Weight of rocket:',t42,f10.2,' lbs' &
          /' ',8x,'Thrust of engines:',t42,f10.2,' lbs' &
          /' ',8x,'Velocity of exhaust:',t42,f10.2,' ft/sec')
   write(*,12) subacl(weight,thrust,velgas,time)
12 format(/' ',8x,'Acceleration of the vehicle:',t42,f10.2,' ft/sec*sec')

end program entomb1

Now we can test the program:

 Input data:
 Weight of rocket (lbs):
10000
 Thrust of the engine (lbs):
80000
 Velocity of exhaust gas (ft/sec):
50000
 Duration thrust (sec):
15
     Output:
         Weight of rocket:                 10000.00 lbs
         Thrust of engines:                80000.00 lbs
         Velocity of exhaust:              50000.00 ft/sec

         Acceleration of the vehicle:        205.74 ft/sec*sec

Also note how format is used here to format the output (the output of the input is somewhat redundant, but I left it here to show the formatting side of things. Especially make note of the / to imply a line break, and t to tab (to column 42).

Coding Fortran: contains (modularity)

Using the contains statement means that the subprograms become part of the main program, i.e. the subprograms are internal, essentially “nested” within the main program.

program subroutine3
implicit none
! solves a quadratic equation, based on the user input.

   real :: p, q, r, root1, root2
   integer :: ifail=0
   logical :: ok=.true.

   call getcoeffs(p,q,r,ok)
   if (ok) then
      call solve(p,q,r,root1,root2,ifail)
      if (ifail == 1) then
         print*,'Complex roots, calculation abandoned'
      else
         print*,'Roots = ',root1,' ',root2
      endif
   else
         print*,'Error in data input.'
   endif

contains

subroutine getcoeffs(a,b,c,ok)
   implicit none
   real, intent(out) :: a,b,c
   logical, intent(out) :: ok
   integer :: io_status=0

   print*,'Enter the coefficients a, b and c'
   read(unit=*,fmt=*,iostat=io_status) a,b,c
   if (io_status == 0) then
      ok = .true.
   else
      ok = .false.
   endif
end subroutine getcoeffs

subroutine solve(x,y,z,root1,root2,ifail)
   implicit none
   real, intent(in) :: x,y,z
   real, intent(out) :: root1,root2
   integer, intent(inout) :: ifail
   real :: term,a2

   term = y*y - 4.*x*z
   a2 = x*2.0
   ! if term < 0, roots are complex
   if (term < 0.0) then
      ifail = 1
   else
      term = sqrt(term)
      root1 = (-y+term)/a2
      root2 = (-y-term)/a2
   endif
end subroutine solve

end program subroutine3

Coding Fortran: Interface blocks (modularity)

Interface blocks provide a means of doing type checking between the calling subprogram, and the called subprogram. Incorrect types between the calling arguments, and the subprogram parameters is a common source of problems. This doesn’t happen as much in languages such as C, because if a real parameter is given as an argument to a subprogram which expects an integer it is merely truncated. Why use an interface block? Largely because it is good programming practice to include explicit interfaces, even if they are not required.

program subroutine2
   implicit none

   interface
   subroutine getcoeffs(a,b,c,ok)
      implicit none
      real, intent(out) :: a,b,c
      logical, intent(out) :: ok
   end subroutine getcoeffs

   subroutine solve(x,y,z,root1,root2,ifail)
      implicit none
      real, intent(in) :: x,y,z
      real, intent(out) :: root1,root2
      integer, intent(inout) :: ifail
   end subroutine solve
   end interface

   real :: p, q, r, root1, root2
   integer :: ifail=0
   logical :: ok=.true.

   call getcoeffs(p,q,r,ok)
   if (ok) then
      call solve(p,q,r,root1,root2,ifail)
      if (ifail == 1) then
         print *,'Complex roots, calculation abandoned.'
      else
         print *,'Roots = ',root1,' ',root2
      endif
   else
      print*,'Error in data input.'
   endif
end program subroutine2

subroutine getcoeffs(a,b,c,ok)
   implicit none
   real, intent(out) :: a,b,c
   logical, intent(out) :: ok
   integer :: io_status=0

   print*,'Enter the coefficients a, b and c'
   read(unit=*,fmt=*,iostat=io_status) a,b,c
   if (io_status == 0) then
      ok=.true.
   else
      ok=.false.
   endif
end subroutine getcoeffs

subroutine solve(x,y,z,root1,root2,ifail)
   implicit none
   real, intent(in) :: x,y,z
   real, intent(out) :: root1,root2
   integer , intent(inout) :: ifail

   real :: term, a2
   term = y*y - 4.*x*z
   a2 = x*2.0
   if(term < 0.0) then
      ifail=1
   else
      term = sqrt(term)
      root1 = (-y+term)/a2
      root2 = (-y-term)/a2
   endif
end subroutine solve

The subroutines contained in the interface can also be placed in an external file.

Coding Fortran: subprograms (modularity)

There are many ways to include subprograms in a Fortran program. The next series of posts deals with these three methods: normal subprograms, interfaces, and use of the contains clause. The program calculates a roots of a quadratic equation, and has two subroutines, getcoeff() and solve(). Here the two subroutines are declared below the main program (in the same file).

program subroutine1
implicit none

   real :: p, q, r, root1, root2
   integer :: ifail=0
   logical :: ok=.true.

   call getcoeffs(p,q,r,ok)
   if (ok) then
      call solve(p,q,r,root1,root2,ifail)
      if (ifail == 1) then
         print*,'Complex roots, calculation abandoned'
      else
         print*,'Roots = ',root1,' ',root2
      endif
   else
         print*,'Error in data input.'
   endif
end program subroutine1

subroutine getcoeffs(a,b,c,ok)
   implicit none
   real, intent(out) :: a,b,c
   logical, intent(out) :: ok
   integer :: io_status=0

   print*,'Enter the coefficients a, b and c'
   read(unit=*,fmt=*,iostat=io_status) a,b,c
   if (io_status == 0) then
      ok = .true.
   else
      ok = .false.
   endif
end subroutine getcoeffs

subroutine solve(x,y,z,root1,root2,ifail)
   implicit none
   real, intent(in) :: x,y,z
   real, intent(out) :: root1,root2
   integer, intent(inout) :: ifail
   real :: term,a2

   term = y*y - 4.*x*z
   a2 = x*2.0
   ! if term < 0, roots are complex
   if (term < 0.0) then
      ifail = 1
   else
      term = sqrt(term)
      root1 = (-y+term)/a2
      root2 = (-y-term)/a2
   endif
end subroutine solve

The only thing to remember is that the types of the arguments and the parameters of a subprogram have to match. If the specifications for the parameters a, b, and c in getcoeffs() were changed to be integers, i.e.

integer, intent(out) :: a,b,c

Then a compiler will pick up on this and produce a compilation error (or most will). For example, compiling with gfortran produces the following error:

    9 |    call getcoeffs(p,q,r,ok)
      |                            1
Error: Type mismatch in argument 'a' at (1); passed REAL(4) to INTEGER(4)

There are two ways to avoid this problem altogether: interface blocks, and contained procedures (and technically modules as well). It is also possible to place the subprograms in a separate file, and associate them using the external clause in the main program.

external :: getcoeffs, solve

Coding Fortran: Validating input

How easy is it to validate user input in Fortran? Super easy. In fact the read() function has a built-in ability to return the validity of the input, using the keyword iostat, and an appropriate integer variable. Here is an example program which reads in an integer. Any other value input will generate the error message.

program validate

   integer :: x, err

   write(*,*) "Integer> "
   read(*,*,iostat=err) x
   if (err /= 0) then
      print *, "wrong input, try again"
   end if

end program validate

What are the values of err?

1. If the value of err is 0, the read() was executed flawlessly and all variables have received their input values.
2. If the value of err is positive, the previous read() has encountered some problem. In general this could mean illegal data, e.g. entering a real number, or a character.
3. If the value of err is negative, it means the end of the input has reached.

Fortran Reengineering: arctan (iv)

The last piece of reengineering deals with the remaining two arithmetic-if statements, which include labels 12, 15 and 16. This piece of code is basically a loop, which tests the condition at the end of the loop (so like a do-while loop in C). The original code snipped is shown below (arctan3.f95). All the arithmetic if statements to is terminate the loop if (term-0.0005) is < 0, AND (-term-0.00005) is < 0. Any other combination jumps the program back up to label 12, effecting the next iteration of the loop.

   12 arct=arct+term
      presxp=prevxp+2.0
      term=-prevxp/presxp*y*term
      prevxp=presxp
      if(term-0.00005) 15,12,12
   15 if(-term-0.00005) 16,12,12
   16 return

So this can be reengineered by adding a generic do loop. The only changes made to the statement in the loop are of an aesthetic kind. The two arithmetic if statements are replaced with a single if statement that decides whether the loop should be exited.

   do
      arct = arct + term
      presxp = prevxp + 2.0
      term = -prevxp / presxp * y * term
      prevxp = presxp
      if (term-0.00005 < 0 .and. -term-0.00005 < 0) then
         exit
      end if
   end do

The program is now reengineered. The only thing remaining is to clean up the code, and add some I/O usability to the main part of the program. Basically add in a user prompt, and format the output, giving it units. The reengineered program (arctan4.f95) is shown below.

program testreeng
   implicit none
   real :: x, y
   write(*,*) "Value of arctan to calculate? "
   read(*,*) x
   y = arctanF(x)
   write(*,10) "Arctan(",x,") = ", y, " rad"
10 format(a7,f7.2,a4,f7.4,a4)
   stop

contains

function arctanF(x) result(arct)
   implicit none
   real, intent(in) :: x
   real :: arct
   real :: term, prevxp, presxp, y

   if (x < 0) then
      return
   end if
   arct = 0.0
   if (x-1.0 > 0) then
      term = -1.0/x
      arct = 1.57079
   else
      term = x
   end if
   prevxp = 1.0
   y = term**2.0
   do
      arct = arct + term
      presxp = prevxp + 2.0
      term = -prevxp / presxp * y * term
      prevxp = presxp
      if (term-0.00005 < 0 .and. -term-0.00005 < 0) then
         exit
      end if
   end do
end function arctanF

end program testreeng

Here is a sample run of the program:

 Value of arctan to calculate?
45
Arctan(  45.00) =  1.5486 rad