#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/vector.hpp>

#include <iostream>


namespace ublas = boost::numeric::ublas;

ublas::vector<double> LTRIS(ublas::matrix<double>& A, ublas::vector<double> b)
{
    ublas::vector <double> x = b;

    ublas::permutation_matrix <std::size_t> P(A.size1());
    int r = ublas::lu_factorize(A, P);

    for(std::size_t i = 1; i <= A.size1(); ++i){
        for(std::size_t j = 1; j < i; j++)
            x(i) = x(i) - A(i,j)*x(j);
        x(i) = x(i) / A(i,i);
    }

    return x;
}

ublas::vector<double> UTRIS(ublas::matrix<double>& A, ublas::vector<double> b)
{
    ublas::vector <double> x = b;

    ublas::permutation_matrix <std::size_t> P(A.size1());
    int r = ublas::lu_factorize(A, P);

    for(std::size_t i = A.size1(); i >= 1; --i){
        for(std::size_t j = i + 1; j < A.size1(); j++)
            x(i) = x(i) - A(i,j)*x(j);
        x(i) = x(i) / A(i,i);
    }

    return x;
}

int main ()
{

}