# 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;
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: mloss = t (F/−V)
• Compute the net mass remaining after t seconds: mt (remaining) = (W/gravity) − mloss
• Compute the weight at Wt: Wt = mt (remaining) gravity
• Compute the unbalanced force that causes acceleration: Ft=F−Wt
• Compute the acceleration: a = Ft / mt

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``````