Performance and Memory Layout

A performance change that was not measured is a guess. Three habits make the measurements trustworthy.

Build for the measurement. Correctness builds and timing builds are different builds. -fcheck=all (bounds, array temporaries, pointer and allocation checks) multiplies the run time of explicit loops, -O0 leaves every variable in memory, and -pg or sanitizers add their own overhead. The numbers below show the size of the effect on the same program (the loop-order program of the next section): the strict development flags are 3.6 times slower than -O2 on the fast loop, and for the explicit stencil loops of A heatflow stencil tuned step by step the strict build is about eleven times slower than -O3. Use -O2 or -O3 and no run-time checks, keep -g (it does not change the generated code) so profilers can name lines, and say which flags were used whenever a number is quoted. Compilers and the flags that matter are on Compilers and Compiler Options.

Time the right thing. Take the timing outside the code under test, repeat it, and keep the best value or the median, because the first repetition also pays for page faults and cold caches and any repetition can be disturbed by another process. Make the result of the computation observable (print a checksum), otherwise an optimizing compiler may remove the work. If the timed routine is pure and called with the same arguments in a loop, the compiler is allowed to call it once; vary the arguments between calls (the vectorization example below shifts its windows for exactly that reason).

Choose the clock. The standard has three timing intrinsics:

The kind of the integer arguments of system_clock matters. All the integer arguments of one reference must have the same kind, and the processor may pick the clock from the types and kinds of the arguments (16.9.202). The standard recommends integer arguments with a decimal exponent range of at least 18, that is 64-bit integers, so that the processor can choose its most accurate clock and the count wraps as rarely as possible. The first example shows both.

Output of the 64-bit and 32-bit queries, then a timed loop (gfortran 16.2, Flang 23 gave the same rates and the same result; the times are indicative):

With default integers both compilers offer a millisecond clock that wraps after 2147483647 ticks, about 24.9 days; with 64-bit integers they offer a nanosecond clock that does not wrap in practice. A 64-bit count therefore avoids both a coarse resolution and a negative difference after a wrap. The kinds come from iso_fortran_env (int64, real64), as in Intrinsic Types and Kinds.

Wall time and CPU time are equal here because the program is serial and uses one core without waiting. They differ when the program waits for I/O or other processes (wall greater than CPU) and when it runs threads (CPU greater than wall). The timing helper that every program on this page uses is the four-line function now(): it converts a 64-bit count and rate to seconds, so a measurement is t = now() before and now() - t after.

Source: J3/23-007r1, clauses 16.9.67 (CPU_TIME), 16.9.69 (DATE_AND_TIME) and 16.9.202 (SYSTEM_CLOCK); GNU Fortran manual: SYSTEM_CLOCK; CPU_TIME

When a program is slow and the reason is not known, measure with a profiler (gprof, perf) instead of timing by hand, and ask the compiler which loops it could not optimize (-fopt-info for gfortran, -Rpass-missed for Flang). Both are demonstrated on the heatflow kernels in Testing, Debugging and Profiling; this page covers what to do with the answer.

Fortran stores the elements of an array in array element order (9.5.3.3): the first subscript varies fastest. For a(4, 4) the elements a(1, 1), a(2, 1), a(3, 1), a(4, 1), a(1, 2) follow each other in memory (see Arrays Fundamentals for the full description and the effect on reshape and on I/O). C, C++, Python and most other languages do the opposite. A processor does not read single elements from memory but cache lines (64 bytes, eight double precision values, on most current processors) and a loop that uses all the values of a line before it is evicted is several times faster than one that uses one value per line.

The figure shows the two loop orders for the same array. The rule that follows is simple: the innermost loop runs over the first subscript, and in a nest over a(i, j, k) the order from outside to inside is k, j, i. The program below sums a 4096 by 4096 double precision array (128 MiB, far larger than any cache) in both orders. Both functions do the same additions; only the loop nest differs.

