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.