Hi,

I've encountered a bug in svdcmp() (NR in C) which I've searched on the web and no one seems to have published a fix for.

In the loop for "Diagonalization of the biadiagonal form", we have the following code:

flag=1;
for (l=k;l>=1;l--) {
nm=l-1;
if ((float)(fabsf(rv1[l])+anorm) == anorm) {
flag=0;
break;
}
if ((float)(fabsf(w[nm])+anorm) == anorm) break;
}


Consider the case where for some reason, the if conditions are never true (either due to precision issues, or NaNs in the variables). Then when l == 1, nm becomes 0 and the access into w[] becomes out-of-bounds.

I tried doing a local fix where I change that line to

if (nm > 0 && (float)(fabsf(w[nm])+anorm) == anorm) break;


However, this doesn't work as code further down the line assumes that l > 1. In an attempt to fix this I've inserted this code prior to the w[nm] access

if (l == 1) { // Added to avoid crash
flag = 0;
break;
}


That seems to avoid all the out-of-bounds accesses (thereby stopping potential crashing). However, does anyone know what the *right* fix is? I've since added some precautions to ensure that the matrix doesn't have NaNs but in general, I think we should have a proper fix to make this code robust. Any ideas?

Thanks,
-Edward

Dear Edward,

This issue is dealt with in the "SVD bad indices" posting
http://www.nr.com/forum/showthread.php?t=605
As you will see there, there is no problem with the code (although the commenting could be more explicit).

Saul Teukolsky

Hi Saul,

I think the reason why I hit it was because I had NaNs then (as previously mentioned). In that case, the if conditions will fail even when rv1[1] is 0 since anorm in my case is likely NaN.

Unfortunately, I'm at home right now and can't double-check. I will try it again tomorrow and report back. Thanks for taking the time to reply!

Cheers,
-Edward

Yep, it's because anorm was NAN. So let's suppose I want to make it robust in face of NAN's (robust in the sense of not crashing). What's the right approach?

Thanks,
-Edward

Hi Edward,

I think the "fixes" you inserted are as good as anything. Once the code has NaN's, the results aren't meaningful anyway.

Saul Teukolsky

I suppose that's true.

Thanks again!
-Edward