possible bug in linalg/balance.c

Aaron Schweiger aschweig@bu.edu
Fri Feb 27 16:00:00 GMT 2004

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;


More information about the Gsl-discuss mailing list