devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

add a bi conjugate gradient stabilized solver

Browse Source
This commit is contained in:
Gael Guennebaud
2011-09-17 10:54:14 +02:00
parent f4122e9f94
commit 9053729d68
5 changed files with 489 additions and 118 deletions
Download Patch File Download Diff File Expand all files Collapse all files

144
unsupported/Eigen/src/IterativeSolvers/ConjugateGradient.h
View File

@@ -83,11 +83,22 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
}
template< typename _MatrixType, int _UpLo=Lower,
typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
class ConjugateGradient;
namespace internal {
template<typename CG, typename Rhs, typename Guess>
class conjugate_gradient_solve_retval_with_guess;
template< typename _MatrixType, int _UpLo, typename _Preconditioner>
struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
{
typedef _MatrixType MatrixType;
typedef _Preconditioner Preconditioner;
};
}
/** \brief A conjugate gradient solver for sparse self-adjoint problems
@@ -137,10 +148,15 @@ class conjugate_gradient_solve_retval_with_guess;
*
* \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
*/
template< typename _MatrixType, int _UpLo=Lower,
typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
class ConjugateGradient
template< typename _MatrixType, int _UpLo, typename _Preconditioner>
class ConjugateGradient : public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
{
typedef IterativeSolverBase<ConjugateGradient> Base;
using Base::mp_matrix;
using Base::m_error;
using Base::m_iterations;
using Base::m_info;
using Base::m_isInitialized;
public:
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
@@ -155,11 +171,7 @@ public:
public:
/** Default constructor. */
ConjugateGradient()
: mp_matrix(0)
{
init();
}
ConjugateGradient() : Base() {}
/** Initialize the solver with matrix \a A for further \c Ax=b solving.
*
@@ -171,90 +183,10 @@ public:
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
ConjugateGradient(const MatrixType& A)
{
init();
compute(A);
}
ConjugateGradient(const MatrixType& A) : Base(A) {}
~ConjugateGradient() {}
/** Initializes the iterative solver with the matrix \a A for further solving \c Ax=b problems.
*
* Currently, this function mostly initialized/compute the preconditioner. In the future
* we might, for instance, implement column reodering for faster matrix vector products.
*
* \warning this class stores a reference to the matrix A as well as some
* precomputed values that depend on it. Therefore, if \a A is changed
* this class becomes invalid. Call compute() to update it with the new
* matrix A, or modify a copy of A.
*/
ConjugateGradient& compute(const MatrixType& A)
{
mp_matrix = &A;
m_preconditioner.compute(A);
m_isInitialized = true;
return *this;
}
/** \internal */
Index rows() const { return mp_matrix->rows(); }
/** \internal */
Index cols() const { return mp_matrix->cols(); }
/** \returns the tolerance threshold used by the stopping criteria */
RealScalar tolerance() const { return m_tolerance; }
/** Sets the tolerance threshold used by the stopping criteria */
ConjugateGradient& setTolerance(RealScalar tolerance)
{
m_tolerance = tolerance;
return *this;
}
/** \returns a read-write reference to the preconditioner for custom configuration. */
Preconditioner& preconditioner() { return m_preconditioner; }
/** \returns a read-only reference to the preconditioner. */
const Preconditioner& preconditioner() const { return m_preconditioner; }
/** \returns the max number of iterations */
int maxIterations() const { return m_maxIterations; }
/** Sets the max number of iterations */
ConjugateGradient& setMaxIterations(int maxIters)
{
m_maxIterations = maxIters;
return *this;
}
/** \returns the number of iterations performed during the last solve */
int iterations() const
{
eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
return m_iterations;
}
/** \returns the tolerance error reached during the last solve */
RealScalar error() const
{
eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
return m_error;
}
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
*
* \sa compute()
*/
template<typename Rhs> inline const internal::solve_retval<ConjugateGradient, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
eigen_assert(rows()==b.rows()
&& "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
return internal::solve_retval<ConjugateGradient, Rhs>(*this, b.derived());
}
/** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
* \a x0 as an initial solution.
*
@@ -265,50 +197,28 @@ public:
solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
{
eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
eigen_assert(rows()==b.rows()
eigen_assert(Base::rows()==b.rows()
&& "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
return internal::conjugate_gradient_solve_retval_with_guess
<ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0);
}
/** \returns Success if the iterations converged, and NoConvergence otherwise. */
ComputationInfo info() const
{
eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
return m_info;
}
/** \internal */
template<typename Rhs,typename Dest>
void _solve(const Rhs& b, Dest& x) const
{
m_iterations = m_maxIterations;
m_error = m_tolerance;
m_iterations = Base::m_maxIterations;
m_error = Base::m_tolerance;
internal::conjugate_gradient(mp_matrix->template selfadjointView<UpLo>(), b, x,
m_preconditioner, m_iterations, m_error);
Base::m_preconditioner, m_iterations, m_error);
m_isInitialized = true;
m_info = m_error <= m_tolerance ? Success : NoConvergence;
m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
}
protected:
void init()
{
m_isInitialized = false;
m_maxIterations = 1000;
m_tolerance = NumTraits<Scalar>::epsilon();
}
const MatrixType* mp_matrix;
Preconditioner m_preconditioner;
int m_maxIterations;
RealScalar m_tolerance;
mutable RealScalar m_error;
mutable int m_iterations;
mutable ComputationInfo m_info;
mutable bool m_isInitialized;
};