Fortran basics

Introduction

Fortran (derived from Formula Translation) is a programming language suitable for numerical computations. The version we will use is Fortran 90/95 which we will refer to as Fortran 90 for brevity. On the one hand, compared to earlier versions (such as Fortran 77) Fortran 90 benefits from many new features. On the other hand, unlike the newer versions (such as Fortran 2003 and 2008) Fortran 90 is supported by most compilers.

See https://en.wikipedia.org/wiki/Fortran for a brief history of Fortran.

Note. This page is in no way a complete description of all the aspects of Fortran. It introduces some basic concepts through a sequence of examples. For more in-depth coverage, refer to the sources listed at the bottom of this page.

Running Fortran programs

A Fortran program must pass through several stages before being executed. We usually write a program in one or more source files that will be converted into object files (with .o extension) that are versions of codes in a machine language. This is done by a compiler. The objects will then be linked together to produce an executable file which can be run in the shell. This is done by a linker. There are several different compilers that can turn a Fortran code into an executable (and hence do both compiling and linking). In this class we will use gfortran, which is an open source compiler. It is part of the GNU Project. Note that there are several commercial compilers, such as the Intel and Portland Group compilers, which sometimes do better optimization and produce faster running executables. They also may have better built-in debugging tools.

For the gfortran compiler, a Fortran code has the extension .f90. Three ways to run a code, say newton.f90, follow:

Way I: We can combine both stages (compiling and linking) in a gfortran command:

$ gfortran newton.f90

This will produce an executable file named a.out. To run the code we would then type:

$ ./a.out

Way II: If you don’t like the strange name a.out, you can specify an output file name using the -o flag with the gfortran command. For example:

$ gfortran newton.f90 -o newton.exe    # or alternatively  $ gfortran -o newton.exe newton.f90
$ ./newton.exe

Way III: We can split up the compile and link steps. This is often useful if there are several source files that need to be compiled and linked. We can do this using the -c flag to compile without linking and then using the -o flag:

$ gfortran -c newton.f90
$ gfortran newton.o -o newton.exe     # or alternatively  $ gfortran -o newton.exe newton.o
$ ./newton.exe

Remark. Note that in the first two ways, there will be no .o file created. By default gfortran removes this file once the executable file is created.

Example 1: a simple code

This simple example assigns a number to a variable, performs a few mathematical operations on the variable, and then prints the output. The comments below the code explain some features.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
! Example 1 for Math 471
program ex1
  implicit none
  real(kind=8) :: a,b,c

  a = 1.d0
  b = exp(A)
  c = a+b
  print *, "The output is = ", c
end program ex1

Notes:

  • In Fortran comments start with an exclamation mark !.
  • A fortran program always starts with the program [name] statement and ends with the end program [name] statement.
  • The implicit none statement on line 3 means that any variable to be used must be explicitly declared. You should always use this. It will give error if a variable is not declared.
  • The variables a,b,c are real numbers stored in 8 bytes (corresponding to double precision). Note that (unlike Perl) Fortran is case insensitive so a and A are the same.
  • The notation 1.0d0 indicates that the number is double precision.
  • Fortran has many built-in functions, such as exp, log, sin, etc.
  • To print to screen the print command can be used. The * indicates that “As much information as possible” should be displayed.

To run the code we can use the third way described above:

$ gfortran -c ex1.f90
$ gfortran -o ex1.x ex1.o
$ ./ex1.x
The output is =    3.7182818284590451

Example 2: do-loop and if-then-else

This example illustrates the use of do-loops and if then else statements. This simple code shows that one has to be careful when performing arithmetic operations with different precision types. You should be careful especially if you are only used to programming in Matlab.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
! Example 2 for Math 471
program ex2
  implicit none
  integer, parameter :: n = 3
  real(kind = 8), parameter :: pi = 3.1415926535897932d0
  integer :: i
  real(kind = 8) :: a
  real :: c = 0.0
  logical :: correct
  do i = 1,n
     if ( i < 2) then
        c = c + pi
        a = abs(dble(i)*pi-c)/(dble(i)*pi)
        correct = .false.
     elseif (i == 2) then
        a = (i/n)*pi
        correct = .false.
     else
        a = (i/n)*pi
        correct = .true.
     endif
     print *, i, a, correct
  end do
end program ex2

