Eigen  3.2.92
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 
15 namespace internal {
16 
17 template<typename MatrixType>
18 struct is_ref_compatible_impl
19 {
20 private:
21  template <typename T0>
22  struct any_conversion
23  {
24  template <typename T> any_conversion(const volatile T&);
25  template <typename T> any_conversion(T&);
26  };
27  struct yes {int a[1];};
28  struct no {int a[2];};
29 
30  template<typename T>
31  static yes test(const Ref<const T>&, int);
32  template<typename T>
33  static no test(any_conversion<T>, ...);
34 
35 public:
36  static MatrixType ms_from;
37  enum { value = sizeof(test<MatrixType>(ms_from, 0))==sizeof(yes) };
38 };
39 
40 template<typename MatrixType>
41 struct is_ref_compatible
42 {
43  enum { value = is_ref_compatible_impl<typename remove_all<MatrixType>::type>::value };
44 };
45 
46 template<typename MatrixType, bool MatrixFree = !internal::is_ref_compatible<MatrixType>::value>
47 class generic_matrix_wrapper;
48 
49 // We have an explicit matrix at hand, compatible with Ref<>
50 template<typename MatrixType>
51 class generic_matrix_wrapper<MatrixType,false>
52 {
53 public:
54  typedef Ref<const MatrixType> ActualMatrixType;
55  template<int UpLo> struct ConstSelfAdjointViewReturnType {
56  typedef typename ActualMatrixType::template ConstSelfAdjointViewReturnType<UpLo>::Type Type;
57  };
58 
59  enum {
60  MatrixFree = false
61  };
62 
63  generic_matrix_wrapper()
64  : m_dummy(0,0), m_matrix(m_dummy)
65  {}
66 
67  template<typename InputType>
68  generic_matrix_wrapper(const InputType &mat)
69  : m_matrix(mat)
70  {}
71 
72  const ActualMatrixType& matrix() const
73  {
74  return m_matrix;
75  }
76 
77  template<typename MatrixDerived>
78  void grab(const EigenBase<MatrixDerived> &mat)
79  {
80  m_matrix.~Ref<const MatrixType>();
81  ::new (&m_matrix) Ref<const MatrixType>(mat.derived());
82  }
83 
84  void grab(const Ref<const MatrixType> &mat)
85  {
86  if(&(mat.derived()) != &m_matrix)
87  {
88  m_matrix.~Ref<const MatrixType>();
89  ::new (&m_matrix) Ref<const MatrixType>(mat);
90  }
91  }
92 
93 protected:
94  MatrixType m_dummy; // used to default initialize the Ref<> object
95  ActualMatrixType m_matrix;
96 };
97 
98 // MatrixType is not compatible with Ref<> -> matrix-free wrapper
99 template<typename MatrixType>
100 class generic_matrix_wrapper<MatrixType,true>
101 {
102 public:
103  typedef MatrixType ActualMatrixType;
104  template<int UpLo> struct ConstSelfAdjointViewReturnType
105  {
106  typedef ActualMatrixType Type;
107  };
108 
109  enum {
110  MatrixFree = true
111  };
112 
113  generic_matrix_wrapper()
114  : mp_matrix(0)
115  {}
116 
117  generic_matrix_wrapper(const MatrixType &mat)
118  : mp_matrix(&mat)
119  {}
120 
121  const ActualMatrixType& matrix() const
122  {
123  return *mp_matrix;
124  }
125 
126  void grab(const MatrixType &mat)
127  {
128  mp_matrix = &mat;
129  }
130 
131 protected:
132  const ActualMatrixType *mp_matrix;
133 };
134 
135 }
136 
142 template< typename Derived>
143 class IterativeSolverBase : public SparseSolverBase<Derived>
144 {
145 protected:
147  using Base::m_isInitialized;
148 
149 public:
150  typedef typename internal::traits<Derived>::MatrixType MatrixType;
151  typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
152  typedef typename MatrixType::Scalar Scalar;
153  typedef typename MatrixType::StorageIndex StorageIndex;
154  typedef typename MatrixType::RealScalar RealScalar;
155 
156  enum {
157  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
158  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
159  };
160 
161 public:
162 
163  using Base::derived;
164 
167  {
168  init();
169  }
170 
181  template<typename MatrixDerived>
183  : m_matrixWrapper(A.derived())
184  {
185  init();
186  compute(matrix());
187  }
188 
189  ~IterativeSolverBase() {}
190 
196  template<typename MatrixDerived>
198  {
199  grab(A.derived());
200  m_preconditioner.analyzePattern(matrix());
201  m_isInitialized = true;
202  m_analysisIsOk = true;
203  m_info = m_preconditioner.info();
204  return derived();
205  }
206 
216  template<typename MatrixDerived>
218  {
219  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
220  grab(A.derived());
221  m_preconditioner.factorize(matrix());
222  m_factorizationIsOk = true;
223  m_info = m_preconditioner.info();
224  return derived();
225  }
226 
237  template<typename MatrixDerived>
239  {
240  grab(A.derived());
241  m_preconditioner.compute(matrix());
242  m_isInitialized = true;
243  m_analysisIsOk = true;
244  m_factorizationIsOk = true;
245  m_info = m_preconditioner.info();
246  return derived();
247  }
248 
250  Index rows() const { return matrix().rows(); }
251 
253  Index cols() const { return matrix().cols(); }
254 
258  RealScalar tolerance() const { return m_tolerance; }
259 
265  Derived& setTolerance(const RealScalar& tolerance)
266  {
267  m_tolerance = tolerance;
268  return derived();
269  }
270 
272  Preconditioner& preconditioner() { return m_preconditioner; }
273 
275  const Preconditioner& preconditioner() const { return m_preconditioner; }
276 
281  Index maxIterations() const
282  {
283  return (m_maxIterations<0) ? 2*matrix().cols() : m_maxIterations;
284  }
285 
289  Derived& setMaxIterations(Index maxIters)
290  {
291  m_maxIterations = maxIters;
292  return derived();
293  }
294 
296  Index iterations() const
297  {
298  eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
299  return m_iterations;
300  }
301 
305  RealScalar error() const
306  {
307  eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
308  return m_error;
309  }
310 
316  template<typename Rhs,typename Guess>
318  solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
319  {
320  eigen_assert(m_isInitialized && "Solver is not initialized.");
321  eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
322  return SolveWithGuess<Derived, Rhs, Guess>(derived(), b.derived(), x0);
323  }
324 
327  {
328  eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
329  return m_info;
330  }
331 
333  template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
334  void _solve_impl(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
335  {
336  eigen_assert(rows()==b.rows());
337 
338  Index rhsCols = b.cols();
339  Index size = b.rows();
342  // We do not directly fill dest because sparse expressions have to be free of aliasing issue.
343  // For non square least-square problems, b and dest might not have the same size whereas they might alias each-other.
345  for(Index k=0; k<rhsCols; ++k)
346  {
347  tb = b.col(k);
348  tx = derived().solve(tb);
349  tmp.col(k) = tx.sparseView(0);
350  }
351  tmp.swap(dest);
352  }
353 
354 protected:
355  void init()
356  {
357  m_isInitialized = false;
358  m_analysisIsOk = false;
359  m_factorizationIsOk = false;
360  m_maxIterations = -1;
361  m_tolerance = NumTraits<Scalar>::epsilon();
362  }
363 
364  typedef internal::generic_matrix_wrapper<MatrixType> MatrixWrapper;
365  typedef typename MatrixWrapper::ActualMatrixType ActualMatrixType;
366 
367  const ActualMatrixType& matrix() const
368  {
369  return m_matrixWrapper.matrix();
370  }
371 
372  template<typename InputType>
373  void grab(const InputType &A)
374  {
375  m_matrixWrapper.grab(A);
376  }
377 
378  MatrixWrapper m_matrixWrapper;
379  Preconditioner m_preconditioner;
380 
381  Index m_maxIterations;
382  RealScalar m_tolerance;
383 
384  mutable RealScalar m_error;
385  mutable Index m_iterations;
386  mutable ComputationInfo m_info;
387  mutable bool m_analysisIsOk, m_factorizationIsOk;
388 };
389 
390 } // end namespace Eigen
391 
392 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H
ColXpr col(Index i)
Definition: DenseBase.h:778
Pseudo expression representing a solving operation.
Definition: SolveWithGuess.h:15
Derived & setTolerance(const RealScalar &tolerance)
Definition: IterativeSolverBase.h:265
A versatible sparse matrix representation.
Definition: SparseMatrix.h:92
IterativeSolverBase(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:182
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
Definition: LDLT.h:16
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
Derived & derived()
Definition: EigenBase.h:44
Derived & factorize(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:217
Index maxIterations() const
Definition: IterativeSolverBase.h:281
Definition: EigenBase.h:28
Derived & analyzePattern(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:197
Derived & compute(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:238
const SolveWithGuess< Derived, Rhs, Guess > solveWithGuess(const MatrixBase< Rhs > &b, const Guess &x0) const
Definition: IterativeSolverBase.h:318
void swap(SparseMatrix &other)
Definition: SparseMatrix.h:722
ComputationInfo info() const
Definition: IterativeSolverBase.h:326
RealScalar tolerance() const
Definition: IterativeSolverBase.h:258
Derived & setMaxIterations(Index maxIters)
Definition: IterativeSolverBase.h:289
Definition: Eigen_Colamd.h:54
ColXpr col(Index i)
Definition: SparseMatrixBase.h:778
Preconditioner & preconditioner()
Definition: IterativeSolverBase.h:272
RealScalar error() const
Definition: IterativeSolverBase.h:305
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
IterativeSolverBase()
Definition: IterativeSolverBase.h:166
ComputationInfo
Definition: Constants.h:430
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:143
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const Preconditioner & preconditioner() const
Definition: IterativeSolverBase.h:275
Index iterations() const
Definition: IterativeSolverBase.h:296