Optimize matrix multiplication cache-friendliness

[frkns]

In the matrix multiplication code, noticed that the last dimension of looping, k, does not match the last dimension of memory access, j:

kactl/content/data-structures/Matrix.h

Lines 18 to 23 in ba85dcd

	M operator*(const M& m) const {
	M a;
	rep(i,0,N) rep(j,0,N)
	rep(k,0,N) a.d[i][j] += d[i][k]*m.d[k][j];
	return a;
	}

I believe this can be made more cache-friendly simply by swapping order of summation so that memory accesses are linear.
For example, here is a one-change that speeds up matrix multiplication by 50% locally:

 	rep(i,0,N) rep(j,0,N) 
 		rep(k,0,N) a.d[i][k] += d[i][j] * m.d[j][k];

[simonlindholm]

Makes sense. Most of the time you're probably going to use this code with modular arithmetic whose runtime is going to dominate that of cache effects (where optimization might involve SIMD Montmomery multiplication, which would also be nice to have in KACTL, and/or some lazy scheme of only modding once in a while), but for the float case, sure!

[bjorn-martinsson]

Makes sense. Most of the time you're probably going to use this code with modular arithmetic whose runtime is going to dominate that of cache effects (where optimization might involve SIMD Montmomery multiplication, which would also be nice to have in KACTL, and/or some lazy scheme of only modding once in a while), but for the float case, sure!

The typical trick to do matrix multiplication with modular arithmetics (mod 10^9 + 7) is to only apply the mod after summing like ~32 terms (the exact number depends on stuff like if you use long long vs unsigned long long). That way you get the mod part practically for free. This is very different to something like NTT where you have to do the mod after each multiplication.