Possible convergence problems in svdcmp, jacobi, tqli, hqr

[Saul Teukolsky]

Recent compiler advances may cause the routines svdcmp, jacobi, tqli and hqr to fail to converge for certain problems. So far we have seen the problem only with the gcc/g++ compiler, version 3.0 or later, on pentium machines only, and only for a single example involving svdcmp. We would like to hear reports involving other compilers/hardware.

The problem arises in convergence tests for negligible x of the form

if (x+y == y)

which can fail when values are kept in 80-bit registers but would be satisfied with 64-bit values for x and y. The fix is to force the compiler to make 64-bit tests by writing the intermediate value to memory:

temp=x+y;
if (temp == y)

As an added precaution, declare temp to be of type "volatile". This prevents the compiler from optimizing away the storage of temp (although the gcc compiler isn't yet smart enough to do this.)

Here is a complete list of changes that need to be implemented:

svdcmp.cpp has 4 changes:

1) After the 7th line of code, insert the declaration

DP volatile temp;

2) Change

if (fabs(rv1[l])+anorm == anorm) {

to

temp=fabs(rv1[l])+anorm;
if (temp == anorm) {

3) Change

if (fabs(w[nm])+anorm == anorm) break;

to

temp=fabs(w[nm])+anorm;
if (temp == anorm) break;

4) Change

if (fabs(f)+anorm == anorm) break;

to

temp=fabs(f)+anorm;
if (temp == anorm) break;

jacobi.cpp has 3 changes:

1) After the 4th line of code from the start of the jacobi routine, insert

DP volatile temp1,temp2;

2) Change

if (i > 4 && (fabs(d[ip])+g) == fabs(d[ip])
&& (fabs(d[iq])+g) == fabs(d[iq]))

to

temp1=fabs(d[ip])+g;
temp2=fabs(d[iq])+g;
if (i > 4 && temp1 == fabs(d[ip]) && temp2 == fabs(d[iq]))

3) Change

if ((fabs(h)+g) == fabs(h))

to

temp1=fabs(h)+g;
if (temp1 == fabs(h))

tqli.cpp has 2 changes:

1) After the 7th line of code, insert

DP volatile temp;

2) Change

if (fabs(e[m])+dd == dd) break;

to

temp=fabs(e[m])+dd;
if (temp == dd) break;

hqr.cpp has 3 changes:

1) After the 8th line of code, insert

DP volatile temp;

2) Change

if (fabs(a[l][l-1]) + s == s) {

to

temp=fabs(a[l][l-1]) + s;
if (temp == s) {

3) Change

if (u+v == v) break;

to

temp=u+v;
if (temp == v) break;

[forall]

wouldnt a more thorough solution be to add some tolerance when testing for equality of floating point numbers?

[Saul Teukolsky]

Yes, it would be better to test against a tolerance. In the next software release we will do so. The idea is something like this:

#include <limits>
...
const DP EPS=numeric_limits<DP>::epsilon();
...
if (fabs(rv1[l]) <= EPS*anorm) {

But the workaround we gave earlier should be fine for most compilers.

[jaydogg]

I have noticed this problem manifested in 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]

Even with the fix as described in this thread added (I used a DOUBLE temp variable), the problem still persisted. I even tried using the epsilon approach, but that did not work either.

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). )

[jpmcphee]

Hi, do these changes have corresponding versions for FORTRAN compilers?

I'm asking because noticed that the results from jacobi differ from those obtained with MATLAB for small eigenvalues.

Thanks,

James

[Saul Teukolsky]

These changes are only necessary to fix convergence problems. They don't affect values obtained.

One reason that small eigenvalues can differ is that the algorithm used only finds the eigenvalues with a certain absolute error, not relative error.

[ernst]

I have observed recently that the svdcmp routine has convergence problems when dealing with matrices which are products of elementary rotation matrices and, therefore, which are orthonormal. Do we have guarantees that the algorithm is indeed well conditionned when dealing with such matrices ?

Damien Ernst

[Andrzej]

Dear all,

A matrix which I placed at the bottom of my homepage
http://www.gub.ruhr-uni-bochum.de/mi...j_niemunis.htm
cannot be decomposed by svdcmp.f90. The recommended
in this forum modifications do not help. The
Koji Mukai's version of the Marquardt routine marq.f90 (improved by Tim Naylor and Lee Howells) works better
but still it cannot SVdecompose my matrix.

It is in contrast to Lapack (the divide-conquer version of SVD) or Mathematica which can handle this matrix ( I have checked multiplying the original by the inverse V . (1/W). U^T and it is ok.)

The diagonal values W are all between 1.0d0
and 2.0d0 so there is no singularity. I would prefer
NR version to the huge Lapack collection of routines
but how to make it work ?

[Saul Teukolsky]

I'm unable to reproduce your problem. svdcmp.f90 works fine on the given matrix on my setup.

Saul Teukolsky

[Andrzej]

....Thank you very much indeed for the prompt answer.

I must apologize for my recent statement about svdcmp.f90 and marq.f90. It was not quite precise. Actually there are matrices which can be solved by the former routine but not by the latter and vice versa.

My yesterday's matrix causes only marq.f90 to fail
and indeed svdcmp.f90 works fine with it.

This time I present TWO matrices. Each matrix
causes the different procedure to fail. I also give
two complete sets of files needed to compile and run both projects so perhaps somebody can check if
the problem is MS-Windows-specific.

I have put everything on my homepage (at the bottom)
http://www.ruhr-uni-bochum.de/gub/mi...j_niemunis.htm
I write explictly which matrix is critical for which routine in the MATRX3.DAT-files.

Aha, please read the matrices column-wise and
not row-wise as the original xsvdcmp.f90 does !!!!!!


Lapack (the divide-conquer version of SVD) and Mathematica can handle both matrices

[Saul Teukolsky]

The problem is taken care of by the fix suggested in post #3 above. In fortran 90, the necessary code is
REAL(DP), PARAMETER :: EPS=epsilon(x)

and then

if (abs(rv1(l)) <= EPS*anorm) exit
if (abs(w(nm)) <= EPS*anorm) then
...
if (abs(f) <= EPS*anorm) exit

at the appropriate places.

[cpgoodri]

Hi,

I am trying to SVdecompose a matrix using the svdcmp method in the c++ v2.10 code with the g++ compiler. For very very small matrices this works fine, but when the matrix gets of the order 20 by 20 I get the message

Numerical Recipes run-time error...
no convergence in 30 svdcmp iterations
...now exiting to system...

I have tried increasing the number of iterations to over 100,000 with no luck. I have also tried including the suggestions from the first post, but that didn't work either. Does anyone know what might be happening?

Thanks so much,
Carl