#ifndef _KLMN_LU_SOLVE_
#define _KLMN_LU_SOLVE_

#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/blas.hpp>
#include <boost/numeric/ublas/lu.hpp>

// boost: разворот вектора на месте (y[0]<->y[n], y[1]<->y[n-1],...)
template <typename T>
void inplace_turn(boost::numeric::ublas::vector<T>& v)
{
	for(size_t i=0, max_i=v.size()/2, j=v.size()-1; i<max_i; i++,j--)
	{
		T tmp = v(i);
		v(i) = v(j);
		v(j) = tmp;
	}
}

// boost: решение системы A*x = b методом LU-разложения за O(n) шагов
// матрица A не сохраняется, решение в b (для оптимизации по скорости и памяти)
// алгоритм требует дополнительную память при копировании A в L и U,
// размер кт. равен памяти, выделяемой под матрицу A.
template <typename T>
void lu_solve (boost::numeric::ublas::matrix<T>& A, boost::numeric::ublas::vector<T>& b)
{
	using namespace boost::numeric::ublas;

	// получаем L и U на месте A
	lu_factorize(A);

	// "вытаскиваем" L из A
	size_t N = A.size1();
	boost::numeric::ublas::triangular_matrix<T> L(N,N);
	for(size_t i=1; i<N; i++)
		for(size_t j=0; j<i; j++)
			L(i,j) = A(i,j);
	for(size_t i=0; i<N; i++)
		L(i,i) = 1;

	// "вытаскиваем" U из A
	boost::numeric::ublas::triangular_matrix<T> U(N,N);
	for(size_t i=0; i<N; i++)
		for(size_t j=i; j<N; j++)
			U(N-1-i,N-1-j) = A(i,j);

	// освобождаем память из-под A
	A.resize(0,0);

	// решаем Ly=Pb
	inplace_solve(L,b,lower_tag());
	inplace_turn(b);
	L.resize(0,0);

	// решаем Ux=y
	inplace_solve(U,b,lower_tag());
	inplace_turn(b);
}

#endif // _KLMN_LU_SOLVE_