Hello.

Chapter 2-6 is a bit unclear on this point. The text reads as if the algorithm depends on the fact that M<=N. But then it says "For an arbitrary matrix A...".

Empirical evidence suggests that the number of rows cannot exceed the number of columns. In that case, the code crashes badly. Note I am not a new user of svdcmp as I have been happily using it on matrices large and small literally millionos of times.

Any clarification of intended support or experience to the contrary would be appreciated. I am not a professional math person so if there is a workaround I would be glad to know it. For example, can I transpose the matrix first?

Thanks for your attention.

Hi David,

svdcmp works on arbitrary matrices. The text describes the format of the output for M >= N in eqn. (2.6.1), and for M < N after eqn. (2.6.4). If you are using the C++ code, make sure you have the 2 bug fixes from v. 2.11 that are listed in the Official Bug Reports in this forum.

Hope this sorts things out for you.

Saul Teukolsky

Thank you for your response. It was helpful to have this point clarified as it forced me to re-examine "my" code. Some time ago, I had based a generic matrix inversion function on the xsvdcmp sample code.

Until now, my matrices have all been square. So some dimension issues have been hidden. I managed to address enough issues to get the first set of rectangular matrices working. But in all cases M>N.

Of course it was only a matter of time before I would need M<N. Consistent failure in all of these cases led to my first post above. As most of you probably expected, the error was not in Numerical Recipes but in "my" code.

Some things I learned in the process...

The number of "valid" singular values seems to be min(M,N).
However you must allocate array w to have max(M,N) elements.
Just looking at the pictures in 2-6, I would have expected w to always be of size N.
If these are indeed true, the xsvdcmp sample code would fail if actual M and N values were used in allocation instead of MP and NP.

When M<N, one of the "valid" singular values is always zero.

When M>=N, Ainverse.A = I.
When M<N, A.Ainverse = I.

Comments welcome.