1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_BASIC_PRECONDITIONERS_H 11 #define EIGEN_BASIC_PRECONDITIONERS_H 12 13 namespace Eigen { 14 15 /** \ingroup IterativeLinearSolvers_Module 16 * \brief A preconditioner based on the digonal entries 17 * 18 * This class allows to approximately solve for A.x = b problems assuming A is a diagonal matrix. 19 * In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for: 20 \code 21 A.diagonal().asDiagonal() . x = b 22 \endcode 23 * 24 * \tparam _Scalar the type of the scalar. 25 * 26 * \implsparsesolverconcept 27 * 28 * This preconditioner is suitable for both selfadjoint and general problems. 29 * The diagonal entries are pre-inverted and stored into a dense vector. 30 * 31 * \note A variant that has yet to be implemented would attempt to preserve the norm of each column. 32 * 33 * \sa class LeastSquareDiagonalPreconditioner, class ConjugateGradient 34 */ 35 template <typename _Scalar> 36 class DiagonalPreconditioner 37 { 38 typedef _Scalar Scalar; 39 typedef Matrix<Scalar,Dynamic,1> Vector; 40 public: 41 typedef typename Vector::StorageIndex StorageIndex; 42 enum { 43 ColsAtCompileTime = Dynamic, 44 MaxColsAtCompileTime = Dynamic 45 }; 46 DiagonalPreconditioner()47 DiagonalPreconditioner() : m_isInitialized(false) {} 48 49 template<typename MatType> DiagonalPreconditioner(const MatType & mat)50 explicit DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols()) 51 { 52 compute(mat); 53 } 54 rows()55 Index rows() const { return m_invdiag.size(); } cols()56 Index cols() const { return m_invdiag.size(); } 57 58 template<typename MatType> analyzePattern(const MatType &)59 DiagonalPreconditioner& analyzePattern(const MatType& ) 60 { 61 return *this; 62 } 63 64 template<typename MatType> factorize(const MatType & mat)65 DiagonalPreconditioner& factorize(const MatType& mat) 66 { 67 m_invdiag.resize(mat.cols()); 68 for(int j=0; j<mat.outerSize(); ++j) 69 { 70 typename MatType::InnerIterator it(mat,j); 71 while(it && it.index()!=j) ++it; 72 if(it && it.index()==j && it.value()!=Scalar(0)) 73 m_invdiag(j) = Scalar(1)/it.value(); 74 else 75 m_invdiag(j) = Scalar(1); 76 } 77 m_isInitialized = true; 78 return *this; 79 } 80 81 template<typename MatType> compute(const MatType & mat)82 DiagonalPreconditioner& compute(const MatType& mat) 83 { 84 return factorize(mat); 85 } 86 87 /** \internal */ 88 template<typename Rhs, typename Dest> _solve_impl(const Rhs & b,Dest & x)89 void _solve_impl(const Rhs& b, Dest& x) const 90 { 91 x = m_invdiag.array() * b.array() ; 92 } 93 94 template<typename Rhs> inline const Solve<DiagonalPreconditioner, Rhs> solve(const MatrixBase<Rhs> & b)95 solve(const MatrixBase<Rhs>& b) const 96 { 97 eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized."); 98 eigen_assert(m_invdiag.size()==b.rows() 99 && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b"); 100 return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived()); 101 } 102 info()103 ComputationInfo info() { return Success; } 104 105 protected: 106 Vector m_invdiag; 107 bool m_isInitialized; 108 }; 109 110 /** \ingroup IterativeLinearSolvers_Module 111 * \brief Jacobi preconditioner for LeastSquaresConjugateGradient 112 * 113 * This class allows to approximately solve for A' A x = A' b problems assuming A' A is a diagonal matrix. 114 * In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for: 115 \code 116 (A.adjoint() * A).diagonal().asDiagonal() * x = b 117 \endcode 118 * 119 * \tparam _Scalar the type of the scalar. 120 * 121 * \implsparsesolverconcept 122 * 123 * The diagonal entries are pre-inverted and stored into a dense vector. 124 * 125 * \sa class LeastSquaresConjugateGradient, class DiagonalPreconditioner 126 */ 127 template <typename _Scalar> 128 class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar> 129 { 130 typedef _Scalar Scalar; 131 typedef typename NumTraits<Scalar>::Real RealScalar; 132 typedef DiagonalPreconditioner<_Scalar> Base; 133 using Base::m_invdiag; 134 public: 135 LeastSquareDiagonalPreconditioner()136 LeastSquareDiagonalPreconditioner() : Base() {} 137 138 template<typename MatType> LeastSquareDiagonalPreconditioner(const MatType & mat)139 explicit LeastSquareDiagonalPreconditioner(const MatType& mat) : Base() 140 { 141 compute(mat); 142 } 143 144 template<typename MatType> analyzePattern(const MatType &)145 LeastSquareDiagonalPreconditioner& analyzePattern(const MatType& ) 146 { 147 return *this; 148 } 149 150 template<typename MatType> factorize(const MatType & mat)151 LeastSquareDiagonalPreconditioner& factorize(const MatType& mat) 152 { 153 // Compute the inverse squared-norm of each column of mat 154 m_invdiag.resize(mat.cols()); 155 for(Index j=0; j<mat.outerSize(); ++j) 156 { 157 RealScalar sum = mat.innerVector(j).squaredNorm(); 158 if(sum>0) 159 m_invdiag(j) = RealScalar(1)/sum; 160 else 161 m_invdiag(j) = RealScalar(1); 162 } 163 Base::m_isInitialized = true; 164 return *this; 165 } 166 167 template<typename MatType> compute(const MatType & mat)168 LeastSquareDiagonalPreconditioner& compute(const MatType& mat) 169 { 170 return factorize(mat); 171 } 172 info()173 ComputationInfo info() { return Success; } 174 175 protected: 176 }; 177 178 /** \ingroup IterativeLinearSolvers_Module 179 * \brief A naive preconditioner which approximates any matrix as the identity matrix 180 * 181 * \implsparsesolverconcept 182 * 183 * \sa class DiagonalPreconditioner 184 */ 185 class IdentityPreconditioner 186 { 187 public: 188 IdentityPreconditioner()189 IdentityPreconditioner() {} 190 191 template<typename MatrixType> IdentityPreconditioner(const MatrixType &)192 explicit IdentityPreconditioner(const MatrixType& ) {} 193 194 template<typename MatrixType> analyzePattern(const MatrixType &)195 IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; } 196 197 template<typename MatrixType> factorize(const MatrixType &)198 IdentityPreconditioner& factorize(const MatrixType& ) { return *this; } 199 200 template<typename MatrixType> compute(const MatrixType &)201 IdentityPreconditioner& compute(const MatrixType& ) { return *this; } 202 203 template<typename Rhs> solve(const Rhs & b)204 inline const Rhs& solve(const Rhs& b) const { return b; } 205 info()206 ComputationInfo info() { return Success; } 207 }; 208 209 } // end namespace Eigen 210 211 #endif // EIGEN_BASIC_PRECONDITIONERS_H 212