Here is a cleaned-up Fortran Discourse post that presents the SYSTEM_CLOCK benchmark and the 1-, 2-, and 4-thread OpenBLAS results clearly.
Title: MATMUL vs cache-friendly loops vs OpenBLAS DGEMM — SYSTEM_CLOCK and thread scaling
Following the earlier discussion about “MATMUL” versus explicit nested loops, I changed the benchmark to use “SYSTEM_CLOCK” rather than “CPU_TIME”, since I wanted to measure elapsed wall-clock time when testing threaded OpenBLAS.
The benchmark compares:
- cache-friendly “j-k-i” nested loops,
- the Fortran intrinsic “MATMUL”,
- OpenBLAS “DGEMM”.
The matrix size is “1000 x 1000”, using “REAL(real64)”, with five timed runs.
I compiled with:
gfortran benchmark2.f90 -O3 -march=native -o benchmark2 -lopenblas
and tested OpenBLAS with 1, 2, and 4 threads:
export OPENBLAS_NUM_THREADS=1
export OMP_NUM_THREADS=1
./benchmark2
export OPENBLAS_NUM_THREADS=2
export OMP_NUM_THREADS=2
./benchmark2
export OPENBLAS_NUM_THREADS=4
export OMP_NUM_THREADS=4
./benchmark2
“SYSTEM_CLOCK” reports:
SYSTEM_CLOCK rate: 1000000000 counts/second
Average wall-clock results
Threads| Nested loops| MATMUL| DGEMM| MATMUL / loops| DGEMM / loops
1| 0.548147 s| 0.167406 s| 0.140996 s| 3.274x| 3.888x
2| 0.548738 s| 0.164382 s| 0.068577 s| 3.338x| 8.002x
4| 0.545909 s| 0.166466 s| 0.146052 s| 3.279x| 3.738x
The corresponding DGEMM/MATMUL ratios were:
1 thread : 1.187x
2 threads: 2.397x
4 threads: 1.140x
An interesting result is that two OpenBLAS threads are by far the fastest configuration on this system.
DGEMM goes from
0.140996 s 1 thread
to
0.068577 s 2 threads
which is a speedup of about 2.06x.
However, increasing the setting to four threads gives:
0.146052 s 4 threads
which is actually slightly slower than the single-thread result.
By contrast, both the explicit loops and “MATMUL” remain almost unchanged:
1 thread 2 threads 4 threads
loops 0.548147 0.548738 0.545909
MATMUL 0.167406 0.164382 0.166466
DGEMM 0.140996 0.068577 0.146052
This also suggests that changing “OPENBLAS_NUM_THREADS” is affecting the explicitly called OpenBLAS “DGEMM”, but not the compiler-generated implementation of “MATMUL” in this particular build.
The explicit loop is deliberately ordered:
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
so that the first Fortran array index varies fastest. This is considerably better than the original “i-j-k” ordering, but “MATMUL” is still about 3.3 times faster.
The unexpected four-thread OpenBLAS result may be particularly interesting. Possible explanations include the number of physical cores actually available to the UserLAnd process, scheduling overhead, CPU frequency/thermal effects, or the characteristics of OpenBLAS on this Android/Linux environment.
Also, “OMP_NUM_THREADS” should not affect the explicit loops here, since the program contains no OpenMP directives and was not compiled with “-fopenmp”. The important variable for this OpenBLAS test is “OPENBLAS_NUM_THREADS”.
Here is the complete program:
program matmul_dgemm_benchmark
use iso_fortran_env, only : real64, int64
implicit none
integer, parameter :: n = 1000
integer, parameter :: nruns = 5
real(real64), allocatable :: A(:,:), B(:,:), C(:,:)
real(real64) :: t_loop(nruns)
real(real64) :: t_matmul(nruns)
real(real64) :: t_dgemm(nruns)
real(real64) :: checksum
integer(int64) :: count1, count2, count_rate
integer :: i, j, k, r
external :: dgemm
allocate(A(n,n), B(n,n), C(n,n))
call random_number(A)
call random_number(B)
call system_clock(count_rate=count_rate)
print *
print *, "Matrix multiplication benchmark"
print *, "Matrix size:", n, "x", n
print *, "Runs:", nruns
print *
print *, "SYSTEM_CLOCK rate:", count_rate, " counts/second"
print *
! Warm up
C = matmul(A,B)
call dgemm('N','N',n,n,n, &
1.0_real64,A,n,B,n, &
0.0_real64,C,n)
! ------------------------------------------------------------
! 1. Cache-friendly explicit loops
! ------------------------------------------------------------
print *, "Nested DO loops"
do r = 1, nruns
C = 0.0_real64
call system_clock(count1)
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
call system_clock(count2)
t_loop(r) = real(count2-count1,real64) / &
real(count_rate,real64)
print '(A,I2,A,F10.6,A)', &
" Run ",r,": ",t_loop(r)," seconds"
end do
checksum = sum(C)
print '(A,ES20.10)', " Checksum = ",checksum
print *
! ------------------------------------------------------------
! 2. Fortran MATMUL
! ------------------------------------------------------------
print *, "Fortran MATMUL"
do r = 1, nruns
call system_clock(count1)
C = matmul(A,B)
call system_clock(count2)
t_matmul(r) = real(count2-count1,real64) / &
real(count_rate,real64)
print '(A,I2,A,F10.6,A)', &
" Run ",r,": ",t_matmul(r)," seconds"
end do
checksum = sum(C)
print '(A,ES20.10)', " Checksum = ",checksum
print *
! ------------------------------------------------------------
! 3. BLAS DGEMM
! C = alpha*A*B + beta*C
! ------------------------------------------------------------
print *, "BLAS DGEMM"
do r = 1, nruns
call system_clock(count1)
call dgemm('N','N',n,n,n, &
1.0_real64,A,n,B,n, &
0.0_real64,C,n)
call system_clock(count2)
t_dgemm(r) = real(count2-count1,real64) / &
real(count_rate,real64)
print '(A,I2,A,F10.6,A)', &
" Run ",r,": ",t_dgemm(r)," seconds"
end do
checksum = sum(C)
print '(A,ES20.10)', " Checksum = ",checksum
print *
! ------------------------------------------------------------
! Summary
! ------------------------------------------------------------
print *, "==============================================="
print *, "Average wall-clock execution times"
print *, "==============================================="
print '(A,F10.6,A)', &
"Nested loops : ",sum(t_loop)/nruns," seconds"
print '(A,F10.6,A)', &
"MATMUL : ",sum(t_matmul)/nruns," seconds"
print '(A,F10.6,A)', &
"DGEMM : ",sum(t_dgemm)/nruns," seconds"
print *
print '(A,F10.3,A)', &
"MATMUL speedup over loops: ", &
sum(t_loop)/sum(t_matmul)," x"
print '(A,F10.3,A)', &
"DGEMM speedup over loops : ", &
sum(t_loop)/sum(t_dgemm)," x"
print '(A,F10.3,A)', &
"DGEMM / MATMUL ratio : ", &
sum(t_matmul)/sum(t_dgemm)," x"
deallocate(A,B,C)
end program matmul_dgemm_benchmark
The main result for me is that the cache-friendly loop ordering reduces the enormous difference seen with the original loop order, while “MATMUL” remains substantially faster. Direct OpenBLAS “DGEMM” is faster still, but its scaling on this particular system peaks at two threads rather than four.
I would be interested to know whether others see a similar difference between 1, 2 and 4 OpenBLAS threads, particularly on ARM systems or under Android/UserLAnd.One useful addition before posting would be the output of nproc and openblas_get_config() (or the OpenBLAS package/version), because the 2-thread DGEMM = 0.0686 s versus 4-thread = 0.1461 s result is the most interesting part and other Fortran Discourse users will probably ask about the available CPU cores and OpenBLAS build.