Run with -O2 (gfortran 16.2; Flang 23 -O2 gave 0.0092 s, 0.2028 s and a factor of 22.1):

The slow nest touches a new cache line on every addition (the stride is 4096 elements, 32 KiB) and, because the array is larger than the cache, reads every line eight times from main memory. It is 21 times slower. The same program under different builds:

Three things are visible. First, the better the compiler optimizes the good loop, the more the bad loop costs in comparison: -ffast-math lets gfortran vectorize the summation (a floating-point reduction cannot be vectorized otherwise, see Procedures, inlining and vectorization) and the good loop becomes twice as fast, while the bad loop cannot be helped by any vectorization because it is waiting for memory. Second, no optimizer rescued the wrong order: neither -O3 nor -O3 -ffast-math -floop-interchange (the GCC loop interchange pass, which is on by default at -O3) changed the time of the second function. The compiler would have to know the extents, the strides of an assumed-shape argument and that the reordering of a floating-point sum is acceptable; gfortran -fopt-info-loop-optimized lists, among its messages, versioned this loop for when certain strides are 1, which is the check that assumed-shape dummies need (see Contiguity and array temporaries). Do not count on loop interchange; write the nest in memory order. Third, an unoptimized build hides the problem (the strict build shows a factor of 3.6 only because everything is slow), so loop-order bugs are found by measuring optimized builds.

Whole-array statements and the array intrinsics (sum(a), a = b + c, maxval) walk the array in memory order for you. The orders to check are the explicit loops, the dim= argument of reductions (sum(a, dim=1) reduces along the fast dimension and is cheap, sum(a, dim=2) strides, but compilers handle it), implied-do loops in I/O and array constructors, and sections: a column a(:, j) is contiguous, a row a(i, :) is a stride-size(a, 1) access (the next section measures it). If an algorithm needs both a row-wise and a column-wise pass over the same data, the usual remedies are to transpose once with blocking so that tiles stay in cache, to store the matrix both ways, or to restructure the algorithm; the matrix product on this page uses the first (matmul, dot_product and BLAS). Passing array sections to other languages and to MPI has the same stride consequences (see MPI with mpi_f08 and Interoperability with C), and do concurrent and GPU Offload explains that the order of the indices in a do concurrent header does not by itself decide the nesting that the compiler generates.

Source: J3/23-007r1, clause 9.5.3.3 (array element order); Fortran Programming Language: Best Practices, Multidimensional arrays

An array is contiguous when its elements occupy consecutive storage in array element order (see Array Arguments: Assumed Shape, Size and Rank for the definitions: the contiguous attribute in 8.5.7, simply contiguous designators in 9.5.4 and is_contiguous, 16.9.115). Whole arrays and columns u(:, j) are contiguous; a row u(i, :), a strided section u(1:n:2) and the interior u(2:nx-1, 2:ny-1) of a grid with a halo are not. Whether it matters depends on the dummy argument:

The program scales every column and every row of a 2048 by 2048 array through the three kinds of dummy arguments:

With -fcheck=all (included in the strict flags) gfortran warns at run time, once per source location, when it creates an array temporary for an actual argument, which is how copy-in/copy-out is found:

and -Warray-temporaries makes the compiler report them at compile time (command gfortran -std=f2023 -O2 -Warray-temporaries -c 04_contiguity.f90):

The first warning is the copy-in, the second the copy-out of line 60; the third is the contiguous dummy at line 72. The call with the assumed-shape dummy on line 66 gets none, and neither do the column calls. The timings (best of five, -O2):

Run-to-run variation of this size (a factor of two on the fastest, memory-resident cases) is normal; the pattern is the same in all of the columns:

