gsl_linalg_SV_decomp goes into infinite loop on reasonable data

[tlb at trevorblackwell dot com]

On my system, FreeBSD 4.5-RELEASE-p2, gsl_linalg_SV_decomp gets stuck
in an infinite loop trying to reduce the matrix to diagonal. The test:

   <linalg/svd.c:103>
        if (gsl_vector_get (&f.vector, b - 1) == 0.0)

repeatedly fails so the loop it's in never terminates. The value it's
getting is 1.83429e-07, noticeably more than DBL_EPS. (I'm using
double throughout.) In the chop_small_elements code, the line which is
supposed to remove these small values:

    <linalg/svdstep.c:14>
      if (fabs (f_i) < GSL_DBL_EPSILON * (fabs (d_i) + fabs (d_ip1)))

has f_i= 1.83429e-07, and d_i and d_ip1 both around 1e-4.

If I hack the test (by changing GSL_DBL_EPSILON to 0.01) it just fails
with an even larger value left at the end of the matrix. So I suspect
that the qrstep function isn't doing its job.

This came up when doing a multifit_linear on some fairly clean data.

You should be able to reproduce with the following test program:


#include <stdlib.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>


main()
{
  gsl_matrix *a,*v;
  gsl_vector *s,*work;
  int i,j;

  double data[25]={
    0, 0.082610689977785284, 0.049070739973785735, -0.066400265730392663, 
    -0.068708124015300312, 0, -0.25203574259670708, -0.18738417047093653, 
    0.2464027866048136, 0.26239284329959578, 0, 0, 0.020510937347991507, 
    -0.024112785427691547, -0.028732791501641798, 0, 0, 0, 
    -0.00014255001561252031, -1.785808972985471e-06, 0, 0, 0, 0, 
    2.1717171634033149e-07
  };

  a=gsl_matrix_alloc(5,5);
  v=gsl_matrix_alloc(5,5);
  s=gsl_vector_alloc(5);
  work=gsl_vector_alloc(5);

  for (i=0; i<5; i++) {
    for (j=0; j<5; j++) {
      gsl_matrix_set(a,i,j,data[i*5+j]);
    }
  }

  gsl_linalg_SV_decomp(a,v,s,work);

}  


--
Trevor Blackwell       tlb@trevorblackwell.com        (650) 776-7870