Chapter 5: Arrays — Fortran's Superpower

"It is better to have 100 functions operate on one data structure than 10 functions on 10 data structures." — Alan J. Perlis, Epigrams on Programming (1982)

Overview

This is the chapter the rest of the book has been waiting for. Everything you have learned so far — types and precision in Chapter 3, loops and decisions in Chapter 4 — has been about single numbers, one at a time. But a climate model does not compute one temperature; it computes a hundred million of them, over and over. A linear solver does not add two numbers; it adds two vectors of a billion. The unit of scientific computing is not the scalar. It is the array, and the reason Fortran is fast — the reason it outlives every language that was supposed to replace it — is that Fortran understands arrays in a way most languages do not.

Perlis's epigram is the whole chapter in one line. An array is a single data structure, and Fortran gives you a hundred operations that act on it as a whole: add two arrays, take their dot product, sum them, multiply matrices, slice out a diagonal — each a single expression, no loop in sight. You will write code that reads almost exactly like the mathematics it implements, and — this is the part that matters — that same code hands the compiler everything it needs to make the computation run at the speed of the hardware. Readable and fast, at the same time, with no trade-off. That is not normal, and by the end of this chapter you will understand exactly why Fortran gets to have both.

We will also meet the single most important performance idea in this entire book, and it is not a compiler flag or a parallel library. It is a fact about memory so simple you could explain it to a child: Fortran stores a 2D array one column at a time. Get your loops to move with that grain and your code flies; move against it and the very same code, giving the very same answer, can run ten times slower. We introduce the fact here and spend all of Part VII cashing it in.

In this chapter, you will learn to:

• Declare arrays of any rank, control their bounds, and speak the vocabulary — rank, shape, extent — that the standard and every compiler error message use.
• Slice arrays with array sections, selecting a row, a column, a block, or every third element as a first-class value you can read, write, and compute with.
• Compute on whole arrays at once — c = a + b, y = sqrt(x) — and build arrays inline with constructors and implied-do loops.
• Reach for the right intrinsic array function — sum, maxval, matmul, dot_product, transpose, and their relatives — instead of writing the loop yourself.
• Size arrays at run time with allocatable arrays that clean up after themselves.
• Understand column-major order and write loop nests that respect it — the habit that separates fast Fortran from slow Fortran.

Learning Paths

How to read this chapter by track. - 🔬 Scientist — this is your chapter; read all of it. §5.3–5.4 (whole-array operations and intrinsics) will change how you write numerical code, and §5.6 (column-major) is why your code will be fast. - 📖 Standard — read straight through; arrays are load-bearing for every chapter that follows. - 🔧 Legacy — modern array syntax (§5.2–5.5) is exactly what FORTRAN 77 lacked; you will meet the old fixed-size, explicit-loop style in Part IV, and this chapter is the modern target you will be converting toward. - ⚡ HPC — §5.6 (column-major layout) is the foundation of everything in Part VII; do not skim it. The matmul note in §5.4 foreshadows why you will call LAPACK rather than roll your own.

5.1 Arrays Are First-Class Objects

Start with a contrast, because it explains everything that follows. In C, when you write double v[5];, the name v is essentially a pointer to the first of five doubles sitting in memory. The language does not remember that there are five of them; it does not remember their shape; v decays to a bare address the moment you pass it to a function, and the function must be told the length separately. An array in C is a convention layered on top of pointer arithmetic.

In Fortran, an array is a genuine object that the language knows things about. When you write

real(dp) :: v(5)

you have declared a real array of five elements, and Fortran remembers that. It knows v has one dimension. It knows that dimension runs from 1 to 5. It knows the array's total size. When you pass v to a procedure, all of that travels with it. The array is not a pointer wearing a costume; it is a first-class value, the same way an integer is. This is the first of the two design decisions from Chapter 1 — arrays are first-class objects — and now we make it concrete.

Three words describe every array, and you will see all three in compiler messages, so learn them now.

Definition (rank, shape, extent). The rank of an array is its number of dimensions: a vector has rank 1, a matrix rank 2, and Fortran allows up to rank 15. The extent along a dimension is the number of elements in it. The shape is the list of all the extents, one per dimension. So a real array declared a(3,4) has rank 2, extents 3 and 4, and shape [3, 4] — three rows and four columns, twelve elements in all.

Declaring and indexing

The simplest declaration gives an extent in parentheses, and the elements are numbered starting at one:

real(dp) :: v(5)          ! v(1), v(2), v(3), v(4), v(5)

That "starting at one" is not a detail; it is a decision with consequences, and it is the single most common thing that trips up programmers arriving from C, Python, or any other zero-based language.

⚠️ Common Pitfall — the first element is v(1), not v(0). Fortran arrays are one-based by default. The first element of v(5) is v(1); the last is v(5); there is no v(0), and touching it is an out-of-bounds error (compile with -fcheck=all and the program will catch it for you at run time). If you are translating a C loop for (i=0; i<n; i++), the Fortran equivalent is do i = 1, n — the bounds shift, and every index inside the loop shifts with them. This will bite you exactly once per language you come from; make the mistake now, on purpose, so you recognize it later.

