BasicPreconditioners.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_BASIC_PRECONDITIONERS_H
26 #define EIGEN_BASIC_PRECONDITIONERS_H
27 
28 namespace Eigen {
29 
47 template <typename _Scalar>
49 {
50  typedef _Scalar Scalar;
52  typedef typename Vector::Index Index;
53 
54  public:
56 
58 
59  template<typename MatrixType>
61  {
62  compute(mat);
63  }
64 
65  Index rows() const { return m_invdiag.size(); }
66  Index cols() const { return m_invdiag.size(); }
67 
68  template<typename MatrixType>
70  {
71  return *this;
72  }
73 
74  template<typename MatrixType>
76  {
77  m_invdiag.resize(mat.cols());
78  for(int j=0; j<mat.outerSize(); ++j)
79  {
80  typename MatrixType::InnerIterator it(mat,j);
81  while(it && it.index()!=j) ++it;
82  if(it && it.index()==j)
83  m_invdiag(j) = Scalar(1)/it.value();
84  else
85  m_invdiag(j) = 0;
86  }
87  m_isInitialized = true;
88  return *this;
89  }
90 
91  template<typename MatrixType>
93  {
94  return factorize(mat);
95  }
96 
97  template<typename Rhs, typename Dest>
98  void _solve(const Rhs& b, Dest& x) const
99  {
100  x = m_invdiag.array() * b.array() ;
101  }
102 
103  template<typename Rhs> inline const internal::solve_retval<DiagonalPreconditioner, Rhs>
104  solve(const MatrixBase<Rhs>& b) const
105  {
106  eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
107  eigen_assert(m_invdiag.size()==b.rows()
108  && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
109  return internal::solve_retval<DiagonalPreconditioner, Rhs>(*this, b.derived());
110  }
111 
112  protected:
115 };
116 
117 namespace internal {
118 
119 template<typename _MatrixType, typename Rhs>
120 struct solve_retval<DiagonalPreconditioner<_MatrixType>, Rhs>
121  : solve_retval_base<DiagonalPreconditioner<_MatrixType>, Rhs>
122 {
124  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
125 
126  template<typename Dest> void evalTo(Dest& dst) const
127  {
128  dec()._solve(rhs(),dst);
129  }
130 };
131 
132 }
133 
140 {
141  public:
142 
144 
145  template<typename MatrixType>
146  IdentityPreconditioner(const MatrixType& ) {}
147 
148  template<typename MatrixType>
149  IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
150 
151  template<typename MatrixType>
152  IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
153 
154  template<typename MatrixType>
155  IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
156 
157  template<typename Rhs>
158  inline const Rhs& solve(const Rhs& b) const { return b; }
159 };
160 
161 } // end namespace Eigen
162 
163 #endif // EIGEN_BASIC_PRECONDITIONERS_H