Numerical Computing in Fortran

Analogy categoryAnalogies

Precision and Kind Selection

Serious numerical work in Fortran starts with choosing floating-point precision deliberately rather than relying on a bare real, whose actual precision is compiler- and platform-dependent. The standard, portable approach is to declare a kind parameter with selected_real_kind, for example integer, parameter :: dp = selected_real_kind(15, 307), which requests at least 15 decimal digits of precision and an exponent range up to 10^307, then declare variables as real(dp) :: x throughout the program; this guarantees double-precision behavior on any standard-conforming compiler, unlike the common but non-portable shortcut of declaring real*8.

🏏

Cricket analogy: Like specifying an exact minimum ball weight and seam height in a match contract rather than just saying 'a normal ball,' so any certified manufacturer's ball meets the same precise standard regardless of country.

Array Operations and Intrinsic Functions

Fortran's whole-array syntax lets numerical code express vector and matrix operations directly without explicit loops: c = a + b adds two same-shaped arrays element-wise, and intrinsic functions like matmul, dot_product, sum, maxval, and transpose operate on entire arrays in one expression. Beyond correctness and readability, this array syntax gives the compiler's optimizer much more freedom than an equivalent hand-written loop, since the compiler can see the whole operation's shape and dependencies at once and is free to vectorize it using SIMD instructions or call an optimized BLAS routine under the hood.

🏏

Cricket analogy: Like a captain declaring a single blanket field placement instruction ('all fielders move two steps in') rather than radioing each fielder individually, letting the whole team reposition in one coordinated motion that's faster to execute.

fortran

program numerical_demo
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: n = 3
  real(dp) :: A(n,n), b(n), c(n)
  integer :: i

  A = reshape([1.0_dp, 2.0_dp, 0.0_dp, &
               0.0_dp, 1.0_dp, 3.0_dp, &
               4.0_dp, 0.0_dp, 1.0_dp], [n,n])
  b = [1.0_dp, 2.0_dp, 3.0_dp]

  c = matmul(A, b)          ! whole-array matrix-vector product

  do i = 1, n
     print '(A,I0,A,F10.6)', 'c(', i, ') = ', c(i)
  end do

  print *, 'sum(c) =', sum(c), '  maxval(c) =', maxval(c)
end program numerical_demo

Interfacing with LAPACK and BLAS

For serious linear algebra — solving dense linear systems, eigenvalue problems, or least-squares fits — production Fortran code almost always calls the battle-tested LAPACK and BLAS libraries rather than hand-writing a solver, since these libraries have been tuned over decades for numerical stability and cache-efficient performance. A typical call such as call dgesv(n, nrhs, A, lda, ipiv, b, ldb, info) solves the dense linear system A*x = b in place, overwriting b with the solution, and checking the returned info argument afterward is mandatory since a nonzero value signals a singular or near-singular matrix rather than a successful solve.

🏏

Cricket analogy: Like a franchise hiring a specialist bowling coach with decades of proven international experience rather than having a rookie player invent bowling technique from scratch, because the specialist's tuned methods are more reliable under match pressure.

Always check the info output argument after calling a LAPACK routine like dgesv or dsyev. A nonzero info does not raise a Fortran runtime error automatically — the routine returns quietly, and ignoring info can let a program silently proceed with garbage results from a singular matrix or a failed eigenvalue decomposition.

Numerical Stability and Floating-Point Pitfalls

Numerical code must account for the limits of finite-precision floating-point arithmetic: subtracting two nearly equal large numbers (catastrophic cancellation) can destroy most of the significant digits in the result, and comparing two reals for exact equality with == is almost always wrong because rounding error makes bitwise equality unreliable even when two values are mathematically identical. Fortran provides the intrinsic epsilon(x) to query the smallest representable relative gap for a given real kind, which is the standard building block for writing a safe tolerance-based comparison such as abs(a - b) <= epsilon(a) * max(abs(a), abs(b)) * tolerance_factor instead of a == b.

🏏

Cricket analogy: Like calculating a batter's exact strike rate by subtracting two very close ball-by-ball run totals late in an innings, where small scoring errors early on can compound into a materially wrong final rate if not tracked carefully.

epsilon(x) returns the machine epsilon appropriate to the kind of x (a much larger value for single precision than for double precision), so tolerance-based comparisons automatically adapt if you later change a kind parameter like dp from single to double precision, without needing to hunt down and update hard-coded tolerance constants throughout the code.

Use selected_real_kind (or the iso_fortran_env kinds real64/real32) to declare precision portably instead of real*8.
Whole-array syntax (a + b, matmul, dot_product) is both more readable and often better-optimized than manual loops.
Production linear algebra should call LAPACK/BLAS routines (e.g. dgesv) rather than hand-rolled solvers.
Always check the info output argument from a LAPACK call; a nonzero value signals failure without raising an error.
Catastrophic cancellation from subtracting nearly equal large numbers can destroy most of a result's significant digits.
Never compare reals with == for equality; use a tolerance built from epsilon(x) instead.
epsilon(x) returns the machine precision for x's specific kind, keeping tolerance-based comparisons portable across kinds.