Why Fortran's Array Slices Still Beat Hand-Rolled Loops for Numerical Work

The problem: loops obscure what the hardware can do

Most engineers reach for nested loops when they need to update a subregion of a 2D field—say, applying a boundary condition to the interior of a fluid grid. The code works, but the compiler sees only scalar operations. It cannot automatically vectorize across the inner dimension, and the memory access pattern often fights the cache. Fortran's array syntax and slice notation were designed to make the intent visible: operate on this whole rectangle at once. The payoff is not just shorter source; it is the difference between a loop the optimizer ignores and a kernel it can turn into aligned SIMD stores.

Slice notation replaces three lines with one

Consider a 1000×1000 real array u where you want to add a scaled neighbor difference to every interior point. In explicit loops:

do j = 2, n-1
   do i = 2, m-1
      u(i,j) = u(i,j) + 0.25 * (u(i+1,j) + u(i-1,j) &
                                + u(i,j+1) + u(i,j-1) - 4.0*u(i,j))
   end do
end do

With array slices the same update becomes:

u(2:n-1, 2:m-1) = u(2:n-1, 2:m-1) + 0.25 * &
   (u(3:n,   2:m-1) + u(1:n-2, 2:m-1) + &
    u(2:n-1, 3:m)   + u(2:n-1, 1:m-2) - 4.0*u(2:n-1, 2:m-1))

The slice u(2:n-1, 2:m-1) denotes a contiguous subarray in column-major order—the first index varies fastest. The compiler emits a single vectorized loop nest with no temporaries, provided the arrays are not aliased. Adding intent(inout) :: u and intent(in) :: ... on dummy arguments in a subroutine makes the no-alias guarantee explicit, letting the optimizer keep values in registers across the whole expression.

WHERE construct handles conditional updates without branches

Suppose the same grid needs a positivity fix: any negative value becomes zero. A scalar loop needs an if inside the inner iteration, which defeats predication on many architectures. Fortran's where construct expresses the same logic as a masked assignment:

where (u < 0.0)
   u = 0.0
elsewhere
   u = u  ! no-op, shown for clarity
end where

The compiler lowers this to a blend instruction (e.g., vblendvpd on AVX2) with zero branch mispredictions. For more complex masks you can nest where or use merge in elemental expressions.

Worked example: verifying the speedup on a small stencil

Compile both versions with gfortran -O3 -march=native -fopt-info-vec (GCC 12+). The vectorization report will show loop vectorized for the slice version and not vectorized: complicated access pattern for the hand-written nest. Run each kernel 100 times on a 2000×2000 double-precision array and measure wall time with cpu_time. Typical results on a Skylake-X core:

• Explicit loops: ~1.8 s
• Array slices: ~0.6 s

The 3× difference comes from 256-bit FMA units processing four doubles per cycle versus scalar code limited by load/store throughput. The slice version also avoids the temporary array that a naïve u = u + ... expression would create; the compiler fuses the arithmetic into a single pass.

Limitation: column-major layout is non-negotiable

Fortran stores arrays column-major. If you write u(i,j) with i as the row index (common in C-interop code), the inner loop strides by the leading dimension, causing a cache miss every iteration. The fix is to declare the array as real, dimension(m,n) :: u and always iterate do j=1,n; do i=1,m so the contiguous dimension is innermost. When calling C libraries that expect row-major, use iso_c_binding and transpose explicitly with transpose(u) or pack into a contiguous buffer—do not pretend the layout matches.

Checklist before you trust the optimizer

1. Add implicit none at the top of every program unit—otherwise i-n default to integer and silent bugs appear.
2. Compile with -fopt-info-vec-missed (GCC) or -qopt-report=5 (Intel) to confirm vectorization.
3. Run a correctness test: compare the slice result against a reference loop on a 4×4 matrix with known values.
4. Profile with perf stat -e cycles,instructions,cache-misses to verify the instruction mix matches expectations.

Closing thought

Array slices and where are not syntactic sugar—they are the contract that lets the compiler generate the code you would write by hand in assembly, but without the maintenance burden. Start by rewriting one hot stencil kernel in slice form, verify the vectorization report, and measure. The performance gain is real, and the source becomes easier to audit for correctness.