Linear System in C++

Linear System Solutions in C++ using Boost::ublas

#include < boost/numeric/ublas/matrix.hpp > #include < boost/numeric/ublas/lu.hpp > #include < boost/numeric/ublas/io.hpp > using namespace boost::numeric::ublas; using namespace std; typedef boost::numeric::ublas::vector< double > ublasVector; typedef boost::numeric::ublas::matrix< double > ublasMatrix; typedef boost::numeric::ublas::matrix_row< ublasMatrix > ublasMatrixRow; typedef boost::numeric::ublas::matrix_range< matrix < double > > ublasMatrixRange;

Gaussian Elimination Method

// Ax=b ==>> transform to Ux=b ublasVector GaussianElimination(ublasMatrix um, ublasVector uv) { std::cout << uv << std::endl; std::cout << um << std::endl; // eliminate and create upper triangle augmented matrix for(int i=0; i<(uv.size()-1); ++i){ for(int j=i+1; j< uv.size(); ++j){ double lem=um(j,i)/um(i,i); ublasMatrixRow(um,j) -= lem*ublasMatrixRow(um,i); uv(j)-=lem*uv(i); } } // back substitution to get result for(int i=uv.size()-1; i>=0; --i){ double in=inner_prod(subrange(uv,i+1,uv.size()),subrange(ublasMatrixRow(um,i),i+1,uv.size())); uv(i)=(uv(i)-in)/um(i,i); } ublasVector retuv(uv); return retuv; }

LU Decomposition

// decompose Ax=b ==>> LUx=b hance tranform A=LU void LUDecomposition(ublasMatrix& um) { std::cout << um << std::endl; for(int i=0; i< (um.size1()-1); ++i){ for(int j=i+1; j< um.size1(); ++j){ double lem=um(j,i)/um(i,i); ublasMatrixRange(um,range(j,j+1),range(i+1,um.size1())) -= lem*ublasMatrixRange(um,range(i,i+1),range(i+1,um.size1())); um(j,i)=lem; } } }

LU Solver using LUDecomposition

// This use LUdecompose matrix to resolve liner equation ublasVector LUSolver(ublasMatrix um, ublasVector uv) { // forward substitution Ly=b to get solution for y for(int i=1; i< uv.size(); ++i){ double in=inner_prod(subrange(ublasMatrixRow(um,i),0,i),subrange(uv,0,i)); uv(i)-=in; } // back substitution to get result Ux=y for(int i=uv.size()-1; i>=0; --i){ double in=inner_prod(subrange(uv,i+1,uv.size()),subrange(ublasMatrixRow(um,i),i+1,uv.size())); uv(i)=(uv(i)-in)/um(i,i); } ublasVector retuv(uv); return retuv; }

Test Programs

void testLUDecompose() { ublasVector uv (3); uv (0) = 6;uv (1) = 3;uv (2) = 7; ublasMatrix um(3,3); um(0,0)=3;um(0,1)=-1;um(0,2)=4; um(1,0)=-2;um(1,1)=0;um(1,2)=5; um(2,0)=7;um(2,1)=2;um(2,2)=-2; std::cout << um << std::endl; LUDecomposition(um); std::cout << um << std::endl; // calculate determinante |A| double dt=1; for(int i=0; i< um.size1();++i) dt*=um(i,i); std::cout << dt << std::endl; ublasVector retuv=LUSolver(um,uv); std::cout << retuv << std::endl; }

Choleski Decomposition Method

// decompose Ax=b ==>> L*trans(L)x=b void CholeskiDecomposition(ublasMatrix& um) { std::cout << um << std::endl; for(int i=0; i< um.size1(); ++i){ try{ ublasVector uv=subrange(ublasMatrixRow(um,i),0,i); double in = inner_prod(uv,uv); if(in < um(i,i)) um(i,i)=sqrt(um(i,i)-inner_prod(uv,uv)); else cout<< "Matirx is not positive definite" << endl; // throw exception here for(int j=i+1; j< um.size1(); ++j){ double in = inner_prod(subrange(ublasMatrixRow(um,i),0,i),subrange(ublasMatrixRow(um,j),0,i)); um(j,i) = (um(j,i) - in)/um(i,i); } } for(int i=1; i< um.size1(); ++i) { ublasMatrixRange(um,range(0,i),range(i,i+1))=zero_matrix(i,1); } }

void testCholeski() { ublasMatrix um1(3,3); um1(0,0)=4;um1(0,1)=-2;um1(0,2)=1; um1(1,0)=-2;um1(1,1)=4;um1(1,2)=-2; um1(2,0)=1;um1(2,1)=-2;um1(2,2)=4; std::cout << um1 << std::endl; CholeskiDecomposition(um1); std::cout << um1 << std::endl; ublasMatrix umt=trans(um1); std::cout << umt << std::endl; std::cout << prod(um1,umt) << std::endl; }