John Knowles
08-20-2008, 09:00 AM
I understand the TTLG principle, having worked over the years with Matlab and the old C version of NR. In applications where this matters, how do people deal with it? I am not a C++ programmer, having just started on the new version of NR two days ago.
Problem 1)
how to translate a matlab matrix, A=[1 2 3;4 5 6;...], that will be used for computations in the C++ side? I need to be able to access the Matlab rows, ideally row i would be A[i], as was done in the old C recipes. Is there a trick in Matlab, say using reshape?
[OK, i see how to do this now, using transpose and redefining the row/colums in C++, pretty easy all of a sudden].
Problem 2)How to do the inverse, so that matrices computed in C come out on the Matlab side with the same rows?
i realize that I could do these operations with a couple of for loops and some new memory allocations, but these are large matrices so I need the most efficient way to do it. In C I would rearrange the elements using pointers without reallocation the memory the data lived in.
davekw7x
08-21-2008, 10:20 AM
...having worked over the years with Matlab and the old C version of NR. ..I am not a C++ programmer
The issue is not Matlab versus Numerical Recipes; it's Fortran-style array accessing versus C/C++ methods.
Under the hood, C and C++ have identical notation and methods of accessing arrays, so whatever worked for C works for C++. The difference between C versions of NR and C++ versions is that the arrays are now zero-based rather than one-based. (So index values in the C++ versions go from zero to N-1 instead of one to N, as they did in the C versions.)
Problem 1)...
[OK, i see how to do this now...]So: problem solved, right?
Problem 2)How to do the inverseDo it the same way: transpose the resultant matrix objects before returning to Matlab.
I gave an example, going both ways, in post number 2 of this thread: http://www.nr.com/forum/showthread.php?t=1506
...In C I would rearrange the elements using pointers without reallocation the memory the data lived in...
Can you give an example of how to perform the notational equivalent of the transposition of an nxm 2-D matrix in C without allocation/reallocation and/or copying? I hate to repeat myself, but whatever works in C also works in C++.
I am assuming that you desire the following:
1. Values stored in a Fortran-style matrix (array) are to be accessed by C++ programs using standard C++ 2-D array notation. In our case, the Fortran-style matrix came from normal usage in Matlab, and rows and columns were defined according to "normal" Matlab operations. See Footnote.
2. The C++ programs themselves (like Numerical Recipes functions) have been written so that an expression like x[r][c] refers to the element in row r and column c.
3. No alterations of the C++ programs are allowed: No new notation, no interchanging the rĂ´les of rows and columns in the source code of the algorithms themselves.
4. In the interests of efficiency, you don't want the program to do a bunch of allocation/reallocation of storage and copying of data elements from the original data structure into your C++ object.
5. You also want to do the reverse operation: Convert from the kind of C++ object used in Numerical Recipes, wherein the data values are stored in a contiguous block of memory in row-major order, back to a Fortran-style data structure with elements in column-major order.
Here's my opinion. It is only an opinion, and I could be wrong. (It wouldn't be the first time.)
It can't be done.
You can't change the way that C++ handles the addressing of something like x[r][c] whether x was declared to be a 2-D array (not very practical in general matrix-handling functions) or whether x is some kind of pointer-to-pointer thing, including a C++ object like an NRmatrix that uses the notation to access an internal data array through a pointer-to-pointer mechanism.
The C++ notation shown in item 2 of the list is, simply, incompatible with accessing a block of memory that contains a 2-D Fortran-style array (values stored in column-major order).
In a C++ program you can hide the work of matrix transposition inside a class that takes care of allocation and copying in its constructors and destructors, but, fundamentally, the program has to re-allocate storage and/or copy the data items. Somehow, somewhere, sometime.
I hate to repeat myself, but please note (and this is very important):
I could be wrong.
It wouldn't be the first time. Not even the first time today.
Now, generally speaking, no one likes to be wrong, but in this case, I would really, really like to be wrong. I would really, really appreciate any insight that someone can offer.
Bottom line: If I am wrong, and all of the desires expressed in the list can be satisfied, then someone will post a way, and I will have learned something valuable.
Unless and until someone shows us the way, I hope that (temporarily) accepting the alternative (doing the transpose thing) will not be unacceptable. I mean, knowing how to do something that can be done is a pretty good fallback position, right?
Additional thoughts (assuming that post-bottom-line ramblings are allowed):
I enumerated the list so that the individual concerns can be addressed in case I have incorrectly deduced your desires. Maybe I extrapolated something you didn't actually express, and your actual requirements are different in one or more details???
You are concerned about efficiency because these are "really large matrices." I think that it's always appropriate to consider efficiency, but here's the way that I see it:
As far as the amount of memory involved: Dynamic allocation of memory used for NRmatrix objects (or other C/C++ dynamically allocated arrays) is limited to something like 2 Gigabytes with compilers that I have for my Windows XP systems and 3 Gigabytes on my 32-bit Linux systems. I can't believe that memory size is an issue. (The ability to access two or three Gigabytes has nothing to do with the amount of physical RAM in the system, although performance may be better with more installed RAM.)
If you are concerned about CPU time, consider that the computational complexity of matrix transposition is Order(n^2)---actually for an nxm matrix, it is Order(n*m). If you are doing meaningful operations on matrices, are they of such an order that the transpose function dominates the CPU time?
For example, Gauss-reduction based matrix inversion routines are Order(n^3), and the most efficient matrix inversion routine that I have seen is something like Order(n^2.376). For larger and larger matrices, the time spent in the transposition function coming from and returning to Matlab becomes less and less important.
Just a thought. Or two.
Regards,
Dave
Footnote: You can access a particular Fortran-style array object in C or C++ without reallocation and copying by using something like the following in your C++ source file:
#define mat(r,c) mat[c][r]
Then, everywhere in your C++ program that you have mat[something][somethingelse], just replace it with mat(something,somethingelse). In my opinion, whether this leads to a meaningful run-time advantage in a practical application is open to question, but it might be fun to play with. Or not.
In my experience, this kind of stuff (use of macros to make C look like Fortran---with or without transposition) drives real C++ programmers bananas. Even dilettantes like yours truly sometime have to "take a moment" to calm down after trying to help people debug programs that incorporate such monstrosities.