Eigen  3.2.91
IterativeSolverBase.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_ITERATIVE_SOLVER_BASE_H
11 #define EIGEN_ITERATIVE_SOLVER_BASE_H
12 
13 namespace Eigen {
14 
20 template< typename Derived>
21 class IterativeSolverBase : public SparseSolverBase<Derived>
22 {
23 protected:
25  using Base::m_isInitialized;
26 
27 public:
28  typedef typename internal::traits<Derived>::MatrixType MatrixType;
29  typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
30  typedef typename MatrixType::Scalar Scalar;
31  typedef typename MatrixType::StorageIndex StorageIndex;
32  typedef typename MatrixType::RealScalar RealScalar;
33 
34 public:
35 
36  using Base::derived;
37 
40  : m_dummy(0,0), mp_matrix(m_dummy)
41  {
42  init();
43  }
44 
55  template<typename MatrixDerived>
57  : mp_matrix(A.derived())
58  {
59  init();
60  compute(mp_matrix);
61  }
62 
64 
70  template<typename MatrixDerived>
72  {
73  grab(A.derived());
74  m_preconditioner.analyzePattern(mp_matrix);
75  m_isInitialized = true;
76  m_analysisIsOk = true;
77  m_info = m_preconditioner.info();
78  return derived();
79  }
80 
90  template<typename MatrixDerived>
92  {
93  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
94  grab(A.derived());
95  m_preconditioner.factorize(mp_matrix);
96  m_factorizationIsOk = true;
97  m_info = m_preconditioner.info();
98  return derived();
99  }
100 
111  template<typename MatrixDerived>
113  {
114  grab(A.derived());
115  m_preconditioner.compute(mp_matrix);
116  m_isInitialized = true;
117  m_analysisIsOk = true;
118  m_factorizationIsOk = true;
119  m_info = m_preconditioner.info();
120  return derived();
121  }
122 
124  Index rows() const { return mp_matrix.rows(); }
125 
127  Index cols() const { return mp_matrix.cols(); }
128 
132  RealScalar tolerance() const { return m_tolerance; }
133 
139  Derived& setTolerance(const RealScalar& tolerance)
140  {
141  m_tolerance = tolerance;
142  return derived();
143  }
144 
146  Preconditioner& preconditioner() { return m_preconditioner; }
147 
149  const Preconditioner& preconditioner() const { return m_preconditioner; }
150 
155  Index maxIterations() const
156  {
157  return (m_maxIterations<0) ? 2*mp_matrix.cols() : m_maxIterations;
158  }
159 
163  Derived& setMaxIterations(Index maxIters)
164  {
165  m_maxIterations = maxIters;
166  return derived();
167  }
168 
170  Index iterations() const
171  {
172  eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
173  return m_iterations;
174  }
175 
179  RealScalar error() const
180  {
181  eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
182  return m_error;
183  }
184 
190  template<typename Rhs,typename Guess>
192  solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
193  {
194  eigen_assert(m_isInitialized && "Solver is not initialized.");
195  eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
196  return SolveWithGuess<Derived, Rhs, Guess>(derived(), b.derived(), x0);
197  }
198 
201  {
202  eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
203  return m_info;
204  }
205 
207  template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
208  void _solve_impl(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
209  {
210  eigen_assert(rows()==b.rows());
211 
212  Index rhsCols = b.cols();
213  Index size = b.rows();
216  // We do not directly fill dest because sparse expressions have to be free of aliasing issue.
217  // For non square least-square problems, b and dest might not have the same size whereas they might alias each-other.
219  for(Index k=0; k<rhsCols; ++k)
220  {
221  tb = b.col(k);
222  tx = derived().solve(tb);
223  tmp.col(k) = tx.sparseView(0);
224  }
225  tmp.swap(dest);
226  }
227 
228 protected:
229  void init()
230  {
231  m_isInitialized = false;
232  m_analysisIsOk = false;
233  m_factorizationIsOk = false;
234  m_maxIterations = -1;
235  m_tolerance = NumTraits<Scalar>::epsilon();
236  }
237 
238  template<typename MatrixDerived>
239  void grab(const EigenBase<MatrixDerived> &A)
240  {
241  mp_matrix.~Ref<const MatrixType>();
242  ::new (&mp_matrix) Ref<const MatrixType>(A.derived());
243  }
244 
245  void grab(const Ref<const MatrixType> &A)
246  {
247  if(&(A.derived()) != &mp_matrix)
248  {
249  mp_matrix.~Ref<const MatrixType>();
250  ::new (&mp_matrix) Ref<const MatrixType>(A);
251  }
252  }
253 
254  MatrixType m_dummy;
255  Ref<const MatrixType> mp_matrix;
256  Preconditioner m_preconditioner;
257 
258  Index m_maxIterations;
259  RealScalar m_tolerance;
260 
261  mutable RealScalar m_error;
262  mutable Index m_iterations;
263  mutable ComputationInfo m_info;
264  mutable bool m_analysisIsOk, m_factorizationIsOk;
265 };
266 
267 } // end namespace Eigen
268 
269 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H
Pseudo expression representing a solving operation.
Definition: SolveWithGuess.h:15
Derived & setTolerance(const RealScalar &tolerance)
Definition: IterativeSolverBase.h:139
A versatible sparse matrix representation.
Definition: SparseMatrix.h:92
IterativeSolverBase(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:56
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
Definition: LDLT.h:16
Derived & derived()
Definition: EigenBase.h:44
Derived & factorize(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:91
Index maxIterations() const
Definition: IterativeSolverBase.h:155
Definition: EigenBase.h:28
Derived & analyzePattern(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:71
Derived & compute(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:112
const SolveWithGuess< Derived, Rhs, Guess > solveWithGuess(const MatrixBase< Rhs > &b, const Guess &x0) const
Definition: IterativeSolverBase.h:192
ComputationInfo info() const
Definition: IterativeSolverBase.h:200
RealScalar tolerance() const
Definition: IterativeSolverBase.h:132
Derived & setMaxIterations(Index maxIters)
Definition: IterativeSolverBase.h:163
Preconditioner & preconditioner()
Definition: IterativeSolverBase.h:146
RealScalar error() const
Definition: IterativeSolverBase.h:179
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
IterativeSolverBase()
Definition: IterativeSolverBase.h:39
ComputationInfo
Definition: Constants.h:422
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:21
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const Preconditioner & preconditioner() const
Definition: IterativeSolverBase.h:149
Index iterations() const
Definition: IterativeSolverBase.h:170