1 // (C) Copyright John Maddock 2006.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 #ifndef BOOST_MATH_TOOLS_SOLVE_HPP
7 #define BOOST_MATH_TOOLS_SOLVE_HPP
8
9 #ifdef _MSC_VER
10 #pragma once
11 #endif
12
13 #include <boost/config.hpp>
14 #include <boost/assert.hpp>
15
16 #ifdef BOOST_MSVC
17 #pragma warning(push)
18 #pragma warning(disable:4996 4267 4244)
19 #endif
20
21 #include <boost/numeric/ublas/lu.hpp>
22 #include <boost/numeric/ublas/matrix.hpp>
23 #include <boost/numeric/ublas/vector.hpp>
24
25 #ifdef BOOST_MSVC
26 #pragma warning(pop)
27 #endif
28
29 namespace boost{ namespace math{ namespace tools{
30
31 //
32 // Find x such that Ax = b
33 //
34 // Caution: this uses undocumented, and untested ublas code,
35 // however short of writing our own LU-decomposition code
36 // it's the only game in town.
37 //
38 template <class T>
solve(const boost::numeric::ublas::matrix<T> & A_,const boost::numeric::ublas::vector<T> & b_)39 boost::numeric::ublas::vector<T> solve(
40 const boost::numeric::ublas::matrix<T>& A_,
41 const boost::numeric::ublas::vector<T>& b_)
42 {
43 //BOOST_ASSERT(A_.size() == b_.size());
44
45 boost::numeric::ublas::matrix<T> A(A_);
46 boost::numeric::ublas::vector<T> b(b_);
47 boost::numeric::ublas::permutation_matrix<> piv(b.size());
48 lu_factorize(A, piv);
49 lu_substitute(A, piv, b);
50 //
51 // iterate to reduce error:
52 //
53 boost::numeric::ublas::vector<T> delta(b.size());
54 for(unsigned k = 0; k < 1; ++k)
55 {
56 noalias(delta) = prod(A_, b);
57 delta -= b_;
58 lu_substitute(A, piv, delta);
59 b -= delta;
60
61 T max_error = 0;
62
63 for(unsigned i = 0; i < delta.size(); ++i)
64 {
65 T err = fabs(delta[i] / b[i]);
66 if(err > max_error)
67 max_error = err;
68 }
69 //std::cout << "Max change in LU error correction: " << max_error << std::endl;
70 }
71
72 return b;
73 }
74
75 }}} // namespaces
76
77 #endif // BOOST_MATH_TOOLS_SOLVE_HPP
78
79
80