How to avoid temporaries in practice: declare dummy arrays assumed-shape (x(:), u(:, :)) and not explicit-shape unless the callee really needs sequence association; store the data so that the pieces that are passed to procedures are columns (the last subscript selects the piece); keep halo cells in separate arrays or accept the stride; add contiguous where the callee is a hot kernel and callers pass whole arrays; and branch with is_contiguous when both cases occur. gfortran also offers -frepack-arrays, which copies a non-contiguous assumed-shape argument into a contiguous block in the procedure prologue (it helps only if the procedure accesses the data many times) and -fstack-arrays, which places temporaries on the stack. Flang 23 has -frepack-arrays, -fstack-repack-arrays and -fstack-arrays too; its help text describes them as options that pack non-contiguous assumed-shape dummies and allocate array temporaries on the stack. Array-valued expressions can create temporaries as well (a function result that is used as an operand, the arguments of cshift); -Warray-temporaries reports those that gfortran creates, as shown for cshift in Array Intrinsics and Reductions.

Pointers add a variant: an array pointer may point at a non-contiguous section, and a contiguous pointer, the contiguous attribute on a dummy pointer or is_contiguous tells the compiler or the program what holds (see Pointers and Targets).

Source: J3/23-007r1, clauses 8.5.7 (CONTIGUOUS attribute), 9.5.4 (simply contiguous array designators) and 16.9.115 (IS_CONTIGUOUS); GNU Fortran manual: Options to request or suppress errors and warnings and Options for code generation conventions; LLVM Flang driver

A derived type puts the fields of one object next to each other: an array of cell_t with four double precision components is an array of structures (AoS), in which one cell occupies 32 bytes and the temperature of the next cell is 32 bytes away. The alternative is a structure of arrays (SoA): one derived type that holds one array per field. A loop that reads and writes two fields of every cell moves 32 bytes per cell in the AoS layout (the other two fields come along in the same cache lines) and 16 in the SoA layout, and the SoA loop is contiguous, so it vectorizes without gather instructions. A loop that uses all fields is about equally good in both, and an algorithm that jumps to a single cell and needs all its fields prefers AoS. The comparison for ten million cells:

The remaining routines (all_fields_aos and all_fields_soa) use the same pattern with all four fields; the complete program times the four routines and checks that the temperatures agree. Output (gfortran 16.2 -O3):

The SoA version is 2.1 to 2.7 times faster for the two-field update; when all four fields are used the gap shrinks to 1.2 to 1.6 times and it does not reverse for this sequential loop. (Without optimization, in the strict build, AoS was as fast or faster, 0.032 s against 0.034 s for two fields and 0.033 s against 0.053 s for four: the layout advantage is something the optimizer realizes, not the language.) Two guidelines follow. Prefer separate arrays, or one array whose last subscript selects the field (u(i, j, field)), for bulk numerical data, and keep the fastest-varying subscript the one the loops run over; an array u(field, i, j) is AoS. Use AoS or a type with allocatable components when the objects are accessed one at a time, are of different lengths, or travel as units (message passing, I/O). The heatflow grid, grid_t with a single allocatable field u(:, :), is an SoA in this sense: when more fields are added (conductivity, source), separate real(dp), allocatable components are a better default than an array of cells. Arrays of derived types with allocatable or pointer components can be slower again, because each element then owns a separately allocated block (pointer chasing); parameterized derived types with len parameters and polymorphic array elements can behave the same way, see Derived Types, Parameterized Derived Types and Object Oriented Programming.

Source: J3/23-007r1, clause 9.5.3.3 (array element order) and clause 7.5 (derived types)

allocate and deallocate are calls to the memory allocator. For a small block with a free block ready, the measurement below shows about ten nanoseconds per allocate and deallocate pair; the cost is higher when the allocator has to obtain fresh pages from the operating system, which these measurements did not show. The rule of thumb is to allocate before the loop that runs many times and to reuse the work arrays inside it. Four places hide allocations in Fortran, and the standard’s semantics for each (10.2.1.3, 9.7) lets a compiler remove only some of them:

The first program calls a three-point smoothing routine with a work array that is allocated on each call, an automatic array, a work array passed in by the caller, and an array-valued function result, for vectors of 8, 256 and 250000 elements (each size does the same amount of arithmetic):

