matrix multiplication speed depends on silly things

Question

I wrote a program for fast matrix multiplication. To use CPU cache to the maximal extent, it divides matrix to 10*10 cells, and multiplies each cell separately (increasing cell size to 20*20 or 50*50 does not significantly change time).

It turns out that the speed significantly depends on whether matrix size and cell size is know in advance or not.

The program is:

#include <cmath>
#include <cstdlib>
#include <iostream>

using namespace std;

#define forall(i,n) for(int i=0; i<(int)(n); i++)

inline void Load(int N, int N2, float* x2, float* x, int iStart, int jStart) {
    int start = iStart * N + jStart;
    forall (i, N2)
        forall (j, N2)
            x2[i * N2 + j] = x[start + i * N + j];
}

inline void Add(int N, int N2, float* x, float* x2, int iStart, int jStart) {
    int start = iStart * N + jStart;
    forall (i, N2)
        forall (j, N2)
            x[start + i * N + j] += x2[i * N2 + j];
}

inline void Mul(int N, float* z, float* x, float* y) {
    forall (i, N)
        forall (j, N) {
            double sum = 0;
            forall (k, N)
                sum += x[i*N+k] * y[k*N+j];
            z[i*N+j] = sum;
        }
}

inline double RandReal() {return random()/((double)RAND_MAX+1);}

int main(int argc, char** argv) {
#if VAR==3
    const int N = atoi(argv[1]), N2 = atoi(argv[2]);
#elif VAR==2
    const int N = 2000, N2 = atoi(argv[2]);
#elif VAR==1
    const int N = atoi(argv[1]), N2 = 10;
#elif VAR==0
    const int N = 2000, N2 = 10;
#else
    #error "Bad VAR"
#endif
    cout << "VAR=" << VAR << " N=" << N << " N2=" << N2 << endl;
    float x[N*N], y[N*N];
    forall (i, N)
        forall (j, N) {
            x[i*N+j] = RandReal();
            y[i*N+j] = RandReal();
        }
    float z[N*N];
    forall (i, N)
        forall (j, N)
            z[i*N+j] = 0;
    for (int i1 = 0; i1 < N; i1 += N2) {
        float x2[N2*N2], y2[N2*N2], z2[N2*N2];
        for (int j1 = 0; j1 < N; j1 += N2) {
            Load(N, N2, x2, x, i1, j1);
            for (int k1 = 0; k1 < N; k1 += N2) {
                Load(N, N2, y2, y, j1, k1);
                Mul(N2, z2, x2, y2);
                Add(N, N2, z, z2, i1, k1);
            }
        }
    }

    double val = 0, val2 = 0;
    forall (i, N)
        forall (j, N)
            val += z[i*N+j], val2 += z[i*N+j]*(i+j);
    cout << "val=" << val << " val2=" << val2 << endl;
}

Now the execution times:

$ for a in 0 1 2 3; do g++ -DVAR=$a -O3 -Wall -o mat mat.cpp; time ./mat 2000 10; done
VAR=0 N=2000 N2=10
val=2.00039e+09 val2=3.99867e+12

real    0m8.127s
user    0m8.108s
sys     0m0.020s
VAR=1 N=2000 N2=10
val=2.00039e+09 val2=3.99867e+12

real    0m3.304s
user    0m3.292s
sys     0m0.012s
VAR=2 N=2000 N2=10
val=2.00039e+09 val2=3.99867e+12

real    0m25.395s
user    0m25.388s
sys     0m0.008s
VAR=3 N=2000 N2=10
val=2.00039e+09 val2=3.99867e+12

real    0m25.515s
user    0m25.495s
sys     0m0.016s

In simple terms:

• Don't know matrix size, know cell size: 3.3 seconds
• Know both matrix size and cell size: 8.1 seconds
• Don't know cell size: 25.5 seconds

Why is it? I use g++ 5.4.0 .

inline does not play any role, if we remove it, the results will be the same.

Answer 1

Introductory note: Much of this post has been re-written, so some of the comments below won't make too much sense anymore. Please look behind the edit for details, if you care. So...

tl;dr

• profile to find hot loops
• use constant counts for these if possible
• try manually unrolling them if not
• OP's results a mystery, nobody else can reproduce anything so extreme

I agree with @user4581301 - the more the compiler knows ahead of time, the more it can do for you in terms of optimising your code.

But you need to profile this code - wall clock times will only take you so far. I don't know a whole lot about profilers for gcc (I have a good one for MSVC) but you could try your luck here.

