## File: Linear-Algebra-Examples.html

package info (click to toggle)
gsl-ref-html 2.3-1
• area: non-free
• in suites: bullseye, buster, sid
• size: 6,876 kB
• ctags: 4,574
• sloc: makefile: 35
 file content (160 lines) | stat: -rw-r--r-- 5,653 bytes parent folder | download
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160` `````` GNU Scientific Library – Reference Manual: Linear Algebra Examples

14.20 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.

``````