Re: [Bug-gsl] Possible bug in gsl_multifit_wlinear_svd

[Brian Gough]

The following patch should fix the problem and will be in the next
release.  Thanks for the bug report.

-- 
Brian Gough
(GSL Maintainer)

Network Theory Ltd,
Publishing the GSL Manual - http://www.network-theory.co.uk/gsl/manual/

Index: gsl/linalg/svdstep.c
diff -u gsl/linalg/svdstep.c:1.9 gsl/linalg/svdstep.c:1.10
--- gsl/linalg/svdstep.c:1.9    Mon Apr 24 19:20:11 2006
+++ gsl/linalg/svdstep.c        Fri Aug 17 10:21:29 2007
@@ -54,7 +54,31 @@
 create_schur (double d0, double f0, double d1, double * c, double * s)
 {
   double apq = 2.0 * d0 * f0;
-  
+
+  if (d0 == 0 || f0 == 0)
+    {
+      *c = 1.0;
+      *s = 0.0;
+      return;
+    }
+
+  /* Check if we need to rescale to avoid underflow/overflow */
+  if (fabs(d0) < GSL_SQRT_DBL_MIN || fabs(d0) > GSL_SQRT_DBL_MAX
+      || fabs(f0) < GSL_SQRT_DBL_MIN || fabs(f0) > GSL_SQRT_DBL_MAX
+      || fabs(d1) < GSL_SQRT_DBL_MIN || fabs(d1) > GSL_SQRT_DBL_MAX)
+    {
+      double scale;
+      int d0_exp, f0_exp;
+      frexp(d0, &d0_exp);
+      frexp(f0, &f0_exp);
+      /* Bring |d0*f0| into the range GSL_DBL_MIN to GSL_DBL_MAX */
+      scale = ldexp(1.0, -(d0_exp + f0_exp)/4);
+      d0 *= scale;
+      f0 *= scale;
+      d1 *= scale;
+      apq = 2.0 * d0 * f0;
+    }
+
   if (apq != 0.0)
     {
       double t;