Fortran Array Slicing and Vectorization: Writing Numerical Code the Compiler Can Optimize

The loop tax in numerical code

If you write Fortran for engineering simulations, you have probably typed a nested DO loop to add two matrices more times than you want to admit. The pattern is mechanical: declare indices, loop over columns, loop over rows, assign the sum. It works, but it spends your attention on index arithmetic instead of the mathematics. Modern Fortran offers a shorter path—whole-array operations and slicing—that can also give the compiler a stronger signal to generate vectorized machine code. The catch is that not every slice is friendly to the CPU cache.

Array syntax as a vectorization hint

In Fortran, an expression like C = A + B is not just a loop in disguise. It is a whole-array operation: the compiler knows every element of C depends only on the corresponding elements of A and B. That independence is exactly what SIMD—Single Instruction, Multiple Data—hardware needs. Instead of loading one float, adding it, and storing it, a vectorized loop loads a register packed with 4, 8, or 16 values and performs the addition in a single instruction.

Explicit loops can also vectorize, but the compiler must prove that iterations do not interfere. Array syntax makes that proof trivial for element-wise expressions. The DO CONCURRENT construct, introduced in Fortran 2008, serves a similar purpose for more complex loop bodies by telling the compiler that iterations are independent and safe to reorder or parallelize.

Worked example: matrix addition and slicing

Consider two n × n matrices. The explicit version requires careful ordering because Fortran stores arrays in column-major order: the first index varies fastest in memory. Looping over the first index in the inner loop keeps memory access contiguous.

program matadd
  implicit none
  integer, parameter :: n = 1024
  real :: a(n,n), b(n,n), c(n,n)
  integer :: i, j

  ! Initialization omitted
  ! a = ... ; b = ...

  ! Explicit nested loop
  do j = 1, n
    do i = 1, n
      c(i,j) = a(i,j) + b(i,j)
    end do
  end do

  ! Array syntax: identical result, one statement
  c = a + b
end program

The array statement c = a + b replaces six lines of loop scaffolding. You can also operate on subarrays without copying them. The following slice adds only the top half of each column:

  c(1:n/2, :) = a(1:n/2, :) + b(1:n/2, :)

Because the sliced dimension is the first index, the elements remain contiguous in memory. The compiler can still emit packed SIMD instructions for each column segment.

Strided slices are different. This statement copies every other column:

  c(:, 1:n:2) = a(:, 1:n:2)

Each column is still a contiguous block, but the blocks are separated by a stride of n elements. The CPU must jump across memory, which can reduce cache and prefetch efficiency. A slice with a stride on the first index, such as a(1:n:2, :), is worse: individual elements are non-contiguous, and the compiler may fall back to scalar loads or slower gather instructions.

Verifying vectorization and avoiding stride penalties

Syntax alone does not guarantee speed. You should verify that the compiler vectorized the hot loop. With GNU Fortran, compile with optimization and ask for a vectorization report. Run this in your project directory where the source file lives:

gfortran -O3 -fopt-info-vec matadd.f90 -o matadd

Expected check: the compiler prints lines such as optimized: loop vectorized using 32 byte vectors. If you see not vectorized: data ref not sufficiently aligned or similar, the slice may be strided or the arrays may need alignment attributes.

To inspect the machine code directly, generate assembly and search for packed instructions:

gfortran -O3 -S matadd.f90 -o matadd.s
grep -E 'addps|addpd|vaddps|vmovups' matadd.s

Expected check: the presence of packed SIMD arithmetic (for example, vaddps on x86-64 with AVX). If only scalar addss instructions appear, vectorization did not happen.

During development, enable bounds checking to catch slice mismatches:

gfortran -O0 -fcheck-all matadd.f90 -o matadd_debug
./matadd_debug

Risk: -fcheck-all adds significant runtime overhead and is usually disabled in production builds. A slice that passes in debug mode can still segfault in an optimized production binary if the shapes are wrong and checking is off.

Array syntax has two well-known limitations. First, complex expressions can create temporary arrays. A statement like C = A + B + D may allocate an intermediate array for A + B, then add D into a second temporary before assigning to C. For very large datasets, that extra memory traffic can dominate runtime. Breaking the expression into two statements—C = A + B followed by C = C + D—often avoids the hidden allocation.

Second, slicing does not eliminate the physics of memory. A non-contiguous slice forces the CPU to fetch cache lines that are only partially used. If your hot loop operates on a(1:n:2, :), you are reading every other element from each column; the memory subsystem still loads entire cache lines, wasting bandwidth. When performance matters, keep the strided dimension as the outer index and the contiguous dimension as the inner one.

Slice expression	Memory pattern	SIMD outlook
a(:, j)	Contiguous column	Excellent
a(i, :)	Strided row	Poor without transpose
a(1:n/2, :)	Contiguous segment	Excellent
a(1:n:2, :)	Every other element	Weak; scalar or gather
a(:, 1:n:2)	Contiguous blocks, large stride	Moderate; vectorizable per block

Practical recommendations

• Use whole-array syntax for element-wise operations; it is shorter and gives the compiler an explicit independence guarantee.
• Prefer slices along the first array index in Fortran to preserve contiguity.
• Compile development builds with -O3 -fopt-info-vec and confirm vectorization at the assembly level for critical kernels.
• Use -fcheck-all in debug builds, but never assume production binaries will catch out-of-bounds slices.
• Profile before and after refactoring loops into array expressions, especially when the working set exceeds last-level cache.

Array slicing in Fortran is a tool for expressing intent, not a magic performance switch. When you respect column-major layout and verify the generated code, it removes boilerplate and lets the hardware do what it does best. When you ignore stride and temporaries, the same syntax can quietly make your kernel memory-bound.