#1




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 80bit registers but would be satisfied with 64bit values for x and y. The fix is to force the compiler to make 64bit 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][l1]) + s == s) { to temp=fabs(a[l][l1]) + s; if (temp == s) { 3) Change if (u+v == v) break; to temp=u+v; if (temp == v) break; Last edited by Saul Teukolsky; 09012003 at 01:21 PM. 
#2




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

#3




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. 
#4




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




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 
#6




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




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 
#8




A resistant matrix
Dear all,
A matrix which I placed at the bottom of my homepage http://www.gub.ruhrunibochum.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 divideconquer 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 ? 
#9




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




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 MSWindowsspecific. I have put everything on my homepage (at the bottom) http://www.ruhrunibochum.de/gub/mi...j_niemunis.htm I write explictly which matrix is critical for which routine in the MATRX3.DATfiles. Aha, please read the matrices columnwise and not rowwise as the original xsvdcmp.f90 does !!!!!! Lapack (the divideconquer version of SVD) and Mathematica can handle both matrices 
#11




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. 
#12




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 runtime 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 
#13




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. 
#14




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. 
#15




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. 
Thread Tools  
Display Modes  

