Eigen  3.3.7
BasicPreconditioners.h
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@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 
35 template <typename _Scalar>
37 {
38  typedef _Scalar Scalar;
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 
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 
127 template <typename _Scalar>
129 {
130  typedef _Scalar Scalar;
131  typedef typename NumTraits<Scalar>::Real RealScalar;
133  using Base::m_invdiag;
134  public:
135 
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.col(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 
201 {
202  public:
203 
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
A preconditioner based on the digonal entries.
Definition: BasicPreconditioners.h:36
Namespace containing all symbols from the Eigen library.
Definition: Core:306
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:279
Jacobi preconditioner for LeastSquaresConjugateGradient.
Definition: BasicPreconditioners.h:128
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:515
Derived & derived()
Definition: EigenBase.h:45
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Definition: Constants.h:432
Index rows() const
Definition: EigenBase.h:59
A naive preconditioner which approximates any matrix as the identity matrix.
Definition: BasicPreconditioners.h:200
const int Dynamic
Definition: Constants.h:21
Pseudo expression representing a solving operation.
Definition: Solve.h:62
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48