| // (C) Copyright John Maddock 2006. |
| // Use, modification and distribution are subject to the |
| // Boost Software License, Version 1.0. (See accompanying file |
| // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
| |
| #ifndef BOOST_MATH_TOOLS_SOLVE_HPP |
| #define BOOST_MATH_TOOLS_SOLVE_HPP |
| |
| #ifdef _MSC_VER |
| #pragma once |
| #endif |
| |
| #include <boost/config.hpp> |
| #include <boost/assert.hpp> |
| |
| #ifdef BOOST_MSVC |
| #pragma warning(push) |
| #pragma warning(disable:4996 4267 4244) |
| #endif |
| |
| #include <boost/numeric/ublas/lu.hpp> |
| #include <boost/numeric/ublas/matrix.hpp> |
| #include <boost/numeric/ublas/vector.hpp> |
| |
| #ifdef BOOST_MSVC |
| #pragma warning(pop) |
| #endif |
| |
| namespace boost{ namespace math{ namespace tools{ |
| |
| // |
| // Find x such that Ax = b |
| // |
| // Caution: this uses undocumented, and untested ublas code, |
| // however short of writing our own LU-decompostion code |
| // it's the only game in town. |
| // |
| template <class T> |
| boost::numeric::ublas::vector<T> solve( |
| const boost::numeric::ublas::matrix<T>& A_, |
| const boost::numeric::ublas::vector<T>& b_) |
| { |
| //BOOST_ASSERT(A_.size() == b_.size()); |
| |
| boost::numeric::ublas::matrix<T> A(A_); |
| boost::numeric::ublas::vector<T> b(b_); |
| boost::numeric::ublas::permutation_matrix<> piv(b.size()); |
| lu_factorize(A, piv); |
| lu_substitute(A, piv, b); |
| // |
| // iterate to reduce error: |
| // |
| boost::numeric::ublas::vector<T> delta(b.size()); |
| for(unsigned i = 0; i < 1; ++i) |
| { |
| noalias(delta) = prod(A_, b); |
| delta -= b_; |
| lu_substitute(A, piv, delta); |
| b -= delta; |
| |
| T max_error = 0; |
| |
| for(unsigned i = 0; i < delta.size(); ++i) |
| { |
| T err = fabs(delta[i] / b[i]); |
| if(err > max_error) |
| max_error = err; |
| } |
| //std::cout << "Max change in LU error correction: " << max_error << std::endl; |
| } |
| |
| return b; |
| } |
| |
| }}} // namespaces |
| |
| #endif // BOOST_MATH_TOOLS_SOLVE_HPP |
| |
| |