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; _______________________________________________ Bug-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/bug-gsl
