Thanks for this excellent explanation. I had questions, but they were answered by a closer reading of your explanation. I recently did some work where, interes

Thanks for this excellent explanation. I had questions, but they were answered by a closer reading of your explanation. 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 Hence in matrix multiplication over binary matrix you very much do want the l.h.s. to be row-major and the r.h.s. to be column major. And, in general, you very much want to read and write matrices in a contiguous manner. Because otherwise, depending on alignment, you are accessing 8, 32, or 64 bits of memory per single bit entry of matrix your are reading or writing. 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.