#ifndef MATRIX_INVERSION_HPP
#define MATRIX_INVERSION_HPP

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

namespace ub = boost::numeric::ublas;

template <class T>
bool invert(const ub::matrix<T>& a, ub::matrix<T>& b) {
	ub::matrix<T> tmp(a);
	ub::permutation_matrix<> pm(tmp.size1());

	if (ub::lu_factorize(tmp, pm) != 0) return false;

	b = ub::identity_matrix<T>(tmp.size1());

	ub::lu_substitute(tmp, pm, b);

	return true;
}

#endif // MATRIX_INVERSION_HPP