Notes:

  • The parameter attribute tells the compiler that the variables that are declared with that attribute will not change.

  • If a variable is declared as real, by default it will become real(kind=4), i.e. single precision. So, real :: c is the same as real(kind = 4) :: c. If you always want to default to real(kind = 8) you can use the compiler flag -fdefault-real-8 (assuming you are using gfortran). This is however not recommended. It is better to try to be consistent in declaring variables to be whatever precision you want them to be. For most of what we’ll do in this class, we will use real numbers with (kind=8).

  • 3.1415926535897932d0 specifies a 8-byte real number. The 4-byte equivalent is 3.1415926535897932e0, with an e instead of d. Be careful to specify constants using the d rather than e notation if you need to use scientific notation.

  • The first if statement illustrates the difference between “single” and “double” precision. The built-in function dble converts its argument into double precision. Hence dble(i)*pi is a double precision computation while c = c + pi is stored in a single precision number. Consequently, the relative difference abs(dble(i)*pi-c)/(dble(i)*pi) is of the order \(10^{-8}\). Note that in Matlab we would get 0, because all computations are carried out with double precision.

  • The second and third parts of the if statement show how integer division “floors” (2/3) to 0, while (3/3) = 1.

  • Note that the control statement can be expressed either with symbols or letters

    (a < b)  same as  (a.lt.b)
(a <= b)          (a.le.b)
(a > b)           (a.gt.b)
(a >= b)          (a.ge.b)
(a == b)          (a.eq.b)

  • We can combine control statements with .and. and .or. to generate more complex control statements. For example ((i>=1) .and. (i<=2)) .or. (i>5) is a boolean variable that is declared with type logical. It takes a logical value: either true or false.

  • Sometimes it is useful to put an exit statement inside a do-loop if a condition is satisfied: if (error <= 1.0d-10) exit.

Executing the above code yields

$ gfortran ex2.f90; ./a.out
  1   2.7827535191837951E-008 F
  2   0.0000000000000000      F
  3   3.1415926535897931      T

The first two computations are wrong. Indeed, we wanted to compute \(a = | \pi - \pi | / | \pi | = 0\) (when \(i=1\)) and \(a = 2 \pi / 3 \neq 0\) (when \(i=2\)).

Example 3: more built-in functions

This example shows more built-in functions.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
! Example 3 for Math 471
program ex3
  implicit none
  real (kind=8) :: pi, x, y
  pi = acos(-1.d0)
  x = sin(pi)
  y = sqrt((log(pi)))**2
  print *, "pi = ", pi
  print *, "x = ", x
  print *, "y = ", y
end program ex3

Notes:

Executing the above code yields

$ gfortran ex3.f90; ./a.out
   pi =    3.1415926535897931
   x =    1.2246467991473532E-016
   y =    1.1447298858494002

Example 4: arrays

This example demonstrate declaring and using arrays.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
! Example 4 for Math 471
program ex4
  implicit none
  integer, parameter :: m = 3, n=2
  real (kind=8), dimension(m,n) :: A
  real (kind=8), dimension(n) :: x
  real (kind=8), dimension(m) :: y
  integer :: i,j
  ! initialize matrix A and vector x
  do j=1,n
      do i=1,m
          A(i,j) = i+j
      enddo
      x(j) = j
  enddo
  ! compute y=A*x
  do i=1,m
      y(i) = 0.
      do j=1,n
          y(i) = y(i) + A(i,j)*x(j)
      enddo
  enddo
  print *, "A = "
  do i=1,m
      print *, A(i,:)   ! i-th row of A
  enddo
  print *, "x = "
  print *, x
  print *, "y = A*x = "
  print *, y
end program ex4

Notes:

  • Arrays are indexed starting at 1 by default, rather than 0 as in Perl. Also note that components of an array are accessed using parentheses, not square brackets.
  • We can also declare arrays and determine their size later, using allocatable arrays. We will see this later.

Executing the above code yields

$ gfortran ex4.f90; ./a.out

A =
  2.0000000000000000        3.0000000000000000
  3.0000000000000000        4.0000000000000000
  4.0000000000000000        5.0000000000000000
x =
  1.0000000000000000        2.0000000000000000
y = A*x =
  8.0000000000000000        11.000000000000000        14.000000000000000

Example 5: array operations