It also pays (as @RetiredNinja said right off the bat) to try to learn some assembler, using Godbolt as a tool, especially if you want to understand a slowdown as dramatic as this.

Now having said all that, your timings make no sense to me so something strange is going on on your system. So I ran some tests of my own, and my results differ markedly from yours. I ran some of these tests on MSVC (because I have such wonderful profiling tools there) and some on gcc on the Mac (although I think that is actually secretly clang under the hood). I have no linux box, sorry.

Firstly, let's deal with the issue of allocating such large objects on the stack. This is obviously not wise, and I can't do it at all on MSVC since it doesn't support variable length arrays, but my tests on the Mac showed that this made no difference to the timings once I had increased ulimit to get it to work at all (see here). So I think we can discount that as a factor, as you yourself say in the comments.

So now let's look at the timings I got on the Mac:

VAR=0 USE_STACK=0 N=2000 (known) N2=10 (known)
user    0m10.813s

VAR=1 USE_STACK=0 N=2000 (unknown) N2=10 (known)
user    0m11.008s

VAR=2 USE_STACK=0 N=2000 (known) N2=10 (unknown)
user    0m12.714s

VAR=3 USE_STACK=0 N=2000 (unknown) N2=10 (unknown)
user    0m12.778s

VAR=0 USE_STACK=1 N=2000 (known) N2=10 (known)
user    0m10.617s

VAR=1 USE_STACK=1 N=2000 (unknown) N2=10 (known)
user    0m10.987s

VAR=2 USE_STACK=1 N=2000 (known) N2=10 (unknown)
user    0m12.653s

VAR=3 USE_STACK=1 N=2000 (unknown) N2=10 (unknown)
user    0m12.673s

Not much to see there; let's move on to what I observed on MSVC (where I can profile):

VAR=0 USE_STACK=0 N=2000 (known) N2=10 (known)
Elapsed: 0:00:06.89

VAR=1 USE_STACK=0 N=2000 (unknown) N2=10 (known)
Elapsed: 0:00:06.86

VAR=2 USE_STACK=0 N=2000 (known) N2=10 (unknown)
Elapsed: 0:00:10.24

VAR=3 USE_STACK=0 N=2000 (unknown) N2=10 (unknown)
Elapsed: 0:00:10.39

Now we have something we can get our teeth into. Like @geza observed, the code takes longer to run when N2 isn't known, which is entirely in line with what one might expect since this is where the hot loops will be, and it's much more likely that the compiler will unroll such a loop when it knows that it's actually comprised of a small, known number of iterations.

So let's get some information from the profiler. It tells me that the hot loop (by quite a long way) is the inner loop in Mul():

inline void Mul(int N, float* z, float* x, float* y) {
    forall (i, N)
        forall (j, N) {
            double sum = 0;

=> forall (k, N) => sum += x[i*N+k] * y[kN+j]; z[iN+j] = sum; } }

Again, I can't say this surprises me much, and when I look at the code I can see that the loop is not unrolled at all (loop setup code omitted for brevity):

$1:
    movss       xmm0,dword ptr [r9+rsi*4]  
    mulss       xmm0,dword ptr [r8+4]  
    movss       xmm1,dword ptr [r9+r15*4]  
    mulss       xmm1,dword ptr [r8]  
    cvtps2pd    xmm2,xmm0  
    cvtps2pd    xmm0,xmm1  
    movss       xmm1,dword ptr [r8+8]  
    mulss       xmm1,dword ptr [r9]  
    addsd       xmm0,xmm3  
    addsd       xmm2,xmm0  
    cvtps2pd    xmm0,xmm1  
    movss       xmm1,dword ptr [r9+r14*4]  
    movaps      xmm3,xmm2  
    mulss       xmm1,dword ptr [r8+0Ch]  
    add         r9,rbp  
    add         r8,10h  
    addsd       xmm3,xmm0  
    cvtps2pd    xmm0,xmm1  
    addsd       xmm3,xmm0  
    sub         rcx,1  
    jne         $1

