# Please, No More Loops (Than Necessary): New Patterns in Fortran 2023

**URL:** <https://fortran-lang.discourse.group/t/please-no-more-loops-than-necessary-new-patterns-in-fortran-2023/10605>\
**Category:** Uncategorized\
**Created:** [January 7, 2026, 8:08pm UTC](https://fortran-lang.discourse.group/t/please-no-more-loops-than-necessary-new-patterns-in-fortran-2023/10605 "2026-01-07T20:08:38Z")\
**Posts on this page:** 1\
**Showing post:** 36

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [January 25, 2026, 2:18am UTC](https://fortran-lang.discourse.group/t/please-no-more-loops-than-necessary-new-patterns-in-fortran-2023/10605/36 "2026-01-25T02:18:14Z")

</div>

> [@zaikunzhang](#):
>
> I also take the viewpoint that loops should be replaced with array (matrix/vector) operations whenever possible.

How would you evaluate the general rank-1 upate A = A + \alpha x y^T, where A is a matrix and x and y are vectors? This is one of the building blocks of the LU factorization algorithm.

The rank-1 update is encapsulated by the BLAS Level 2 call:

```fortran
call sger(m, n, alpha, x, 1, y, 1, a, lda)

```

It is possible to express this operation with array intrinsics:

```fortran
a(1:m,1:n) = a(1:m,1:n) + spread(x,2,n) * spread(y,1,m)

```

(This statement was used in the data parallel [maspar BLAS](https://netlib.org/maspar/f90/sblas/sger.f), presumably named after the [MasPar](https://en.wikipedia.org/wiki/MasPar) parallel computer.)

I used an LLM to generate a micro-benchmark measuring 7 different variants (loops, `matmul` and `spread`) and an 8th variant which just forwards the call to BLAS:

> **dger\_abstraction\_penalty.f90**
>
> ```fortran
> ! dger_abstraction_penalty.f90 -- Evaluate speed of SGER operation
> !
> ! Most of this example has been generated using Google Gemini.
> ! Ivan Pribec, 25/01/2026
> 
> module sger_kernels
> implicit none
> contains
> ! --- KERNELS ---
> 
> subroutine sger_nested(m, n, alpha, x, y, a)
> integer, intent(in) :: m, n
> real, intent(in) :: alpha, x(m), y(n)
> real, intent(inout) :: a(m, n)
> integer :: i, j
> do j = 1, n
> do i = 1, m
> a(i, j) = a(i, j) + alpha * x(i) * y(j)
> end do
> end do
> end subroutine
> 
> subroutine sger_col(m, n, alpha, x, y, a)
> integer, intent(in) :: m, n
> real, intent(in) :: alpha, x(m), y(n)
> real, intent(inout) :: a(m, n)
> integer :: j
> do j = 1, n
> a(:, j) = a(:, j) + (alpha * y(j)) * x
> end do
> end subroutine
> 
> subroutine sger_concurrent(m, n, alpha, x, y, a)
> integer, intent(in) :: m, n
> real, intent(in) :: alpha, x(m), y(n)
> real, intent(inout) :: a(m, n)
> integer :: i, j
> do concurrent (j=1:n, i=1:m)
> a(i, j) = a(i, j) + alpha * x(i) * y(j)
> end do
> end subroutine
> 
> subroutine sger_ptr(m, n, alpha, x, y, a)
> integer, intent(in) :: m, n
> real, intent(in) :: alpha
> real, intent(in), target :: x(m), y(n)
> real, intent(inout) :: a(m, n)
> real, pointer :: xm(:,:), ym(:,:)
> xm(1:m, 1:1) => x
> ym(1:1, 1:n) => y
> a = a + alpha * matmul(xm, ym)
> end subroutine
> 
> subroutine sger_ptr_wrap(m, n, alpha, x, y, a)
> integer, intent(in) :: m, n
> real, intent(in) :: alpha, x(m), y(n)
> real, intent(inout) :: a(m, n)
> call sger_ptr(m, n, alpha, x, y, a)
> end subroutine
> 
> subroutine sger_reshape(m, n, alpha, x, y, a)
> integer, intent(in) :: m, n
> real, intent(in) :: alpha, x(m), y(n)
> real, intent(inout) :: a(m, n)
> a = a + alpha * matmul(reshape(x, [m, 1]), reshape(y, [1, n]))
> end subroutine
> 
> subroutine sger_spread(m, n, alpha, x, y, a)
> integer, intent(in) :: m, n
> real, intent(in) :: alpha, x(m), y(n)
> real, intent(inout) :: a(m, n)
> a = a + alpha * spread(x, 2, n) * spread(y, 1, m)
> end subroutine
> 
> elemental function ef(al, val, xi, yj)
> real, intent(in) :: al, val, xi, yj
> real :: ef
> ef = val + al * xi * yj
> end function
> 
> subroutine sger_elem(m, n, alpha, x, y, a)
> integer, intent(in) :: m, n
> real, intent(in) :: alpha, x(m), y(n)
> real, intent(inout) :: a(m, n)
> a = ef(alpha, a, spread(x, 2, n), spread(y, 1, m))
> end subroutine
> 
> subroutine sger_blas(m, n, alpha, x, y, a)
> integer, intent(in) :: m, n
> real, intent(in) :: alpha, x(m), y(n)
> real, intent(inout) :: a(m, n)
> external :: sger
> call sger(m,n,alpha,x,1,y,1,a,m)
> end subroutine
> 
> end module
> 
> program abstraction_penalty
> use sger_kernels
> implicit none
> integer, parameter :: M = 100, N = 100
> integer, parameter :: TRIALS = 10
> integer, parameter :: NUM_VARIANTS = 8
> 
> real(8), parameter :: G_OPS = (2.0d0 * M * N) / 1.0d9
> 
> real, allocatable, target :: A(:,:), x(:), y(:)
> real :: alpha, ref_time
> real :: times(NUM_VARIANTS), ratios(NUM_VARIANTS), gflops(NUM_VARIANTS)
> character(len=25) :: names(NUM_VARIANTS)
> 
> character(len=*), parameter :: fmt = '(A25, F12.4, F12.3)'
> 
> allocate(A(M, N), x(M), y(N))
> call random_number(x); call random_number(y); alpha = 1.1
> 
> names = [character(len=25) :: "Classic Nested Loops", &
> "Column-wise (1-loop)", &
> "Do Concurrent", &
> "Pointer Remap + Matmul", &
> "Reshape + Matmul", &
> "Spread (High-Level)", &
> "Elemental (Mapped)",&
> "BLAS"]
> 
> print '(A, I0, A, I0, A)', "Benchmarking Matrix ", M, " x ", N
> print '(A)', "-----------------------------------------------------------------------"
> print '(A25, A12, A12, A12)', "Variant", "Avg Time(s)", "Speedup (x)", "GFLOPS"
> 
> ! Reference run
> ref_time = run_kernel(sger_nested)
> 
> call record_stats(1, ref_time)
> call record_stats(2, run_kernel(sger_col))
> call record_stats(3, run_kernel(sger_concurrent))
> call record_stats(4, run_kernel(sger_ptr_wrap))
> call record_stats(5, run_kernel(sger_reshape))
> call record_stats(6, run_kernel(sger_spread))
> call record_stats(7, run_kernel(sger_elem))
> call record_stats(8, run_kernel(sger_blas))
> 
> print '(A)', "-----------------------------------------------------------------------"
> print '(A, F12.4)', "Global Abstraction Penalty (Log-Avg Speedup): ", &
> exp(sum(log(ratios)) / NUM_VARIANTS)
> 
> contains
> 
> subroutine record_stats(idx, t_avg)
> integer, intent(in) :: idx
> real, intent(in) :: t_avg
> times(idx) = t_avg
> ratios(idx) = ref_time / t_avg
> gflops(idx) = real(G_OPS / t_avg)
> print '(A25, E12.4, F12.3, G12.3)', &
> names(idx), times(idx), ratios(idx), gflops(idx)
> end subroutine
> 
> function run_kernel(kernel) result(avg_t)
> interface
> subroutine kernel(m, n, alpha, x, y, a)
> integer, intent(in) :: m, n
> real, intent(in) :: alpha, x(m), y(n)
> real, intent(inout) :: a(m, n)
> end subroutine
> end interface
> integer :: current_trials, t
> real(kind(1.0d0)) :: tstart, tend, total
> real :: avg_t
> 
> current_trials = TRIALS
> do
> A = 0.0
> call cpu_time(tstart)
> do t = 1, current_trials
> call kernel(M, N, alpha, x, y, A)
> end do
> call cpu_time(tend)
> total = tend - tstart
> if (total >= 0.1) exit
> if (current_trials > 1000) exit
> current_trials = current_trials * 2
> end do
> avg_t = total / current_trials
> end function
> 
> end program
> 
> ```

The results I got on an Apple M2 processor (using BLAS from Apple Accelerate):

```txt
$ gfortran -O3 -mcpu=native dger_abstraction_penalty.F90 -lblas && ./a.out
Benchmarking Matrix 100 x 100
-----------------------------------------------------------------------
                  Variant Avg Time(s) Speedup (x) GFLOPS
Classic Nested Loops 0.8633E-06 1.000 23.2
Column-wise (1-loop) 0.8570E-06 1.007 23.3
Do Concurrent 0.8641E-06 0.999 23.1
Pointer Remap + Matmul 0.2793E-05 0.309 7.16
Reshape + Matmul 0.4941E-05 0.175 4.05
Spread (High-Level) 0.1045E-04 0.083 1.91
Elemental (Mapped) 0.1044E-04 0.083 1.91
BLAS 0.5812E-06 1.485 34.4
-----------------------------------------------------------------------
Global Abstraction Penalty (Log-Avg Speedup): 0.3914

```

```txt
$ flang -O3 -mcpu=native dger_abstraction_penalty.F90 -lblas && ./a.out
Benchmarking Matrix 100 x 100
-----------------------------------------------------------------------
                  Variant Avg Time(s) Speedup (x) GFLOPS
Classic Nested Loops 0.4844E-06 1.000 41.3
Column-wise (1-loop) 0.4742E-06 1.021 42.2
Do Concurrent 0.4734E-06 1.023 42.2
Pointer Remap + Matmul 0.2844E-05 0.170 7.03
Reshape + Matmul 0.2849E-05 0.170 7.02
Spread (High-Level) 0.1385E-03 0.003 0.144
Elemental (Mapped) 0.1404E-03 0.003 0.142
BLAS 0.5578E-06 0.868 35.9
-----------------------------------------------------------------------
Global Abstraction Penalty (Log-Avg Speedup): 0.1540

```

For this particular example, using `spread` or `matmul` is slower than using a `do` loop. The BLAS result on the other hand is consistent when switching compilers.

I suppose the `matmul` results could be improved by using a special kernel for the scenario `(m x 1) * (1 x n)` (inner dimension of product is 1). The same goes for `spread() * spread()`. But it is a lot extra work to add such corner cases to a compiler.

(Note: I’ve brought this example up previously - [`target` attribute seemingly affecting performances - #25 by ivanpribec](https://fortran-lang.discourse.group/t/target-attribute-seemingly-affecting-performances/9137/25))

---

_[View the full topic](https://fortran-lang.discourse.group/t/please-no-more-loops-than-necessary-new-patterns-in-fortran-2023/10605)._