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

using namespace std;
using namespace boost::numeric::ublas;

void PrintMatrix(const string &tag, const matrix<float> &M);

int main()
{
	matrix<float>  A(2, 2);
	
	A(0,0) = 4; A(1,0) = 3;
	A(0,1) = 6; A(1,1) = 3;
	
	PrintMatrix("A", A);
	
	permutation_matrix<> pm(A.size1());
	
	int r = lu_factorize< matrix<float> >(A);
//	int r = lu_factorize< matrix<float>, permutation_matrix<> >(A, pm);
	
	PrintMatrix("lu_factorize A", A);
	
	matrix<float>  L(2, 2), U(2, 2);
	
	L(0,0) = 1;      L(1,0) = 0;
	L(0,1) = A(0,1); L(1,1) = 1;
	U(0,0) = A(0,0); U(1,0) = A(1,0);
	U(0,1) = 0;      U(1,1) = A(1,1);
	
	PrintMatrix("L", L);
	PrintMatrix("U", U);
	
	// L and U are unique if I impose the diagonal elements of L to be equal to 1.
	
	matrix<float>  inverse(A.size1(), A.size2());
	
	inverse.assign(identity_matrix<float>(A.size1()));
	
	lu_substitute(A, pm, inverse);
	
	PrintMatrix("inverse", inverse);
}


void PrintMatrix(const string &tag, const matrix<float> &M)
{
	cout << tag << "[i,j]:  " << M << endl;
	for (int j = 0; j < M.size2(); ++j) {
		for (int i = 0; i < M.size1(); ++i)
			cout << "\t[i=" << i << ",j=" << j << "] " << M(i, j) << "\t";
		cout << endl;
	}
	cout << endl;
}