Fortran DO CONCURRENT: The Smallest Suitable Design for Data-Parallel Loops

The Problem: Expressing Parallelism Without the Machinery

You have an array computation—element-wise arithmetic, a stencil update, a reduction—that runs too slowly in a plain DO loop. You could add OpenMP directives, but that introduces a dependency, a runtime, and a mental model of threads. You could rewrite for coarrays, but that changes the execution model entirely. What you need is a way to tell the compiler: these iterations don't depend on each other; you figure out the rest. That is exactly what DO CONCURRENT does.

Requirements: When DO CONCURRENT Is the Right Tool

The construct fits when:

• The loop body performs independent work per iteration (no cross-iteration reads/writes to the same memory).
• The iteration space is a simple integer range or a FORALL-style triplet.
• You want the compiler to choose vectorization, threading, or both—without you specifying which.
• You accept that floating-point results may differ bitwise from sequential execution due to reordering.

If you need ordered reductions, shared mutable state, or deterministic bitwise reproducibility, DO CONCURRENT is not the smallest suitable design—it's the wrong one.

The Construct: Syntax and Locality

Standardized in Fortran 2008 and refined in 2018, the basic form is:

do concurrent (i = 1:n)
    a(i) = b(i) * c(i) + d(i)
end do

Fortran 2018 added locality specifiers that make the data boundary explicit:

do concurrent (i = 1:n) local(x) shared(y) default(none)
    x = b(i) * 2.0
    a(i) = x + y
end do

• LOCAL(var) — each iteration gets its own private copy.
• SHARED(var) — one instance shared across iterations (read-only unless you synchronize externally).
• DEFAULT(NONE) — forces declaration of every variable used in the loop; the compiler errors on undeclared ones. Adopt this as a cheap audit mechanism.

Without locality specifiers, the compiler infers locality, but inference can be wrong when pointers, module variables, or assumed-shape arrays hide aliasing.

Trust Boundary: The Programmer Guarantees Independence

DO CONCURRENT is an assertion, not a request. The compiler is permitted to trust it. Violating the independence guarantee—through hidden aliasing via pointers, TARGET attributes, or module-scope variables—produces silent data races and nondeterministic wrong results. No compile error, no runtime warning.

Example of a violation that compiles cleanly but breaks under optimization:

module m
    real, target :: shared_val
end module m

subroutine broken(n, a, b)
    use m
    integer, intent(in) :: n
    real, intent(out) :: a(n)
    real, intent(in) :: b(n)
    real, pointer :: p(:)
    p => a  ! aliasing hidden from the loop
    do concurrent (i = 1:n)
        shared_val = b(i) * 2.0  ! cross-iteration write
        a(i) = shared_val + p(i) ! read through alias
    end do
end subroutine

Under -O2 or -O3, iterations may be reordered or vectorized, causing shared_val to be overwritten before other iterations read it. The result is nondeterministic.

Operational Checks: Verify Before You Trust

Before committing a DO CONCURRENT rewrite, run these checks:

1. Numerical equivalence test. Compile both the sequential DO and DO CONCURRENT versions. Run on representative inputs and compare results within a stated tolerance (e.g., relative error < 1e-12). Bitwise equality is not expected.
2. Multi-compiler, multi-optimization test. Build with at least two compilers (e.g., gfortran and ifx/nvfortran) at -O2 and -O3. Compare outputs and runtimes. A loop that vectorizes under one compiler but runs sequentially under another reveals a portability gap.
3. Locality audit with DEFAULT(NONE). Add default(none) and declare every variable. The compiler will flag any undeclared reference—forcing you to confront hidden dependencies.
4. Inspect optimization reports or assembly. Look for vectorization messages (-fopt-info-vec for GCC, -qopt-report for Intel) or scan assembly for packed SIMD instructions (vmulps, vaddpd, etc.). If the loop runs scalar, the construct bought you nothing.

Failure Modes: What Goes Wrong in Practice