One-based is the default, but Fortran lets you choose any bounds you like, which is a genuine convenience when your indices carry meaning:

real(dp) :: decades(1960:2020)    ! index BY the year: decades(1994) is legal
integer  :: stencil(-1:1)         ! index from -1 to 1, centered on zero

Here is a complete first program that declares arrays, fills them with a loop, and asks the language about their shape:

program array_basics
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: v(5)                 ! indices 1..5 (Fortran counts from 1)
  real(dp) :: decades(1960:1964)   ! indices you choose: 1960..1964
  integer  :: i

  do i = 1, 5
    v(i) = real(i, dp)**2          ! v = [1, 4, 9, 16, 25]
  end do

  do i = 1960, 1964
    decades(i) = real(i, dp)
  end do

  print '(a, f6.1)', 'first element v(1) = ', v(1)
  print '(a, f6.1)', 'last  element v(5) = ', v(5)
  print '(a, i0)',   'size(v)            = ', size(v)
  print '(a, f8.1)', 'decades(1962)      = ', decades(1962)
end program array_basics

$ gfortran -std=f2018 -Wall array_basics.f90 -o basics && ./basics
first element v(1) =    1.0
last  element v(5) =   25.0
size(v)            = 5
decades(1962)      =   1962.0

The size intrinsic reports the total number of elements — you never have to track it in a separate variable the way C forces you to, because the array carries its own size. That single fact eliminates an entire category of bug.

📜 From History. One-based indexing and column-major storage (§5.6) both date to FORTRAN I in 1957, and both were chosen to match how scientists already wrote mathematics on paper. A matrix entry is $a_{ij}$ with $i$ and $j$ starting at 1, and a matrix column is the natural unit of linear algebra. Other languages later chose zero-based indexing because it makes pointer arithmetic cleaner for systems programming — a perfectly good reason for their domain. Fortran optimized for the mathematician, not the pointer, and seventy years of numerical code has been written comfortably in that notation ever since.

5.2 Multidimensional Arrays and Array Sections

Scientific data is rarely a simple list. A temperature field on a plate is a grid — a 2D array. A stack of images, or a field evolving over time, is a 3D array. Fortran declares these by giving one extent per dimension:

real(dp) :: plate(100, 100)       ! rank 2: 100 x 100 grid, 10 000 points
real(dp) :: movie(100, 100, 500)  ! rank 3: the grid at 500 time steps

By near-universal convention we read a(i, j) as row i, column j, so the first subscript is the row and the second is the column. The array answers questions about itself through intrinsics you will use constantly:

integer :: a(3, 4)
! rank(a)   is 2        -- two dimensions
! size(a)   is 12       -- total number of elements
! size(a,1) is 3        -- extent of dimension 1 (rows)
! size(a,2) is 4        -- extent of dimension 2 (columns)
! shape(a)  is [3, 4]   -- the whole shape, as a rank-1 array

The superpower: array sections

Here is where Fortran pulls ahead of the pack. You do not have to work with a whole array or a single element; you can name a rectangular piece of it and treat that piece as an array in its own right.

Definition (array section). An array section is a subarray selected with a subscript triplet of the form lower:upper:stride in one or more dimensions. v(2:5) is the section of v from index 2 through 5; v(1:10:2) is every second element from 1 to 10; a lone colon : means the whole extent of that dimension. A section is a first-class array value — you can read it, assign to it, pass it to a procedure, and use it in any array expression, exactly as if it were a freshly declared array of that shape.

The cleanest way to see sections is to fill a small matrix so that every entry announces its own address — a(i,j) = 10*i + j, so the value 23 sits at row 2, column 3 — and then slice it every way:

program array_sections
  implicit none
  integer :: a(3, 4)            ! 3 rows, 4 columns
  integer :: row, col

  do row = 1, 3
    do col = 1, 4
      a(row, col) = 10*row + col     ! a(2,3) = 23, and so on
    end do
  end do

  print '(a, i0)',        'rank(a)   = ', rank(a)
  print '(a, i0)',        'size(a)   = ', size(a)
  print '(a, i0, a, i0)', 'size dims = ', size(a,1), ' x ', size(a,2)

  print '(a)',   'row 2      a(2, :)     :'
  print '(4i5)', a(2, :)                    ! a whole row
  print '(a)',   'col 3      a(:, 3)     :'
  print '(3i5)', a(:, 3)                    ! a whole column
  print '(a)',   'block      a(1:2, 2:3) :'
  print '(2i5)', a(1:2, 2:3)                ! a 2x2 sub-block
  print '(a)',   'stride     a(1, 1:4:2) :'
  print '(2i5)', a(1, 1:4:2)                ! row 1, every other column
end program array_sections

$ gfortran -std=f2018 -Wall array_sections.f90 -o sections && ./sections
rank(a)   = 2
size(a)   = 12
size dims = 3 x 4
row 2      a(2, :)     :
   21   22   23   24
col 3      a(:, 3)     :
   13   23   33
block      a(1:2, 2:3) :
   12   22
   13   23
stride     a(1, 1:4:2) :
   11   13

