====================================
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:
.. code-block:: none
$ gfortran newton.f90
This will produce an executable file named ``a.out``. To run the code we would then type:
.. code-block:: none
$ ./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:
.. code-block:: none
$ 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:
.. code-block:: none
$ 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.
.. code-block:: fortran
:linenos:
! 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:
.. code-block:: none
$ 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.
.. code-block:: fortran
:linenos:
! 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
:math:`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
.. code-block:: fortran
(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
.. code-block:: none
$ 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
:math:`a = | \pi - \pi | / | \pi | = 0` (when :math:`i=1`) and :math:`a = 2 \pi / 3 \neq 0` (when :math:`i=2`).
Example 3: more built-in functions
+++++++++++++++++++++++++++++++++++++++++++++++++++
This example shows more built-in functions.
.. code-block:: fortran
:linenos:
! 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:
* We can use the built-in function arc-cosine of -1 to comput
pi. Note that we need to write ``-1.d0`` for full precision.
* For a list of built-in functions see http://www.nsc.liu.se/~boein/f77to90/a5.html.
Executing the above code yields
.. code-block:: none
$ 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.
.. code-block:: fortran
:linenos:
! 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
.. code-block:: none
$ 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.
.. code-block:: fortran
:linenos:
! 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
.. code-block:: none
$ 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.
.. code-block:: fortran
:linenos:
! 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
.. code-block:: none
$ 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``.
.. code-block:: fortran
:linenos:
! 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.
.. code-block:: fortran
:linenos:
! 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.
.. code-block:: none
$ 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:
.. code-block:: fortran
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:
.. code-block:: fortran
:linenos:
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:
.. code-block:: fortran
:linenos:
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:
.. code-block:: fortran
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:
.. code-block:: none
$ 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:
.. code-block:: fortran
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:
.. code-block:: fortran
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 :math:`n \times n`
matrix-matrix multiplication. Try :math:`n=100`, :math:`n=200`, and :math:`n=400`, and
check if the CPU time agrees with the theoretical :math:`{\mathcal
O}(n^3)` flops. What happens if you set :math:`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:
.. code-block:: fortran
module <MODULENAME>
! declare variables
contains
! define subroutines and functions
end module <MODULENAME>
In the main program, a module can be used as:
.. code-block:: fortran
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:
.. code-block:: none
$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:
.. code-block:: fortran
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
__ https://en.wikipedia.org/wiki/Fortran
__ https://en.wikipedia.org/wiki/Compiler
__ https://en.wikipedia.org/wiki/Linker_%28computing%29
__ https://gcc.gnu.org/fortran/
__ https://www.gnu.org/
__ http://math.unm.edu/~motamed/Teaching/Fall20/HPSC/hw4.html