Failure Mode	Symptom	Mitigation
Floating-point non-associativity	Results differ from sequential by more than rounding error	Define acceptable tolerance before adoption; use REPRODUCIBLE pragma if vendor supports it
Hidden aliasing	Correct at -O0, wrong at -O2/-O3	Prefer allocatable/explicit-shape arrays; avoid pointers and module variables in loop body; use DEFAULT(NONE)
Performance regression	Parallel overhead exceeds compute work	Ensure loop body has enough work (rule of thumb: >100 FLOPs/iteration); consider IF (n > threshold) guard
Compiler ignores parallelism	No speedup, no vectorization report	Check vendor docs for required flags (e.g., -floop-parallelize-all for GCC, -parallel for Intel); consider OpenMP as fallback

Conditions That Would Change the Design

Switch away from DO CONCURRENT when:

• Deterministic bitwise reproducibility is required (e.g., regression tests that compare binary output). Use a plain DO loop or a vendor-specific reproducible mode.
• Cross-iteration communication is needed (e.g., a prefix sum, a stencil with halo exchange). Use coarrays for distributed memory or OpenMP with explicit reduction/atomic clauses for shared memory.
• GPU offload is the target. DO CONCURRENT can offload with vendor support (NVIDIA's nvfortran supports -gpu=autopar on DO CONCURRENT), but OpenMP target directives give more control over data mapping and kernel launch.
• Loop-carried dependencies exist but are benign (e.g., a reduction). OpenMP REDUCTION handles this correctly; DO CONCURRENT does not have a standard reduction clause (though some vendors extend it).

Concrete Example: Stencil Update with Locality Audit

A 5-point 2D stencil on a 1000×1000 grid:

subroutine stencil5(n, m, u, unew)
    integer, intent(in) :: n, m
    real, intent(in) :: u(0:n+1, 0:m+1)
    real, intent(out) :: unew(1:n, 1:m)
    integer :: i, j
    do concurrent (i = 1:n, j = 1:m) default(none) shared(u, unew)
        unew(i,j) = 0.25 * (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1))
    end do
end subroutine

Verification steps (run on Linux with gfortran ≥12 and ifx):

# 1. Compile sequential reference
gfortran -O2 -o stencil_seq stencil.f90
# 2. Compile DO CONCURRENT version
gfortran -O2 -floop-parallelize-all -o stencil_concurrent stencil.f90
# 3. Run both, compare max absolute difference
./stencil_seq > out_seq.txt
./stencil_concurrent > out_concurrent.txt
python3 -c "
import numpy as np
seq = np.loadtxt('out_seq.txt')
con = np.loadtxt('out_concurrent.txt')
print('max abs diff:', np.max(np.abs(seq - con)))
print('rel diff:', np.max(np.abs(seq - con)) / np.max(np.abs(seq)))
"
# 4. Check vectorization report
gfortran -O2 -floop-parallelize-all -fopt-info-vec stencil.f90 2>&1 | grep -i stencil

Expected checks: Relative difference < 1e-12 (typical for double precision). Vectorization report shows loop was vectorized. If gfortran emits no vectorization, try -floop-nest-optimize or switch to a vendor compiler.

Risks: The halo cells u(0,:), u(n+1,:), etc. must be valid—this is a caller responsibility. The DEFAULT(NONE) clause catches any accidental reference to a module variable or undeclared scalar.

Limitations and Practical Verification

• Compiler support varies. Syntax acceptance ≠ parallel execution. GCC 12+ parallelizes with -floop-parallelize-all; Intel ifx does so with -parallel; NVIDIA nvfortran offloads with -gpu=autopar. Test each target compiler.
• Assumed-shape and pointer arrays hide aliasing. Prefer explicit-shape or allocatable arrays with clear ownership. If you must use assumed-shape, ensure the actual argument has no overlapping regions.
• No standard reduction clause. Vendor extensions exist (e.g., Intel's REDUCTION on DO CONCURRENT), but they're non-portable.
• Bitwise differences are expected. Define your numerical tolerance before adopting. A test that asserts bitwise equality will fail.

To verify a specific compiler actually parallelizes: compile a known-independent loop (e.g., a(i) = b(i) + c(i)) at -O3, run with a large n (10^7), and measure speedup over sequential. If speedup ≈ 1×, the compiler didn't parallelize—check flags and reports.