Read those outputs against the 10*i + j rule and every one confirms itself: row 2 is 21 22 23 24; column 3 is 13 23 33; the block a(1:2, 2:3) is the four corner-values 12, 13, 22, 23; the strided a(1, 1:4:2) picks columns 1 and 3 of row 1, giving 11 13.

There is one subtlety hiding in the block, and it is worth pausing on because it is a preview of the most important idea in the chapter. The 2×2 block printed as 12 22 then 13 23 — that is, down the first column, then down the second, not left-to-right across the rows. When Fortran streams a multidimensional array out (or through memory), it walks the first index fastest. Hold that thought; §5.6 is entirely about it.

🧩 Try It Yourself. Before you run it, predict the output of print '(2i5)', a(2:3, 3:4) for the same matrix. Write down the four values and their order. (Answer: the section is rows 2–3, columns 3–4, i.e. 23, 33, 24, 34 — streamed first-index-fastest as 23 33 then 24 34.) Compile it and check. Getting the order right, not just the set, is the skill.

🔄 Check Your Understanding. 1. What are the rank, shape, and size of an array declared real(dp) :: g(0:9, 5)? 2. Write the section that selects the last column of a rank-2 array m(n, n) without hard-coding n. 3. True or false: an array section can appear on the left of an assignment.

Answers

1. Rank 2; shape [10, 5] (the first dimension has extent 10, from index 0 to 9); size 50. 2. m(:, n) — or, size-agnostic, m(:, size(m,2)). 3. True. a(1, :) = 0.0_dp zeroes the first row; sections are assignable, first-class values.

5.3 Whole-Array Operations and Constructors

Now we cash in the fact that an array is a single object. If two arrays have the same shape, you can add them with a single +, and Fortran applies it element by element with no loop from you.

Definition (whole-array operation). A whole-array operation is an operation written on entire arrays rather than on individual elements. If a and b are conformable (same shape), then a + b, a * b, and a - b produce arrays of that shape, computed elementwise; a scalar combines with every element (2.0_dp * a); and an elemental intrinsic such as sqrt or sin, applied to an array, is applied to each element. Two arrays are conformable when they have the same shape, or when one is a scalar.

The distinction to keep straight from the first day: a * b on two arrays is the elementwise product, a(i)*b(i) for each i — it is not a dot product and not a matrix product. Those have their own intrinsics (§5.4). Fortran's * is elementwise, always.

Definition (array constructor). An array constructor builds a rank-1 array value inline from a list of elements between square brackets: [1.0_dp, 2.0_dp, 3.0_dp]. An implied-do inside a constructor generates elements with a loop-like expression: [(i*i, i = 1, 5)] builds [1, 4, 9, 16, 25]. Constructors make only rank-1 arrays; to build higher-rank arrays you pour a constructor into reshape (§5.4).

Putting whole-array operations and constructors together:

program whole_array_ops
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer  :: i
  real(dp) :: a(4), b(4), c(4)
  integer  :: squares(5)

  a = [2.0_dp, 4.0_dp, 6.0_dp, 8.0_dp]   ! array constructor
  b = [1.0_dp, 1.0_dp, 2.0_dp, 2.0_dp]

  c = a + b               ! whole-array add: c(i) = a(i) + b(i), all i, no loop
  print '(a, 4f7.2)', 'a + b   = ', c

  c = a * b               ! ELEMENTWISE product (not a dot or matrix product)
  print '(a, 4f7.2)', 'a * b   = ', c

  c = 0.5_dp * a          ! a scalar broadcasts across every element
  print '(a, 4f7.2)', '0.5*a   = ', c

  c = sqrt(a)             ! an elemental intrinsic, applied element by element
  print '(a, 4f7.2)', 'sqrt(a) = ', c

  squares = [ (i*i, i = 1, 5) ]           ! implied-do: [1, 4, 9, 16, 25]
  print '(a, 5i5)', 'squares = ', squares
end program whole_array_ops

$ gfortran -std=f2018 -Wall whole_array_ops.f90 -o waops && ./waops
a + b   =    3.00   5.00   8.00  10.00
a * b   =    2.00   4.00  12.00  16.00
0.5*a   =    1.00   2.00   3.00   4.00
sqrt(a) =    1.41   2.00   2.45   2.83
squares =     1    4    9   16   25

Every result is checkable in your head: a + b is [3, 5, 8, 10]; the elementwise a * b is [2, 4, 12, 16]; halving a gives [1, 2, 3, 4]; and sqrt([2,4,6,8]) is [1.41, 2.00, 2.45, 2.83] to two places. The squares line was generated by an implied-do with no separate loop at all.