Now it doesn't look that there would be any savings to be made there by unrolling that loop since looping around is going to be cheap compared to executing all the rest of the code in there, but if you look at the disassembly of the same loop when N2 is known, you get a surprise:

    movss       xmm0,dword ptr [rax-8]  
    mulss       xmm0,dword ptr [rcx-50h]  
    cvtps2pd    xmm2,xmm0  
    movss       xmm0,dword ptr [rcx-28h]  
    mulss       xmm0,dword ptr [rax-4]  
    addsd       xmm2,xmm7  
    cvtps2pd    xmm1,xmm0  
    movss       xmm0,dword ptr [rcx]  
    mulss       xmm0,dword ptr [rax]  
    addsd       xmm2,xmm1  
    cvtps2pd    xmm1,xmm0  
    movss       xmm0,dword ptr [rcx+28h]  
    mulss       xmm0,dword ptr [rax+4]  
    addsd       xmm2,xmm1  
    cvtps2pd    xmm1,xmm0  
    movss       xmm0,dword ptr [rcx+50h]  
    mulss       xmm0,dword ptr [rax+8]  
    addsd       xmm2,xmm1  
    cvtps2pd    xmm1,xmm0  
    movss       xmm0,dword ptr [rcx+78h]  
    mulss       xmm0,dword ptr [rax+0Ch]  
    addsd       xmm2,xmm1  
    cvtps2pd    xmm1,xmm0  
    movss       xmm0,dword ptr [rcx+0A0h]  
    mulss       xmm0,dword ptr [rax+10h]  
    addsd       xmm2,xmm1  
    cvtps2pd    xmm1,xmm0  
    movss       xmm0,dword ptr [rcx+0C8h]  
    mulss       xmm0,dword ptr [rax+14h]  
    addsd       xmm2,xmm1  
    cvtps2pd    xmm1,xmm0  
    movss       xmm0,dword ptr [rcx+0F0h]  
    mulss       xmm0,dword ptr [rax+18h]  
    addsd       xmm2,xmm1  
    cvtps2pd    xmm1,xmm0  
    movss       xmm0,dword ptr [rcx+118h]  
    mulss       xmm0,dword ptr [rax+1Ch]  
    addsd       xmm2,xmm1  
    cvtps2pd    xmm1,xmm0  
    addsd       xmm2,xmm1  
    cvtpd2ps    xmm0,xmm2  
    movss       dword ptr [rdx+rcx],xmm0  

Now there is no loop, and the number of instructions that will be executed overall is clearly reduced. Perhaps MS are not such a bunch of dumb clucks after all.

So finally, as an exercise, let's just quickly unroll that loop manually and see what timings we get (I didn't look at the generated code):

#define UNROLL_LOOP 1

inline void Mul(int N, float* z, float* x, float* y) {
    forall (i, N)
        forall (j, N) {
            double sum = 0;
#if UNROLL_LOOP
            assert (N == 10);
            sum += x[i*N] * y[0*N+j];
            sum += x[i*N+1] * y[1*N+j];
            sum += x[i*N+2] * y[2*N+j];
            sum += x[i*N+3] * y[3*N+j];
            sum += x[i*N+4] * y[4*N+j];
            sum += x[i*N+5] * y[5*N+j];
            sum += x[i*N+6] * y[6*N+j];
            sum += x[i*N+7] * y[7*N+j];
            sum += x[i*N+8] * y[8*N+j];
            sum += x[i*N+9] * y[9*N+j];
#else
            forall (k, N)
                sum += x[i*N+k] * y[k*N+j];
#endif
            z[i*N+j] = sum;
        }
}

And when I did that, I got:

VAR=3 USE_STACK=0 N=2000 (unknown) N2=10 (unknown)
Elapsed: 0:00:07.48 (compared with 10.39 / 6.86, not bad, more may be possible).

So that's the process you need to go through to analyse performance issues like this, and you need good tools. I don't know what's happening in your case, because I can't reproduce the issue, but loop unrolling is (as expected) the primary factor on MSVC when (small) loop counts are unknown.

The test code I used is here in case anyone wants to refer to it. I think you owe me a vote, OP.

Edit:

Played around a little bit at Wandbox with gcc 9.0.0. Timings (these are rather slower and a bit more imprecise since we are running on a shared box, or, more likely, in a VM):

VAR=0 USE_STACK=0 N=2000 (known) N2=10 (known) Elapsed time = ~8sec

VAR=3 USE_STACK=0 N=2000 (unknown) N2=10 (unknown) Elapsed time = ~15.5sec

VAR=3 USE_STACK=0 N=2000 (unknown) N2=10 (unknown), loop unrolled Elapsed time = ~ 13.5sec

So that needs a bit more investigation - with a profiler and by looking at the generated code - and is still a million miles away from what the OP is getting.

Live demo - you can play around with that yourself if you want to try different things, OP.