How to achieve ultimate summing speed?
Jean-Francois
I'm stuck with the following problem:
I have a large matrix R, say of size 2000 by 500, and a threshold value k. For each row, I seek the numbers of elements that are below k. Hence, the output would be a 2000 by 1 vector.
So far, I've used "Q = SUM( R < k , 2 )", which proved to be much faster than everything I could think of, e.g. HISTC or (obviously) SORT+DIFF. A single calculation is in the order of milliseconds, so it wouldn't seem like a big deal. However, in the context of a Monte Carlo simulation, that summation turns out to be a huge bottleneck (> 60%); hence, even the slightest increase in efficiency would be very valuable.
Thanks for your suggestion.
us
before we get started: two important questions...
- does this run in a loop
- does your K change, eg, in the loop or according to some rule...
us
Walter Roberson
If you are looking for minor increases in efficiency, then consider
building your R matrix sideways. Summing along columns is slightly
faster than summing along rows, because Matlab data is stored down
columns rather than being stored across rows.
Rune Allnor
1) Make sure your vectors are stored sequentially in memory
2) MEX the routine
3) Make sure to optimize the compiler
Rune
*Very* rudimentary - no safeguards; not at all tested - C++ code:
#include "mex.h"
void mexFunction(
int nlhs,
mxArray *plhs[],
int nrhs,
const mxArray *prhs[]
)
{
const size_t M(mxGetM(prhs[0]));
const size_t N(mxGetN(prhs[0]));
const double* const p(mxGetPr(prhs[0]));
const double k((mxGetPr(prhs[1]))[0]);
plhs[0] = mxCreateNumericMatrix(N,1,mxUINT32_CLASS,mxREAL);
size_t* d((size_t*)mxGetData(plhs[0]));
for (size_t m = 0; m<M; ++m)
{
size_t n = 0;
double K = 0;
const double* const q(p+m*M);
while((n < N) && (K < k))
{
K+=q[n];
++n;
}
d[m]=--n;
}
}
James Tursa
The R < k calculation returns a 2000 x 500 matrix, which then must be accessed again to do the SUM calculation. This overhead could be eliminated with a simple C-mex routine.
// File sum2.c
// Implements sum(R<k,2), called as sum2(R,k)
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwSize i, j, m, n;
double *pr, *pr0, *R;
double k;
m = mxGetM(prhs[0]);
n = mxGetN(prhs[0]);
R = mxGetPr(prhs[0]);
k = mxGetScalar(prhs[1]);
plhs[0] = mxCreateDoubleMatrix(m, 1, mxREAL);
pr0 = mxGetPr(plhs[0]);
for( j=0; j<n; j++ ) {
pr = pr0;
for( i=0; i<m; i++ ) {
*pr++ += (*R++) < k;
}
}
}
My tests with the lcc compiler show about a 15% - 20% improvement in speed over the sum(R<k,2) method.
James Tursa
James Tursa
Here is a slightly faster version that does the summing as integers instead of floating point. My tests show a 20% - 25% speed improvement for this version.
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwSize i, j, m, n;
mwSize *s;
double *pr, *R;
double k;
m = mxGetM(prhs[0]);
n = mxGetN(prhs[0]);
R = mxGetPr(prhs[0]);
k = mxGetScalar(prhs[1]);
plhs[0] = mxCreateDoubleMatrix(m, 1, mxREAL);
pr = mxGetPr(plhs[0]);
s = mxCalloc(m, sizeof(*s));
for( j=0; j<n; j++ ) {
for( i=0; i<m; i++ ) {
if( (*R++) < k ) {
++s[i];
}
}
}
for( i=0; i<m; i++ ) {
pr[i] = s[i];
}
mxFree(s);
}
James Tursa
Rune Allnor
<aclassyguy_with_a_k_not_...@hotmail.com> wrote:
> s = mxCalloc(m, sizeof(*s));
> for( i=0; i<m; i++ ) {
> pr[i] = s[i];
> }
> mxFree(s);
Those calloc/free calls are expensive. What about
letting s be a local scalar and move the result to
pr[i] after treating each row?
I wouldn't be surprised if you end up saving >> 50%
by doing that.
Rune
James Tursa
I can try that. My thought was to traverse the R array only once, hence the s allocation. If I use only one scalar then I have to traverse the R array m times, which I expect to be slower than allocating s, but I haven't actually tried it yet. I will give it a shot ...
James Tursa
James Tursa
Here is the result, about 70% slower than my previous post using an allocated s. This is about what I would have expected given the multiple traverses of R involved. All that redundant memory access just kills the running times.
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwSize i, j, m, n;
mwSize s;
double *pr, *R, *R0;
double k;
m = mxGetM(prhs[0]);
n = mxGetN(prhs[0]);
R0 = mxGetPr(prhs[0]);
k = mxGetScalar(prhs[1]);
plhs[0] = mxCreateDoubleMatrix(m, 1, mxREAL);
pr = mxGetPr(plhs[0]);
for( i=0; i<m; i++ ) {
R = R0++;
s = 0;
for( j=0; j<n; j++ ) {
if( *R < k ) {
++s;
}
if( j < n-1 ) R += m;
}
pr[i] = s;
}
}
James Tursa
Rune Allnor
<aclassyguy_with_a_k_not_...@hotmail.com> wrote:
> "James Tursa" <aclassyguy_with_a_k_not_...@hotmail.com> wrote in message <hmrtu8$go...@fred.mathworks.com>...
> > Rune Allnor <all...@tele.ntnu.no> wrote in message <44ae8909-03ae-4219-a6b7-14d9329eb...@z35g2000yqd.googlegroups.com>...
> > > On 5 Mar, 20:13, "James Tursa"
> > > <aclassyguy_with_a_k_not_...@hotmail.com> wrote:
>
> > > > s = mxCalloc(m, sizeof(*s));
> > > > for( i=0; i<m; i++ ) {
> > > > pr[i] = s[i];
> > > > }
> > > > mxFree(s);
>
> > > Those calloc/free calls are expensive. What about
> > > letting s be a local scalar and move the result to
> > > pr[i] after treating each row?
>
> > > I wouldn't be surprised if you end up saving >> 50%
> > > by doing that.
>
> > > Rune
>
> > I can try that. My thought was to traverse the R array only once, hence the s allocation. If I use only one scalar then I have to traverse the R array m times, which I expect to be slower than allocating s, but I haven't actually tried it yet. I will give it a shot ...
>
> > James Tursa
>
> Here is the result, about 70% slower than my previous post using an allocated s. This is about what I would have expected given the multiple traverses of R involved. All that redundant memory access just kills the running times.
You traverse the rows? I didn't see that. I assumed
you traversed the columns. Traversing the columns one
could use the scalar local variable at the same time
there would be no need to traverse the array more than
once.
Rune
James Tursa
Yes. That issue was brought up earlier by another poster. If OP could rearrange his data as the transpose, then one could do what you suggest.
James Tursa
Jean-Francois
Yes, it does. In fact, the whole code is as follows:
***************************************
function d = trshld(R,bins,k,p)
[ndraw n] = size(R); Nb = length(bins)-1;
Bins_l = repmat(bins(1:end-1),[ndraw 1]);
Bins_h = repmat(bins(2:end),[ndraw 1]);
d = 0;
for j = 1:length(k)
Q = sum((R < k(j)),2)/n;
Q = repmat(Q,[1 Nb]);
Ql = (Q > Bins_l);
Qh = ~(Q > Bins_h);
d = d + p(j)*sum(Ql.*Qh,1);
end
d = d/ndraw;
***************************************
Basically, the code receives a vector 'k'. For each element of 'k', it first calculates for each row of 'R' the fraction of elements that fall below k. Then it assigns the resulting fractions to pre-specified bins. 'p' is just a weight vector and is innocuous. There is a huge bottleneck at 'Q = sum((R < k(j)),2)/n', and there is also one, somehow smaller, at 'Q = repmat(Q,[1 Nb])'. Should the whole code go to a mex file, or parts of it? Thanks for your answer.
James Tursa
>
> Should the whole code go to a mex file, or parts of it? Thanks for your answer.
Have you tried the provided mex code yet?
James Tursa
Jean-Francois
With R a 2000 x 200 matrix and 14,000 runs, I got the following times:
1. sum(R<k,2) --> 17.515 s
2. sum(R<k,1) --> 12.545 s
3. mex code from post 6 --> 15.810 s
4. mex code from post 9 --> 20.341 s
I am not familiar with C, so I am still figuring out how to traverse columns as was suggested previously. Also, I'm not sure I understood well, but does traversing rows vs columns make a difference, as it does in Matlab? Thanks all for the help.
Walter Roberson
> I am not familiar with C, so I am still figuring out how to traverse
> columns as was suggested previously. Also, I'm not sure I understood
> well, but does traversing rows vs columns make a difference, as it does
> in Matlab? Thanks all for the help.
Yes, in C (and in Matlab) it can make an enormous difference in the timings!!
But whether it is a major difference or a minor difference depends on the
quality of the compiler and the details of the processor instruction set and
upon the exact size of the matrix as compared to the exact size of the
processor's primary cache.
In general, processors are designed under the assumption that values that are
close together in the address space will be accessed fairly close together in
time in the program, so they are generally designed so that when a value is
read from (relatively slow) main memory to the processor, that a group of
adjacent values are also read in at the same time and are held in very fast
memory in or "very near" the CPU. When the CPU then goes to access that nearby
value it is then read from that very fast memory instead of going through the
slow memory interface to the main memory subsystem.
But there is a limit to how much fast memory is available, so when you go to
read in a new value, the fast-memory subsystem ("primary cache") needs to
discard (or write out to main memory if needed) something to make room. There
are different techniques to determine what gets discarded. One of the commonly
used techniques is to divide main memory into columns (you could describe it
as; "cache lines" is the technical term), and when a new item is to be read
in, what is discarded or written out is the old item that happened to be from
the same column of memory. Thus if you are handling an item at address R, and
the processor then needs an item from address R + k*S where k is an integer
and S is a fixed size dependent on the system configuration, then the item at
R is discarded and replaced by the item at R + k*S . Which is all well and
fine unless your array size happens to be a multiple of S. If A(i,j) happens
to be a multiple of S away from A(i+1,j) then you want to avoid accesses that
go along the first dimension, and if A(i,j) happens to be a multiple of S away
from A(i,j+1) then you want to avoid access that go along the second dimension.
In Matlab, A(i+1,j) is always immediately adjacent to A(i,j); in C, A(i,j+1)
is always immediately adjacent to A(i,j). So in C, you want to iterate along
the columns, and in Matlab you want to iterate along the rows -- because if
you mix up the orders, there is a chance that the hardware architecture will
make it necessary to throw away A(i,j) from fast memory in order to access the
new value.
James Tursa
Traversing rows vs columns may or may not make a big difference ... it depends on how the algorithm is coded and potentially on machine architecture. For the C mex routines, there will likely be very little difference between the timing of the post 6 sum2.c code which sums along the rows and the sum1.c code posted below which sums along the columns. That is because they are both coded in such a way that they do not create a large intermediate array (like MATLAB does for the R<k calculation) and they only traverse the R array once. The only important difference is the allocation of the s array and the multiple traversing of the s array in the sum2.c code, whereas in the sum1.c code a single scalar s is sufficient. That will in all likelihood only make a small timing increase in sum2.c compared to sum1.c. I would expect that they both would outperform the MATLAB code you have shown
because of the overhead savings ... the amount of the overhead savings dependent on computer, MATLAB version, etc. etc. If the MATLAB behind the scenes algorithm traverses the R array multiple times you could easily see a doubling of the run times in comparison to the C code, but if MATLAB uses good algorithms you may only see a small (e.g., 10% - 15%) improvement. I don't think you will be able to improve this calculation speed any more than these C routines give you. Of course, a better C compiler can sometimes make use of registers etc. more efficiently, particularly for looping code such as this. I have seen cases where Microsoft Visual C/C++ outperformed the supplied lcc compiler by a factor of 2-3 in speed. If you do any C work like this at all, IMO it is worth it to invest in a good C compiler.
James Tursa
// sum1.c, calculates the equivalent of sum(R<k,1)
// called as sum1(R,k)
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwSize i, j, m, n;
register mwSize s;
double *pr, *R;
double k;
m = mxGetM(prhs[0]);
n = mxGetN(prhs[0]);
R = mxGetPr(prhs[0]);
k = mxGetScalar(prhs[1]);
plhs[0] = mxCreateDoubleMatrix(1, n, mxREAL);
pr = mxGetPr(plhs[0]);
for( j=0; j<n; j++ ) {
s = 0;
for( i=0; i<m; i++ ) {
if( (*R++) < k ) {
++s;
}
}
*pr++ = s;
}
}
Jean-Francois
> // called as sum1(R,k)
> #include "mex.h"
> void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
> {
> mwSize i, j, m, n;
> register mwSize s;
> double *pr, *R;
> double k;
>
> m = mxGetM(prhs[0]);
> n = mxGetN(prhs[0]);
> R = mxGetPr(prhs[0]);
> k = mxGetScalar(prhs[1]);
> plhs[0] = mxCreateDoubleMatrix(1, n, mxREAL);
> pr = mxGetPr(plhs[0]);
> for( j=0; j<n; j++ ) {
> s = 0;
> for( i=0; i<m; i++ ) {
> if( (*R++) < k ) {
> ++s;
> }
> }
> *pr++ = s;
> }
> }
I tried using sum1.c, with no significantly improved performance over using Matlab's sum. Let me point out that in performing sum(R<k), ~66% of the time is spent on (R<k), and the rest on the summation itself, which confirms what you said above about the overhead. Is there any way to make the comparison between R and k more efficient, either with Matlab or a mex code?
Jan Simon
> > // sum1.c, calculates the equivalent of sum(R<k,1)
> > // called as sum1(R,k)
> I tried using sum1.c, with no significantly improved performance over using Matlab's sum. Let me point out that in performing sum(R<k), ~66% of the time is spent on (R<k), and the rest on the summation itself, ...
I'm confused. sum1.c includes the the R<k comparison already, so if 66% of the time of SUM(R<k,1) is spent for R<k, then I expect sum1.c to be faster. Can you show us your speed comparison with some details?
Which compiler do you use?
You can distribute the calculations for the different columns to different threads, if your computer has multiple cores. E.g.:
http://www.mathworks.com/matlabcentral/fileexchange/21233
shows how this can be programmed.
Kind regards, Jan
Oliver Woodford
> plhs[0] = mxCreateDoubleMatrix(1, n, mxREAL);
> pr = mxGetPr(plhs[0]);
> for( j=0; j<n; j++ ) {
> s = 0;
> for( i=0; i<m; i++ ) {
> if( (*R++) < k ) {
> ++s;
> }
> }
> *pr++ = s;
> }
James, I'm surprised a speed demon like you is using mxCreateDoubleMatrix like that. Since it sets the matrix to zero first it drags the whole matrix through the cache once, before you even write to it. Since you then set every entry later you don't need to do that. See a discussion on the subject, and solution, here:
http://wwwuser.gwdg.de/~mleuten/MATLABToolbox/CmexWrapper.html
Oliver
Oliver Woodford
> You can distribute the calculations for the different columns to different threads, if your computer has multiple cores. E.g.:
> http://www.mathworks.com/matlabcentral/fileexchange/21233
> shows how this can be programmed.
For this simple iteration over each column, with no dependencies between each one, the easiest way (IMHO) to parallelize this is with an OpenMP pragma before the first (outer) for loop. And it's completely platform independent.
Oliver
James Tursa
Yes, I am aware of that technique but hadn't thought of it for this application. Thanks for pointing it out. I don't know how much difference it will make, however, but I will try to run some tests later.
I have made previous attempts at speed improvement tests like this using mxMalloc (I was attempting to create a fast C mex preallocation routine that didn't zero out the memory) but have noticed that mxMalloc (or a previous mxFree) *always* zeroes out the memory. I made this conclusion after allocating very large blocks, setting all the elements to non-zero, freeing the memory with mxFree, and then immediately using mxMalloc again to grab the exact same block of memory and noting that all the memory is suddenly zeroed out. Don't know if it is mxMalloc or mxFree that is doing it. Bottom line is using mxCreateDoubleMatrix with 1 x n vs creating 0 x 0 and then mxMalloc for the memory may *still* drag the memory through the cache once to get the memory set to zero and you may not have any control over this. This may be a security design in the MATLAB memory manager, or something else. I will
quite frankly admit that I don't fully understand this behavior yet & who is doing the zeroing, and don't know all the implications with regards to cache memory, etc.
James Tursa
Jan Simon
You need a compiler with OpenMP support. E.g. the Express version of MSVC allows this after some cumbersome extra downloads only (my trials to follow the instructions from the net have not been successful).
Kind regards, Jan
Oliver Woodford
Shame. I was able to get the Express version of MSVC working admirably with OpenMP. If I recall correctly, there were two issues: Firstly downloading the extra libraries, which are in the Windows SDK. This went straighfowardly for me. I installed MSVC first, then the SDK. Secondly you need to set the OpenMP flag in the mexopts.bat file. Add the extra line:
set COMPFLAGS=%COMPFLAGS% /openmp
On linux, gcc has support for OpenMP built in, but you again need to add the flag in the mexopts.sh file. The extra line in this case is:
CFLAGS="$CFLAGS -fopenmp"
HTH,
Oliver
Jan Simon
> > plhs[0] = mxCreateDoubleMatrix(1, n, mxREAL);
> James, I'm surprised a speed demon like you is using mxCreateDoubleMatrix like that. Since it sets the matrix to zero first it drags the whole matrix through the cache once, before you even write to it. Since you then set every entry later you don't need to do that. See a discussion on the subject, and solution, here:
> http://wwwuser.gwdg.de/~mleuten/MATLABToolbox/CmexWrapper.html
Although your arguments are convincing and I've expected a remarkable acceleration also, I cannot reproduce advantages for your method for Matlab 6.5, 7.7 and 7.8 using different compilers, e.g. OpenWatcom 1.8 and MSVC 2008 Express.
Is it necessary to start Matlab with a specific memory manager?
Kind regards, Jan
Jan Simon
> > You need a compiler with OpenMP support. E.g. the Express version of MSVC allows this after some cumbersome extra downloads only (my trials to follow the instructions from the net have not been successful).
> >
> > Kind regards, Jan
>
> Shame. I was able to get the Express version of MSVC working admirably with OpenMP. If I recall correctly, there were two issues: Firstly downloading the extra libraries, which are in the Windows SDK. This went straighfowardly for me. I installed MSVC first, then the SDK. Secondly you need to set the OpenMP flag in the mexopts.bat file. Add the extra line:
> set COMPFLAGS=%COMPFLAGS% /openmp
Shame - on the compiler or the user?
I'll try it again. I assume you are using MSVC 2008?
Do you have experiences with the GCC on Windows+MinGW or Cygwin and Mex compilations?
Kind regards, Jan
Oliver Woodford
> Shame - on the compiler or the user?
> I'll try it again. I assume you are using MSVC 2008?
Yes (Express edition).
> Do you have experiences with the GCC on Windows+MinGW or Cygwin and Mex compilations?
No.
Oliver
Oliver Woodford
> Although your arguments are convincing and I've expected a remarkable acceleration also, I cannot reproduce advantages for your method for Matlab 6.5, 7.7 and 7.8 using different compilers, e.g. OpenWatcom 1.8 and MSVC 2008 Express.
> Is it necessary to start Matlab with a specific memory manager?
I wouldn't call them my argument or approach, rather Marcel Leutenegger's. I must confess I haven't tested the theory too much - I always just took it as read. Perhaps more fool me, given your and James' experience. Which does beg the question why MATLAB insists on setting the buffer to zeros even when you don't need it to - essentially (based on what James has said) there seems to be little practical difference between mxMalloc and mxCalloc.
Oliver
Jean-Francois
*************************
I'm using MSVC C++ Express. Which compiler do you recommend?
By the way, aren't Matlab basic functions (SUM, TIMES, etc) already using a system's multithreading capabilities?
So, I split SUM(R<k,1) into Q=(R<k) and SUM(Q,1) to see if the bottleneck arises from the comparison or from the summation. For R a 2000x200 matrix, I ran 122000 repetitions, and total times were:
- Q=(R<k) --> 73.259s
- SUM(Q,1) --> 37.812s
Hopefully, I answered your question. Do you have any suggestion at this stage? Thanks.
James Tursa
I ran some spot checks using the raw mxMalloc method and I get no discernible timing difference compared to the original method. Used various R2006b - R2009b with lcc and MSVC8. So, at least based on these timings I don't see any advantage. It *was* interesting to see how the sum calculation has changed from R2006b ro R2009b, however. Using R=rand(3000) and k=0.5 the sum(R<k,1) calculation took about 0.20 seconds in R2006b. As the versions of MATLAB increase, the timing of this calculation drops significantly. Partly because of moving to multi-threading, and partly because of more efficient coding/compiling. At R2009b the same calculation takes 0.04 seconds ... quite an improvement over the R2006b runtime of 0.20 seconds. The equivalent calculation in a C mex routine using MSVC8 was consistently 0.04 seconds for all MATLAB versions I tested.
The other thing I did was re-test the mxMalloc zeroing stuff. Re-confirmed that something was zeroing out the memory between mxFree and subsequent mxMalloc calls. Then I did the same test with malloc and free and got the same results ... something was zeroing out the memory. Then I ran similar test code for malloc and free on an Alpha mainframe and a Sun workstation and those machines did *not* zero out the memory between calls. So maybe this is a Windows security thing and not related to MATLAB or C memory managers at all.
James Tursa
Jean-Francois
******************
So now I use a parfor loop instead of cellfun to call the function TRSHLD from post 12, with a 50% increase in speed. On top of that, using sum1.c instead of Matlab's sum in TRSHLD gives an additional 33% gain in speed, for an overall gain of 66%. On the other hand, when I used cellfun, there was no noticeable difference between using sum1.c or Matlab's sum. Is there a scientific explanation?
Also, what is meant in Rune's post (#4) by 'storing vectors sequentially' and 'optimizing the compiler'? I use MSVC++, I guess with defaults settings. What should I do?
Jan Simon
> Also, what is meant in Rune's post (#4) by 'storing vectors sequentially' and 'optimizing the compiler'? I use MSVC++, I guess with defaults settings. What should I do?
Storing vectors sequentially has the advantage, that the function can access neighbouring elements. Neighbouring elements can be accessed faster, because it is "easier" to get them from the RAM. E.g. the summation over rows of a matrix is slower than over columns:
x = rand(1000);
tic; for i=1:100; v = sum(x, 1); end; toc % fast
tic; for i=1:100; v = sum(x, 2); end; toc % slow
Transposing the input is usually not recommended, because the TRANSPOSE itself suffers from accessing elements, which are not neighbouring:
tic; for i=1:100; v = sum(transpose(x), 1); end; toc % slow also
So it is recommended to create the arrays such, that the orientation allows a sequential access (here the example is ridiculous, but it is just a demonstration):
x = transposes(x); tic; for i=1:100; v = sum(x, 1); end; toc % fast
You can try to set the /arch:SSE2 flag in the mexopts.bat file in the folder PREFDIR. This sometimes help.
Kind regards, Jan
Jean-Francois
**********************
Hi Jan,
1. How do I set the /arch:SSE2 flag in the mexopts.bat file? There is no reference to /arch.
2. Also, is using --> mex -O "filename" <-- supposed to significantly optimize the code? I tried it, but with no noticeable difference.
3. Are there exceptions to programming column wise in C that you can think of? I have just begun low-level programming, so rules of thumb would make the game relatively easier at this stage.
Thanks for your help.
Jan Simon
> 1. How do I set the /arch:SSE2 flag in the mexopts.bat file? There is no reference to /arch.
Create the mexopts.bat file with "mex -setup"
cd(prefdir)
edit mexopts.bat
==> insert "/arch:SSE2" in the line:
set OPTIMFLAGS= ...
You can try this also: "/fp:fast"
In both cases you have to check the (accuracy of) the results and the speed can be in- or decreased!
> 2. Also, is using --> mex -O "filename" <-- supposed to significantly optimize the code? I tried it, but with no noticeable difference.
This can happen.
> 3. Are there exceptions to programming column wise in C that you can think of? I have just begun low-level programming, so rules of thumb would make the game relatively easier at this stage.
Accessing neighbouring elements is faster in every programming language, because this allows an faster caching by the processor. The really important and general exception is, that columnwise storing of values can be confusing and the time to debug the "optimal" algorithm can be 1000 times longer than the total processing time. Never optimize a part of a program, which is not frequently used - just concentrate on the bottlenecks.
Kind regards, Jan