[Help-gsl] gsl_linalg_SV_decomp segfaults, but I think I've called it correctly

[Tom Vajzovic]

Hi!

I'm trying to solve a system of linear equations.

When I call the function gsl_linalg_SV_decomp, my program segfaults.

I compiled with  

$ gcc -lgsl -lgslcblas -lm test.c -o test  

I'd appreciate any help. 

Please ask if I've left out any details that are relavant.

Thanks!

Tom

Code follows:

/* test.c */

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

int main(void){

  gsl_vector *work, *b, *x;
  gsl_matrix *A, *V, *S;

  gsl_set_error_handler(NULL);

  work= gsl_vector_alloc( 3 );
  b= gsl_vector_alloc( 3 );
  x= gsl_vector_alloc( 3 );
  A= gsl_matrix_alloc( 3, 3 );
  V= gsl_matrix_alloc( 3, 3 );
  S= gsl_matrix_alloc( 3, 3 );

  if( !work || !b || !x || !A || !V || !S ){
    fprintf( stderr, "fail\n" );
    gsl_vector_free( work );
    gsl_vector_free( b );
    gsl_vector_free( x );
    gsl_matrix_free( A );
    gsl_matrix_free( V );
    gsl_matrix_free( S );
    exit( EXIT_FAILURE );
  }

  gsl_vector_set( b, 0, 82.0 );
  gsl_vector_set( b, 1, 15.0 );
  gsl_vector_set( b, 2, 21.0 );

  gsl_matrix_set( A, 0, 0, 2.0 );
  gsl_matrix_set( A, 0, 1, 3.0 );
  gsl_matrix_set( A, 0, 2, 9.0 );
  gsl_matrix_set( A, 1, 0, 1.0 );
  gsl_matrix_set( A, 1, 1, 1.0 );
  gsl_matrix_set( A, 1, 2, 1.0 );
  gsl_matrix_set( A, 2, 0, 0.0 );
  gsl_matrix_set( A, 2, 1, 0.0 );
  gsl_matrix_set( A, 2, 2, 3.0 );

  fprintf( stderr, "got to line %d\n", __LINE__ );

  if(  gsl_linalg_SV_decomp( A, V, S, work )
    || (fprintf( stderr, "got here\n"),0)
    || gsl_linalg_SV_solve ( A, V, S, b, x )
     ){

    fprintf( stderr, "fail\n" );
    gsl_vector_free( work );
    gsl_vector_free( b );
    gsl_vector_free( x );
    gsl_matrix_free( A );
    gsl_matrix_free( V );
    gsl_matrix_free( S );
    exit( EXIT_FAILURE );
  }

  fprintf( stderr, "got to line %d\n", __LINE__ );

  fprintf( stdout, " x= %.1f \n y= %.1f \n z= %.1f \n",
                 gsl_vector_get ( x, 0 ),
                 gsl_vector_get ( x, 1 ),
                 gsl_vector_get ( x, 2 ) );

  gsl_vector_free( work );
  gsl_vector_free( b );
  gsl_vector_free( x );
  gsl_matrix_free( A );
  gsl_matrix_free( V );
  gsl_matrix_free( S );

    exit( EXIT_SUCCESS );
}
/*

This should solve:

 Ax = b

[ 2 3 9 ] [ x ]   [ 82 ]
[ 1 1 1 ] [ y ] = [ 15 ]
[ 0 0 3 ] [ z ]   [ 21 ]

which has solution:

[ x ]   [ 5 ]
[ y ] = [ 3 ]
[ z ]   [ 7 ]

But instead it segfaults on gsl_linalg_SV_decomp()

*/