Last Update: Feb. 9, 2006
#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;
}
$ ./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
#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;
}
$ ./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