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

int main()
{
	typedef double value_type;
    typedef boost::numeric::ublas::matrix<value_type> matrix_type;
    typedef boost::numeric::ublas::permutation_matrix<value_type> pmatrix_type;

    matrix_type A(4,4);

    A(0,0) = 0.555950; A(0,1) = 0.274690; A(0,2) = 0.540605; A(0,3) = 0.798938;
    A(1,0) = 0.108929; A(1,1) = 0.830123; A(1,2) = 0.891726; A(1,3) = 0.895283;
    A(2,0) = 0.948014; A(2,1) = 0.973234; A(2,2) = 0.216504; A(2,3) = 0.883152;
    A(3,0) = 0.023787; A(3,1) = 0.675382; A(3,2) = 0.231751; A(3,3) = 0.450332;

	int res;

	// LU without pivoting
	matrix_type A1(A);

	res = boost::numeric::ublas::lu_factorize(A1);

	// LU with partial pivoting
	matrix_type A2(A);

	pmatrix_type P(A2.size1());

	res = boost::numeric::ublas::lu_factorize(A2, P);

	std::cout << "Input matrix" << std::endl;
	std::cout << "\tmatrix: " << A << std::endl;
	std::cout << "LU *without* permutation matrix" << std::endl;
	std::cout << "\tres: " << res << std::endl;
	std::cout << "\tmatrix: " << A1 << std::endl;
	std::cout << "LU *with* permutation matrix" << std::endl;
	std::cout << "\tres: " << res << std::endl;
	std::cout << "\tmatrix: " << A2 << std::endl;
}