24 #ifndef _LINEARSOLVER_HPP_ 25 #define _LINEARSOLVER_HPP_ 29 #include <boost/numeric/ublas/vector.hpp> 30 #include <boost/numeric/ublas/matrix.hpp> 31 #include <boost/numeric/ublas/lu.hpp> 32 #include <boost/numeric/ublas/triangular.hpp> 33 #include <boost/numeric/ublas/vector_proxy.hpp> 34 #include <boost/numeric/ublas/io.hpp> 50 template <LinearSolverType solver_type = LinearSolverType_LU>
61 template<
typename matrix_type,
typename vector_type>
62 vector_type
solve(
const matrix_type&
A,
67 matrix_type A_factorized =
A;
68 ublas::permutation_matrix<size_t> pm(y.size());
70 int singular = lu_factorize(A_factorized, pm);
71 if (singular)
throw std::runtime_error(
"[LinearSolver<LU>::solve()] A is singular.");
73 vector_type result(y);
74 lu_substitute(A_factorized, pm, result);
88 template<
typename matrix_type,
typename vector_type>
89 vector_type
solve(
const matrix_type&
A,
const vector_type&
y)
91 typedef typename matrix_type::size_type size_type;
92 typedef typename matrix_type::value_type value_type;
96 matrix_type Q(A.size1(), A.size2()), R(A.size1(), A.size2());
100 vector_type b = prod(trans(Q), y);
103 if (R.size1() > R.size2())
105 size_type min = (R.size1() < R.size2() ? R.size1() : R.size2());
107 result = ublas::solve(subrange(R, 0, min, 0, min),
113 result = ublas::solve(R, b, ublas::upper_tag());
123 #endif // _LINEARSOLVER_HPP_
vector_type solve(const matrix_type &A, const vector_type &y)
solve system of linear equations Ax = y using boost::ublas; note: extra copying inefficiencies for ea...
void qr(const matrix_type &A, matrix_type &Q, matrix_type &R)
vector_type solve(const matrix_type &A, const vector_type &y)
solve system of linear equations Ax = y using boost::ublas; note: extra copying inefficiencies for ea...
KernelTraitsBase< Kernel >::space_type::ordinate_type y