- publishing free software manuals
 GNU Scientific Library Reference Manual - Third Edition (v1.12) by M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, P. Alken, M. Booth, F. RossiPaperback (6"x9"), 592 pages, 60 figuresISBN 0954612078RRP £24.95 (\$39.95)

## 13.15 Examples

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

```[ 0.18 0.60 0.57 0.96 ] [x0]   [1.0]
[ 0.41 0.24 0.99 0.58 ] [x1] = [2.0]
[ 0.14 0.30 0.97 0.66 ] [x2]   [3.0]
[ 0.51 0.13 0.19 0.85 ] [x3]   [4.0]
```

and the solution is found using LU decomposition of the matrix 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);
gsl_vector_free (x);
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 x by the original matrix 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, b, in accordance with the equation A x = b.

 ISBN 0954612078 GNU Scientific Library Reference Manual - Third Edition (v1.12) See the print edition