Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Numerical Recipes Official Announcements > Old Bug Reports (Obsolete Editions of NR)

Reply
 
Thread Tools Display Modes
  #1  
Old 09-01-2003, 10:52 AM
Saul Teukolsky Saul Teukolsky is online now
Numerical Recipes Author
 
Join Date: Dec 2001
Posts: 211
Possible convergence problems in svdcmp, jacobi, tqli, hqr

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;

Last edited by Saul Teukolsky; 09-01-2003 at 01:21 PM.
Reply With Quote
  #2  
Old 03-03-2005, 09:40 PM
forall forall is offline
Registered User
 
Join Date: Jul 2002
Posts: 10
wouldnt a more thorough solution be to add some tolerance when testing for equality of floating point numbers?
Reply With Quote
  #3  
Old 03-04-2005, 07:12 AM
Saul Teukolsky Saul Teukolsky is online now
Numerical Recipes Author
 
Join Date: Dec 2001
Posts: 211
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.
Reply With Quote
  #4  
Old 04-07-2005, 06:37 PM
jaydogg jaydogg is offline
Member
 
Join Date: Apr 2005
Posts: 3
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). )
Reply With Quote
  #5  
Old 07-20-2005, 12:58 PM
jpmcphee jpmcphee is offline
Member
 
Join Date: Jul 2005
Posts: 2
fortran?

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
Reply With Quote
  #6  
Old 07-20-2005, 06:36 PM
Saul Teukolsky Saul Teukolsky is online now
Numerical Recipes Author
 
Join Date: Dec 2001
Posts: 211
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.
Reply With Quote
  #7  
Old 08-05-2005, 10:13 AM
ernst ernst is offline
Member
 
Join Date: Aug 2005
Location: Belgium
Posts: 1
Bug in svdcmp

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
Reply With Quote
  #8  
Old 05-23-2006, 03:29 PM
Andrzej Andrzej is offline
Registered User
 
Join Date: May 2006
Posts: 2
A resistant matrix

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 ?
Reply With Quote
  #9  
Old 05-23-2006, 08:08 PM
Saul Teukolsky Saul Teukolsky is online now
Numerical Recipes Author
 
Join Date: Dec 2001
Posts: 211
I'm unable to reproduce your problem. svdcmp.f90 works fine on the given matrix on my setup.

Saul Teukolsky
Reply With Quote
  #10  
Old 05-24-2006, 09:00 AM
Andrzej Andrzej is offline
Registered User
 
Join Date: May 2006
Posts: 2
a resistant matrix 2

....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
Reply With Quote
  #11  
Old 08-14-2006, 01:59 PM
Saul Teukolsky Saul Teukolsky is online now
Numerical Recipes Author
 
Join Date: Dec 2001
Posts: 211
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.
Reply With Quote
  #12  
Old 06-14-2011, 08:47 AM
cpgoodri cpgoodri is offline
Registered User
 
Join Date: Jun 2011
Posts: 1
Convergence error in svdcmp

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
Reply With Quote
  #13  
Old 07-01-2013, 06:01 PM
mayes88 mayes88 is offline
Registered User
 
Join Date: Jul 2013
Posts: 2
Urgent problem with svdfit.f, svdcmp.f,svbksb.f f77 bugs?

Hello,
Using the three subroutines above as prescribed in NR fortran, the results are not converging.
I am using the f77 code downloaded 09 June 2013. Has this code been modified with tolerance checks as the f90 mentioned previously?
If there are mods, please reply in f77, full subroutines if possible.
I can also provide data if necessary.
Thank you.
Reply With Quote
  #14  
Old 07-03-2013, 04:09 PM
mayes88 mayes88 is offline
Registered User
 
Join Date: Jul 2013
Posts: 2
SVD routines not working

Is this forum still active?

Please see my previous post. Please also provide debugged version of svdcmp.f ( f77 ) with the new tolerance checks inserted. The new checks referenced in earlier posts refer to svdcmp? Name of subroutine and line number for insertion would be helpful, but I'd rather have the corrected routines, purchased June 2013.
Reply With Quote
  #15  
Old 07-08-2013, 10:22 PM
Bill Press Bill Press is offline
Numerical Recipes Author
 
Join Date: Dec 2001
Posts: 227
Hello. This is a forum where, with a few exceptions, questions are asked and answers provided by NR users, not by the authors.

The NR 2nd edition Fortran code dates from 1992 and has not been a supported product since 2007. So it would not have any changes since then. It might not even use the same logic as the 3rd edition code.

Hope this helps.
Reply With Quote
Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off

Forum Jump


All times are GMT -5. The time now is 12:50 PM.


Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.