dev/cpp/lu_solve/lu_solve.hpp

62 lines
1.8 KiB
C++
Raw Normal View History

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