RE: possible bug in linalg/balance.c

["Aaron Schweiger" <aschweig at bu dot edu>]

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;
}

------------------------------