## 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}.