Eigen  3.2.91
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 
33 template <typename _Scalar>
35 {
36  typedef _Scalar Scalar;
38  public:
39  typedef typename Vector::StorageIndex StorageIndex;
40  // this typedef is only to export the scalar type and compile-time dimensions to solve_retval
42 
43  DiagonalPreconditioner() : m_isInitialized(false) {}
44 
45  template<typename MatType>
46  explicit 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 && it.value()!=Scalar(0))
69  m_invdiag(j) = Scalar(1)/it.value();
70  else
71  m_invdiag(j) = Scalar(1);
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 
84  template<typename Rhs, typename Dest>
85  void _solve_impl(const Rhs& b, Dest& x) const
86  {
87  x = m_invdiag.array() * b.array() ;
88  }
89 
90  template<typename Rhs> inline const Solve<DiagonalPreconditioner, Rhs>
91  solve(const MatrixBase<Rhs>& b) const
92  {
93  eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
94  eigen_assert(m_invdiag.size()==b.rows()
95  && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
96  return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived());
97  }
98 
99  ComputationInfo info() { return Success; }
100 
101  protected:
102  Vector m_invdiag;
103  bool m_isInitialized;
104 };
105 
121 template <typename _Scalar>
123 {
124  typedef _Scalar Scalar;
125  typedef typename NumTraits<Scalar>::Real RealScalar;
127  using Base::m_invdiag;
128  public:
129 
131 
132  template<typename MatType>
133  explicit LeastSquareDiagonalPreconditioner(const MatType& mat) : Base()
134  {
135  compute(mat);
136  }
137 
138  template<typename MatType>
139  LeastSquareDiagonalPreconditioner& analyzePattern(const MatType& )
140  {
141  return *this;
142  }
143 
144  template<typename MatType>
145  LeastSquareDiagonalPreconditioner& factorize(const MatType& mat)
146  {
147  // Compute the inverse squared-norm of each column of mat
148  m_invdiag.resize(mat.cols());
149  for(Index j=0; j<mat.outerSize(); ++j)
150  {
151  RealScalar sum = mat.innerVector(j).squaredNorm();
152  if(sum>0)
153  m_invdiag(j) = RealScalar(1)/sum;
154  else
155  m_invdiag(j) = RealScalar(1);
156  }
157  Base::m_isInitialized = true;
158  return *this;
159  }
160 
161  template<typename MatType>
162  LeastSquareDiagonalPreconditioner& compute(const MatType& mat)
163  {
164  return factorize(mat);
165  }
166 
167  ComputationInfo info() { return Success; }
168 
169  protected:
170 };
171 
178 {
179  public:
180 
182 
183  template<typename MatrixType>
184  explicit IdentityPreconditioner(const MatrixType& ) {}
185 
186  template<typename MatrixType>
187  IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
188 
189  template<typename MatrixType>
190  IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
191 
192  template<typename MatrixType>
193  IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
194 
195  template<typename Rhs>
196  inline const Rhs& solve(const Rhs& b) const { return b; }
197 
198  ComputationInfo info() { return Success; }
199 };
200 
201 } // end namespace Eigen
202 
203 #endif // EIGEN_BASIC_PRECONDITIONERS_H
A preconditioner based on the digonal entries.
Definition: BasicPreconditioners.h:34
Definition: LDLT.h:16
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:252
Jacobi preconditioner for LeastSquaresConjugateGradient.
Definition: BasicPreconditioners.h:122
Definition: Constants.h:424
A naive preconditioner which approximates any matrix as the identity matrix.
Definition: BasicPreconditioners.h:177
Pseudo expression representing a solving operation.
Definition: Solve.h:63
ComputationInfo
Definition: Constants.h:422
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48