1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2011 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 * This preconditioner is suitable for both selfadjoint and general problems. 27 * The diagonal entries are pre-inverted and stored into a dense vector. 28 * 29 * \note A variant that has yet to be implemented would attempt to preserve the norm of each column. 30 * 31 */ 32 template <typename _Scalar> 33 class DiagonalPreconditioner 34 { 35 typedef _Scalar Scalar; 36 typedef Matrix<Scalar,Dynamic,1> Vector; 37 typedef typename Vector::Index Index; 38 39 public: 40 // this typedef is only to export the scalar type and compile-time dimensions to solve_retval 41 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType; 42 DiagonalPreconditioner()43 DiagonalPreconditioner() : m_isInitialized(false) {} 44 45 template<typename MatType> DiagonalPreconditioner(const MatType & mat)46 DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols()) 47 { 48 compute(mat); 49 } 50 rows()51 Index rows() const { return m_invdiag.size(); } cols()52 Index cols() const { return m_invdiag.size(); } 53 54 template<typename MatType> analyzePattern(const MatType &)55 DiagonalPreconditioner& analyzePattern(const MatType& ) 56 { 57 return *this; 58 } 59 60 template<typename MatType> factorize(const MatType & mat)61 DiagonalPreconditioner& factorize(const MatType& mat) 62 { 63 m_invdiag.resize(mat.cols()); 64 for(int j=0; j<mat.outerSize(); ++j) 65 { 66 typename MatType::InnerIterator it(mat,j); 67 while(it && it.index()!=j) ++it; 68 if(it && it.index()==j) 69 m_invdiag(j) = Scalar(1)/it.value(); 70 else 71 m_invdiag(j) = 0; 72 } 73 m_isInitialized = true; 74 return *this; 75 } 76 77 template<typename MatType> compute(const MatType & mat)78 DiagonalPreconditioner& compute(const MatType& mat) 79 { 80 return factorize(mat); 81 } 82 83 template<typename Rhs, typename Dest> _solve(const Rhs & b,Dest & x)84 void _solve(const Rhs& b, Dest& x) const 85 { 86 x = m_invdiag.array() * b.array() ; 87 } 88 89 template<typename Rhs> inline const internal::solve_retval<DiagonalPreconditioner, Rhs> solve(const MatrixBase<Rhs> & b)90 solve(const MatrixBase<Rhs>& b) const 91 { 92 eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized."); 93 eigen_assert(m_invdiag.size()==b.rows() 94 && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b"); 95 return internal::solve_retval<DiagonalPreconditioner, Rhs>(*this, b.derived()); 96 } 97 98 protected: 99 Vector m_invdiag; 100 bool m_isInitialized; 101 }; 102 103 namespace internal { 104 105 template<typename _MatrixType, typename Rhs> 106 struct solve_retval<DiagonalPreconditioner<_MatrixType>, Rhs> 107 : solve_retval_base<DiagonalPreconditioner<_MatrixType>, Rhs> 108 { 109 typedef DiagonalPreconditioner<_MatrixType> Dec; 110 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 111 112 template<typename Dest> void evalTo(Dest& dst) const 113 { 114 dec()._solve(rhs(),dst); 115 } 116 }; 117 118 } 119 120 /** \ingroup IterativeLinearSolvers_Module 121 * \brief A naive preconditioner which approximates any matrix as the identity matrix 122 * 123 * \sa class DiagonalPreconditioner 124 */ 125 class IdentityPreconditioner 126 { 127 public: 128 129 IdentityPreconditioner() {} 130 131 template<typename MatrixType> 132 IdentityPreconditioner(const MatrixType& ) {} 133 134 template<typename MatrixType> 135 IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; } 136 137 template<typename MatrixType> 138 IdentityPreconditioner& factorize(const MatrixType& ) { return *this; } 139 140 template<typename MatrixType> 141 IdentityPreconditioner& compute(const MatrixType& ) { return *this; } 142 143 template<typename Rhs> 144 inline const Rhs& solve(const Rhs& b) const { return b; } 145 }; 146 147 } // end namespace Eigen 148 149 #endif // EIGEN_BASIC_PRECONDITIONERS_H 150