Why does this matter beyond saving keystrokes? Because c = a + b tells the compiler the entire operation at once. It sees a single array add over a known shape, with no possibility (Fortran's no-aliasing rule from Chapter 1) that a, b, and c secretly overlap. That is precisely the information it needs to emit vectorized machine instructions — the ones that add eight numbers in a single step — and, later, to spread the work across cores. The readable form and the fast form are the same form. We prove this with measurements in Chapter 27.

🐍 Python Comparison. In NumPy you would write c = a + b too, and it would be fast — because NumPy's + calls down into a compiled C loop, the very kind of code Fortran emits natively. So for operations you can express as a few big array statements, NumPy and Fortran are in the same league. The cliff appears the moment your algorithm cannot be phrased that way — when each element depends on the last, as in the time-stepping heat solver you are building, and you are forced back into an explicit Python for loop over elements. There, pure Python is commonly 50–100× slower than Fortran, because every iteration pays the interpreter's overhead. The professional pattern, which we build in Chapter 15, is to keep the orchestration in Python and hand the hot elementwise loop to Fortran.

🚪 Threshold Concept — an array is a single object you compute with. The mental shift this chapter is really about is this: stop seeing an array as a bag of elements you visit one at a time with a loop, and start seeing it as one value — a vector, a matrix, a field — that you add, scale, slice, and transform whole. Once the shift happens, a great deal of numerical code collapses from a page of nested loops into three lines that read like the mathematics. And it is not only prettier: the whole-array form is exactly what lets the compiler go fast, because it hands over structure instead of hiding it inside a loop the compiler has to reverse-engineer. Readability and speed stop being a trade-off. That is the superpower.

5.4 The Intrinsic Array Functions

Fortran ships with a toolbox of intrinsic functions that act on whole arrays. Knowing them is the difference between writing a loop and writing a word. Here are the ones you will use daily.

Reductions and inquiries

A reduction collapses an array to a scalar (or to a lower-rank array). The workhorses:

Watch them work on the first eight digits of $\pi$:

program array_reductions
  implicit none
  integer :: v(8) = [3, 1, 4, 1, 5, 9, 2, 6]

  print '(a, i0)', 'size(v)          = ', size(v)
  print '(a, i0)', 'sum(v)           = ', sum(v)
  print '(a, i0)', 'product(v)       = ', product(v)
  print '(a, i0)', 'maxval(v)        = ', maxval(v)
  print '(a, i0)', 'minval(v)        = ', minval(v)
  print '(a, i0)', 'maxloc(v, dim=1) = ', maxloc(v, dim=1)
  print '(a, i0)', 'minloc(v, dim=1) = ', minloc(v, dim=1)
  print '(a, i0)', 'count(v > 3)     = ', count(v > 3)
end program array_reductions

$ gfortran -std=f2018 -Wall array_reductions.f90 -o reduce && ./reduce
size(v)          = 8
sum(v)           = 31
product(v)       = 6480
maxval(v)        = 9
minval(v)        = 1
maxloc(v, dim=1) = 6
minloc(v, dim=1) = 2
count(v > 3)     = 4

Each answer is worth verifying by hand once, so you trust the tool forever: the eight values sum to 31 and multiply to 6480; the maximum 9 sits in position 6; the minimum 1 first appears in position 2; and four of the eight values (4, 5, 9, 6) exceed 3. Notice maxloc and minloc return a position, not a value — with dim=1 on a rank-1 array they hand back a plain integer, which is what you usually want. And notice count(v > 3): the expression v > 3 is itself a whole-array operation, producing a logical array of the same shape, which count then tallies. We put that idea to work in §5.7.

Linear algebra: matmul, dot_product, transpose

Three intrinsics cover the linear-algebra operations that elementwise * deliberately does not:

• dot_product(x, y) — the scalar $\sum_i x_i y_i$ of two rank-1 arrays.
• matmul(A, B) — the true matrix product, $C_{ij} = \sum_k A_{ik} B_{kj}$.
• transpose(A) — the matrix with rows and columns swapped.

program array_intrinsics
  implicit none
  integer :: A(2,2), B(2,2), C(2,2), T(2,2)
  integer :: p(3) = [1, 2, 3]
  integer :: q(3) = [4, 5, 6]
  integer :: x(2) = [1, 1], Ax(2)
  integer :: row

  ! order=[2,1] lets us list entries row by row (see Section 5.6 for the "why").
  A = reshape([1, 2,  3, 4], [2,2], order=[2,1])   ! A = [[1, 2], [3, 4]]
  B = reshape([5, 6,  7, 8], [2,2], order=[2,1])   ! B = [[5, 6], [7, 8]]

  C  = matmul(A, B)        ! matrix * matrix
  T  = transpose(A)        ! rows <-> columns
  Ax = matmul(A, x)        ! matrix * vector

  print '(a)', 'C = matmul(A, B):'
  do row = 1, 2
    print '(2i5)', C(row, :)
  end do
  print '(a)', 'T = transpose(A):'
  do row = 1, 2
    print '(2i5)', T(row, :)
  end do
  print '(a, 2i4)', 'matmul(A, [1,1]) = ', Ax
  print '(a, i0)',  'dot_product(p,q) = ', dot_product(p, q)
end program array_intrinsics

$ gfortran -std=f2018 -Wall array_intrinsics.f90 -o intr && ./intr
C = matmul(A, B):
   19   22
   43   50
T = transpose(A):
    1    3
    2    4
matmul(A, [1,1]) =    3   7
dot_product(p,q) = 32

Trace the matrix product once, so you own it: $C_{11} = 1\cdot5 + 2\cdot7 = 19$, $C_{12} = 1\cdot6 + 2\cdot8 = 22$, $C_{21} = 3\cdot5 + 4\cdot7 = 43$, $C_{22} = 3\cdot6 + 4\cdot8 = 50$. The transpose flips [[1,2], [3,4]] into [[1,3],[2,4]]. The matrix-vector product A·[1,1] sums each row: [3, 7]. And dot_product([1,2,3],[4,5,6]) = 4 + 10 + 18 = 32. (We print each matrix a row at a time, with C(row, :), precisely to display it in the order you would write it — a habit worth keeping, for reasons §5.6 will make clear.)

⚠️ Common Pitfall — A * B is not matmul(A, B). For two rank-2 arrays, A * B is the elementwise product A(i,j)*B(i,j) and requires identical shapes; matmul(A, B) is the matrix product and requires the inner dimensions to agree (A is $m\times k$, B is $k\times n$). They give completely different answers and are among the most common array bugs. When you mean linear algebra, name it: matmul, dot_product, transpose.

🔗 Connection — matmul foreshadows LAPACK. matmul is perfect for small matrices and for prototyping, and the compiler implements it competently. But when your matrices are large, or you need to solve $A\mathbf{x} = \mathbf{b}$ rather than merely multiply, you do not write the algorithm yourself — you call LAPACK and BLAS, the Fortran libraries that sit underneath NumPy, MATLAB, and R. They are named in Chapter 16 and put to work in Chapter 21, where you call dgesv to solve a linear system. The reason a decades-tuned library beats your hand-rolled loop is a story we return to in Chapter 29; for now, know that matmul is the on-ramp and LAPACK is the highway.

5.5 Allocatable Arrays: Sizing at Run Time

Every array so far has had a size fixed when you wrote the program. But real programs learn their sizes at run time: the grid resolution comes from an input file, the number of particles from a command-line argument, the length of a dataset from however many rows the file happens to have. For that, Fortran gives you the allocatable array.

Definition (allocatable array). An allocatable array is declared with the allocatable attribute and a deferred shape — a colon for each dimension, a(:) or grid(:,:) — which reserves the rank but not the size. You give it a size at run time with allocate, release its memory with deallocate, and test its state with the allocated intrinsic. Crucially, a local allocatable array is automatically deallocated when the procedure it lives in returns — the language cleans up for you.

program allocatable_demo
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp), allocatable :: samples(:)     ! deferred shape: rank known, size not
  integer :: n, i

  n = 6                                    ! in real code, from a file or user input
  allocate(samples(n))                     ! size chosen now, at run time

  print '(a, l1)', 'allocated after allocate? ', allocated(samples)
  print '(a, i0)', 'size(samples)            = ', size(samples)

  do i = 1, n
    samples(i) = real(i, dp) / 2.0_dp      ! 0.5, 1.0, 1.5, 2.0, 2.5, 3.0
  end do

  print '(a, f6.2)', 'mean = sum/size          = ', &
        sum(samples) / real(size(samples), dp)

  deallocate(samples)                      ! explicit release (optional here)
  print '(a, l1)', 'allocated after dealloc?  ', allocated(samples)
