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