# Solving a linear system of equation with LU decomposition

Last Update: Feb. 9, 2006

## MPFR

### Sample code

lu_mpf.c

#include <stdio.h>
#include <math.h>

/* Header file of BNCpack */
#include "bnc.h"

/* LU decomposion */

/* Main function */
int main()
{
long int ierror;
MPFMatrix dmat;
MPFVector dv, dans;

/* Set default precision in decimal digits */
set_bnc_default_prec_decimal(50);

/* Initialize */
dmat = init_mpfmatrix(3, 3);
dv   = init_mpfvector(3);
dans = init_mpfvector(3);

/* Set elements of matrix */
/* A = [  1  1  1 ] */
/*     [ 20 24 30 ] */
/*     [ 12 10  8 ] */
set_mpfmatrix_ij_ui(dmat, 0, 0, 1);
set_mpfmatrix_ij_ui(dmat, 0, 1, 1);
set_mpfmatrix_ij_ui(dmat, 0, 2, 1);

set_mpfmatrix_ij_ui(dmat, 1, 0, 20);
set_mpfmatrix_ij_ui(dmat, 1, 1, 24);
set_mpfmatrix_ij_ui(dmat, 1, 2, 30);

set_mpfmatrix_ij_ui(dmat, 2, 0, 12);
set_mpfmatrix_ij_ui(dmat, 2, 1, 10);
set_mpfmatrix_ij_ui(dmat, 2, 2,  8);

/* Print matrix */
printf("Matrix:\n");
print_mpfmatrix(dmat);

/* call LU decomposion */
ierror = MPFLUdecomp(
dmat    /* Original matrix and its LU decomposion */
);

/* Print LU decomposition of dmat */
printf("LU decomp (Error status %d): \n", ierror);
print_mpfmatrix(dmat);

/* Set dv */
/* b = [3 74 30]^T */
set_mpfvector_i_ui(dv, 0,  3);
set_mpfvector_i_ui(dv, 1, 74);
set_mpfvector_i_ui(dv, 2, 30);

/* Print vector */
printf("Constant Vector:\n");
print_mpfvector(dv);

/* Solve linear system of equation Ax = b */
SolveMPFLS(dans, dmat, dv);

/* Print solution x */
printf("Solution:\n");
print_mpfvector(dans);

/* Clear */
free_mpfmatrix(dmat);
free_mpfvector(dv);
free_mpfvector(dans);

/* Exit */
return 0;
}

### Example execution

\$ ./lu_mpf
-------------------------------------------------------------------------------
BNC Default Precision    : 167 bits(50.3 decimal digits)
BNC Default Rounding Mode: Round to Nearest
-------------------------------------------------------------------------------
Matrix:
0 1.000000000000000000000000000000000000000000000000000 1.000000000000000000000000000000000000000000000000000 1.000000000000000000000000000000000000000000000000000
1 2.000000000000000000000000000000000000000000000000000e1 2.400000000000000000000000000000000000000000000000000e1 3.000000000000000000000000000000000000000000000000000e1
2 1.200000000000000000000000000000000000000000000000000e1 1.000000000000000000000000000000000000000000000000000e1 8.000000000000000000000000000000000000000000000000000
LU decomp (Error status 0):
0 1.000000000000000000000000000000000000000000000000000 1.000000000000000000000000000000000000000000000000000 1.000000000000000000000000000000000000000000000000000
1 2.000000000000000000000000000000000000000000000000000e1 4.000000000000000000000000000000000000000000000000000 1.000000000000000000000000000000000000000000000000000e1
2 1.200000000000000000000000000000000000000000000000000e1 -5.000000000000000000000000000000000000000000000000000e-1 1.000000000000000000000000000000000000000000000000000
Constant Vector:
0 3.000000000000000000000000000000000000000000000000000
1 7.400000000000000000000000000000000000000000000000000e1
2 3.000000000000000000000000000000000000000000000000000e1
Solution:
0 1.000000000000000000000000000000000000000000000000000
1 1.000000000000000000000000000000000000000000000000000
2 1.000000000000000000000000000000000000000000000000000

## IEEE754 double

### Sample code

lu_d.c

#include <stdio.h>
#include <math.h>

/* Header file of BNCpack */
#include "bnc.h"

/* LU decomposion */

/* Main function */
int main()
{
long int ierror;
DMatrix dmat;
DVector dv, dans;

/* Initialize */
dmat = init_dmatrix(3, 3);
dv   = init_dvector(3);
dans = init_dvector(3);

/* Set elements of matrix */
/* A = [  1  1  1 ] */
/*     [ 20 24 30 ] */
/*     [ 12 10  8 ] */
set_dmatrix_ij(dmat, 0, 0, 1.0);
set_dmatrix_ij(dmat, 0, 1, 1.0);
set_dmatrix_ij(dmat, 0, 2, 1.0);

set_dmatrix_ij(dmat, 1, 0, 20.0);
set_dmatrix_ij(dmat, 1, 1, 24.0);
set_dmatrix_ij(dmat, 1, 2, 30.0);

set_dmatrix_ij(dmat, 2, 0, 12.0);
set_dmatrix_ij(dmat, 2, 1, 10.0);
set_dmatrix_ij(dmat, 2, 2,  8.0);

/* Print matrix */
printf("Matrix:\n");
print_dmatrix(dmat);

/* call LU decomposion */
ierror = DLUdecomp(
dmat    /* Original matrix and its LU decomposion */
);

/* Print LU decomposition of dmat */
printf("LU decomp (Error status %d): \n", ierror);
print_dmatrix(dmat);

/* Set dv */
/* b = [3 74 30]^T */
set_dvector_i(dv, 0,  3.0);
set_dvector_i(dv, 1, 74.0);
set_dvector_i(dv, 2, 30.0);

/* Print vector */
printf("Constant Vector:\n");
print_dvector(dv);

/* Solve linear system of equation Ax = b */
SolveDLS(dans, dmat, dv);

/* Print solution x */
printf("Solution:\n");
print_dvector(dans);

/* Clear */
free_dmatrix(dmat);
free_dvector(dv);
free_dvector(dans);

/* Exit */
return 0;
}

### Example execution

\$ ./lu_d
Matrix:
0   1.00000000000000000e+00   1.00000000000000000e+00   1.00000000000000000e+00
1   2.00000000000000000e+01   2.40000000000000000e+01   3.00000000000000000e+01
2   1.20000000000000000e+01   1.00000000000000000e+01   8.00000000000000000e+00
LU decomp (Error status 0):
0   1.00000000000000000e+00   1.00000000000000000e+00   1.00000000000000000e+00
1   2.00000000000000000e+01   4.00000000000000000e+00   1.00000000000000000e+01
2   1.20000000000000000e+01  -5.00000000000000000e-01   1.00000000000000000e+00
Constant Vector:
0   3.00000000000000000e+00
1   7.40000000000000000e+01
2   3.00000000000000000e+01
Solution:
0   1.00000000000000000e+00
1   1.00000000000000000e+00
2   1.00000000000000000e+00

<- Go back