Skip to content
Numerical Linear Algebra
step 1/5

Reading — step 1 of 5

Learn

~2 min readOOP, Parallel, Numerical

Fortran's killer apps are scientific computing and numerical linear algebra. The standard library (intrinsic functions + LAPACK/BLAS) covers most of what scientists need.

Built-in array operations

real :: A(3, 3), B(3, 3), C(3, 3), v(3), x(3)

! Matrix multiply
C = matmul(A, B)

! Transpose
B = transpose(A)

! Dot product
x = dot_product(v, v)

! Norm
real :: r
r = sqrt(sum(v**2))

! Outer product (use spread)
C = spread(v, 2, 3) * spread(v, 1, 3)

Most array ops are vectorized — single instruction generates SIMD code.

LAPACK and BLAS

For solving linear systems, eigenvalues, decompositions — Fortran has the mature standard:

! LAPACK example: solve Ax = b for x
real :: A(3, 3), b(3)
integer :: ipiv(3), info

call sgesv(3, 1, A, 3, ipiv, b, 3, info)
if (info /= 0) print *, "singular!"
! b now contains x

sgesv solves a general dense system. There are similar routines for:

  • Symmetric: ssysv
  • Triangular: strsv
  • Banded: sgbsv
  • Sparse: sssm and friends

The naming convention: 1st letter: precision (s = single, d = double, c = complex, z = double complex) 2nd-3rd: matrix type (ge = general, sy = symmetric, tr = triangular) Rest: operation

Writing for performance

! Bad — Fortran is column-major; this is the fast direction:
do j = 1, n
    do i = 1, n
        A(i, j) = ...
    end do
end do

! Reversed — slow due to cache misses:
do i = 1, n
    do j = 1, n
        A(i, j) = ...
    end do
end do

C is row-major (last index changes fastest); Fortran is column-major (first index changes fastest). This affects loop order for cache.

Numerical types

  • real(kind=4) — single precision (32-bit)
  • real(kind=8) or double precision — double precision (64-bit)
  • real(kind=16) — quad precision (128-bit, gfortran supports)
  • complex(kind=8) — double precision complex

Random numbers

real :: x
call random_number(x)           ! uniform [0, 1)

real :: arr(100)
call random_number(arr)         ! fills the whole array

call random_seed()              ! seed from time

For reproducible results: call random_seed(put=fixed_seed_array).

Where Fortran shines

  • Climate and atmospheric models
  • Computational fluid dynamics (CFD)
  • Computational chemistry (DFT, ab-initio)
  • Lattice QCD (particle physics)
  • Many ML / numerical libraries' inner loops (PyTorch, NumPy use BLAS underneath)

Modern alternatives

For new scientific code:

  • Julia — modern, similar performance, much better ergonomics
  • C++ with Eigen / Kokkos — performance + flexibility, more complex
  • Python + NumPy — productivity, slower without Numba

Fortran isn't dying — it's the best tool for its niche, and rewriting decades of validated code is risky. New code may be elsewhere; existing Fortran will live for decades.

Discussion

Ask a question, share an insight, or help someone who’s stuck.

Sign in to post a comment or reply.

Loading…