Examples

The following program solves the linear system @math{A x = b}. The system to be solved is,

and the solution is found using LU decomposition of the matrix @math{A}.

#include <stdio.h>
#include <gsl/gsl_linalg.h>

int
main (void)
{
  double a_data[] = { 0.18, 0.60, 0.57, 0.96,
                      0.41, 0.24, 0.99, 0.58,
                      0.14, 0.30, 0.97, 0.66,
                      0.51, 0.13, 0.19, 0.85 };

  double b_data[] = { 1.0, 2.0, 3.0, 4.0 };

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

  gsl_vector_view b
    = gsl_vector_view_array(b_data, 4);

  gsl_vector *x = gsl_vector_alloc (4);
  
  int s;

  gsl_permutation * p = gsl_permutation_alloc (4);

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

Here is the output from the program,

x = -4.05205
-12.6056
1.66091
8.69377

This can be verified by multiplying the solution @math{x} by the original matrix @math{A} using GNU OCTAVE,

octave> A = [ 0.18, 0.60, 0.57, 0.96;
              0.41, 0.24, 0.99, 0.58; 
              0.14, 0.30, 0.97, 0.66; 
              0.51, 0.13, 0.19, 0.85 ];

octave> x = [ -4.05205; -12.6056; 1.66091; 8.69377];

octave> A * x
ans =

  1.0000
  2.0000
  3.0000
  4.0000

This reproduces the original right-hand side vector, @math{b}, in accordance with the equation @math{A x = b}.