Output with gfortran 16.2 -O2 (the time of each variant for about the same amount of arithmetic; indicative):

The honest summary of this measurement is that the allocation is not the dominant cost of a routine that does a pass over its data: the variants are within a factor of two of each other at every size, and -fstack-arrays changed nothing beyond the run-to-run variation (gfortran -O2 -fstack-arrays: 0.151, 0.119, 0.128 and 0.096 s for 8 elements, the same pattern as without it). The function result was the fastest on gfortran because the other variants copy the input into the work array first, an extra pass over the data that the function formulation does not need. On Flang 23 -O2 the same program gave 0.092, 0.076, 0.031 and 0.058 s for 256 elements (0.090, 0.073, 0.030 and 0.057 s with -fstack-arrays): there the caller-provided work array was three times faster than allocating on every call. The cost of an allocation becomes visible when the work per allocation is tiny. The second program creates a three-element work vector in the innermost loop, and assigns a four-element allocatable array 20 million times with and without the section form:

The results, in seconds for 20 million repetitions (gfortran -O2 and -O2 -fno-realloc-lhs gave the same figures; -O2, indicative):

gfortran removed the allocate/deallocate pair of the small, constant-size array at -O2 (the loop is as fast as the one with the fixed array); Flang 23 did not, and is 15 times slower than its own fixed-size loop. Without optimization (the strict -fcheck=all build) the same loop costs 0.188 s against 0.065 s for the fixed array in gfortran. Whether a compiler can prove that an allocation is unnecessary is not something the language promises, so for small, fixed-size work vectors use a fixed-size array (a named constant size, real(dp) :: w(3)) and for variable-size work arrays allocate once outside the loop. The realloc-on-assignment check was not measurable here: v = a*v + 1 and v(:) = a*v + 1 took the same time on both compilers, and -Wrealloc-lhs printed only the assignment of the constructor in sum_allocating (Code for reallocating the allocatable array at (1) will be added). Keep v(:) = ... anyway in code that must also be fast on compilers that do not remove the check.

The stack has a size limit (commonly 8 MiB for the main thread on Linux, set with ulimit -s; OpenMP threads have their own stack size, set with OMP_STACKSIZE). Put the big arrays in allocatable variables, which live on the heap, and use -fstack-arrays only when the arrays are known to be small enough. -fmax-stack-var-size=n sets the size above which gfortran moves a local array of known size from the stack to static memory. Memory that is allocated is not necessarily memory that is available: on Linux a large allocate may succeed and the process may be killed when the pages are touched, as noted on Allocatable Arrays and Dynamic Memory.

Source: J3/23-007r1, clauses 9.7 (dynamic association) and 10.2.1.3 (interpretation of intrinsic assignments); GNU Fortran manual: Options for code generation conventions; Fortran Programming Language: Best Practices, Allocatable arrays

Two names are aliased when they refer to the same storage. The standard gives the compiler the guarantee that makes Fortran loops easy to optimize: while an entity is associated with a dummy argument, the procedure shall not modify it through any other name, and if it is modified through the dummy argument it shall be referenced only through that dummy (15.5.2.14, with exceptions for the pointer and target attributes). Two ordinary array dummies therefore never overlap, a loop over them can be reordered and vectorized without a run-time check, and restrict is not needed. The guarantee is a promise by the programmer, which the compiler does not check. This program breaks it, on purpose: the dummies x and y of add_to are v(1:11) and v(2:12), which overlap.

The program is not conforming and its output depends on how the loop was compiled. The same source with different compilers and optimization levels (gfortran 16.2 -std=f2023, Flang 23):