This example illustrates the use of array operations and array intrinsic functions.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
! Example 5 for Math 471
program ex5
  implicit none
  real(kind=8), dimension(3) :: x, y, z
  real(kind=8), dimension(3,2) :: a
  real(kind=8), dimension(2,3) :: b
  real(kind=8), dimension(3,3) :: c
  integer i

  x = (/1.,2.,3./)
  y = (/1.,4.,9./)
  a = reshape((/1,2,3,4,5,6/), (/3,2/))

  print *, "x**2 + x*y = "
  print *, x**2 + x*y

  print *, "sqrt(y) = "
  print *, sqrt(y)

  print *, "dot_product(x,y) = "
  print *, dot_product(x,y)

  print *, "a = "
  do i=1,3
      print *, a(i,:)
  enddo

  b = transpose(a)
  print *, "b = transpose(a) = "
  do i=1,2
      print *, b(i,:)
  enddo

  c = matmul(a,b)
  print *, "c = a*b = "
  do i=1,3
      print *, c(i,:)
  enddo

  z = matmul(c,x)
  print *, "z = c*x = ",z
end program ex5

Notes:

  • Array operations (such as +, *, **) and the intrinsic functions (such as sqrt) apply component-wise.
  • For matrix multiplication, we can use the intrinsic function matmul. Other useful intrinsic functions include reshape and transpose.

Executing the above code yields

$ gfortran ex5.f90; ./a.out
x**2 + x*y =
  2.0000000000000000        12.000000000000000        36.000000000000000
sqrt(y) =
  1.0000000000000000        2.0000000000000000        3.0000000000000000
dot_product(x,y) =
  36.000000000000000
a =
  1.0000000000000000        4.0000000000000000
  2.0000000000000000        5.0000000000000000
  3.0000000000000000        6.0000000000000000
b = transpose(a) =
  1.0000000000000000        2.0000000000000000        3.0000000000000000
  4.0000000000000000        5.0000000000000000        6.0000000000000000
c = a*b =
  17.000000000000000        22.000000000000000        27.000000000000000
  22.000000000000000        29.000000000000000        36.000000000000000
  27.000000000000000        36.000000000000000        45.000000000000000
z = c*x =    142.00000000000000        188.00000000000000        234.00000000000000

Example 6: functions

There are two kinds of procedures in Fortran, functions and subroutines. A function returns a single scalar variable of some designated type but can take multiple inputs. The code below shows how to compute the square of a double precision number x using a function.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
! Example 6 for Math 471
program ex6
  implicit none
  real(kind = 8) :: x,y
  real(kind = 8), external  :: myfun
  x = 3.d0
  y = myfun(x)
  write(*,*) "x = ",x, "and y = ",y
end program ex6

real(kind = 8) function myfun(x)
  real(kind = 8), intent(in) :: x
  myfun = x*x
end function myfun

Notes:

  • Functions must be declared as external in the main program. External tells the compiler that it is a function defined elsewhere rather than a variable.
  • Use intent(in) to help the compiler know that the variable x is passed into the function and won’t be changed inside the function.
  • If you need to modify/change the input arguments to a function, use intent(inout). This is however not very common. In such cases, subrutines are usually used.
  • The Write statement can also be used to output the results (instead of print). The first argument takes a file identification number or if you want to output to std out use * as is done here. The second argument is for specifying the format of the output. Again * means “as much as possible”.
  • Here the main program and the function are in the same file. We can also break things up so that functions or subroutines are in separate files.

Executing the above code yields

$ gfortran ex6.f90; ./a.out
        x =    3.0000000000000000      and y =    9.0000000000000000

Example 7: subroutines

Subroutines are used much more frequently than functions. While functions produce a single output variable, subroutines can have any number of inputs and outputs. In particular they may have not any output variable. This example shows how we can use a subroutine for computing the square of x.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
! Example 7 for Math 471
program ex7
  implicit none
  real(kind = 8) :: x,y
  x = 3.d0
  call mysub(x,y)
  write(*,*) "x = ",x, "and y = ",y
end program ex7

subroutine mysub(x,y)
  implicit none
  real(kind = 8), intent(in)  :: x
  real(kind = 8), intent(out) :: y
  y =  x*x
end subroutine mysub

Notes:

  • Subroutines may not explicitly return anything. They typically change some of its arguments.
  • Subroutines do not have to be declared.
  • The intent specification helps the compiler (and the programmer) to recognize what is to be changed inside the subroutine and what is not.

Example 8: subroutines, cont.

We need to be careful with different precision types when calling functions and subroutines. This example shows what happens when a subroutine intended to be called with 8 byte numbers is called with a real(kind = 4) number and an integer.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
! Example 8 for Math 471
program ex8
  implicit none
  real(kind = 8) :: x,y
  real :: s
  integer :: k
  x = 3.d0 ; s = 3.0 ; k = 3
  call mysub(x,y)
  write(*,*) "x = ",x, "and y = ",y
  call mysub(s,y)
  write(*,*) "x = ",s, "and y = ",y
  call mysub(k,y)
  write(*,*) "x = ",k, "and y = ",y
