This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
RE: possible bug in linalg/balance.c
- From: "Aaron Schweiger" <aschweig at bu dot edu>
- To: "Brian Gough" <bjg at network-theory dot co dot uk>, "Gsl-Discuss" <gsl-discuss at sources dot redhat dot com>
- Date: Fri, 27 Feb 2004 11:00:42 -0500
- Subject: RE: possible bug in linalg/balance.c
Dear Brian,
I modified the multilinear fitting example to create this test problem that
illustrates the hang that occurs when one of the input data is infinity.
(GSL 1.4, GCC 3.2.3)
Note, my suggested fix to linalg/balance.c is to check whether s is
infinite, and if so, invoke GSL_ERROR. However - I don't know if this is
the 'correct' behavior of the balance algorithm - as I am not familiar with
it. Thanks again,
Aaron Schweiger
----------------------------------
/*
lstest.c - demonstrates bug in linalg/balance.c
this code should hang GSL 1.4 due to loops in linalg/balance.c that
try to halve infinity.
*/
#include <stdio.h>
#include <gsl/gsl_multifit.h>
int main (void)
{
int i, n = 4;
double chisq;
gsl_matrix *X, *cov;
gsl_vector *y, *w, *c;
puts("Allocating vectors and matrix.");
X = gsl_matrix_alloc (n, 3);
y = gsl_vector_alloc (n);
w = gsl_vector_alloc (n);
c = gsl_vector_alloc (3);
cov = gsl_matrix_alloc (3, 3);
puts("Setting up values, with spurious infinity.");
for (i = 0; i < n; i++)
{
gsl_matrix_set (X, i, 0, 1.0);
gsl_matrix_set (X, i, 1, (double) i );
gsl_matrix_set (X, i, 2, 1.0/((double) i));
gsl_vector_set (y, i, sqrt((double) i));
gsl_vector_set (w, i, 1.0);
}
puts("Before Regression block, setting up workspace.");
{
gsl_multifit_linear_workspace * work
= gsl_multifit_linear_alloc (n, 3);
puts("Running regression.");
gsl_multifit_wlinear (X, w, y, c, cov,
&chisq, work);
puts("De-allocating workspace.");
gsl_multifit_linear_free (work);
}
puts("Reporting results.");
#define C(i) (gsl_vector_get(c,(i)))
printf ("# best fit: Y = %g + %g X + %g X^2\n",
C(0), C(1), C(2));
return 0;
}
------------------------------