#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_