end program ex8

subroutine mysub(x,y)
  implicit none
  real(kind = 8), intent(in)  :: x
  real(kind = 8), intent(out) :: y
  y =  x*x
end subroutine mysub

Below is the output. Again, not what you expect if you are used to working with Matlab.

$ gfortran ex8.f90; ./a.out
 x =    3.0000000000000000      and y =    9.0000000000000000
 x =    3.00000000     and y =    0.0000000000000000
 x =            3 and y =    2.8033264416021267E+306

More on arrays

  1. To generate a 1D array x of n evenly spaced points between a and b, after declaring and initializing the variables, you can type:

    h=(b-a)/(n-1)
x=a+h*(/(i,i=0,n-1)/)

    This will act similar to x=linspace(a,b,n) in Matlab.

    You can write a Fortran function, named linspace, and use it whenever needed. Since such a function returns an array, which is not a single value, you will need to include an interface (line 3-11 in the code below) wherever you call such a function:

     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
    program TESTlinspace

  interface linspace
        function linspace(a,b,n)
          real(kind=8), intent(in) :: a,b
          integer, intent(in) :: n
          real(kind=8) :: linspace(n)
          real(kind=8) :: h
          integer :: i
        end function linspace
  end interface linspace

  real(kind=8) :: a,b
  integer, parameter :: n=5
  real(kind=8), dimension(n) :: x
  a=1.d0
  b=2.d0
  x = linspace(a,b,n)
  write(*,*) "x = ",x

end program TESTlinspace


function linspace(a,b,n)
  real(kind=8), intent(in) :: a,b
  integer, intent(in) :: n
  real(kind=8) :: linspace(n)
  real(kind=8) :: h
  integer :: i
  h=(b-a)/(n-1)
  linspace=a+h*(/(i,i=0,n-1)/)
end function linspace

    Remark: The function input areguments (here a, b, and n) are called dummy variables or dummy arguments, and linspace is the result variable. Note that the intent attribute can be specified for dummy variables only.

  2. The following code illustrates how we can work with fixed-size and allocatable arrays:

     1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
    integer, parameter :: n = 3
integer :: i
real(kind = 8) :: A(1:3,-2:0), B(n,n)
real(kind = 8), dimension(:,:), allocatable :: C

allocate(C(0:n-1,0:n-1))

A = 1.d0
B = 2.d0
B(3,3) = 0.d0
C = A+B
do i = 0,n-1
   write(*,*) C(i,:)
end do
deallocate(C)

    Remark 1: The arrays A and B are fixed-size, while C is allocatable.

    Remark 2: We can index arrays as we please. Here, for example, the rows of A are indexed from 1 to 3, while its columns are indexed from -2 to 0. The array B is indexed starting at 1 by default. The array C is indexed from 0 to 2.

    Remark3: In line 6, we can alternatively write allocate(C(n,n)). If we do so, then we will need to use the default indexing from 1 to n in line 12.

    Remark4: Do not forgoet to deallocate the memory (as done here in line 15). You also need to deallocate an array before making another allocation to that array. In general, whenever you need to make multiple allocations, before any allocation a good practice is to deallocate the array in case it is already allocated:

    if (allocated(C)) deallocate(C)
allocate(C(5,5))

Useful gfortran flags

gfortran has a variety of flags (or command line options) that control what the compiler does. Here, I list a few useful ones. For more information type the shell command:

