I am using the single precision version of the C (not C++) routine svdcmp.c -

I am using Solaris 5.7 with gcc 2.95.2.

I have generated an example where the u matrix returned by svdcmp() will have NaN's in it's last column. The input matrix I used was 6 rows by 6 columns, with every row equal to the following:

[55 52 88 30 58 58]

Can anyone else reproduce these results so I know it's not just me?

( I have not yet noticed this problem in the DOUBLE precision version of svdcmp (dsvdcmp). )

Wow that is weird - I got the same thing. I am not using Solaris, but I am using some version of Linux. It must be something wrong with the algorithm because I downloaded this straight from NR. Some numerical problems or something....

Sorry guys, I am unable to reproduce your problem. I get singular values of 355.817 and 5 zeros. No NaN's.

Saul Teukolsky

I think I may have found the problem -

In the very beginning of the svd algorithm there is a Householder reduction. Within this reduction there is a check performed on each absolute column sum of the 'a' matrix to see if the sum is more than zero. The code looks like:

for (k=i;k<=m;k++) scale += fabs(a[k][i]);
if (scale) {
.....

The problem was that this column sum was very close to zero (on the order of 10e-30), but not quite zero. Consequently, the w[i] variable was set to a small value:

w[i] = scale*g;

I won't illustrate the path, but this small w[i] used in subsequent calculations eventually led to an INF value, which in turn led to a INF/zero calculation, and INF/zero = NaN.

Looking at the description of the Householder algorithm in the NR book Section 11.2 (page 473), the criteria for "skipping the transformation" is that the absolute column sum is zero "to machine precision".

In other words, the "if(scale)" condition above should really be something like "if(scale>machine_zero)", where

machine_zero = (smallest_positive number representable) / (machine precision)

When I made the above modification, the output did not have NaN's anymore, and the results were accurate.