Intrinsic Fortran function matmul has speed advantage over do loops

1

Caution, skipping hostkey check for localhost

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

userland@localhost:~$

2

As ChatGPT says,

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

 matmul time:                   0.03125000   seconds
 original i-j-k loops time:     1.26562500   seconds
 cache-friendly j-k-i time:     0.10937500   seconds
 blocked loops time:            0.09375000   seconds

so matmul is still faster than loops, but not by 2 orders of magnitude.

1 Like
3

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.

1 Like
4

A simple change is to reorder the loops to achieve better use of memory

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

To emphasise this, the following array syntax may help

do j = 1, N
    do k = 1, N
        C(:,j) = C(:,j) + A(:,k) * B(k,j)
    end do
end do
1 Like