end program allocatable_demo

$ gfortran -std=f2018 -Wall allocatable_demo.f90 -o alloc && ./alloc
allocated after allocate? T
size(samples)            = 6
mean = sum/size          =   1.75
allocated after dealloc?  F

The six samples are 0.5, 1.0, 1.5, 2.0, 2.5, 3.0; they sum to 10.5, and the mean is $10.5/6 = 1.75$. The allocated intrinsic reports T after allocate and F after deallocate, which is how defensive code checks an array's state before touching it.

That automatic deallocation is a bigger deal than it looks. In C you malloc and you must free, by hand, on every path out of the function — and forgetting, on even one error branch, is a memory leak. Fortran's local allocatables free themselves when they go out of scope, which quietly removes an entire class of bug that has haunted systems programming for fifty years. You can (and, for clarity in long-lived arrays, should) still deallocate explicitly; but you will not leak if you forget.

🐛 Find the Bug. This program will not compile. Why?

fortran real(dp) :: a(3), b(4) a = [1.0_dp, 2.0_dp, 3.0_dp] b = a ! <-- the compiler stops here

Diagnosis

a has shape [3] and b has shape [4]; they are not conformable, so the whole-array assignment b = a is illegal, and gfortran rejects it at compile time with a message about different shapes for array assignment. Whole-array operations demand matching shapes — that strictness is a feature, catching a mismatch the moment you write it rather than corrupting memory at run time as a C memcpy of the wrong length would. The fix is to make the shapes agree, or to assign a section: b(1:3) = a.

When an allocate can fail — a grid so large it exhausts memory — you will want to check for it rather than crash. Fortran's allocate accepts a stat= result and an errmsg= message for exactly that; we fold it into the broader error-handling discipline in Chapter 13. And when an array becomes a natural bundle with its dimensions and spacing — as the heat field will — you will wrap it in a derived type in Chapter 9. For now, one array, allocated when you know its size, is exactly the tool.

