wait, so a row major stored matrix should be multiplied by extracting the columns first? Or what?
Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage?
Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage?
No, the other way around. In row-major storage, all the elements of the same row are contiguous, so it is faster to process the array row by row.
> Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage?
The explanations here are decent enough. All you need is the figure at the top, really: https://en.m.wikipedia.org/wiki/Row-_and_column-major_order .
On a computer, such operations should almost never be implemented with scalar products, because the speed of a scalar product is limited by the latency of the FMA operation, not by its throughput. It is possible to accelerate the computation of a scalar product by ordering the elementary operations in a tree, but that complicates the program and it still does not reach the maximum speed of the hardware.
Except for vector-vector operations (i.e. for level 2 BLAS or greater levels), it is possible to change the order of the loops so that the scalar products are replaced by AXPY operations (whose result vectors should be kept in registers, to not be limited by the memory transfer throughput), which are not limited by the latency of the FMA, like the scalar products.
Except for vector-vector operations and matrix-vector operations (i.e. for level 3 BLAS or greater levels), it is possible to change the order of the loops so that the scalar products are replaced by rank-one matrix updates.
Most of the computation of a rank-one matrix update consists of the tensor product of 2 vectors.
Therefore, a matrix-matrix multiplication can be computed either by scalar products of the rows of the 1st matrix with the columns of the 2nd matrix, or by the tensor products of the columns of the 1st matrix with the rows of the 2nd matrix. The second method is much faster. The result of the tensor product must be kept in registers for the duration of the computation, so large matrices are partitioned in small blocks, i.e. submatrices, which are multiplied directly (an additional complication that increases the number of nested loops is that the small blocks that can be multiplied directly must be grouped in larger blocks that can be kept in the level-2 cache memory and reused for computations).
For matrix-matrix multiplications it is less important whether they are stored in column major order or in row major order, because when submatrices are brought into the L2 cache, the matrix elements will be gathered or scattered to be in the best order, before being loaded into registers repeatedly.
For matrix-vector multiplications, it is better to have the matrix in column major order, in order to compute the product by AXPY operations between columns and elements of the vector, and not by scalar products between rows and the vector (computing AXPY operations with shorter parts of a column at a time, so that the result vector can be kept in registers).
I recently did some work where, interestingly enough, this did not apply. I was doing linear algebra modulo 2 (i.e. binary entries, adding is xor, multiplication is and). In order to have any kind of memory (throughput) efficiency, you want to pack your matrix entries into bits of integers so a single 64 int can represent 64 entries.
This packing (when done row-wise or column-wise) gives a very quick primitive for computing inner products. <a, b> is just (in pseudo-code)
sum(popcnt(a[i] & b[i]) for i in range len(a)) % 2
I wonder if a similar very fast scalar product primitive for floats (or ints?) could work. If not, I imagine that is due to FMA having some inherent efficiencies that cannot be gained by a very quick accumulation step. One of the beauties in the binary case is the popcnt instruction doing a 64-input accumulation very quickly.