Hi all

I have started implementing some wrappers to numerical recipes code in matlab and there is a bug in the NRMatrix class template in the NR3matlab.h header file where it is defined. The method ncols returns the number of rows and nrows returns the number of columns, instead of what they should return. Internally nn is used to store the number of columns and mm the number of rows and it uses the worng one. There may be other problems elsewhere with row column confusion but I have not found them yet.

I am editing the template in my version to fix this so I can move forward with the coding but thought I should bring it to your attention, so that it can be fixed for future users to download.

Will Handler

... there is a bug...Internally nn is used to store the number of columns and mm the number of rows I don't see it that way. The value of nn is the number of rows, and the value of mm is the number of columns. See Footnote.

If you are thinking that the functions that create NR3 MatDoub objects from Matlab arguments have the wrong number of rows and columns, then I will say this:

It's not a "bug"; it's a "feature," as explained in the Tutorial on using NR3 code from within MATLAB (http://www.nr.com/nr3_matlab.html)


"Matrices: "Through the Looking Glass" (TTLG)
.
.
.
when a Matlab matrix crosses the interface into C++, what appears on the other side as its transpose."
(I think that's supposed to be "other side is its transpose.")


Look carefully at the results from NRmatrixDemo.cpp in the tutorial.

If you want to use unmodified NR stuff, then transpose all input and output matrices:


/*
* xgauss.cpp
*
* Demo of using gaussj from GNU Octave (Should be OK
* with Matlab, I think, but I have no way of testing.)
* See Footnote.
*
* Note that a "real" application would check the
* nature and validity of inputs. This is just a
* simple demo.
*
* davekw7x
*/
#include "nr3matlab.h"
#include "../code/gaussj.h"

MatDoub transpose(const MatDoub & x)
{
MatDoub ret(x.ncols(), x.nrows());
for (int i = 0; i < ret.nrows(); i++) {
for (int j = 0; j < ret.ncols(); j++) {
ret[i][j] = x[j][i];
}
}
return ret;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
MatDoub a(prhs[0]);
MatDoub b(prhs[1]);
MatDoub c(b.nrows(), b.ncols(), plhs[0]);

// The NR version of the matrices
MatDoub nra = transpose(a);
MatDoub nrb = transpose(b);

// Do the deed with NR3
gaussj(nra, nrb);

//Convert back to Matlab/Octave matrix form
c = transpose(nrb);
}

/*

Here's a simple test script (the answer is 1, 2, 3):

%testgauss.m
printf('Will solve ax = b using NR3 gaussj function\n')
a = [1 2 3;2 3 4;1 2 4]
b = [14; 20; 17]
x = xgauss(a,b)


Footnote:
I don't have access to Matlab. I am using Octave version 3.0.1
on my Linux systems. Using "mex xxxx.cpp" from an Octave command
line works OK, or from a Linux command line,
"mkoctfile --mex xxxx.cpp"
also does the trick.
*/


I think that's about the simplest way without modifying any NR stuff. See Footnote.

Here's what GNU Octave shows me for this test case:


octave:1> mex xgauss.cpp
octave:2> testgauss
Will solve ax = b using NR3 gaussj function
a =

1 2 3
2 3 4
1 2 4

b =

14
20
17

x =

1
2
3


In Matlab, b has one column and three rows. We need the corresponding MatDoub in xgauss.cpp to have one column and three rows in order to use it in the nr3 gaussj() function.

If you put in print statements that show the shapes of b and nrb, you can see something like:

b.nrows() = 1, b.ncols() = 3 <---From prhs(1)
nrb.nrows() = 3, nrb.ncols() = 1 <---From transpose(b)



Of course, Octave and Matlab have built-in functions/operators for solving systems of equations, so you don't need NR3 gaussj(), but I just wanted to show a little exercise for others who might be interested in getting started without any modifications to the distributed code. (And to point out that the published code seems to work as advertised.)

At least it works for me. As always: Your Mileage May Vary.

Regards,

Dave

Footnote:

The NRMatrix definitions in nr3matlab.h are pretty much the same as in code/nr3.h of the distribution, except that there is some new stuff having to do with Matlab mxArrays and other Matlab-specific stuff. If you change anything about non-Matlab NRMatrix constructors or functions, all of the nr3 functions that use NRMatrix objects are down the drain. (There are a couple of other changes in the definitions, but they don't interfere with functionality of the original nr3 code in the distribution.)

Furthermore, just interchanging the rôles of rows and columns in the NRMatrix definitions will not make it work with Matlab matrices. Matlab matrices are contiguous blocks of memory stored in column-major order (Ã* la Fortran). Even though the NRMatrix internal data values are stored in a contiguous block of memory, the mechanism to access them, C++ style, is defined as pointer-to-pointer for each row. So the internal storage is in row-major order, as it always is in C and C++. It is simply incompatible with the storage scheme of the Fortran-style matrices of Matlab. (At least, that's my opinion. I would really, really (really) like to see anything that people have to offer that wouldn't screw up the rest of the nr3 code.)

You have to copy them an element at a time. The simple "transpose" function in my example does that. You can hide the transpose functionality in a constructor and/or overloaded assignment operator for the rhs arguments of mexFunction. Handling the lhs argument is a little more interesting, I'm thinking.

If you are going to use NRMatrix objects in unmodified nr3 functions, somehow, sometime, some way, the program still has to do the work of allocating storage and copying the matrices an element at a time coming into and going out of mexFunction, regardless of how you dress it up.

Unless, of course you make the Matlab part of the code use transposes of all matrices as I/O data. For me, that's simply unacceptably ugly for "simple" stuff, but might be worthwhile if performance is an issue.

But if performance is really an issue, I (probably) won't be calling C++ functions from Matlab in the first place. If the functionality isn't available or is not efficient enough in Matlab implementations, I will (probably) just use C++ (or Fortran or whatever other language that I feel is appropriate) natively. Again: Just a vague generality, expressed as a personal opinion, YMMV.

I completely missed the transpose mention on the examples page, so thanks for pointing that that out. I will have to code my wrappers so that the transpose is taken into account and see how that works. Its nice that you have another confusing element in an already confusing mixture of code.

I too would normally code something in a native language when I needed it. I have about 20 people though who code in matlab, where we sometimes do not have the functionality we need and by adding some nifty numerical recipes codes to matlab they can do a lot more, like fitting functions without me paying 1000 dollars a year for a toolbox license as one example.

Matlab is easy to code in and easy to learn and except for a few slow spots is an excellent platform for modeling, simulation and analysis, made better with some extra mex functions.

Will Handler

...by adding some nifty numerical recipes codes to matlab...

I like the sound of that. I hope that you can share with us here...

Regards,

Dave