🔄 Check Your Understanding. 1. What does allocated(x) return for an allocatable x that has been declared but not yet allocated? 2. Why does a local allocatable array not leak memory even if you never call deallocate? 3. Declare a 2D allocatable real array grid, then allocate it as nx by ny.

Answers

1. .false. — declaration reserves the name and rank but allocates no storage; allocated is true only between allocate and deallocate. 2. Because a local allocatable is automatically deallocated when the procedure returns; the language releases it for you. (This is not true of pointer arrays — a distinction Chapter 11 makes.) 3. real(dp), allocatable :: grid(:,:) then allocate(grid(nx, ny)).

5.6 Column-Major Order: The Most Important Idea in This Book

We arrive at the fact that, more than any other single thing in this book, separates fast Fortran from slow Fortran. It is not subtle and it is not advanced. It is about how a two-dimensional array is laid out in a one-dimensional memory.

Computer memory is a single long line of addresses, 0, 1, 2, 3, …. A 2D array is conceptually a grid, but it must be flattened onto that line somehow, and there are two natural choices: store it row by row, or store it column by column. Fortran, since 1957, stores it column by column.

Definition (column-major order). In column-major order, the elements of a multidimensional array are laid out in memory so that the first index varies fastest. For a 2D array a, the storage order is a(1,1), a(2,1), a(3,1), … — all of column 1 first, contiguously, then all of column 2, and so on. C, C++, and NumPy (by default) use the opposite, row-major order, where the last index varies fastest. This single difference is the source of most Fortran-versus-C performance folklore and most Fortran↔Python array bugs.

Here is the layout, drawn out for a 3×3 array. The picture is worth memorizing:

A 3 x 3 array  a(i,j)             Linear memory (column-major):
                                  the FIRST index varies fastest.

        j=1    j=2    j=3         addr   holds
      +------+------+------+       0    a(1,1)  \
 i=1  | a11  | a12  | a13  |       1    a(2,1)   |  column 1  (contiguous)
      +------+------+------+       2    a(3,1)  /
 i=2  | a21  | a22  | a23  |       3    a(1,2)  \
      +------+------+------+       4    a(2,2)   |  column 2  (contiguous)
 i=3  | a31  | a32  | a33  |       5    a(3,2)  /
      +------+------+------+       6    a(1,3)  \
                                   7    a(2,3)   |  column 3  (contiguous)
   stepping DOWN a column          8    a(3,3)  /
   walks ADJACENT memory

Now the consequence, which is everything. Modern processors do not fetch one number from memory at a time; they fetch a whole cache line — typically 64 bytes, eight real(dp) values — in a single go, betting that if you wanted one you will soon want its neighbors. When your loop reads memory in the order it is stored, that bet pays off: each fetched cache line is used to the last byte before the next is loaded. When your loop reads across the grain of storage — jumping a whole column's length between consecutive accesses — most of every fetched line is wasted, and the processor stalls waiting on memory it will not use. Same array, same arithmetic, same answer; wildly different speed.

For Fortran, "with the grain" means: the innermost loop should run over the first index. Sweep down a column on the inside, move across columns on the outside:

program column_major
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: n = 3
  real(dp) :: a(n, n)
  integer  :: i, j

  ! FAST: inner loop over the first index i -> consecutive iterations touch
  ! consecutive memory (down a column). This is column-major loop order.
  do j = 1, n                 ! outer: over columns
    do i = 1, n               ! inner: over rows  <-- first index varies fastest
      a(i, j) = real(10*i + j, dp)
    end do
  end do

  print '(a)', 'a(i,j) = 10*i + j, printed row by row:'
  do i = 1, n
    print '(3f7.1)', a(i, :)
  end do
end program column_major

$ gfortran -std=f2018 -Wall column_major.f90 -o colmaj && ./colmaj
a(i,j) = 10*i + j, printed row by row:
   11.0   12.0   13.0
   21.0   22.0   23.0
   31.0   32.0   33.0

The loop nest above is correct. So is the version with the two do lines swapped — do i outside, do j inside — which produces the identical grid. But the swapped version walks memory in strides of n, touching one element of each cache line and discarding the other seven, and for a large n it can run several times to ten times slower for no other reason. The answer is the same; only the speed is destroyed. Getting this loop order right, every time, by reflex, is the most valuable habit this book will give you.

🚪 Threshold Concept — respect the memory order, and the speed follows. Most programmers think of a 2D array as a grid and never ask how it sits in memory — and in a slow language it does not much matter. In Fortran it matters enormously, because the language is fast enough that memory traffic, not arithmetic, is usually the bottleneck. Once you internalize first index fastest, inner loop over the first index, you will read a loop nest and instantly see whether it flies or crawls — and you will write the fast one without thinking. This one idea, applied everywhere, is a large part of what "knowing Fortran" actually means.

⚡ Performance Note. The 10× figure is an illustrative order of magnitude, not a promise — the true factor depends on the array size, the cache, and the compiler. What is not negotiable is the direction: first-index-inner is with the grain, last-index-inner is against it. We measure the real gap on real hardware in Chapter 27, and turn it into an optimization discipline (loop ordering, cache blocking, vectorization) in Chapter 29. Everything there rests on the fact you learned here.