The serial, sequential meaning (a running sum, 1 2 3 ... 12) is what -O0 and Flang produce; gfortran -O2 and -O3 vectorize the loop under the no-overlap assumption and produce different, equally wrong, answers. -fcheck=all did not detect the aliasing and printed the sequential result, as the strict run shows. Nothing in the language or the tools detects the violation reliably, so the rule is the discipline: never pass overlapping pieces of the same array to two dummy arguments of a procedure that modifies one of them; copy one of them, or give the dummies the target or pointer attribute so that the compiler assumes aliasing is possible. Remember that the guarantee is also what makes intent(in) arguments useful to the optimizer, and that the target and pointer attributes, and volatile and asynchronous, take it away for the objects that have them. The pointer rules, including why a target dummy and an actual argument that is not contiguous interact, are on Pointers and Targets.

Source: J3/23-007r1, clause 15.5.2.14 (restrictions on entities associated with dummy arguments)

matmul (16.9.134) and dot_product (16.9.73) are intrinsic, but how they are implemented is the compiler’s business. gfortran inlines small products (-finline-matmul-limit=n, default 30, active when front-end optimization is on, which every -O option except -O0 and -Og enables, per the manual) and otherwise calls a library routine in libgfortran; with -fexternal-blas it can call the gemm of an external BLAS for matrices larger than -fblas-matmul-limit=n. Flang lowers matmul to its own runtime. The examples below multiply two 1000 by 1000 matrices with the intrinsic and with three hand-written loop nests: the dot-product order (the innermost loop over k, which walks a along a row), the column order (j, k, i, the innermost loop over the first subscript) and the column order with blocking over k and i so that a block of a stays in cache.

The blocked version (product_blocked, block size 64) is in the complete program. Results for n = 1000, the 2 n^3 floating-point operations (the three loop versions agree with matmul to 3.4e-13 in gfortran, to 0 in Flang; indicative):

The numbers say several things. The loop order alone is a factor of 4.6 at -O3 (4.9 against 22.8 GFlop/s), and the compiler’s vectorization turns the correct order into a further factor of 3 from -O2 to -O3; the blocking is worth only 9 per cent at this size, because the matrices are small enough that the working set of the column order is partly cached, and its value grows with n. The intrinsic is 1.7 to 1.9 times faster than the best simple loops in gfortran at -O3 (and 6 times at -O2), and twice as fast again when it is routed to OpenBLAS with -fexternal-blas (with -fblas-matmul-limit=1 it reached 102 GFlop/s). In Flang 23 matmul was about four times slower than gfortran’s on this machine and only on par with a good hand-written loop. These figures are for one machine and one size; on a different processor the ratios differ, as the comparison with reference BLAS, OpenBLAS and dgemm on BLAS, LAPACK and Numerical Libraries shows. The practical advice: use matmul for small and medium products where the code’s clarity matters, call a tuned BLAS (dgemm) or LAPACK for large dense linear algebra, and write loops only when the operation is not one that a library provides, in memory order. The matmul of a transposed argument and the other temporaries it can create are discussed on Array Intrinsics and Reductions.

Source: J3/23-007r1, clauses 16.9.73 (DOT_PRODUCT) and 16.9.134 (MATMUL); GNU Fortran manual: MATMUL; Options for code generation conventions; BLAS (netlib)

A stencil sweep does about four floating-point operations per updated point and moves 16 bytes through main memory (one read and one write of a double precision value, when the neighbours come from the cache), so it is limited by memory bandwidth long before it is limited by arithmetic. More threads help only until the memory system is saturated, and cpu_time adds the time of all threads. The program below runs the same sweep serially and with an OpenMP parallel do over the columns (the static schedule gives each thread a block of whole columns, so each thread runs the stride-1 inner loop), and initializes the arrays with the same schedule so that each thread’s pages are placed in memory close to that thread (first touch, relevant on machines with several memory nodes).

The main program times 40 sweeps over two 4096 by 4096 arrays (256 MiB in total) with system_clock and cpu_time. Built with -O3 -fopenmp and run with OMP_NUM_THREADS set (indicative):

