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.