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
13namespace 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  */
32template <typename _Scalar>
33class 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
43    DiagonalPreconditioner() : m_isInitialized(false) {}
44
45    template<typename MatType>
46    DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols())
47    {
48      compute(mat);
49    }
50
51    Index rows() const { return m_invdiag.size(); }
52    Index cols() const { return m_invdiag.size(); }
53
54    template<typename MatType>
55    DiagonalPreconditioner& analyzePattern(const MatType& )
56    {
57      return *this;
58    }
59
60    template<typename MatType>
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>
78    DiagonalPreconditioner& compute(const MatType& mat)
79    {
80      return factorize(mat);
81    }
82
83    template<typename Rhs, typename Dest>
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>
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
103namespace internal {
104
105template<typename _MatrixType, typename Rhs>
106struct 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  */
125class 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