Fortran linear algebra
last modified November 14, 2026
A matrix is a rectangular table of numbers, and linear algebra is the set of operations on such tables: multiply them, transpose them, solve the system of equations they describe, invert them, find the numbers they scale. Fortran was designed for this work, and its array syntax makes the short operations one line each. The long ones, Gaussian elimination and eigenvalue iteration, are still loops, but short loops that take a page each.
This article multiplies matrices with matmul and verifies the
result against explicit loops, solves a three-by-three system by elimination
with partial pivoting, computes determinants and inverses, calls the LAPACK
routine dgesv for the same system, and estimates the dominant
eigenvalue of a matrix by repeated multiplication. Along the way it explains why
the numbers of a matrix fill a Fortran array column by column.
Matrices, shapes and matmul
A matrix is a two-dimensional array, and the intrinsic matmul
multiplies two of them when their inner dimensions agree: an
m-by-n times an n-by-p gives
an m-by-p. transpose swaps the rows and
columns, shape returns the two extents, and
dot_product reduces two vectors to one number. A matrix-vector
product is a matmul in which the second argument is
one-dimensional.
program j1
use iso_fortran_env, only: real64
implicit none
real(real64) :: a(2, 3), b(3, 2), v(3), w(2)
real(real64) :: product(2, 2), vector(2), diagonal
integer :: i
a = reshape([1.0_real64, 2.0_real64, 3.0_real64, 4.0_real64, &
5.0_real64, 6.0_real64], [2, 3])
b = reshape([1.0_real64, 0.0_real64, 0.0_real64, 1.0_real64, &
2.0_real64, 1.0_real64], [3, 2])
v = [1.0_real64, 2.0_real64, 3.0_real64]
w = [4.0_real64, 0.0_real64]
print '(a, 2i3)', 'shape of a: ', shape(a)
print '(a, i0, a, i0)', 'size: ', size(a), ' rows: ', size(a, 1)
call show('a', a)
call show('transpose(a)', transpose(a))
product = matmul(a, b)
call show('matmul(a, b)', product)
vector = matmul(a, v)
print '(a, 2f8.2)', 'matmul(a, v): ', vector
print '(a, f8.2)', 'dot_product(v, v): ', dot_product(v, v)
print '(a, f8.2)', 'dot_product(v, [2,2,2]): ', &
dot_product(v, [2.0_real64, 2.0_real64, 2.0_real64])
print '(a, 3f8.2)', 'a(1, :) is ', a(1, :)
print '(a, 2f8.2)', 'a(:, 2) is ', a(:, 2)
diagonal = 0.0_real64
do i = 1, min(size(a, 1), size(a, 2))
diagonal = diagonal + a(i, i)
end do
print '(a, f8.2)', 'the sum of the diagonal: ', diagonal
print '(a, f8.2)', 'a(1,1) plus a(2,2): ', a(1, 1) + a(2, 2)
contains
subroutine show(label, m)
character(len=*), intent(in) :: label
real(real64), intent(in) :: m(:, :)
integer :: row
print '(a, 2i3)', label//':', shape(m)
do row = 1, size(m, 1)
print '(8f8.2)', m(row, :)
end do
end subroutine show
end program j1
The array constructor is the part that surprises newcomers. Fortran stores a
matrix column by column, and reshape fills it in that order, so the
six values 1 to 6 become the rows [1, 3, 5] and
[2, 4, 6] rather than two rows of three consecutive numbers.
Column-major order is also why a slice such as a(1, :) is a row
across memory and a(:, 2) is a contiguous column. The product with
a two-by-two matrix gives a two-by-two result, the product with the vector
v gives a vector, and the dot products reduce to single numbers.
$ gfortran j1.f90 -o j1
$ ./j1
shape of a: 2 3
size: 6 rows: 2
a: 2 3
1.00 3.00 5.00
2.00 4.00 6.00
transpose(a): 3 2
1.00 2.00
3.00 4.00
5.00 6.00
matmul(a, b): 2 2
1.00 12.00
2.00 16.00
matmul(a, v): 22.00 28.00
dot_product(v, v): 14.00
dot_product(v, [2,2,2]): 12.00
a(1, :) is 1.00 3.00 5.00
a(:, 2) is 3.00 4.00
the sum of the diagonal: 5.00
a(1,1) plus a(2,2): 5.00
Checking the built-in product
A library intrinsic is only as trustworthy as the test that surrounds it. The
definition of the matrix product is a triple loop, c(i,j) = sum over k of
a(i,k)*b(k,j), and writing that loop once by hand confirms that
matmul computes what the definition says, that multiplying by the
identity changes nothing, and that the transpose of a product is the product of
the transposes in the opposite order.
program j2
use iso_fortran_env, only: real64
implicit none
real(real64) :: a(3, 2), b(2, 4), c(3, 4), manual(3, 4), eye(2, 2)
real(real64) :: checked(3, 2), transposed(4, 3)
integer :: i, j, k
a = reshape([1.0_real64, 4.0_real64, -2.0_real64, 3.0_real64, &
0.0_real64, 5.0_real64], [3, 2])
b = reshape([2.0_real64, 1.0_real64, 0.0_real64, -1.0_real64, &
3.0_real64, 2.0_real64, 1.0_real64, 4.0_real64], [2, 4])
c = matmul(a, b)
do j = 1, 4
do i = 1, 3
manual(i, j) = 0.0_real64
do k = 1, 2
manual(i, j) = manual(i, j) + a(i, k) * b(k, j)
end do
end do
end do
print '(a)', 'c = a * b:'
call show(c)
print '(a, es10.2)', 'the largest difference from the loops: ', &
maxval(abs(c - manual))
eye = 0.0_real64
eye(1, 1) = 1.0_real64
eye(2, 2) = 1.0_real64
checked = matmul(a, eye)
print '(a, es10.2)', 'a times the identity differs by: ', &
maxval(abs(checked - a))
transposed = transpose(matmul(a, b))
print '(a, es10.2)', '(ab)^T against b^T a^T: ', &
maxval(abs(transposed - matmul(transpose(b), transpose(a))))
print '(a, es10.2)', 'a(ab)^T entry check at (3,2): ', &
abs(transposed(3, 2) - c(2, 3))
contains
subroutine show(m)
real(real64), intent(in) :: m(:, :)
integer :: row
do row = 1, size(m, 1)
print '(8f8.2)', m(row, :)
end do
end subroutine show
end program j2
Every difference printed is exactly zero, because all the entries are small integers and every product and sum fits a double without rounding. With fractional entries the differences would be of the order of 1e-16 instead of zero, and the test would need a tolerance rather than an equality. The three identities checked here are the ones worth keeping in a test suite: agreement with the definition, the identity matrix, and the transpose rule.
$ gfortran j2.f90 -o j2
$ ./j2
c = a * b:
5.00 -3.00 9.00 13.00
8.00 0.00 12.00 4.00
1.00 -5.00 4.00 18.00
the largest difference from the loops: 0.00E+00
a times the identity differs by: 0.00E+00
(ab)^T against b^T a^T: 0.00E+00
a(ab)^T entry check at (3,2): 0.00E+00
Gaussian elimination with pivoting
Solving A x = b means finding the vector that the matrix maps to
b. Gaussian elimination subtracts multiples of one equation
from the others until the matrix is triangular, and then solves the triangle
from the bottom up by back substitution. Partial pivoting
swaps in the row with the largest entry in the current column before dividing,
because dividing by a small number amplifies rounding errors; the price is a
record of the row order, and the benefit is a method that survives most
well-conditioned systems.
module linalg
use iso_fortran_env, only: real64
implicit none
private
public :: solve, determinant, inverse
real(real64), parameter :: tiny = 1.0e-12_real64
contains
subroutine solve(a, b, x, ok)
real(real64), intent(in) :: a(:, :), b(:)
real(real64), intent(out) :: x(:)
logical, intent(out) :: ok
real(real64) :: m(size(a, 1), size(a, 2)), rhs(size(b))
real(real64) :: factor, total, biggest, swap
integer :: n, i, j, k, pivot
n = size(a, 1)
m = a
rhs = b
ok = .true.
do k = 1, n - 1
pivot = k
biggest = abs(m(k, k))
do i = k + 1, n
if (abs(m(i, k)) > biggest) then
biggest = abs(m(i, k))
pivot = i
end if
end do
if (biggest < tiny) then
ok = .false.
x = 0.0_real64
return
end if
if (pivot /= k) then
do j = 1, n
swap = m(k, j)
m(k, j) = m(pivot, j)
m(pivot, j) = swap
end do
swap = rhs(k)
rhs(k) = rhs(pivot)
rhs(pivot) = swap
end if
do i = k + 1, n
factor = m(i, k) / m(k, k)
m(i, k:n) = m(i, k:n) - factor * m(k, k:n)
rhs(i) = rhs(i) - factor * rhs(k)
end do
end do
if (abs(m(n, n)) < tiny) then
ok = .false.
x = 0.0_real64
return
end if
do i = n, 1, -1
total = rhs(i)
do j = i + 1, n
total = total - m(i, j) * x(j)
end do
x(i) = total / m(i, i)
end do
end subroutine solve
real(real64) function determinant(a) result(det)
real(real64), intent(in) :: a(:, :)
real(real64) :: m(size(a, 1), size(a, 2)), factor, swap
integer :: n, i, j, k, pivot, flips
n = size(a, 1)
m = a
flips = 0
det = 0.0_real64
do k = 1, n - 1
pivot = k
do i = k + 1, n
if (abs(m(i, k)) > abs(m(pivot, k))) pivot = i
end do
if (abs(m(pivot, k)) < tiny) return
if (pivot /= k) then
do j = 1, n
swap = m(k, j)
m(k, j) = m(pivot, j)
m(pivot, j) = swap
end do
flips = flips + 1
end if
do i = k + 1, n
factor = m(i, k) / m(k, k)
m(i, k:n) = m(i, k:n) - factor * m(k, k:n)
end do
end do
if (abs(m(n, n)) < tiny) return
det = 1.0_real64
do i = 1, n
det = det * m(i, i)
end do
if (mod(flips, 2) == 1) det = -det
end function determinant
subroutine inverse(a, inv, ok)
real(real64), intent(in) :: a(:, :)
real(real64), intent(out) :: inv(:, :)
logical, intent(out) :: ok
real(real64) :: m(size(a, 1), size(a, 2)), unit(size(a, 1), size(a, 2))
real(real64) :: factor, swap, biggest
integer :: n, i, j, k, pivot
n = size(a, 1)
m = a
unit = 0.0_real64
do i = 1, n
unit(i, i) = 1.0_real64
end do
ok = .true.
do k = 1, n
pivot = k
biggest = abs(m(k, k))
do i = k + 1, n
if (abs(m(i, k)) > biggest) then
biggest = abs(m(i, k))
pivot = i
end if
end do
if (biggest < tiny) then
ok = .false.
inv = 0.0_real64
return
end if
if (pivot /= k) then
do j = 1, n
swap = m(k, j)
m(k, j) = m(pivot, j)
m(pivot, j) = swap
swap = unit(k, j)
unit(k, j) = unit(pivot, j)
unit(pivot, j) = swap
end do
end if
factor = m(k, k)
m(k, :) = m(k, :) / factor
unit(k, :) = unit(k, :) / factor
do i = 1, n
if (i /= k) then
factor = m(i, k)
m(i, :) = m(i, :) - factor * m(k, :)
unit(i, :) = unit(i, :) - factor * unit(k, :)
end if
end do
end do
inv = unit
end subroutine inverse
end module linalg
The subroutine copies its arguments, because an intent(in) array
must not be modified by the elimination. The outer loop runs over the columns,
finds the pivot row, swaps it into place, and subtracts the scaled pivot row
from every row below. A pivot smaller than 1e-12 in absolute value is treated as
a zero column, and the routine returns ok = false instead of
dividing by it. The three public routines, solve,
determinant and inverse, share the same test.
program j3
use iso_fortran_env, only: real64
use linalg
implicit none
real(real64) :: a(3, 3), b(3), x(3), residual(3)
logical :: ok
a = reshape([2.0_real64, -3.0_real64, -2.0_real64, 1.0_real64, &
-1.0_real64, 1.0_real64, -1.0_real64, 2.0_real64, 2.0_real64], &
[3, 3])
b = [8.0_real64, -11.0_real64, -3.0_real64]
call solve(a, b, x, ok)
print '(a, l1)', 'solved: ', ok
print '(a, 3f10.6)', 'the solution: ', x
residual = matmul(a, x) - b
print '(a, 3es11.2)', 'the residual: ', residual
print '(a, es10.2)', 'the largest part: ', maxval(abs(residual))
end program j3
The three equations have the solution 2, 3 and -1, which the routine returns to
six decimals. The residual, computed as A x - b with
matmul, is zero in two components and 8.9e-16 in the third, which
is the size of a rounding error in double precision for numbers of this
magnitude. A residual near zero does not prove that the answer is accurate, but
a large one proves that it is not, which makes this a cheap check to keep in
any program that solves a system.
$ gfortran linalg.f90 j3.f90 -o j3 $ ./j3 solved: T the solution: 2.000000 3.000000 -1.000000 the residual: 0.00E+00 0.00E+00 8.88E-16 the largest part: 8.88E-16
The determinant
The determinant is a single number attached to a square matrix. It is zero exactly when the matrix is singular, that is when its equations have no unique solution, and its absolute value is the factor by which the matrix stretches areas or volumes. Elimination computes it almost for free: the determinant of a triangular matrix is the product of its diagonal, and each row swap flips the sign, so the determinant of the original matrix is that product times minus one to the number of swaps.
program j4
use iso_fortran_env, only: real64
use linalg
implicit none
real(real64) :: a(3, 3), flat(2, 2), singular(2, 2), identity(3, 3)
integer :: i
a = reshape([2.0_real64, -3.0_real64, -2.0_real64, 1.0_real64, &
-1.0_real64, 1.0_real64, -1.0_real64, 2.0_real64, 2.0_real64], &
[3, 3])
flat = reshape([3.0_real64, 5.0_real64, 2.0_real64, 4.0_real64], [2, 2])
singular = reshape([1.0_real64, 2.0_real64, 2.0_real64, 4.0_real64], &
[2, 2])
identity = 0.0_real64
do i = 1, 3
identity(i, i) = 1.0_real64
end do
print '(a, f9.4)', 'the 3 by 3 system matrix: ', determinant(a)
print '(a, f9.4)', 'the 2 by 2: ', determinant(flat)
print '(a, f9.4)', 'by the formula: ', 3.0_real64 * 4.0_real64 - &
2.0_real64 * 5.0_real64
print '(a, f9.4)', 'the singular matrix: ', determinant(singular)
print '(a, f9.4)', 'the identity: ', determinant(identity)
print '(a, f9.4)', 'doubling every entry: ', &
determinant(2.0_real64 * flat)
print '(a, f9.4)', 'four times the original: ', 4.0_real64 * &
determinant(flat)
print '(a, f9.4)', 'the transpose: ', &
determinant(transpose(flat))
print '(a, f9.4)', 'swapping two rows: ', &
determinant(reshape([2.0_real64, 4.0_real64, 3.0_real64, &
5.0_real64], [2, 2]))
end program j4
The three-by-three system matrix has determinant -1, the two-by-two matrix has 2,
the same value as the formula ad - bc, and the matrix with a row
proportional to the other has 0. Doubling every entry of a two-by-two multiplies
the determinant by four, because each of the two dimensions is stretched by
two, and the transpose leaves it unchanged. The last line swaps the two rows and
the sign follows, from 2 to -2.
$ gfortran linalg.f90 j4.f90 -o j4 $ ./j4 the 3 by 3 system matrix: -1.0000 the 2 by 2: 2.0000 by the formula: 2.0000 the singular matrix: 0.0000 the identity: 1.0000 doubling every entry: 8.0000 four times the original: 8.0000 the transpose: 2.0000 swapping two rows: -2.0000
The inverse
The inverse of a matrix is the matrix that undoes it, so that A
times its inverse is the identity. Gauss-Jordan elimination computes it
by applying the same row operations to the matrix and to the identity at once:
when the left side has become the identity, the right side has become the
inverse. The work is about twice that of a single solution, and the routine
refuses a singular matrix in the same way as the solver.
program j5
use iso_fortran_env, only: real64
use linalg
implicit none
real(real64) :: a(3, 3), inv(3, 3), product(3, 3), singular(2, 2)
real(real64) :: bad(2, 2)
logical :: ok
a = reshape([2.0_real64, -3.0_real64, -2.0_real64, 1.0_real64, &
-1.0_real64, 1.0_real64, -1.0_real64, 2.0_real64, 2.0_real64], &
[3, 3])
call inverse(a, inv, ok)
print '(a, l1)', 'inverted: ', ok
print '(a)', 'the inverse:'
call show(inv)
product = matmul(a, inv)
print '(a)', 'a times its inverse:'
call show(product)
print '(a, es10.2)', 'the largest departure from one: ', &
maxval(abs(product - identity(3)))
singular = reshape([1.0_real64, 2.0_real64, 2.0_real64, 4.0_real64], &
[2, 2])
call inverse(singular, bad, ok)
print '(a, l1)', 'a singular matrix: ', ok
print '(a, es10.2)', 'and its inverse comes out as ', maxval(abs(bad))
contains
subroutine show(m)
real(real64), intent(in) :: m(:, :)
integer :: row
do row = 1, size(m, 1)
print '(8f10.3)', m(row, :)
end do
end subroutine show
function identity(n) result(unit)
integer, intent(in) :: n
real(real64) :: unit(n, n)
integer :: k
unit = 0.0_real64
do k = 1, n
unit(k, k) = 1.0_real64
end do
end function identity
end program j5
The inverse is printed to three decimals and the product with the original is the identity to within 1.8e-15, which is the rounding error of the whole elimination. The smallest entry of the product prints as -0.000 because the computed value is a tiny negative number, not exactly zero; a program that compares such an entry with zero would be wrong, and a program that compares its absolute value with a tolerance is right. The singular matrix is refused, and the output array is filled with zeros rather than with meaningless numbers.
$ gfortran linalg.f90 j5.f90 -o j5
$ ./j5
inverted: T
the inverse:
4.000 3.000 -1.000
-2.000 -2.000 1.000
5.000 4.000 -1.000
a times its inverse:
1.000 0.000 0.000
0.000 1.000 -0.000
0.000 0.000 1.000
the largest departure from one: 1.78E-15
a singular matrix: F
and its inverse comes out as 0.00E+00
Calling LAPACK
Nobody writes production linear algebra by hand: the reference is
LAPACK, a library of Fortran routines that has been tuned for decades,
and it is installed on most Linux systems together with its BLAS
building blocks. The routine for a general system in double precision is
dgesv, which takes the matrix, the right-hand side and an integer
array for the pivot rows, and overwrites both arrays with the factors and the
solution. Its name encodes the types: d for double precision,
ge for a general matrix, sv for solve.
program j6
use iso_fortran_env, only: real64
implicit none
real(real64) :: a(3, 3), original(3, 3), b(3, 1), rhs(3)
integer :: pivots(3), info
interface
subroutine dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)
import :: real64
integer, intent(in) :: n, nrhs, lda, ldb
integer, intent(out) :: info, ipiv(*)
real(real64), intent(inout) :: a(lda, *), b(ldb, *)
end subroutine dgesv
end interface
a = reshape([2.0_real64, -3.0_real64, -2.0_real64, 1.0_real64, &
-1.0_real64, 1.0_real64, -1.0_real64, 2.0_real64, 2.0_real64], &
[3, 3])
original = a
rhs = [8.0_real64, -11.0_real64, -3.0_real64]
b(:, 1) = rhs
call dgesv(3, 1, a, 3, pivots, b, 3, info)
print '(a, i0)', 'info: ', info
print '(a, 3f10.6)', 'the solution: ', b(:, 1)
print '(a, 3es11.2)', 'the residual: ', matmul(original, b(:, 1)) - rhs
print '(a, 3i4)', 'the pivot rows: ', pivots
print '(a)', 'the lu factors left in the array:'
call show(a)
contains
subroutine show(m)
real(real64), intent(in) :: m(:, :)
integer :: row
do row = 1, size(m, 1)
print '(3f10.4)', m(row, :)
end do
end subroutine show
end program j6
The interface block is required, because the compiler must know the types of
every argument of an external routine whose source it cannot see. A zero
info means success; a positive value would name the column that
turned out to be singular. The solution matches the hand-written solver, the
residual is again at the 1e-16 level, and the pivot rows 2, 3, 3 record the
swaps that the routine made. The array a no longer holds the
original matrix: it holds the lower and upper triangular factors, which is why
the residual is computed from the copy saved before the call.
$ gfortran j6.f90 -o j6 -llapack -lblas
$ ./j6
info: 0
the solution: 2.000000 3.000000 -1.000000
the residual: -8.88E-16 0.00E+00 8.88E-16
the pivot rows: 2 3 3
the lu factors left in the array:
-3.0000 -1.0000 2.0000
0.6667 1.6667 0.6667
-0.6667 0.2000 0.2000
Eigenvalues by power iteration
An eigenvector of a matrix is a vector that the matrix only scales, and the factor is its eigenvalue. Multiplying any vector repeatedly by the matrix stretches it more along the direction of the largest eigenvalue than along any other, so after a few steps the direction of the product is dominated by that eigenvector. Power iteration normalises the product at every step to keep the numbers in range and estimates the eigenvalue with the Rayleigh quotient, the dot product of the vector with its image divided by the squared length.
program j7
use iso_fortran_env, only: real64
implicit none
real(real64) :: a(2, 2), v(2), w(2), lambda, previous, exact(2)
integer :: k
a = reshape([4.0_real64, 2.0_real64, 1.0_real64, 3.0_real64], [2, 2])
v = [1.0_real64, 0.0_real64]
previous = 0.0_real64
print '(a)', ' k rayleigh change vector'
do k = 1, 9
w = matmul(a, v)
v = w / norm2(w)
lambda = dot_product(v, matmul(a, v))
print '(i3, f13.9, es11.3, 2f9.5)', k, lambda, abs(lambda - previous), v
previous = lambda
end do
exact(1) = 0.5_real64 * (7.0_real64 + sqrt(7.0_real64**2 - 40.0_real64))
exact(2) = 0.5_real64 * (7.0_real64 - sqrt(7.0_real64**2 - 40.0_real64))
print '(a, 2f8.4)', 'the exact eigenvalues: ', exact
print '(a, f8.4)', 'the trace: ', a(1, 1) + a(2, 2)
print '(a, f8.4)', 'the determinant: ', determinant(a)
contains
real(real64) function determinant(m) result(d)
real(real64), intent(in) :: m(2, 2)
d = m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)
end function determinant
end program j7
The matrix has eigenvalues 5 and 2, which the trace and the determinant confirm: their sum is 7 and their product is 10. Starting from the vector (1, 0), the estimate climbs from 5.08 to 5.0002 in nine steps and the vector approaches (0.7071, 0.7071), the normalised eigenvector. The change falls by a factor that approaches 2.5, the ratio of the larger eigenvalue to the smaller, which is why the method converges quickly when the eigenvalues are separated and painfully slowly when they are close.
$ gfortran j7.f90 -o j7 $ ./j7 k rayleigh change vector 1 5.000000000 5.000E+00 0.89443 0.44721 2 5.076923077 7.692E-02 0.78935 0.61394 3 5.041543027 3.538E-02 0.74072 0.67181 4 5.018196664 2.335E-02 0.72064 0.69331 5 5.007521409 1.068E-02 0.71253 0.70164 6 5.003046751 4.475E-03 0.70928 0.70493 7 5.001224768 1.822E-03 0.70798 0.70624 8 5.000490875 7.339E-04 0.70745 0.70676 9 5.000196505 2.944E-04 0.70725 0.70697 the exact eigenvalues: 5.0000 2.0000 the trace: 7.0000 the determinant: 10.0000
Next steps
The arrays article explains the column-major
order and the sections used here, and the
array intrinsics article covers
sum, maxval and reshape. The
array arguments article describes the
assumed-shape arguments m(:, :) of the printing routine, the
modules article explains the separate
compilation of linalg.f90, and the
numerical integration article
applies the same care about tolerances and convergence to a different problem.
The floating-point article explains the
1e-16 residuals that appear in place of exact zeros.
Source
Fortran 2018 standard (J3/18-007r1)
LAPACK: Linear Algebra PACKage
In this article we have built matrices with reshape and seen the
column-major order, multiplied them with matmul and verified the
result against triple loops, solved a system by Gaussian elimination with
partial pivoting and checked the residual, computed determinants and their
properties, inverted a matrix with Gauss-Jordan elimination and detected a
singular one, called LAPACK's dgesv and read its pivot array, and
found the dominant eigenvalue by power iteration. The previous article,
Fortran numerical integration and root
finding, uses the same convergence tests on scalar problems.
Author
List all Fortran tutorials.