🔄 Check Your Understanding. 1. In memory, which element immediately follows a(3,1) in a real(dp) :: a(3,3)? 2. You must sum a large a(n, n) with a nested loop. Which index belongs on the inner loop, and why? 3. NumPy defaults to row-major. When you pass a NumPy array to Fortran, what must you do about layout?

Answers

1. a(1,2) — after the last element of column 1 comes the first element of column 2 (first index fastest). 2. The first index, i — so consecutive inner iterations a(i,j) walk adjacent memory down a column. 3. Store or transpose it into column-major (order='F') so Fortran reads it correctly; this is exactly the NumPy↔Fortran issue handled in Chapter 15.

5.7 where and Masks on Real Data

Often you want to operate not on a whole array but on the elements that satisfy some condition — clamp the negatives to zero, cap the outliers, average only the valid readings. You met where in Chapter 4 as masked whole-array assignment; here we point it at realistic data and pair it with masked reductions.

A where construct applies an assignment only where a logical mask — itself a whole-array comparison — is true, with an optional elsewhere for the rest:

program where_masks
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: reading(6) = [ -1.0_dp, 3.5_dp, 20.0_dp, -0.2_dp, 9.0_dp, 14.5_dp ]
  real(dp) :: clean(6)

  clean = reading

  ! Replace invalid (negative) readings with 0 -- a masked assignment.
  where (clean < 0.0_dp)
    clean = 0.0_dp
  end where

  ! Cap anything above 10 at 10 (a saturation limit); the single-line form.
  where (clean > 10.0_dp) clean = 10.0_dp

  print '(a, 6f6.1)', 'raw    = ', reading
  print '(a, 6f6.1)', 'clean  = ', clean

  ! A MASKED reduction: work with only the originally-valid readings.
  print '(a, i0)',   'valid count   = ', count(reading >= 0.0_dp)
  print '(a, f6.2)', 'mean of valid = ', &
        sum(reading, mask = reading >= 0.0_dp) / real(count(reading >= 0.0_dp), dp)
end program where_masks

$ gfortran -std=f2018 -Wall where_masks.f90 -o wmask && ./wmask
raw    =   -1.0   3.5  20.0  -0.2   9.0  14.5
clean  =    0.0   3.5  10.0   0.0   9.0  10.0
valid count   = 4
mean of valid =  11.75

Follow the two masks in turn. The first, clean < 0.0_dp, is true for the two negative readings, so those become 0. The second, clean > 10.0_dp, is true for 20.0 and 14.5, so those saturate to 10. The result clean is [0.0, 3.5, 10.0, 0.0, 9.0, 10.0]. Then the masked reduction: count(reading >= 0.0_dp) finds the four non-negative readings, and sum(reading, mask = reading >= 0.0_dp) adds only those four — 3.5 + 20.0 + 9.0 + 14.5 = 47.0 — so their mean is $47.0/4 = 11.75$, computed ignoring the invalid entries entirely.

This pattern is everywhere in real scientific code. A climate array masks land cells against ocean cells. A sensor pipeline masks readings flagged as missing (often stored as a sentinel like -999). An image filter masks the pixels above a threshold. where and the masked sum/count/maxval let you express all of it without a single explicit loop or if — the condition rides along as an array, and the reduction respects it. (The elemental merge(t_source, f_source, mask), which picks elementwise between two arrays, is a handy cousin worth knowing; reach for it when you want an expression rather than a construct.)

Project Checkpoint

Time to make the heat solver real in the one way this chapter enables: the temperature field becomes a 2D allocatable array, and its Laplacian becomes a single whole-array statement built from array sections.

Recall the plan from Chapter 1: a square plate whose interior temperature evolves over time. The plate is a grid — exactly a rank-2 array u(nx, ny) — and the physics (the heat equation) needs, at every interior point, a quantity called the Laplacian: each point's four neighbors, added up, minus four times the point itself. Written point by point that is a double loop; written with array sections it is one line, and it is the whole interior at once:

lap(2:n-1, 2:n-1) =           u(1:n-2, 2:n-1) + u(3:n,   2:n-1)   &   ! up + down
                   +          u(2:n-1, 1:n-2) + u(2:n-1, 3:n  )   &   ! left + right
                   - 4.0_dp * u(2:n-1, 2:n-1)                         ! - 4 * centre

Each term is a section of the same 2×2 (for our tiny grid) interior, just shifted one step up, down, left, or right. The complete checkpoint program allocates the field, fills it with a stand-in pattern chosen so you can check the answer by hand, and computes the Laplacian in that single statement:

