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:
call sger(m, n, alpha, x, 1, y, 1, a, lda)
It is possible to express this operation with array intrinsics:
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, presumably named after the 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
! 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):
$ 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
$ 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)