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 (at) 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 47 DiagonalPreconditioner() : m_isInitialized(false) {} 48 49 template<typename MatType> 50 explicit DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols()) 51 { 52 compute(mat); 53 } 54 55 Index rows() const { return m_invdiag.size(); } 56 Index cols() const { return m_invdiag.size(); } 57 58 template<typename MatType> 59 DiagonalPreconditioner& analyzePattern(const MatType& ) 60 { 61 return *this; 62 } 63 64 template<typename MatType> 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> 82 DiagonalPreconditioner& compute(const MatType& mat) 83 { 84 return factorize(mat); 85 } 86 87 /** \internal */ 88 template<typename Rhs, typename Dest> 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> 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 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 136 LeastSquareDiagonalPreconditioner() : Base() {} 137 138 template<typename MatType> 139 explicit LeastSquareDiagonalPreconditioner(const MatType& mat) : Base() 140 { 141 compute(mat); 142 } 143 144 template<typename MatType> 145 LeastSquareDiagonalPreconditioner& analyzePattern(const MatType& ) 146 { 147 return *this; 148 } 149 150 template<typename MatType> 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 if(MatType::IsRowMajor) 156 { 157 m_invdiag.setZero(); 158 for(Index j=0; j<mat.outerSize(); ++j) 159 { 160 for(typename MatType::InnerIterator it(mat,j); it; ++it) 161 m_invdiag(it.index()) += numext::abs2(it.value()); 162 } 163 for(Index j=0; j<mat.cols(); ++j) 164 if(numext::real(m_invdiag(j))>RealScalar(0)) 165 m_invdiag(j) = RealScalar(1)/numext::real(m_invdiag(j)); 166 } 167 else 168 { 169 for(Index j=0; j<mat.outerSize(); ++j) 170 { 171 RealScalar sum = mat.innerVector(j).squaredNorm(); 172 if(sum>RealScalar(0)) 173 m_invdiag(j) = RealScalar(1)/sum; 174 else 175 m_invdiag(j) = RealScalar(1); 176 } 177 } 178 Base::m_isInitialized = true; 179 return *this; 180 } 181 182 template<typename MatType> 183 LeastSquareDiagonalPreconditioner& compute(const MatType& mat) 184 { 185 return factorize(mat); 186 } 187 188 ComputationInfo info() { return Success; } 189 190 protected: 191 }; 192 193 /** \ingroup IterativeLinearSolvers_Module 194 * \brief A naive preconditioner which approximates any matrix as the identity matrix 195 * 196 * \implsparsesolverconcept 197 * 198 * \sa class DiagonalPreconditioner 199 */ 200 class IdentityPreconditioner 201 { 202 public: 203 204 IdentityPreconditioner() {} 205 206 template<typename MatrixType> 207 explicit IdentityPreconditioner(const MatrixType& ) {} 208 209 template<typename MatrixType> 210 IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; } 211 212 template<typename MatrixType> 213 IdentityPreconditioner& factorize(const MatrixType& ) { return *this; } 214 215 template<typename MatrixType> 216 IdentityPreconditioner& compute(const MatrixType& ) { return *this; } 217 218 template<typename Rhs> 219 inline const Rhs& solve(const Rhs& b) const { return b; } 220 221 ComputationInfo info() { return Success; } 222 }; 223 224 } // end namespace Eigen 225 226 #endif // EIGEN_BASIC_PRECONDITIONERS_H 227
