coding Fortran

A recursive function for reading lines (from a file)

Sometimes all we want to do is read a line of text from a file. But let’s make it more interesting and use recursion. The following Fortran module readlines reads lines from a text file. The read statement uses non-advancing I/O to read a small piece of the! current line in a file. If the end-of-line is reached, the file pointer is moved to the start of the next line. It uses a recursive subroutine piecemeal_readline() to do all the hard graft. In this case it reads chunks of text, 10 characters in length, and appends them to the line. When the end-of-line is encountered, the line is returned. 

module readlines
   use iso_fortran_env
   implicit none

contains

subroutine readline(lun, line, success)
   integer, intent(in) :: lun
   character(len=:), allocatable, intent(out) :: line
   logical, intent(out) :: success

   character(len=0) :: newline

   success = .true.
   call piecemeal_readline(newline)

contains

   recursive subroutine piecemeal_readline(newline)
      character(len=*), intent(in) :: newline

      character(len=10) :: piece
      integer :: ierr, sz

      read(lun, '(a)', advance='no', size=sz, iostat=ierr) piece
      if (ierr /= 0 .and. ierr /= iostat_eor) then
         allocate(character(len=len(newline)):: line)
         line = newline
         success = .false.
         return
      endif

      ! End of line?
      if (sz >= len(piece)) then
         call piecemeal_readline(newline // piece)
      else
         allocate(character(len=len(newline)+sz):: line)
         line = newline // piece(1:sz)
         success = .true.
      endif
   end subroutine piecemeal_readline
end subroutine readline

end module readlines

Here is a passage to test what is happening (from Sense and Sensibility, Jane Austen):

I am excessively fond of a cottage; there is always so much comfort, so much
elegance about them. And I protest, if I had any money to spare, I should buy
a little land and build one myself, within a short distance of London, where I
might drive myself down at any time, and collect a few friends about me and be
happy. I advise everybody who is going to build, to build a cottage.

Basically let’s look at what happens at Line 34, which determines whether the end of a line has been encountered. If we process just the first line of the above text, using string chunks 10 in length, we get the following if the if statement is true:

 I am exces
 I am excessively fon
 I am excessively fond of a cot
 I am excessively fond of a cottage; ther
 I am excessively fond of a cottage; there is alway
 I am excessively fond of a cottage; there is always so much
 I am excessively fond of a cottage; there is always so much comfort, s
 I am excessively fond of a cottage; there is always so much comfort, so much

Basically, because it is not end-of-line, the current piece of text is appended to newline, which is then used as input to the recursive call of piecemeal_readline(). Each recursive call adds a chuck onto the line of text. When the end of line is encountered, the other part of the if statement is processed, returning the line as a string to the calling program. Here is a sample program that uses the module readlines.

program dynamo
   use readlines
   implicit none

   character(*), parameter :: filename = 'passage2.txt'

   integer :: lun
   logical :: ios
   character (len=:), allocatable :: str

   open(unit=lun, file=filename, status='old', action='read')

   ios = .true.
   do while (ios .eqv. .true.)
      call readline(lun,str,ios)
      if (ios .eqv. .true.) then
         print*, str
         print*, "Length = ", len(str)
      endif
   enddo
   close(lun)

end program dynamo

Here is the program running with the passage from above (in passage4.txt):

 I am excessively fond of a cottage; there is always so much comfort, so much
 Length =           77
 elegance about them. And I protest, if I had any money to spare, I should buy
 Length =           78
 a little land and build one myself, within a short distance of London, where I
 Length =           79
 might drive myself down at any time, and collect a few friends about me and be
 Length =           79
 happy. I advise everybody who is going to build, to build a cottage.
 Length =           68

Some Fortran string tips

it can get a little confused with strings in Fortran. For example, consider the following piece of code:

program str
   implicit none
   character (len=20) :: s1, s2

   s1 = "Obi-Wan"
   s2 = 'Kenobi'

   print*, s1, " ", s2

end program str

The string s1 is created with the first seven elements taking up the characters “Obi-Wan“, and the remaining 13 elements filled with trailing spaces. Here either single or double quotes can be used to designate a string – although the best way is to choose one and stick with it, just for the sake of consistency. Functionally there is no difference.

Now if you want to derive the length of a string, you can use a function called len_trim(), which calculates the length of a string minus trailing blanks. So the length of s1, len_trim(s1), would be 7. If you just use plain old len(), then the value returned would be 20, the size of the string.

String concatenation can be done using the // operator. For example (assuming s3 is the same string as s1 and s2):

s3 = s1 // s2

The problem here is that the answer would just be “Obi-Wan”. Why? Because s3 is only 20 characters in length, and s1 alone is 20 characters, so there is no room for s2. The solution is to trim away the trailing blanks:

s3 = trim(s1) // trim(s2)

Then the solution is Obi-WanKenobi. Finally, let’s consider strings that are built character by character. For example the string s4, declared just like s1. Now a loop is used to read four characters from standard input, and assign them to the first four elements of the string. 

 do i = 1,4
    read*, s4(i:i)
 end do
 print*, s4, len_trim(s4)

Now when this code is run, and the characters l, e, n, and s are input, the following output is produced:

lens M�l��          20

Obviously there is a problem here, and it relates to the fact that the trailing characters have not been filled with spaces, but rather garbage. But the solution is easy, just add, s4="" before the loop. Note that this is not an issue if the string is input whole.

More 2D strings in Fortran

Arrays of strings are easy to do. For example the following declaration produces an array of strings, where there are 5 strings, each of 10 characters in length:

character(len=12), dimension(4) :: a

This could also be written as:

character(len=12) :: a(4)

However the first declaration is the one that should be used, the second declaration is associated with old versions of Fortran, and should be avoided. Here is a program which deals with some strings, and prints out various things. 

program strarr
   implicit none

   character(len=12), dimension(4) :: a
   a(1) = "Mandalorian"
   a(2) = "BobaFett"
   a(3) = "CaraDune"
   a(4) = "IG-11"

   print*, a
   print*, a(2)(5:8)
end program strarr

You will notice there are two types of indices: one for the index of the array, and another for the position of the character. These two indices are specified separately, with the index first, and the position of the character in the string second. So in the example a(2)(5:8) means array element 2, string indices 5 to 8. Note that because this is a string and not a character array, to access a single element you have to use a range in the for x:x. Here is the output from the program above.

 Mandalorian BobaFett    CaraDune    IG-11
 Fett

The strings in the array could also be allocated, so that the storage for each string in the array is variable length. Here is a program that uses an array of allocatable strings.

program arrvarstr
   implicit none

   type :: varl
      character(len=:), allocatable :: name
   end type varl

   type(varl), dimension(4) :: a(4)
   integer :: i

   allocate(character(len=11) :: a(1)%name)
   a(1)%name = "Mandalorian"
   allocate(character(len=8) :: a(2)%name)
   a(2)%name = "BobaFett"
   allocate(character(len=8) :: a(3)%name)
   a(3)%name = "CaraDune"
   allocate(character(len=5) :: a(4)%name)
   a(4)%name = "IG-11"

   do i = 1,4
      print*, a(i)%name
   end do
   print*, a(2)%name(5:8)

end program arrvarstr

Here is the output:

 Mandalorian
 BobaFett
 CaraDune
 IG-11
 Fett

Finally, there a variable length array of (static) strings. Here is an example of a variable length array of strings of size 20. The length of the array of strings is allocated to 40.

program arrvarlen

   implicit none

   character(len=20), dimension(:), allocatable :: c
   allocate(c(40))

   c(1) = "do"
   c(2) = "or"
   c(3) = "do not"

   c = [c, "there is no try"]

   print*, c
   deallocate(c)

end program arrvarlen

Note there is also an interesting piece of code on Line 12. Here a new string is appended to the next position in the array (in this case element 4). Here is the output from the program:

 do                  or                  do not              there is no try

Calculating wind chill in Fortran

As temperatures fall and the wind howls, we begin hearing about the danger of “wind chill”. The “wind chill” factor W is reported by meteorologists during winter. W is an equivalent temperature that accounts for the increased chilling effects of wind on a human body. The wind chill factor combines the temperature and wind speed to tell you how cold the wind makes it “feel”. The coldest wind chill recorded in Canada was at Pelly Bay, Nunavut, on January 13, 1975, when 56 km/h winds (a wind chill factor of 3,357 watts/m²) made the temperature of -51°C feel more like -92°C.

The first wind chill formula and associated tables were based on research conducted by scientists Paul A. Siple and Charles F. Passel in Antarctica in the 1940s. The research found that the rate at which water freezes depends on three factors: how warm it was to begin with, the outside temperature and the wind speed. This Wind Chill Index was a 3-4 digit number with units of kilocalories per square metre per hour (kcal/m2/hr). For example 1400 was the threshold for frostbite. The original Siple-Passel equation for wind chill factor (°C) is:

C = (10.0√V - V + 10.5)(33-Ta)

where V is the wind velocity (m/s), and Ta was the air temperature (°C). One of the problems with this formula was that the units were not very user-friendly for the general public. In 2001 a new formula was derived , based on a model of how fast a human face loses heat and incorporates modern heat transfer theory, that is, the theory of how much heat is lost by the body to its surroundings during cold and windy days. (The face is the part of the body most often exposed to severe winter weather). It must be noted that although the wind chill factor is expressed on a temperature scale (the Celsius scale in Canada), it is not a temperature: it only expresses a human sensation. There are two formulas used by Environment Canada:

W = 13.12 + 0.6215Tair - 11.37V(10m)0.16 + 0.3965Tair V(10m)0.16
W = Tair + [(-1.59 + 0.1345Tair)/5]V(10m)

where W is the wind chill factor based on the Celsius temperature scale, Tair is the air temperature in degrees celsius, and V10m is the wind speed in km/h at 10 meters. The first equation is used when the temperature of the air is ≤ 0°C and the wind speed is ≥ 5km/h. Then second equation is used when the temperature of the air is ≤ 0°C and the wind speed is > 0km/h, but < 5km/h.

Here is the Fortran program to perform the calculation:

program wind_chill

    real :: airTemp, windS, windCF

    ! Obtain the user input
    write (*,*) 'Air temperature (Celsius): '
    read (*,*) airTemp
    write (*,*) 'Wind speed (km/hr): '
    read (*,*) windS

    if (airTemp <= 0.0) then
        if (windS >= 5.0) then
            windCF = 13.12 + 0.6215*airTemp - 11.37*windS**0.16 + &
                     0.3965*airTemp*windS**0.16;
        else if ((windS > 0.0) .and. (windS < 5.0)) then
            windCF = airTemp + ((-1.59+0.1345*airTemp)/5.0)*windS;
        else
            write (*,*) 'There is no wind!'
        end if
        write (*,100) windCF
        100 format ('The temperature feels like ', F7.2, ' degrees Celsius')
    else
        write (*,*) 'Unable to calculate - the air temperature is too high'
    end if

end program wind_chill

The program is inherently easy to write in Fortran. Here’s the program running:

Air temperature (Celsius):
-3
Wind speed (km/hr):
15
The temperature feels like   -8.12 degrees Celsius

Coding Fortran: A dynamic array of ragged strings

One of the tricks with Fortran is dealing with strings. Easy strings are easy, but more, shall we say challenging forms of strings can be tricky. Imagine if we wanted a dynamic array of ragged strings, i.e. strings that are all of different sizes.

One way of doing this is to create some derived data types. First we create one for a ragged string, ragged_str:

type :: ragged_str
   character(len=:), allocatable :: word
end type ragged_str

It has one member, word. This specifies that word is a string of deferred length, i.e. allocatable. Dynamic arrays are created using the attribute allocatable. Now we can create a second derived data type, strList, which creates a dynamic length array comprised of ragged strings. It’s member, list is a deferred length “array” of type ragged_str.

type :: strList
   type (ragged_str), dimension(:), allocatable :: list
end type strList

Now we can create a variable of this type:

type(strList) :: d

Now that all the data types have been created, they need to be allocated before being used. In the first case, we need to allocate the length of the list. Let’s say we need a list 10 words in length.

allocate(d%list(10))

Note that the percentage symbol % is used to access the members of a derived type. Now we need to allocate the size for each of the 10 ragged string elements. For the first string, 7 characters in length, this can be done in the following manner:

allocate(character(len=7) :: d%list(1)%word)

As an example of how to use this to create a structure contains n words input by the user:

character (len=20) :: text
integer :: n, i, wl
print*, "How many words? "
read(*,*) n

allocate(d%list(n))

do i = 1,n
   read(*,*) text
   wl = len(trim(text))
   allocate(character(len=wl) :: d%list(i)%word)
   d%list(i)%word = trim(text)
end do

Here is what happens in Lines 10-12:

  1. The length of the word, wl, is calculated (after using the function trim() to remove trailing blank characters from the string text).
  2. Allocate the size of the member word, based on wl.
  3. Assign the word text to the member d%list(i)%word.

NB: The allocatable arrays are always measured in number of array elements, not the number of bytes or other storage units. When they are allocated, the number of array elements in each dimension are specified, not the total amount of storage the array will consume.

Memories of Fortran II

What about if you want the array to hang around? Add the save attribute to the array. For example, consider the following subroutine, alloc_arr(), where save is used to maintain the array arr.

subroutine alloc_arr(n)
   integer, intent(in) :: n
   integer, dimension(:), allocatable, save :: arr
   if (.not. allocated(arr)) then
      allocate(arr(n))
      do i=1,n
         arr(i) = i
      end do
   else
      do i=1,n
         arr(i) = arr(i) + 2
      end do
   end if
   print *, arr
end subroutine alloc_arr

Here a check is made to make sure the array arr has not already been allocated, if not, it is allocated and assigned the values 1..n, where n is the size of the array to create. If the subroutine is called again, and arr already exists, the value of its elements are incremented by 2. For example if alloc_arr() is called three times, the output looks like this:

           1           2           3           4           5           6
           3           4           5           6           7           8
           5           6           7           8           9          10

The magical save attribute can also be applied to other variables, for example to store the state of a counter. Here is another program which counts the number of calls of a recursive function.

recursive function factorial(n) result(r)
   integer, intent(in) :: n
   integer :: r
   integer, save :: nr = 0

   if (n == 1) then
      print*, "No. calls =", nr
      r = 1
   else
      nr = nr + 1
      r = n * factorial(n-1)
   end if

end function factorial

Here there is one saved variables, nr, which counts the number of function calls. The counter nr will increment each time there is a recursive call (the first call to factorial() is not recursive). The variable nr will be maintained, even when the function terminates.

Memories of Fortran

Memory management is a horrible thing. It’s one of the things novice programmers hate most about languages like C and C++. In C, poor memory management eventually leads to “memory leaks”, when you create memory that you loose track of somewhere along the way. I mean it’s not hard to do. Of course it’s also possible to use up all the memory because someone forgot to deallocate it after it has been used. Other languages like Java have “garbage collectors” that periodically clean up things (likely they should be re-termed “memory reuse managers”). C++ uses constructors and destructors (sounds like a transformer of sorts).

Fortran handles thing a little differently, i.e. it automatically deallocates any pointer to memory that goes out of scope. For example here is a piece of code that creates a “dynamic” array, x, using allocate(), which is created in heap memory.

program mem2
   integer, dimension(:), allocatable :: x
   integer :: i, n
   write(*,*) "How large an array? "
   read(*,*) n
   allocate(x(n))
   do i = 1,n
      x(i) = i*i
   end do
end program mem2

Notice though that x is not deallocated. In C similar code would be a disaster, but not in Fortran. Why? Because the Fortran standard requires that when an allocatable array goes out of scope it should be deallocated. You can use deallocate to explicitly manage the lifetime of a piece of memory. If in a function, the function will allocate the array as requested and then deallocate it again as soon as the function is over. Note that x can be used like a normal array once allocated. In order to initialize and array, the keyword source can be used. For example:

allocate(x(n), source=0)

Fortran 2003 also allows array arguments to a subprogram to be allocatable. For example the subroutine print_arr() shown below prints out a 1D array, and could be used in the above example. The parameter arr is an allocatable array. To make sure it is actually initialized, the function allocated() can be used. 

subroutine print_arr(arr)
   integer, dimension(:), allocatable, intent(in) :: arr
   if (allocated(arr)) then
      do i = 1,size(arr)
         write(*,'(i3,x)',advance="no") x(i)
      end do
   end if
end subroutine print_arr

It is also possible to move allocations can be moved between different arrays with the move_alloc() subroutine. The subroutine shown below, resize(), resizes an allocatable array, to a new length, n. If the array arr is already allocated, it moves the contents to tmp using move_alloc()(line 9), and in the process deallocates arr. Then arr is allocated to size n (line 12), and if tmp is allocated, transfers the relevant number of elements back to arr (lines 13-16).

subroutine resize(arr, n)
   integer, dimension(:), allocatable, intent(inout) :: arr
   integer, intent(in) :: n
   integer, dimension(:), allocatable :: tmp
   integer :: arsz,nsz

   if (allocated(arr)) then
      arsz = size(arr,1)
      call move_alloc(arr,tmp)
   end if

   allocate(arr(n))
   if (allocated(tmp)) then
      nsz = min(arsz,n)
      arr(:nsz) = tmp(:nsz)
   end if
end subroutine resize

Fortran – the greatest thing (for arrays) since bread came sliced

Many people likely find the idea of programming in Fortran offensive. It’s from a bygone era. Modern languages are cooler. But are they? I have written many a blog post on why Fortran is a better language for novice programmers. Better than C anyway. I know some people will disagree, but frankly unless you have spent 20 years trying to teach C as an introductory programming language, you have no real context. Even for certain tasks, Fortran may be a better choice, and by this I mean crunching numbers.

Fortran has an easy-to-use array syntax, and that’s not surprising considering it was developed within a culture of high-performance computing – crunching numbers. Nobody will really ever build as OS, or indeed a compiler using Fortran, for its forte is numbers (and truly, C was designed as a systems language). This is apparent in the fact that most algorithms implemented in Fortran aren’t much slower than those implement in C, and in some cases they are actually faster.

C sucks at dealing with arrays, at least in a less-than-complicated way. For example when you teach someone about how to implement an algorithm to process images, you don’t really want to have to dive deeply into the processes of memory management. The good thing about Fortran is that the system handles all the dirty work associated with memory allocation, and so the focus can be on the algorithm. Let’s face it, Fortran introduced the original array syntax. The intuitive use of arrays in Fortran is one of the reasons that it is still heavily used in the physics community. Consider the following array operations in Fortran:

integer, dimension(1000) :: a,b,c
integer :: s
...
a = b
a = 1.73*b
c = a * b
c = matmul(a,b)
s = sum(c)

Here’s what each line of code does:

  1. Array b is copied to array a.
  2. Array b is multiplied by the scalar 1.73, and the result copied to array a.
  3. The element-by-element multiplication of arrays a and b (assuming they are the same size).
  4. The matrix multiplication of arrays a and b.
  5. Sum the values in array c, and assign the value calculated to s.

Almost all of the intrinsic math functions in Fortran can take arrays as input. To perform many of these tasks in C requires cycling through all the elements of the array using a for loop.

In addition, Fortran arrays can be indexed using any starting value, making an the language conform to the algorithm, rather than the other way around (and if you really like arrays that start at zero, you can use them). For anyone looking for a good example of an algorithm that uses negative array indices, look no further than Niklaus Wirth’s classic 8-Queens problem [1]. For a multidimensional array in C, one has to use the syntax x[i][j], whereas in Fortran it is simply x(i,j). Here are some other cool features of arrays:

integer, dimension(1000) :: x, z
integer, dimension(40) :: y
x = (/ (i, i = 1,1000) /)
y = x(1:1000:25)
z = 4
z(300:339) = y

Here’s what each line of code does:

  1. The array x is initialized using an implicit do loop, called an array constructor.
  2. The array y is created from every 25th element in array x.
  3. The array z is initialized to the value 4 in every element.
  4. The array y is copied into array z, in elements 300 to 339.

It is also possible to selectively process an array using the where-elsewhere construct. For example in the code below, only the values of the array x that are not 0 are processed by the log() function, all zero values are set to 0.1 in the output array logx.

real, dimension(10) :: x, logx
x = (/ (i, i = 1,10) /)
x(5) = 0.0
where (x .ne. 0.0)
   logx = log(x)
elsewhere
   logx = 0.1
end where

Calling Fortran old is like calling C old. Yes, there are new languages like Julia which handle arrays nicely too (Python does not even have real arrays you need to add numpy), but Julia has become a somewhat bloated language with no backwards compatibility (and it still suffers from compiler messages that *really* suck).

Fortran excels at processing arrays, just as C excels at being a systems language, Ada excels at being a real-time language, and Cobol excels at being a data processing language. Each does their own thing, and there is no need for one to supplant another. There will always be new languages, but time has shown that sometimes sticking with a known entity for a particular task is the most optimal approach.

  1. Wirth, N., Algorithms + Data Structures = Programs, Prentice-Hall (1976).

Coding Fortran: external vs. module vs. interface

People seem to get these two concepts mixed up quite often. There are a number of ways of accessing external subprograms in a Fortran program. First there are external subprograms, which are defined in a separate file. They are defined using the external keyword. Here is a simple example, using a function for a recursive factorial. The main program, factext.f03, is:

program factorialR
   implicit none
   integer :: n, fact
   integer, external :: factorial

   write(*,*) "n? "
   read(*,*) n
   fact = factorial(n,1)
   write(*,'("Factorial ", i2, " = ", i20)') n, fact
end program factorialR

The external file, factR.f03, looks like this:

recursive function factorial(n,acc) result(r)
   integer, intent(in) :: n, acc
   integer :: r
   if (n == 0) then
      r = acc
   else
      r = factorial(n-1, n*acc)
   end if
end function factorial

The program is compiled in the following way:

gfortran factR.f95 -c
gfortran factext.f95 factR.o

In reality though, external is not the best way to approach the problem of external files. A better way is using the interface block. This allows an external function to be used by listing it an interface block along with the declaration of its arguments and their types and the type of the function value. For example, the above factR.f03 could be modified to include an interface block (keeping factR.f03 the same).

program factorialR
   implicit none
   integer :: n, fact, acc
   interface
      recursive function factorial(n,acc) result(r)
         integer :: r
         integer, intent(in) :: n, acc
      end function factorial
   end interface

   write(*,*) "n? "
   read(*,*) n
   acc = 1
   fact = factorial(n,acc)
   write(*,'("Factorial ", i2, " = ", i20)') n, fact
end program factorialR

Finally, there are modules. A module helps to modularize code. The subprograms used in modules can be used in any other program unit, as long as a use statement is added. Modules are used to group subprograms and data structures (but not data). Constants can be defined in a module, but avoid defining other variables. Here is the factorial function added to a module, factRm.f95.

module factRm
   implicit none
contains
   recursive function factorial(n,acc) result(r)
      integer, intent(in) :: n, acc
      integer :: r
      if (n == 0) then
         r = acc
      else
         r = factorial(n-1, n*acc)
      end if
   end function factorial

end module factRm

Here is the program, factmod.f95, which uses the module factRm:

program factorialR
   use factRm
   implicit none
   integer :: n, fact, acc

   write(*,*) "n? "
   read(*,*) n
   acc = 1
   fact = factorial(n,acc)
   write(*,'("Factorial ", i2, " = ", i20)') n, fact
end program factorialR

It can be compiled in the following manner:

fortran factRm.f95 factmod.f95

Note that the module must precede the calling program. During compilation, a module file – factrm.mod – is created.