$ man gfortran
  1. Output flags: We have already seen the two output flags -c and -o that control what kind of output gfortran generates:
    • -c compiles (multiple source files) to an object file. It is particularly useful if the source code is split into multiple files.
    • -o specifies the name of the output file.
  2. Warning flags: They tell gfortran to warn you about questionable sections of code. Warnings will often identify bugs:
    • -Wall warns about all. It tells gfortran to warn about many common sources of bugs, for example, when you have a function/subroutine with the same name as a built-in function, or when you pass the same variable both as an intent(in) and an intent(out) argument of the same subroutine.
    • -Wextra warns about more potential problems, for example, when a subroutine argument is never used.
  3. Debugging flags: They tell the compiler to include information inside the compiled program to help find bugs:
    • -fbacktrace produces a backtrace, showing what functions or subroutines were being called at the time of the error, if the program crashes.
    • -fbounds-check checks that the array index is within the bounds of the array every time an array element is accessed. It is useful for finding bugs related to arrays.
    • -ffpe-trap=LIST traps the listed floating point exceptions (fpe). It will send a signal and interrupts the program, producing a core file useful for debugging. LIST is a comma-separated list of the following IEEE exceptions:
      • invalid: invalid floating point operation. Exmples include sqrt(-1.0) or log(-1.0) or (0.0/0.0).
      • zero: division by zero. The code will die if you divide by zero. Otherwise, it may set the result to +INFINITY and continue.
      • overflow: overflow in a floating point operation. The code will halt if the value of a variable exceeds the largest value that can be represented by the corresponding data type.
      • underflow: underflow in a floating point operation. the code will halt if you compute a number that is too small because the exponent is a very large negative number. If you don’t trap underflows, such numbers will just be set to 0, which is generally the correct thing to do. But computing such small numbers may indicate a bug of some sort in the program, so you might want to trap them.
  4. Optimization flags: They control how the compiler optimizes your code:
    • -O0: No optimization (the default).
    • -O1, O2, O3: Higher levels usually produce faster code but take longer to compile.

Timing Fortran codes

There are two useful built-in subroutines in Fortran: cpu_time and system_clock.

cpu_time returns a REAL value representing the elapsed CPU time in seconds. We can use it to determine execution time for segments of code. For example:

program test_cpu_time
  real(kind=8) :: start, finish
  real(kind=8) :: elapsed_time
  call cpu_time(start)

  ! <code segment to be timed>

  call cpu_time(finish)
  elapsed_time = finish-start
  print '("Time = ",f6.3," seconds.")', elapsed_time
end program test_cpu_time

system_clock can also be used to measure the elapsed time between two successive calls. For example:

program test_system_clock
  integer(kind=8) :: tclock1, tclock2, clock_rate
  real(kind=8) :: elapsed_time
  call system_clock(tclock1)

  ! <code segment to be timed>

  call system_clock(tclock2, clock_rate)
  elapsed_time = float(tclock2 - tclock1) / float(clock_rate)
  print '("Time = ",f6.3," seconds.")', elapsed_time
end program test_system_clock

Exercise. Write a Fortran code and measure the CPU time for \(n \times n\) matrix-matrix multiplication. Try \(n=100\), \(n=200\), and \(n=400\), and check if the CPU time agrees with the theoretical \({\mathcal O}(n^3)\) flops. What happens if you set \(n=800\)? Repeat with adding the -O3 optimization flag and see how much your codes gets speed up.

You can find the code for this exercise (called TESTtiming.f90) in the course repository under Labs/lab2. Try to write your own code and then compare it with my code.

Modules

A module can be viewed as a library containing different things, such as a list of functions and subroutines that you often use, or a list of constants and variables that are common to many routines (global variables) and need to be shared by all these routines. With modules we can define “global variables” in one program unit and use them in another, without the need to pass the variable as a subroutine or function argument. Modules are most useful for organizing and structuring a program and linking routines across many different programs.

The structure of a module:

module <MODULENAME>
  ! declare variables

contains
  ! define subroutines and functions

end module <MODULENAME>

In the main program, a module can be used as:

program <PROGRAMNAME>
  use <MODULENAME>
  ! declare variables

  ! put the body of the program

end program <PROGRAMNAME>

Modules must be compiled before any program units that use the module. Compiling a module will create both a .o file and a .mod file to be used in executing the program. Here is how we can compile and execute a program that uses a module:

$gfortran -c modulename.f90 programname.f90
$gfortran -o main.x modulename.o programname.o
$./main.x

When declaring module variables, the statement save can be used to indicate that variables defined in the module should have values saved between one use of the module to another. In general, you should use this when declaring module variables.

As an example, in Homework4 we have a module which contains a variables and two functions:

module xycoord
  real(kind = 8), parameter :: pi = acos(-1.d0)
  save

  contains

  real(kind=8) function x_coord(r,s)
    implicit none
    real(kind=8) r,s
    x_coord = (2.d0+r+0.2*sin(5.d0*pi*s))*cos(0.5d0*pi*s)
  end function x_coord

  real(kind=8) function y_coord(r,s)
    implicit none
    real(kind=8) r,s
    y_coord = (2.d0+r+0.2*sin(5.d0*pi*s))*sin(0.5d0*pi*s)
  end function y_coord

end module xycoord