I'm measuring three approaches to matrix multiplication performance: a naive blocked OpenMP implementation, Eigen, and SGEMM from MKL 2021.4.0. For simplicity all matrices are square, type float, size n x n, aligned at 64-bytes. The compiler is GCC 8.3.1 with compilation flags -msse4.2 -O3 -fopenmp. OS is CentOS 7
I don't understand why MKL SGEMM is the slowest. Why is a naive OpenMP implementation faster than a fancy-optimized library?
Blocked OpenMP (BS = n / 64):
#pragma omp for collapse(2)
for(int i=0; i<n; i++)
for(int j=0; j<n; j++)
C[i*n+j] *= beta;
#pragma omp parallel for schedule(dynamic)
for (int i = 0; i < n; i+=BS)
for (int k = 0; k < n; k+=BS)
for (int j = 0; j < n; j+=BS)
for (int ii = i; ii < i+BS; ii++)
for (int kk = k; kk < k+BS; kk++)
for (int jj = j; jj < j+BS; jj++)
C[ii*n+jj] += alpha*A[ii*n+kk]*B[kk*n+jj];
Eigen
Eigen::Map<const Eigen::MatrixXf> AM(A, n, n);
Eigen::Map<const Eigen::MatrixXf> BM(B, n, n);
Eigen::Map<Eigen::MatrixXf> CM(C, n, n);
CM.noalias() = beta*CM + alpha*(BM * AM); // fortran order!
MKL SGEMM
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
n, n, n, alpha, A, n, B, n, beta, C, n);
The google benchmark results for Intel Xeon Silver 4114, 2 sockets, 2 NUMA nodes:
Benchmark Time CPU Iterations
-------------------------------------------------------------------------------
MatMul/OmpBlk/4096/64/real_time 1132 ms 1038 ms 1
MatMul/OmpBlk/16384/64/real_time 83668 ms 80612 ms 1
MatMul/OmpBlk/32768/64/real_time 1562980 ms 1492184 ms 1
MatMul/Eigen/4096/real_time 878 ms 867 ms 1
MatMul/Eigen/16384/real_time 36140 ms 31629 ms 1
MatMul/Eigen/32768/real_time 259762 ms 246788 ms 1
MatMul/Blas/4096/real_time 4091 ms 3719 ms 1
MatMul/Blas/16384/real_time 219940 ms 219581 ms 1
MatMul/Blas/32768/real_time 1773874 ms 1750015 ms 1
Simple average time of three runs (no google benchmark, one warm-up not included)
OmpBlk/4096: 1452 ms
OmpBlk/16384: 87494 ms
Eigen/4096: 818 ms
Eigen/16384: 34719 ms
Blas/4096: 4060 ms
Blas/16384: 225647 ms
ldd snippet:
libmkl_intel_ilp64.so.1 => /opt/intel/oneapi/mkl/2021.4.0/lib/intel64/libmkl_intel_ilp64.so.1
libmkl_core.so.1 => /opt/intel/oneapi/mkl/2021.4.0/lib/intel64/libmkl_core.so.1
libmkl_intel_thread.so.1 => /opt/intel/oneapi/mkl/2021.4.0/lib/intel64/libmkl_intel_thread.so.1
libiomp5.so => /opt/intel/oneapi/compiler/latest/linux/compiler/lib/intel64/libiomp5.so