Fortran playground

Fortran — the FORmula TRANslator — came out of John Backus’s team at IBM in 1957 and was the first high-level language that produced code as fast as hand-written assembly. It has never stopped being the language of big numerical codes: weather and climate models, computational chemistry, and the LAPACK routines underneath NumPy and R are Fortran. Modern Fortran (90 onwards) is free form, has modules, and above all treats arrays as values — a = b + c adds whole arrays.

These blocks are compiled by LFortran, a new LLVM-based Fortran compiler, itself compiled to WebAssembly (about 18 MB on the first Run). LFortran is still in beta, so a few corners of the language are missing, but everything here compiles and runs. The stdin button supplies input for read *.

Arrays are values: whole-array arithmetic and intrinsics

program arrays
  implicit none
  real :: temps(7) = [14.2, 17.8, 21.5, 19.0, 12.6, 9.3, 23.1]
  real :: fahrenheit(7)
  integer :: i

  fahrenheit = temps * 9.0 / 5.0 + 32.0          ! every element at once, no loop
  print '(a, 7f6.1)', 'F:    ', fahrenheit
  print '(a, f6.2)',  'mean: ', sum(temps) / size(temps)
  print '(a, f6.1, a, i0)', 'max:  ', maxval(temps), ' on day ', maxloc(temps, 1)
  print '(a, i0)',    'days above 18: ', count(temps > 18.0)

  print '(a, 3f6.1)', 'mid-week: ', temps(3:5)                ! array sections
  print '(a, 4f6.1)', 'odd days: ', temps(1:7:2)              ! with a stride
  print '(a, 7f6.1)', 'reversed: ', temps(7:1:-1)
  print '(a, *(i4))', 'squares:  ', [(i * i, i = 1, 8)]       ! an implied-do constructor
end program arrays

Matrices: matmul, transpose, and reshape

program matrices
  implicit none
  real :: a(3, 3), b(3, 3), v(3)
  integer :: i

  ! Fortran is column-major: reshape fills down the columns first
  a = reshape([2., 1., 0., 1., 3., 1., 0., 1., 4.], [3, 3])
  v = [1., 2., 3.]

  print *, 'A:'
  do i = 1, 3
    print '(3f7.2)', a(i, :)                     ! row i: every column
  end do

  print '(a, 3f8.2)', 'A·v    = ', matmul(a, v)
  print '(a, f8.2)',  'v·v    = ', dot_product(v, v)
  b = matmul(a, transpose(a))
  print '(a, 3f8.2)', 'diag(AAᵀ) = ', (b(i, i), i = 1, 3)
  print '(a, l2)',    'A symmetric? ', all(a == transpose(a))
end program matrices

WHERE: masked assignment without a loop

program masks
  implicit none
  real :: x(10), y(10)
  integer :: i

  x = [(real(i) - 5.5, i = 1, 10)]
  where (x > 0.0)
    y = sqrt(x)
  elsewhere
    y = -1.0                                     ! mark the values sqrt cannot take
  end where

  print '(10f6.2)', x
  print '(10f6.2)', y
  print '(a, f6.3)', 'sum of the roots: ', sum(y, mask=y >= 0.0)
end program masks

Modules, derived types, and elemental functions

module geometry
  implicit none
  real, parameter :: pi = 4.0 * atan(1.0)

  type :: circle
    real :: x, y, r
  end type circle

contains

  elemental real function area(c)                ! elemental: works on one circle or an array of them
    type(circle), intent(in) :: c
    area = pi * c%r**2
  end function area

  logical function overlaps(a, b)
    type(circle), intent(in) :: a, b
    overlaps = hypot(a%x - b%x, a%y - b%y) < a%r + b%r
  end function overlaps

end module geometry

program shapes
  use geometry
  implicit none
  type(circle) :: cs(3)

  cs(1) = circle(0.0, 0.0, 1.0)
  cs(2) = circle(1.5, 0.0, 1.0)
  cs(3) = circle(5.0, 5.0, 2.0)

  print '(a, 3f8.3)', 'areas: ', area(cs)
  print '(a, f8.3)',  'total: ', sum(area(cs))
  print '(a, l2)',    '1 overlaps 2? ', overlaps(cs(1), cs(2))
  print '(a, l2)',    '1 overlaps 3? ', overlaps(cs(1), cs(3))
end program shapes

A numerical kernel: the heat equation on a rod

! A metal rod, 1 m long, starts at 20 °C; its left end is held at 100 °C and its right end at 0 °C.
! Explicit finite differences: each step, every interior point moves toward the mean of its neighbours.
program heat
  implicit none
  integer, parameter :: n = 21, steps = 2000
  real, parameter :: alpha = 1.0e-4, dx = 1.0 / (n - 1), dt = 0.4 * dx**2 / alpha
  real :: u(n)
  integer :: step

  u = 20.0
  u(1) = 100.0
  u(n) = 0.0

  do step = 1, steps
    u(2:n-1) = u(2:n-1) + alpha * dt / dx**2 * (u(1:n-2) - 2.0 * u(2:n-1) + u(3:n))
    if (step == 50 .or. step == 400 .or. step == steps) then
      print '(a, i5, a, f8.0, a)', 'step', step, '  (t =', step * dt, ' s)'
      print '(11f6.1)', u(1:n:2)
    end if
  end do

  ! At steady state the profile is a straight line from 100 to 0
  print '(a, es9.2)', 'max deviation from a straight line: ', &
        maxval(abs(u - [(100.0 * (1.0 - real(step - 1) / (n - 1)), step = 1, n)]))
end program heat

The whole update is one line: u(2:n-1) = … u(1:n-2) … u(3:n) shifts the array left and right and combines all three at once — the array syntax that lets Fortran compilers vectorise loops without being asked. The step size dt is chosen at 0.4 of the stability limit; try 0.6 and watch the solution blow up.

Reading input: Newton’s method for a square root

program newton
  implicit none
  real(8) :: a, x
  integer :: i

  read *, a                                      ! press stdin and type a number first
  x = a
  do i = 1, 8
    x = 0.5d0 * (x + a / x)
    print '(i2, f22.15)', i, x
  end do
  print '(a, es10.2)', 'error vs sqrt(): ', abs(x - sqrt(a))
end program newton