80 lines
1.9 KiB
C++
80 lines
1.9 KiB
C++
// (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
|
|
|
|
|