62 lines
1.8 KiB
C++
62 lines
1.8 KiB
C++
#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_
|
|
|