The speed-up stops at about two on this machine however many threads there are (it has ten cores, and a sweep is memory-bound), cpu_time grows with the number of threads (it is the sum over the threads in both compilers), and the serial sweep was already limited by memory, so tuning the arithmetic of a memory-bound loop gains nothing: reduce the memory traffic (fewer arrays, fewer sweeps, blocking in time, a smaller type) before adding threads. With one thread the parallel version is slightly slower than the serial one because of the OpenMP runtime. The strict -fcheck=all build ran the serial sweep in 3.078 s and the parallel one in 1.397 s with two threads. do concurrent, coarrays and MPI have the same memory-system limits; their usage is on OpenMP and OpenACC, do concurrent and GPU Offload, Coarrays, Images and Synchronization and MPI with mpi_f08. Parallel execution of small problems costs more than it gains because of thread start-up and synchronization; measure with the problem size that matters. Avoid threads writing to neighbouring elements of the same cache line in the same loop (false sharing): when each thread owns whole columns, as here, it does not happen.

Source: OpenMP specifications; GNU Fortran manual: CPU_TIME

Formatted I/O converts every value to text and back; unformatted stream I/O moves the bytes. For a checkpoint of a large grid the difference is two orders of magnitude, and the stream form is exact by construction (the formatted form with es24.16, 17 significant digits, agreed to 1e-15 here, at a much higher cost). The program writes and reads two million doubles through a scratch file with each method:

(gfortran 16.2 -O2; Flang 23 -O2 gave 0.268 s / 0.492 s and 0.001 s / 0.003 s. The unformatted figures are largely the page cache: the 16 MB were never written to a disk.) Use access='stream' with form='unformatted' for data that only Fortran programs read (the heatflow write_stream and read_stream do), write whole arrays in one statement (not one element at a time), and use a format with enough digits (es24.16 for double precision) when text is required. Details of units, files, formats and non-advancing I/O are on Input/Output: Units and Files, Format Specifications and Edit Descriptors and Derived-Type and Asynchronous I/O.

Source: J3/23-007r1, clause 12 (input/output statements)

The Jacobi kernel of heatflow (jacobi_sweep in heatflow_kernels) is the example for the rest of the page. In the project the solver calls it, copies the new interior back into the grid, and applies the boundary condition on every step. The program below builds that loop in five steps, always updating a fixed-boundary grid with the same arithmetic, and measures each step on two grid sizes with the same total number of cell updates (about 10^8 per measurement): 256 by 256 (two arrays of 512 KiB) and 2048 by 2048 (two arrays of 32 MiB). The steps:

The kernels that carry the steps:

The drivers run_wrong_order, run_baseline and run_lazy_residual, and the timing program, are in the same file. Output with gfortran 16.2 -O3 (rates in millions of cell updates per second; max diff vs 1 compares the final grids with step 1, and a value of zero means the results are identical to the last bit):

The same steps with the other builds (Mcell/s, grid 2048 by 2048; the 256 by 256 values are in parentheses):

What each step teaches:

Applying the steps to the project means: keep the stencil loops in memory order (step 1, already done), avoid the copy-back by alternating two arrays (step 2), compute the residual every k-th step (step 3), build with -O3 and no run-time checks, and parallelize the outer loop with OpenMP or do concurrent once the serial kernel is bandwidth-bound. The result of each step is bit-identical to the previous one here, which is a property to preserve when tuning: tuning that changes the numbers (reassociation with -ffast-math, a different order of the additions) needs a tolerance and a test, as described in Testing, Debugging and Profiling.

The options are described in the manuals linked below; the build systems that carry them (fpm profiles, CMake build types) are on Fortran Package Manager (fpm) and Building with Make, CMake and Meson. A plain fpm build passes no warning flags and uses its default (debug) profile; use fpm build --profile release for timing.

The two Fortran books consulted for this section predate the standard of 2023 and most of today’s compiler technology. For performance and memory layout the differences are:

Not checked here, because no compiler for them is installed: Intel ifx, NAG, HPE Cray and NVIDIA HPC compilers. See History, Standards and Compilers for the per-compiler support table.

See also Best Practices and Style and Common Gotchas and Pitfalls.