[Bug-gsl] possible matrix solver bug

[Andrew Wang]

Hi,

I tried the linear algebra example with slight modification (do
decompose and solve one more time) and got different results from the
same matrix and vector.
Did I miss something here?

Andrew Wang
-----------------------------------------------------------------------------------------
#include <stdio.h>
#include <gsl/gsl_linalg.h>

int
main (void)
{
  double a_data[] = { 1, 1, 1,
                      1, -1, 1,
                      4, 2, 1 };

  double b_data[] = { 6, 8, 11 };

  gsl_matrix_view m
    = gsl_matrix_view_array (a_data, 3, 3);

  gsl_vector_view b
    = gsl_vector_view_array (b_data, 3);

  gsl_vector *x = gsl_vector_alloc (3);

  int s;

  gsl_permutation * p = gsl_permutation_alloc (3);

  gsl_linalg_LU_decomp (&m.matrix, p, &s);

  gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);

  printf ("x = \n");
  gsl_vector_fprintf (stdout, x, "%g");

// one more time, same matrix and vector

  gsl_linalg_LU_decomp (&m.matrix, p, &s);

  gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);

  printf ("x = \n");

  gsl_vector_fprintf (stdout, x, "%g");

  gsl_permutation_free (p);

  return 0;
}


x =
2
-1
5
x =
-1.41379
0.12931
11.3966