How to write a matrix matrix product that can compete with Eigen? : <p>Below is the C++ implementation comparing the time taken by Eigen and For Loop to perform matrix-matrix products. The For loop has been optimised to minimise cache misses. The for loop is faster than Eigen initially but then eventually becomes slower (upto a factor of 2 for 500 by 500 matrices). What else should I do to compete with Eigen? Is blocking the reason for the better Eigen performance? If so, how should I go about adding  blocking to the for loop?</p>

<pre><code>#include&lt;iostream&gt;
#include&lt;Eigen/Dense&gt;
#include&lt;ctime&gt;

int main(int argc, char* argv[]) {
    srand(time(NULL));
    //  Input the size of the matrix from the user
    int N   =   atoi(argv[1]);

    int M   =   N*N;

    //  The matrices stored as row-wise vectors
    double a[M];
    double b[M];
    double c[M];

    //  Initializing Eigen Matrices
    Eigen::MatrixXd a_E =   Eigen::MatrixXd::Random(N,N);
    Eigen::MatrixXd b_E =   Eigen::MatrixXd::Random(N,N);
    Eigen::MatrixXd c_E(N,N);

    double CPS  =   CLOCKS_PER_SEC;

    clock_t start, end;


    //  Matrix vector product by Eigen
    start   =   clock();
    c_E     =   a_E*b_E;
    end     =   clock();

    std::cout &lt;&lt; "\nTime taken by Eigen is: " &lt;&lt; (end-start)/CPS &lt;&lt; "\n";

    //  Initializing For loop vectors
    int count   =   0;
    for (int j=0; j&lt;N; ++j) {
        for (int k=0; k&lt;N; ++k) {
            a[count]    =   a_E(j,k);
            b[count]    =   b_E(j,k);
            ++count;
        }
    }

    //  Matrix vector product by For loop
    start   =   clock();
    count   =   0;
    int count1, count2;
    for (int j=0; j&lt;N; ++j) {
        count1  =   j*N;
        for (int k=0; k&lt;N; ++k) {
            c[count]    =   a[count1]*b[k];
            ++count;
        }
    }

    for (int j=0; j&lt;N; ++j) {
        count2  =   N;
        for (int l=1; l&lt;N; ++l) {
            count   =   j*N;
            count1  =   count+l;
            for (int k=0; k&lt;N; ++k) {
                c[count]+=a[count1]*b[count2];
                ++count;
                ++count2;
            }
        }
    }
    end     =   clock();

    std::cout &lt;&lt; "\nTime taken by for-loop is: " &lt;&lt; (end-start)/CPS &lt;&lt; "\n";
}
</code></pre>
