# GSL sample: Linear System of Equations

Last Update: Oct. 4, 2007

### Sample code

linear_system.c

```/***********************************************/
/* T.Kouya's GSL sample program collection     */
/*          Solving Linear System of Equations */
/*                            LU decomposition */
/*                   Written by Tomonori Kouya */
/*                                             */
/* Version 0.01: 2007-10-02                    */
/***********************************************/
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

/* Dimension of Matrix and Vectors */
#define DIM 5

int main(void)
{
int i, j, lotkin_signum, frank_signum;
gsl_matrix *lotkin_a, *frank_a;
gsl_vector *x, *lotkin_b, *frank_b, *lotkin_x, *frank_x;
gsl_permutation *lotkin_perm, *frank_perm;

/* allocate a, x, b */
lotkin_a = gsl_matrix_alloc(DIM, DIM);
frank_a = gsl_matrix_alloc(DIM, DIM);
x = gsl_vector_alloc(DIM);
lotkin_b = gsl_vector_alloc(DIM);
frank_b = gsl_vector_alloc(DIM);
lotkin_x = gsl_vector_alloc(DIM);
frank_x = gsl_vector_alloc(DIM);

/* set x = [1 2 ... DIM] */
for(i = 0; i < DIM; i++)
gsl_vector_set(x, i, (double)i);

/* set Lotkin matrix		      */
/* a_ij = 1 (i = 1) or 1/(i+j-1) (i != 1) */
for(i = 0; i < DIM; i++)
gsl_matrix_set(lotkin_a, 0, i, 1.0);
for(i = 1; i < DIM; i++)
for(j = 0; j < DIM; j++)
gsl_matrix_set(lotkin_a, i, j, 1.0 / (double)(i + j + 1));

/* set Frank matrix       */
/* a_ij = DIM - min(i,j) + 1 */
for(i = 0; i < DIM; i++)
for(j = 0; j < DIM; j++)
gsl_matrix_set(frank_a, i, j, (double)DIM - (double)GSL_MAX(i, j) );

/* Print matrix */
printf("Lotkin Matrix(DIM = %d)\n", DIM);
for(i = 0; i < DIM; i++)
{
printf("%3d: ", i);
for(j = 0; j < DIM; j++)
printf("%g ", gsl_matrix_get(lotkin_a, i, j));
printf("\n");
}
printf("\n");
printf("Frank Matrix(DIM = %d)\n", DIM);
for(i = 0; i < DIM; i++)
{
printf("%3d: ", i);
for(j = 0; j < DIM; j++)
printf("%g ", gsl_matrix_get(frank_a, i, j));
printf("\n");
}
printf("\n");

/* b = A * x */
gsl_blas_dgemv(CblasNoTrans, 1.0, lotkin_a, x, 0.0, lotkin_b);
gsl_blas_dgemv(CblasNoTrans, 1.0, frank_a, x, 0.0, frank_b);

/* LU decomposition and forward&backward substition */
lotkin_perm = gsl_permutation_alloc(DIM);
gsl_linalg_LU_decomp(lotkin_a, lotkin_perm, &lotkin_signum);
gsl_linalg_LU_solve(lotkin_a, lotkin_perm, lotkin_b, lotkin_x);
gsl_permutation_free(lotkin_perm);

frank_perm = gsl_permutation_alloc(DIM);
gsl_linalg_LU_decomp(frank_a, frank_perm, &frank_signum);
gsl_linalg_LU_solve(frank_a, frank_perm, frank_b, frank_x);
gsl_permutation_free(frank_perm);

/* print */
printf("index	  Lotkin		  Frank			 True\n");
for(i = 0; i < DIM; i++)
printf("%3d %25.17e %25.17e %25.17e\n", i, gsl_vector_get(lotkin_x, i), gsl_vector_get(frank_x, i), gsl_vector_get(x, i));

/* free a, x, b */
gsl_matrix_free(lotkin_a);
gsl_matrix_free(frank_a);
gsl_vector_free(x);
gsl_vector_free(lotkin_b);
gsl_vector_free(frank_b);
gsl_vector_free(lotkin_x);
gsl_vector_free(frank_x);

return 0;
}```

### Example execution

```\$ ./linear_system
Lotkin Matrix(DIM = 5)
0: 1 1 1 1 1
1: 0.5 0.333333 0.25 0.2 0.166667
2: 0.333333 0.25 0.2 0.166667 0.142857
3: 0.25 0.2 0.166667 0.142857 0.125
4: 0.2 0.166667 0.142857 0.125 0.111111

Frank Matrix(DIM = 5)
0: 5 4 3 2 1
1: 4 4 3 2 1
2: 3 3 3 2 1
3: 2 2 2 2 1
4: 1 1 1 1 1

index          Lotkin                  Frank                         True
0  -2.82440737464639824e-13  -1.42108547152020045e-15   0.00000000000000000e+00
1   1.00000000000282596e+00   1.00000000000000355e+00   1.00000000000000000e+00
2   1.99999999999155187e+00   1.99999999999999645e+00   2.00000000000000000e+00
3   3.00000000000982148e+00   3.00000000000000178e+00   3.00000000000000000e+00
4   3.99999999999608313e+00   4.00000000000000000e+00   4.00000000000000000e+00```

<- Go back