Welcome to Ubuntu in UserLAnd!
userland@localhost:~$ nano arraytest.f
userland@localhost:~$ gfortran arraytest.f -ffree-form
userland@localhost:~$ ./a.out
matmul time: 8.40620026E-02 seconds
do loops time: 7.21938992 seconds
userland@localhost:~$ cat arraytest.f
program matmul_speed_test
implicit none
integer, parameter :: N = 1000
real :: A(N,N), B(N,N), C(N,N)
integer :: i, j, k
real :: start, finish
! Initialize matrices
call random_number(A)
call random_number(B)
! matmul
call cpu_time(start)
C = matmul(A, B)
call cpu_time(finish)
print *, "matmul time: ", finish - start, " seconds"
! explicit do loops
C = 0.0
call cpu_time(start)
do i = 1, N
do j = 1, N
do k = 1, N
C(i,j) = C(i,j) + A(i,k) * B(k,j)
end do
end do
end do
call cpu_time(finish)
print *, "do loops time: ", finish - start, " seconds"
end program matmul_speed_test
The main problem is loop ordering. Fortran stores arrays in column-major order, so the first index should vary fastest. Your original inner loop varies k, causing A(i,k) to be accessed with a stride of n.
Compiling the program it gave with gfortran xmatmul_speed.f90 -march=native -O3 and running gives on my PC
Here is an interesting exercise for people who are experimenting with this for the first time. There are six different ways to order the three do loops (3!=6). It is not too much effort to write a code that times all six versions. Then after you run the code and look at the timings, try to understand why each of the loop orders behaves the way it does. You will see that the innermost loop is either a dot product or an axpy operation, and that is important. Does your code compile to use fused multiply-add operations? Does it use vector instructions? There are also issues involving array strides and cache use. Consider the same thing for the transpose(A) times B operation; surprisingly that is very different. What about subblocking the matrices; does that help or hurt performance? Understanding all of that information will prove useful for the rest of your programming career.
Intrinsic Fortran “MATMUL” versus nested “DO” loops
I performed a simple benchmark comparing Fortran’s intrinsic “MATMUL” function with a direct implementation using three nested “DO” loops.
My purpose is to understand:
Why “MATMUL” is much faster in this example.
Whether the compiler is replacing “MATMUL” with an optimized matrix-multiplication routine.
How much of the difference is caused by the ordering of the explicit loops.
What constitutes a fair benchmark between an intrinsic operation and manually written loops.
Whether explicit loops can approach the performance of “MATMUL” after loop reordering, blocking, vectorization, or parallelization.
I compiled the program using:
gfortran arraytest.f -ffree-form -O3
The result on my system was:
matmul time: 8.05520043E-02 seconds
do loops time: 2.58868003 seconds
Thus, in this particular run, “MATMUL” completed in approximately 0.081 seconds, while the explicit loops required approximately 2.59 seconds. This makes the explicit implementation about 32 times slower.
The test program was:
program matmul_speed_test
implicit none
integer, parameter :: n = 1000
real :: a(n,n), b(n,n), c(n,n)
integer :: i, j, k
real :: start, finish
! Initialize matrices.
call random_number(a)
call random_number(b)
! Intrinsic MATMUL.
call cpu_time(start)
c = matmul(a, b)
call cpu_time(finish)
print *, "matmul time: ", finish - start, " seconds"
! Explicit nested loops.
c = 0.0
call cpu_time(start)
do i = 1, n
do j = 1, n
do k = 1, n
c(i,j) = c(i,j) + a(i,k) * b(k,j)
end do
end do
end do
call cpu_time(finish)
print *, "do loops time:", finish - start, " seconds"
end program matmul_speed_test
I understand that Fortran arrays are stored in column-major order. The loop ordering above may therefore give poor memory-access behaviour, particularly when accessing “b(k,j)” and repeatedly updating a single element of “c”.
A potentially better explicit ordering may be:
c = 0.0
do j = 1, n
do k = 1, n
do i = 1, n
c(i,j) = c(i,j) + a(i,k) * b(k,j)
end do
end do
end do
Here, “i” is the innermost loop, so a(i,k) and c(i,j) are accessed consecutively in memory.
For a more reliable comparison, I also intend to:
run each version several times;
include a warm-up calculation;
verify that both versions produce equivalent results;
prevent the compiler from eliminating an unused result;
report the compiler version and processor;
examine compiler optimization reports;
compare different loop orderings;
test compilation with and without an optimized BLAS library.
My main questions are:
Is the original comparison meaningful, or is the explicit loop ordering too inefficient to provide a fair test?
Does gfortran normally implement “MATMUL” inline, call an internal optimized routine, or use an external BLAS library?
Which compiler options are useful for determining whether the loops were vectorized?
What explicit-loop implementation would provide the fairest comparison with “MATMUL”?
Should elapsed wall-clock time be used instead of “CPU_TIME”, particularly when threaded libraries may be involved?
I would be interested in results from other compilers, optimization settings, processors, and BLAS implementations.One important correction is that the original i-j-k loop order is particularly unfavourable for Fortran. The revised post acknowledges that issue rather than presenting the measured 32-times difference as a general result.
I will try to answer some of your questions, as I have investigated why the Gfortran MATMUL is more efficient than the basic DO loop approaches I have tested.
I have discussed this with Thomas Koenig, who was involved in its development.
Basically, they use a 4x4 submatrix multiplication where the sub-matrices are sized to be stored in L2 cache. This is a further refinement that the do loop order suggested above, which tries to address L3 cache optimisation.
It is my “understanding” of AVX instructions that full efficiency is only achieved when the vectors are available in L2 cache, which is assisted if the extended vectors are available in L3 cache. There are no Fortran instructions to manage this, only strategies mainly of sequential memory addressing.
Addressing L2 cache usage can be a much more complex approach(due to smaller size), while the L3 cache approach can be as simple as ordering do loops and ordering of array subscripts.
A further “L3cache” improvement to matrix multiplication can be to calculate [C] = [A]transpose x [B] using dot_product inner loop, so storing [A] transpose (or recognising A symmetry) can improve cache usage.
Is the original comparison meaningful, or is the explicit loop ordering too inefficient to provide a fair test? Loop ordering is a useful approach, by optimising sequential memory addressing in the inner loop.
Does gfortran normally implement “MATMUL” inline, call an internal optimized routine, or use an external BLAS library? I don’t know, but give the complexity of MATMUL code this is unlikely. I assume optimised BLAS libraries would address L2 efficiency.
Which compiler options are useful for determining whether the loops were vectorized? I don’t know.
What explicit-loop implementation would provide the fairest comparison with “MATMUL”? Perhaps you could investigate a sub-matrix approach, which becomes more difficult for a general value of N. Multi-threaded approach is also of use for large N, although cache efficiency and memory access bandwidth are not easily solved with Fortran code.
Should elapsed wall-clock time be used instead of “CPU_TIME”, particularly when threaded libraries may be involved? Definitely, SYSTEM_CLOCK is a much better measure of efficiency.
Although cache access efficiency and memory access bandwidth are significant issues to address when looking at MATMUL (and OpenMP), the idea of controlling what is in the cache from a Fortran program would be a terrible capability to have. It is best to try to develop strategies and leave those controls to the OS.
The concept of sequential access of memory has been key for many version of computer hardware and OS. I first learnt this when listening to virtual memory paging systems reducing large sets of linear equations, where you could write a new solver and test it before the old solver finished. My first introduction to this problem was when my boss ran a FEA problem using “BNDSOL” on a Pr1me 300 that ran for 5 days. It did however produce an answer !
An (untested) example of using sub-matrices, without unrolling the inner loop is:
! for sub-matrix size s, this could be achieved with
do j = 1, N
do k2 = 1,n,s
k2 = min (k1+s-1,n)
do i1 = 1,n,s
i2 = min(i1+s-1,N)
do k = k1, k2
! C(i1:i2,j) = C(i1:i2,j) + A(i1:i2,k) * B(k,j)
do i = i1, i2
C(i,j) = C(i,j) + A(i,k) * B(k,j)
end do
end do
end do
end do
end do
program matmul_speed_test
implicit none
integer, parameter :: N = 1000
real*8 :: A(N,N), B(N,N), C(N,N), D(N,N)
integer :: i, j, k, si=8, sk=8, k1,k2,i1,i2
real :: start, finish
! Initialize matrices
start = elapse_sec ()
call random_number(A)
call random_number(B)
finish = elapse_sec () - start
print *, "initialise time: ", finish, " seconds"
! matmul
start = elapse_sec ()
D = matmul (A, B)
finish = elapse_sec () - start
print *, "matmul time: ", finish, " seconds"
! explicit do loops
C = 0.0
start = elapse_sec ()
do i = 1, N
do j = 1, N
do k = 1, N
C(i,j) = C(i,j) + A(i,k) * B(k,j)
end do
end do
end do
finish = elapse_sec () - start
print *, "ijk loops time: ", finish, " seconds : err=", max_err ( c, d )
! better do loops
C = 0.0
start = elapse_sec ()
do j = 1, N
do k = 1, N
do i = 1, N
C(i,j) = C(i,j) + A(i,k) * B(k,j)
end do
end do
end do
finish = elapse_sec () - start
print *, "jki loops time: ", finish, " seconds : err=", max_err ( c, d )
! array do loops
C = 0.0
start = elapse_sec ()
do j = 1, N
do k = 1, N
C(:,j) = C(:,j) + A(:,k) * B(k,j)
end do
end do
finish = elapse_sec () - start
print *, "jk: loops time: ", finish, " seconds : err=", max_err ( c, d )
! for sub-matrix size s, this could be achieved with
C = 0.0
start = elapse_sec ()
do j = 1, N
do k1 = 1,N,sk
k2 = min (k1+sk-1,N)
do i1 = 1,N,si
i2 = min(i1+si-1,N)
do k = k1, k2
C(i1:i2,j) = C(i1:i2,j) + A(i1:i2,k) * B(k,j)
end do
end do
end do
end do
finish = elapse_sec () - start
print *, "sub loops time: ", finish, " seconds : err=", max_err ( c, d )
write (*,*) 'si =',si
write (*,*) 'sk =',sk
contains
real function elapse_sec ()
integer*8 :: tick, rate, last = -1
real :: sec
call system_clock ( tick, rate)
if ( last < 0 ) last = tick
sec = dble (tick-last) / dble (rate)
last = tick
elapse_sec = sec
end function elapse_sec
real*8 function max_err ( c, d )
real*8 :: c(:,:), d(:,:)
integer :: i, j
real*8 x
x = 0
do i = 1,size(c,2)
do j = 1,size(c,1)
x = max ( x, abs(c(j,i)-d(j,i)) )
end do
end do
max_err = x
end function max_err
end program matmul_speed_test
Dos batch file
set options=%1.f90 -O3 -march=native -ffast-math -o %1.exe
gfortran %options%
%1 >> %1.log
type %1.log
Certainly the loop order does improve with -O3 -ffast-math
si=8 for 64-bit reals is best on my AMD hardware, but could vary with other L2 cache sizes.
Surprised that sub-matrix did not approach the jki loop order, but that could be due to Gfortran not fully optimising " C(i1:i2,j) = C(i1:i2,j) + A(i1:i2,k) * B(k,j)" ?
The difference between “jk:” and “sub” for Gfortran vs Flang is very interesting.
This reflects the different optimisation targets of the two compilers.
Can anyone test in ifx/ifort ?
On AMD Ryzen™ AI 7 PRO 350, so results might not be optimal (due to ifort historically changing its behavior on non-Intel CPUs).
The -fast flag caused an ICE due to the implied -static flag (the libc.a on Debian 13 causes a mismatch during linking), so I set the other “fast” flags explicitly.
(ins)$ ifort -V -ipo -O3 -no-prec-div -fp-model fast=2 -xHost -diag-disable=10448 -static-intel mat.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.13.1 Build 20240703_000000
Copyright (C) 1985-2024 Intel Corporation. All rights reserved.
Intel(R) Fortran 2021.13.1-1693
GNU ld (GNU Binutils for Debian) 2.44
(ins)$ time ./a.out
initialise time: 2.0700000E-02 seconds
matmul time: 0.4088450 seconds
ijk loops time: 0.4068690 seconds : err= 0.000000000000000E+000
jki loops time: 0.4059490 seconds : err= 0.000000000000000E+000
jk: loops time: 0.4081930 seconds : err= 0.000000000000000E+000
sub loops time: 0.4896010 seconds : err= 0.000000000000000E+000
si = 8
sk = 8
real 0m2.158s
user 0m2.150s
sys 0m0.008s
For ifx:
Debian 13 VM
The -fast flag refused to work since it doesn’t detect a supported CPU, so, again, the “fast” flags are explicit.
$ ifx -V -ipo -O3 -no-prec-div -fp-model fast=2 -xHost -static-intel mat.f90
Intel(R) Fortran Compiler for applications running on Intel(R) 64, Version 2026.0.0 Build 20260331
Copyright (C) 1985-2026 Intel Corporation. All rights reserved.
Intel(R) Fortran 26.0-1156
GNU ld (GNU Binutils for Debian) 2.44
$ time ./a.out
initialise time: 1.1572000E-02 seconds
matmul time: 3.5103999E-02 seconds
ijk loops time: 2.7476000E-02 seconds : err= 0.000000000000000E+000
jki loops time: 2.7339000E-02 seconds : err= 0.000000000000000E+000
jk: loops time: 2.6957000E-02 seconds : err= 0.000000000000000E+000
sub loops time: 3.5011001E-02 seconds : err= 0.000000000000000E+000
si = 8
sk = 8
real 0m0.178s
user 0m0.165s
sys 0m0.012s
I noticed that the operations are still only computing one column at a time of the result. That is, the matrix-matrix product is being computed as a series of matrix-vector products. I thought there might be a way to improve the times a little because of that. I thought I would try subblocking all three matrices. When you do this in the code, you end up going from three nested loops to six nested loops, so there are now not just 3!=6 loop orders, but (3!)^2=36 loop orders (assuming all the subblock loops are grouped together). So things are now much more complicated. Here is a modified version of the above code with four of those possibilities.
program matmul_speed_test
use, intrinsic :: iso_fortran_env, only: int64, wp => real64
implicit none
integer, parameter :: N = 1000
real(wp) :: A(N,N), B(N,N), C(N,N), D(N,N)
integer :: i, j, k, i1, i2, j1, j2, k1, k2
integer, parameter :: si = 8, sk = 8, SB = 4
real :: start, finish
real(wp) :: AB(SB,SB), BB(SB,SB), CB(SB,SB) ! local subblock arrays.
character(*), parameter :: cfmt = '(a,t18,f8.4,a,es0.1e2,*(g0:1x))'
if ( mod(N,SB) .ne. 0 ) stop 'SB must divide N'
! Initialize matrices
start = elapse_sec ()
call random_number(A)
call random_number(B)
finish = elapse_sec () - start
print cfmt, "initialise time: ", finish, " seconds"
! matmul. This result is used below for error reference.
start = elapse_sec ()
D = matmul (A, B)
finish = elapse_sec () - start
print cfmt, "matmul time: ", finish, " seconds"
! explicit do loops, dot product inner loop, noncontiguous array order.
C = 0.0
start = elapse_sec ()
do i = 1, N
do j = 1, N
do k = 1, N
C(i,j) = C(i,j) + A(i,k) * B(k,j)
end do
end do
end do
finish = elapse_sec () - start
print cfmt, "ijk loops time: ", finish, " seconds : err=", max_err ( C, D )
! explicit do loops, axpy inner loop, worst possible array access.
C = 0.0
start = elapse_sec ()
do i = 1, N
do k = 1, N
do j = 1, N
C(i,j) = C(i,j) + A(i,k) * B(k,j)
end do
end do
end do
finish = elapse_sec () - start
print cfmt, "ikj loops time: ", finish, " seconds : err=", max_err ( C, D )
! axpy inner loop, contiguous array order.
C = 0.0
start = elapse_sec ()
do j = 1, N
do k = 1, N
do i = 1, N
C(i,j) = C(i,j) + A(i,k) * B(k,j)
end do
end do
end do
finish = elapse_sec () - start
print cfmt, "jki loops time: ", finish, " seconds : err=", max_err ( C, D )
! array do loops, axpy inner loop, contiguous array order.
C = 0.0
start = elapse_sec ()
do j = 1, N
do k = 1, N
C(:,j) = C(:,j) + A(:,k) * B(k,j)
end do
end do
finish = elapse_sec () - start
print cfmt, "jk: loops time: ", finish, " seconds : err=", max_err ( C, D )
! for sub-matrix size s, this could be achieved with
C = 0.0
start = elapse_sec ()
do j = 1, N
do k1 = 1,N,sk
k2 = min (k1+sk-1,N)
do i1 = 1,N,si
i2 = min(i1+si-1,N)
do k = k1, k2
C(i1:i2,j) = C(i1:i2,j) + A(i1:i2,k) * B(k,j)
end do
end do
end do
end do
finish = elapse_sec () - start
print cfmt, "sub loops time: ", finish, " seconds : err=", max_err ( C, D ), ' si =',si,' sk =',sk
! matrix subblock version with axpy inner loop, contiguous *B array order.
C = 0.0
start = elapse_sec ()
do j1 = 1, N, SB
j2 = j1 + SB - 1
do k1 = 1,N,SB
k2 = k1 + SB - 1
BB = B(k1:k2,j1:j2)
do i1 = 1,N,SB
i2 = i1 + SB - 1
AB = A(i1:i2,k1:k2)
CB = 0.0
do j = 1, SB
do k = 1, SB
do i = 1, SB
CB(i,j) = CB(i,j) + AB(i,k) * BB(k,j)
enddo
enddo
enddo
C(i1:i2,j1:j2) = C(i1:i2,j1:j2) + CB
enddo
enddo
end do
finish = elapse_sec () - start
print cfmt, "AXPY SB time: ", finish, " seconds : err=", max_err ( C, D ), ' SB =', SB
! accumulate matrix subblock version with axpy inner loop, contiguous *B array access.
start = elapse_sec ()
do j1 = 1, N, SB
j2 = j1 + SB - 1
do i1 = 1, N, SB
i2 = i1 + SB - 1
CB = 0.0
do k1 = 1,N,SB
k2 = k1 + SB - 1
BB = B(k1:k2,j1:j2)
AB = A(i1:i2,k1:k2)
do j = 1, SB
do k = 1, SB
do i = 1, SB
CB(i,j) = CB(i,j) + AB(i,k) * BB(k,j)
enddo
enddo
enddo
enddo
C(i1:i2,j1:j2) = CB
enddo
end do
finish = elapse_sec () - start
print cfmt, "Acc AXPY SB time: ", finish, " seconds : err=", max_err ( C, D ), ' SB =', SB
! matrix subblock version with dot product inner loop, contiguous *B array order.
C = 0.0
start = elapse_sec ()
do j1 = 1, N, SB
j2 = j1 + SB - 1
do k1 = 1,N,SB
k2 = k1 + SB - 1
BB = B(k1:k2,j1:j2)
do i1 = 1,N,SB
i2 = i1 + SB - 1
AB = transpose(A(i1:i2,k1:k2))
CB = 0.0
do j = 1, SB
do i = 1, SB
do k = 1, SB
CB(i,j) = CB(i,j) + AB(k,i) * BB(k,j)
enddo
enddo
enddo
C(i1:i2,j1:j2) = C(i1:i2,j1:j2) + CB
enddo
enddo
end do
finish = elapse_sec () - start
print cfmt, "DOT SB time: ", finish, " seconds : err=", max_err ( C, D ), ' SB =', SB
! accumulate matrix subblock version with dot product inner loop, contiguous *B array access.
start = elapse_sec ()
do j1 = 1, N, SB
j2 = j1 + SB - 1
do i1 = 1, N, SB
i2 = i1 + SB - 1
CB = 0.0
do k1 = 1,N,SB
k2 = k1 + SB - 1
BB = B(k1:k2,j1:j2)
AB = TRANSPOSE( A(i1:i2,k1:k2) )
do j = 1, SB
do i = 1, SB
do k = 1, SB
CB(i,j) = CB(i,j) + AB(k,i) * BB(k,j)
enddo
enddo
enddo
enddo
C(i1:i2,j1:j2) = CB
enddo
end do
finish = elapse_sec () - start
print cfmt, "Acc DOT SB time: ", finish, " seconds : err=", max_err ( C, D ), ' SB =', SB
contains
real function elapse_sec ()
integer(int64) :: tick, rate, last = -1
real :: sec
call system_clock ( tick, rate)
if ( last < 0 ) last = tick
sec = dble (tick-last) / dble (rate)
last = tick
elapse_sec = sec
end function elapse_sec
real(wp) function max_err ( C, D )
real(wp) :: c(:,:), d(:,:)
integer :: i, j
real(wp) :: x
x = 0
do i = 1,size(c,2)
do j = 1,size(c,1)
x = max ( x, abs(c(j,i)-d(j,i)) )
end do
end do
max_err = x
end function max_err
end program matmul_speed_test
This is with Apple M2 arm64 hardware. There are a couple of things that surprised me. First, I added another do loop order denoted “ikj” in the above listing. I did not add this to improve the times, I added it because I thought that was the worst possible loop order, so I wanted to see just how bad things could get by going through ALL of the arrays in the wrong loop order. The surprising thing is that this was not the slowest time, the original “ijk” loop order always seemed to be the worst. I think the compiler must recognize just how bad this “ikj” code is, and it rearranges the operations somehow. Even with -O0 that supposedly worst loop order was not the slowest. I’m not sure what is happening here, it seems like that code should perform even worse than it does.
Another odd thing is related to the transpose() operations in the code. I first wrote the subblock version without those just to make sure I was getting the correct output. Of course the inner loops do not access the *B arrays contiguously in those cases. When I added the transpose() in order to get contiguous array element references, I was expecting improved timings. However, that did not occur, the timings remained the same. I left the transpose operations in the above code, but I wanted to note that for some reason that I do not understand, they are not necessary to get good timings. Maybe the compiler optimizer is helping out the code in these cases too? Maybe this is an L1 cache feature?
Another thing that surprised me was the timings for the two accumulate loop orders. The idea here was that the output array C(:,:) needs to be referenced less often, but at the expense of more references to the other arrays. Those other array references are just copies, whereas the C(:,:) references are both fetch and store, so I thought the accumulation loop orders would be a little faster. The actual timings show that with -O0 that prediction is correct, the accumulation timings are about 10% faster, but with -O3 the accumulation timings are about 2x slower! I again assume that the compiler optimizer is doing some magic here, but I don’t really know what it is.
You might wonder, if the compiler can generate code that is some 30% faster than the intrinsic MATMUL(), why isn’t the intrinsic faster? I think the answer is that the above fortran code is cheating a little bit. Note that there is a test at the beginning of the code to make sure that MOD(N,SB)==0. This means that the fortran code doesn’t need to test for border cases. The general MATMUL() intrinsic must handle all possible matrix dimensions, so it must test for those cases. Also, the intrinsic MATMUL() must also test for the matrix-vector, vector-matrix, and vector-vector product special cases. So I think that must be what is going on.
The fortran programmer does not have direct control of the cache, but there is some indirect control. That is what the AB(:,:),BB(:,:), and CB(:,:) arrays do in the above code. By copying the original noncontiguous array subblocks into those local arrays, and then referencing only those copies in the innermost loops, the fortran programmer can indirectly put those array subblocks into L1 cache.
Finally one other comment. If you add the option -fexternal-blas -framework accelerate, then you get the following:
Here the matmul() intrinsic blows the doors off any of the fortran codes by linking eventually to the dgemm() BLAS-3 subroutine. That is a multithreaded code that uses the vector instructions. Also there are now some 0.0 errors computed, but I think that is just coincidence; all these errors are in the least significant bits, so they depend on the order that the roundoff errors occur as the results are computed. The only way I know to consistently beat this matmul() time is to directly access dgemm() from the fortran code.