81 lines · plain
1// (C) Copyright John Maddock 2006.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_TOOLS_SOLVE_HPP7#define BOOST_MATH_TOOLS_SOLVE_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/config.hpp>14#include <boost/math/tools/assert.hpp>15 16#ifdef _MSC_VER17#pragma warning(push)18#pragma warning(disable:4996 4267 4244)19#endif20 21#include <boost/numeric/ublas/lu.hpp>22#include <boost/numeric/ublas/matrix.hpp>23#include <boost/numeric/ublas/vector.hpp>24 25#ifdef _MSC_VER26#pragma warning(pop)27#endif28 29namespace boost{ namespace math{ namespace tools{30 31//32// Find x such that Ax = b33//34// Caution: this uses undocumented, and untested ublas code,35// however short of writing our own LU-decomposition code36// it's the only game in town.37//38template <class T>39boost::numeric::ublas::vector<T> solve(40 const boost::numeric::ublas::matrix<T>& A_,41 const boost::numeric::ublas::vector<T>& b_)42{43 //BOOST_MATH_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 using std::abs;66 T err = abs(delta[i] / b[i]);67 if(err > max_error)68 max_error = err;69 }70 //std::cout << "Max change in LU error correction: " << max_error << std::endl;71 }72 73 return b;74}75 76}}} // namespaces77 78#endif // BOOST_MATH_TOOLS_SOLVE_HPP79 80 81