From 35f7829af10c61e33dd2e2a7a015058e11a11ea0 Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sat, 25 Mar 2017 14:17:07 +0100 Subject: update --- .../IterativeLinearSolvers/BasicPreconditioners.h | 114 ++++++++++++++++----- 1 file changed, 88 insertions(+), 26 deletions(-) (limited to 'eigen/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h') diff --git a/eigen/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h b/eigen/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h index 1f3c060..358444a 100644 --- a/eigen/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h +++ b/eigen/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 +// Copyright (C) 2011-2014 Gael Guennebaud // // 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 class DiagonalPreconditioner { typedef _Scalar Scalar; typedef Matrix 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 MatrixType; + typedef typename Vector::StorageIndex StorageIndex; + enum { + ColsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic + }; DiagonalPreconditioner() : m_isInitialized(false) {} template - DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols()) + explicit DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols()) { compute(mat); } @@ -80,46 +84,102 @@ class DiagonalPreconditioner return factorize(mat); } + /** \internal */ template - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const Rhs& b, Dest& x) const { x = m_invdiag.array() * b.array() ; } - template inline const internal::solve_retval + template inline const Solve solve(const MatrixBase& 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(*this, b.derived()); + return Solve(*this, b.derived()); } + + ComputationInfo info() { return Success; } protected: Vector m_invdiag; bool m_isInitialized; }; -namespace internal { - -template -struct solve_retval, Rhs> - : solve_retval_base, 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 +class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar> { - typedef DiagonalPreconditioner<_MatrixType> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) + typedef _Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef DiagonalPreconditioner<_Scalar> Base; + using Base::m_invdiag; + public: - template void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; + LeastSquareDiagonalPreconditioner() : Base() {} + + template + explicit LeastSquareDiagonalPreconditioner(const MatType& mat) : Base() + { + compute(mat); + } + + template + LeastSquareDiagonalPreconditioner& analyzePattern(const MatType& ) + { + return *this; + } + + template + LeastSquareDiagonalPreconditioner& factorize(const MatType& mat) + { + // Compute the inverse squared-norm of each column of mat + m_invdiag.resize(mat.cols()); + for(Index j=0; j0) + m_invdiag(j) = RealScalar(1)/sum; + else + m_invdiag(j) = RealScalar(1); + } + Base::m_isInitialized = true; + return *this; + } + + template + 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 - IdentityPreconditioner(const MatrixType& ) {} + explicit IdentityPreconditioner(const MatrixType& ) {} template IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; } @@ -142,6 +202,8 @@ class IdentityPreconditioner template inline const Rhs& solve(const Rhs& b) const { return b; } + + ComputationInfo info() { return Success; } }; } // end namespace Eigen -- cgit v1.2.3