Matrix multiplication speed BLAS/Eigen and Numpy

I am benchmarking matrix multiplication for different libraries as I am thinking of rewriting some cython code to native c++. However the simple tests seem to imply that numpy is faster than BLAS or eigen for simple matrix multiplications.

I have written the following files:

#!test_blas.cpp
#include <random>
#include <cstdio>
#include <stdlib.h>
#include <iostream>
#include <cblas.h>

int main ( int argc, char* argv[] ) {

    // Random numbers
    std::mt19937_64 rnd;
    std::uniform_real_distribution<double> doubleDist(0, 1);

    // Create arrays that represent the matrices A,B,C
    const int n = 2000;
    double*  A = new double[n*n];
    double*  B = new double[n*n];
    double*  C = new double[n*n];

    // Fill A and B with random numbers
    for(uint i =0; i <n; i++){
        for(uint j=0; j<n; j++){
            A[i*n+j] = doubleDist(rnd);
            B[i*n+j] = doubleDist(rnd);
        }
    }

    // Calculate A*B=C
    clock_t start = clock();
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1.0, A, n, B, n, 0.0, C, n);
    clock_t end = clock();
    double time = double(end - start)/ CLOCKS_PER_SEC;
    std::cout<< "Time taken : " << time << std::endl; 
    // Clean up
    delete[] A;
    delete[] B;
    delete[] C;

    return 0;
}

#!test_eigen.cpp
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{

    int n_a_rows = 2000;
    int n_a_cols = 2000;
    int n_b_rows = n_a_cols;
    int n_b_cols = 2000;

    MatrixXd a(n_a_rows, n_a_cols);

    for (int i = 0; i < n_a_rows; ++ i)
        for (int j = 0; j < n_a_cols; ++ j)
            a (i, j) = n_a_cols * i + j;

    MatrixXd b (n_b_rows, n_b_cols);
    for (int i = 0; i < n_b_rows; ++ i)
        for (int j = 0; j < n_b_cols; ++ j)
            b (i, j) = n_b_cols * i + j;

    MatrixXd d (n_a_rows, n_b_cols);

    clock_t begin = clock();

    d = a * b;
    clock_t end = clock();
    double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
    std::cout << "Time taken : " << elapsed_secs << std::endl;

}

#!test_numpy.py

import numpy as np
import time
N = 2000

a = np.random.rand(N, N)
b = np.random.rand(N, N)
start = time.time()
c = a.dot(b)
print(f"Time taken : {time.time() - start}")

Finally I created the following test file

#!test_matrix.sh
c++ -O2 -march=native -std=c++11 -I /usr/include/eigen3 test_eigen.cpp -o eigen
c++ -O2 -march=native -std=c++11 test_blas.cpp  -o blas -lcblas

echo "testing BLAS"
./blas
echo "testing Eigen"
./eigen
echo "testing numpy"
python test_numpy.py

which yields the output

testing BLAS
Time taken : 1.63807
testing Eigen
Time taken : 0.795115
testing numpy
Time taken : 0.28397703170776367

Now my question is, how come numpy is the fastest of these tests? Am I missing something with regards to optimizations?

One thing could be that numpy uses threading to compute the matrix product. Adding the compiler flag -fopenmp however yields worse performance for eigen and BLAS.

I am using g++ version 9.0.3-1. Numpy is version 1.18.1 using python 3.8.2. Thanks in advance.