Re: [Bug-gsl] gsl_eigen_symmv

[Brian Gough]

At Tue, 10 Nov 2009 08:02:38 +0100,
Katrin Wolff wrote:
> I sent the wrong version of the program before, which will not reproduce 
> the error. I attach the correctly faulty one to this mail.

I have been able to reproduce the problem now, which occurs because of
the presence of very small non-zero values that prevent convergence.
The patch below should avoid the problem.  Thanks for the bug report.

-- 
Brian Gough

GNU Scientific Library -
http://www.gnu.org/software/gsl/

diff --git a/eigen/qrstep.c b/eigen/qrstep.c
index 26705a2..a8b51c4 100644
--- a/eigen/qrstep.c
+++ b/eigen/qrstep.c
@@ -104,6 +104,16 @@ qrstep (const size_t n, double d[], double sd[], double 
gc[], double gs[])
 
   double mu = trailing_eigenvalue (n, d, sd);
 
+  /* If mu is large relative to d_0 and sd_0 then the Givens rotation
+     will have no effect, leading to an infinite loop.  
+
+     We set mu to zero in this case, which at least diagonalises the
+     submatrix [d_0, sd_0 ; sd_0, d_0] and allows further progress. */
+
+  if (GSL_DBL_EPSILON * fabs(mu) > (fabs(d[0]) + fabs(sd[0]))) { 
+    mu = 0;
+  }
+
   x = d[0] - mu;
   z = sd[0];