program heat_field
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter   :: n = 4          ! a tiny 4x4 grid, checkable by hand
  real(dp), allocatable :: u(:,:), lap(:,:)
  integer               :: i, j

  allocate(u(n, n), lap(n, n))
  lap = 0.0_dp                            ! interior is overwritten; edges stay 0

  do j = 1, n
    do i = 1, n                           ! inner loop over the first index (Section 5.6)
      u(i, j) = real(i, dp)**3            ! stand-in field u = i**3
    end do
  end do

  lap(2:n-1, 2:n-1) =           u(1:n-2, 2:n-1) + u(3:n,   2:n-1)   &
                     +          u(2:n-1, 1:n-2) + u(2:n-1, 3:n  )   &
                     - 4.0_dp * u(2:n-1, 2:n-1)

  print '(a)', 'interior Laplacian lap(2:3, 2:3), row by row:'
  do i = 2, n-1
    print '(2f8.1)', lap(i, 2:n-1)
  end do

  deallocate(u, lap)
end program heat_field

$ gfortran -std=f2018 -Wall project-checkpoint.f90 -o checkpoint && ./checkpoint
interior Laplacian lap(2:3, 2:3), row by row:
    12.0    12.0
    18.0    18.0

Here is the hand check, because you should never trust a stencil you have not verified once. With u(i,j) = i**3 the field is constant along each row: column values are 1, 8, 27, 64 for i = 1..4. At interior point (2,2) the Laplacian is u(1,2) + u(3,2) + u(2,1) + u(2,3) - 4*u(2,2) = 1 + 27 + 8 + 8 - 32 = 12; at (3,3) it is 8 + 64 + 27 + 27 - 108 = 18. Both interior rows come out constant because the field varies only with i, and the discrete second difference of i**3 is 6i — which is 12 at i=2 and 18 at i=3, exactly what printed.

Keep the numerics informal on purpose. This checkpoint is about the array machinery — a field that sizes itself, and a stencil expressed as sections — not about physics. The real finite-difference method: the $1/\Delta x^2$ scaling that turns this neighbor sum into an actual second derivative, the boundary conditions that hold the plate's edges fixed, the timestep and the stability limit that keep the simulation from blowing up — all of that is the heart of Chapter 24, and we defer it there deliberately. What you have built now is the data structure and the core spatial operation that Chapter 24 will wrap in time-stepping and the Chapter 38 capstone will make parallel. Save this as heat-solver/heat_solver.f90; it is the first real numerical kernel of your solver.

Summary

Arrays are why Fortran exists and why it is fast. This chapter turned "an array" from a list you loop over into a value you compute with.

The two things to memorize. First: c = a + b is not shorthand for a loop — it is more information than a loop, and that is why it is fast. Second, and above all: Fortran is column-major, so the inner loop runs over the first index. Everything in the performance half of this book is a consequence of that one sentence.

Spaced Review

Retrieval practice on the chapters that built up to this one. Answer before peeking.

1. (Ch. 3) What is the value of the integer expression 7 / 2 in Fortran, and how do you get 3.5 instead?

   Answer

   7 / 2 is 3 — integer division truncates toward zero. To get 3.5 you must make at least one operand real: 7.0_dp / 2.0_dp, or real(7, dp) / 2. This is the integer-division trap from Chapter 3; it appears inside array code too, e.g. an average written sum(v) / size(v) on an integer array truncates.

2. (Ch. 3) Why do we declare reals as real(dp) with dp = selected_real_kind(15, 307) rather than plain real?

   Answer

   Plain real is single precision (~7 digits) on most compilers; real(dp) requests at least 15 significant digits and an exponent range to $10^{307}$ — double precision — portably, without depending on a compiler default. Numerical code needs the extra precision to control round-off, so dp is our single precision kind everywhere.

3. (Ch. 4) What is the difference between exit and cycle inside a do loop?
   

   Answer

   `exit` leaves the loop entirely (jumps to the statement after `end do`); `cycle` skips the rest of the current iteration and goes on to the next one. With named loops you can target an outer loop: `exit outer` / `cycle outer`. Both are from [Chapter 4](../chapter-04-control-flow/index.md).

4. (Ch. 4) You met where in Chapter 4 as masked assignment. In §5.7 we also used count(reading >= 0.0_dp). What kind of value is the expression reading >= 0.0_dp on its own?
   

   Answer

   It is a whole-array operation producing a **logical array** of the same shape as `reading`, with `.true.` wherever the element is non-negative. That logical array is the mask that `where`, `count`, and a masked `sum` all consume — the same array-thinking from Chapter 4, now on real data.

5. (Ch. 3 + 4) Predict the output: integer :: k(4) = [1,2,3,4]; print *, sum(k) / 2.
   

   Answer

   `sum(k)` is `10` (an integer), and `10 / 2` is integer division giving `5`. Had `k` been `real(dp)`, `sum(k)/2.0_dp` would give `5.0`. The trap survives the jump to arrays — watch the type of a reduction before you divide.

What's Next

You can now build arrays, slice them, compute on them whole, size them at run time, and — most importantly — loop over them the fast way. But every example so far has lived inside one program. Real code is organized into reusable pieces, and those pieces pass arrays back and forth safely. Chapter 6 introduces procedures — subroutines and functions — and the feature that makes passing arrays around both safe and fast: intent, which lets you promise the compiler whether an argument is read, written, or both, and assumed-shape array arguments, which let a procedure accept an array of any size and ask it its own shape. That is where your solver's Laplacian becomes a reusable subroutine step — and where the arrays of this chapter start to travel.