tl;dr Should I store a matrix of Nx3 elements row-wise or column-wise, if I'm going to access 3 elements at a time? Does it matter?I'm setting up a numerical simulation with matrices of x,y,z positions in 3 dimensions, and I wonder whether it's better to store them row-wise or column-wise in R, Matlab, and Armadillo (C++). Maybe it doesn't matter at all.The simulation considers N particles, say N=5000, and each is described by 3 coordinates, 3 angles, and 3 sizes. So I've 3 matrices of Nx3 elements. In the code (setting up a linear system of 3N equations), I will need to access positions, sizes, and angles one particle at a time. Therefore, I'm tempted to store those as columns, i.e. those matrices would have 3 rows, and N columns. This is because R, Matlab and Armadillo all have a column-major ordering, so my thinking is that accessing a column will be (ever so slightly) faster. On the other hand, maybe there is a penalty of having many columns (memory addresses)? I do not know.By-the-way, I'm mentioning three programming languages because I'll be writing the code in all three (it's not much effort, since the syntax is pretty much the same).
matrix storage, many rows or many columns?
matlab;r
Which storage scheme is optimal depends very strongly on the accesspattern. You already said that you're accessing them one particle at atime, but are you going to access the particles sequentially in order,or randomly? Another conventional name for this issue isstruct-of-arrays vs array-of-structs.The key thing here is to know how the memory hierarchy works on a CPU(see, e.g.,https://www.scss.tcd.ie/Jeremy.Jones/CS3021/5%20caches.pdf, as well as the many questions on this topic on Stack Overflow). Readsfrom main memory have a very large latency, which will lead to longstalls in the execution, unless this is mitigated by prefetching. Eachread from main memory will read one cache line worth of data, placingit in cache, and the reads from cache are comparatively free.So for example if you have a matrix $(a_{ij})_{1\leq i,j\leq 3}$,and it is stored in column-major order in memory as $$ a_{11}, a_{21}, a_{31}, a_{12}, a_{22}, a_{32}, a_{13}, a_{23},a_{33}, $$then reading $a_{11}$ (or any of $a_{11},\ldots,a_{23}$) will place all the eight elements from $a_{11}$ to$a_{23}$ in cache (assuming correct alignment and a typical 64-bytecache line and an element size of 8 bytes), making them free to access.If you have several different arrays, a struct-of-arrays, like$$ x_{1}, \ldots, \qquad y_{1}, \ldots, \qquad z_{1}, \ldots $$and you read $x_1, y_1, z_1$, you will effectively read the firsteight elements of each array. If you end up using the next sevenelements immediately, that's fine. If not, then that wasted 7/8 ofyour memory bandwidth by reading 24 elements instead of 3, which is inefficient.If instead you store them as array-of-structs,$$ x_1, y_1, z_1, x_2, y_2, z_2, \ldots $$then reading $x_1,y_1,z_1$ will read the first eight elements. If you use those nextelements, there is no great difference, otherwise you waste 5/8 (reading 8 instead of 3) of your memory bandwidth - note this is still better than reading 24 elements as in the previous example.With nine 8-byte elements, if you know you'll access themconsecutively there should be no big difference either way, as above. If you willaccess them in random order, then reading one set of nine elementswill lead to either $9\times 8 = 72$ elements being read from memory (ifstored as struct-of-arrays) or $\lceil \frac{9}{8}\rceil\times 8 = 16$elements (if stored as array-of-structs) - this could be avery large performance gain.The difference between column-major and row-major, I feel, is a bit ofred herring here - the key is knowing the access pattern, and makingsure the data layout matches the access pattern. The next mostimportant thing would be prefetching, which is a technique forhiding the relatively large memory access latency.A separate related issue to consider is how you access data that you already read previously, and whether it is likely to still remain in the cache or whether it was purged to make way for other data, causing cache misses - an array of $5000\times 9\times 8$ bytes is $360\,\mathrm{KB}$ and will likely fit entirely in either the L2 or L3 cache. It gets even more complicated if your program involves multiple threads running on different cores, with potentially separate lower-level caches per core.Note, incidentally, that for non-sequential access patterns, having nine elements per particle is somewhat inefficient in itself if a cache line is 8 elements long as reading nine elements would be essentially equivalent to reading 16 elements - one technique is to split your data into often-used versus rarely-used parts so that the often-used parts fit better.There is no penalty for having either many rows or columns. Not sure what you mean by a penalty for many memory addresses, but there is no such penalty.Some of the things I said depend on exact properties of the architecture your program will run on. Using an interpreted language like R or Matlab can also affect this severely - for C++ it should hold.In general, it is best to use the array-of-structs pattern, unless you have some specific reasons not to.