diff options
Diffstat (limited to 'Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h')
| -rw-r--r-- | Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h | 118 |
1 files changed, 90 insertions, 28 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h index 73ca9bfde..358444aff 100644 --- a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h +++ b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -17,33 +17,37 @@ namespace Eigen { * * This class allows to approximately solve for A.x = b problems assuming A is a diagonal matrix. * In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for: - * \code - * A.diagonal().asDiagonal() . x = b - * \endcode + \code + A.diagonal().asDiagonal() . x = b + \endcode * * \tparam _Scalar the type of the scalar. * + * \implsparsesolverconcept + * * This preconditioner is suitable for both selfadjoint and general problems. * The diagonal entries are pre-inverted and stored into a dense vector. * * \note A variant that has yet to be implemented would attempt to preserve the norm of each column. * + * \sa class LeastSquareDiagonalPreconditioner, class ConjugateGradient */ template <typename _Scalar> class DiagonalPreconditioner { typedef _Scalar Scalar; typedef Matrix<Scalar,Dynamic,1> Vector; - typedef typename Vector::Index Index; - public: - // this typedef is only to export the scalar type and compile-time dimensions to solve_retval - typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; + typedef typename Vector::StorageIndex StorageIndex; + enum { + ColsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic + }; DiagonalPreconditioner() : m_isInitialized(false) {} template<typename MatType> - DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols()) + explicit DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols()) { compute(mat); } @@ -65,10 +69,10 @@ class DiagonalPreconditioner { typename MatType::InnerIterator it(mat,j); while(it && it.index()!=j) ++it; - if(it && it.index()==j) + if(it && it.index()==j && it.value()!=Scalar(0)) m_invdiag(j) = Scalar(1)/it.value(); else - m_invdiag(j) = 0; + m_invdiag(j) = Scalar(1); } m_isInitialized = true; return *this; @@ -80,46 +84,102 @@ class DiagonalPreconditioner return factorize(mat); } + /** \internal */ template<typename Rhs, typename Dest> - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const Rhs& b, Dest& x) const { x = m_invdiag.array() * b.array() ; } - template<typename Rhs> inline const internal::solve_retval<DiagonalPreconditioner, Rhs> + template<typename Rhs> inline const Solve<DiagonalPreconditioner, Rhs> solve(const MatrixBase<Rhs>& b) const { eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized."); eigen_assert(m_invdiag.size()==b.rows() && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval<DiagonalPreconditioner, Rhs>(*this, b.derived()); + return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived()); } + + ComputationInfo info() { return Success; } protected: Vector m_invdiag; bool m_isInitialized; }; -namespace internal { - -template<typename _MatrixType, typename Rhs> -struct solve_retval<DiagonalPreconditioner<_MatrixType>, Rhs> - : solve_retval_base<DiagonalPreconditioner<_MatrixType>, Rhs> +/** \ingroup IterativeLinearSolvers_Module + * \brief Jacobi preconditioner for LeastSquaresConjugateGradient + * + * This class allows to approximately solve for A' A x = A' b problems assuming A' A is a diagonal matrix. + * In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for: + \code + (A.adjoint() * A).diagonal().asDiagonal() * x = b + \endcode + * + * \tparam _Scalar the type of the scalar. + * + * \implsparsesolverconcept + * + * The diagonal entries are pre-inverted and stored into a dense vector. + * + * \sa class LeastSquaresConjugateGradient, class DiagonalPreconditioner + */ +template <typename _Scalar> +class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar> { - typedef DiagonalPreconditioner<_MatrixType> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) + typedef _Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef DiagonalPreconditioner<_Scalar> Base; + using Base::m_invdiag; + public: - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; + LeastSquareDiagonalPreconditioner() : Base() {} + + template<typename MatType> + explicit LeastSquareDiagonalPreconditioner(const MatType& mat) : Base() + { + compute(mat); + } + + template<typename MatType> + LeastSquareDiagonalPreconditioner& analyzePattern(const MatType& ) + { + return *this; + } + + template<typename MatType> + LeastSquareDiagonalPreconditioner& factorize(const MatType& mat) + { + // Compute the inverse squared-norm of each column of mat + m_invdiag.resize(mat.cols()); + for(Index j=0; j<mat.outerSize(); ++j) + { + RealScalar sum = mat.innerVector(j).squaredNorm(); + if(sum>0) + m_invdiag(j) = RealScalar(1)/sum; + else + m_invdiag(j) = RealScalar(1); + } + Base::m_isInitialized = true; + return *this; + } + + template<typename MatType> + LeastSquareDiagonalPreconditioner& compute(const MatType& mat) + { + return factorize(mat); + } + + ComputationInfo info() { return Success; } -} + protected: +}; /** \ingroup IterativeLinearSolvers_Module * \brief A naive preconditioner which approximates any matrix as the identity matrix * + * \implsparsesolverconcept + * * \sa class DiagonalPreconditioner */ class IdentityPreconditioner @@ -129,7 +189,7 @@ class IdentityPreconditioner IdentityPreconditioner() {} template<typename MatrixType> - IdentityPreconditioner(const MatrixType& ) {} + explicit IdentityPreconditioner(const MatrixType& ) {} template<typename MatrixType> IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; } @@ -142,6 +202,8 @@ class IdentityPreconditioner template<typename Rhs> inline const Rhs& solve(const Rhs& b) const { return b; } + + ComputationInfo info() { return Success; } }; } // end namespace Eigen |