Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
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 #include "./InternalHeaderCheck.h"
14 
15 namespace Eigen {
16 
17 namespace internal {
18 
19 template<typename MatrixType>
20 struct is_ref_compatible_impl
21 {
22 private:
23  template <typename T0>
24  struct any_conversion
25  {
26  template <typename T> any_conversion(const volatile T&);
27  template <typename T> any_conversion(T&);
28  };
29  struct yes {int a[1];};
30  struct no {int a[2];};
31 
32  template<typename T>
33  static yes test(const Ref<const T>&, int);
34  template<typename T>
35  static no test(any_conversion<T>, ...);
36 
37 public:
38  static MatrixType ms_from;
39  enum { value = sizeof(test<MatrixType>(ms_from, 0))==sizeof(yes) };
40 };
41 
42 template<typename MatrixType>
43 struct is_ref_compatible
44 {
45  enum { value = is_ref_compatible_impl<remove_all_t<MatrixType>>::value };
46 };
47 
48 template<typename MatrixType, bool MatrixFree = !internal::is_ref_compatible<MatrixType>::value>
49 class generic_matrix_wrapper;
50 
51 // We have an explicit matrix at hand, compatible with Ref<>
52 template<typename MatrixType>
53 class generic_matrix_wrapper<MatrixType,false>
54 {
55 public:
56  typedef Ref<const MatrixType> ActualMatrixType;
57  template<int UpLo> struct ConstSelfAdjointViewReturnType {
58  typedef typename ActualMatrixType::template ConstSelfAdjointViewReturnType<UpLo>::Type Type;
59  };
60 
61  enum {
62  MatrixFree = false
63  };
64 
65  generic_matrix_wrapper()
66  : m_dummy(0,0), m_matrix(m_dummy)
67  {}
68 
69  template<typename InputType>
70  generic_matrix_wrapper(const InputType &mat)
71  : m_matrix(mat)
72  {}
73 
74  const ActualMatrixType& matrix() const
75  {
76  return m_matrix;
77  }
78 
79  template<typename MatrixDerived>
80  void grab(const EigenBase<MatrixDerived> &mat)
81  {
82  internal::destroy_at(&m_matrix);
83  internal::construct_at(&m_matrix, mat.derived());
84  }
85 
86  void grab(const Ref<const MatrixType> &mat)
87  {
88  if(&(mat.derived()) != &m_matrix)
89  {
90  internal::destroy_at(&m_matrix);
91  internal::construct_at(&m_matrix, mat);
92  }
93  }
94 
95 protected:
96  MatrixType m_dummy; // used to default initialize the Ref<> object
97  ActualMatrixType m_matrix;
98 };
99 
100 // MatrixType is not compatible with Ref<> -> matrix-free wrapper
101 template<typename MatrixType>
102 class generic_matrix_wrapper<MatrixType,true>
103 {
104 public:
105  typedef MatrixType ActualMatrixType;
106  template<int UpLo> struct ConstSelfAdjointViewReturnType
107  {
108  typedef ActualMatrixType Type;
109  };
110 
111  enum {
112  MatrixFree = true
113  };
114 
115  generic_matrix_wrapper()
116  : mp_matrix(0)
117  {}
118 
119  generic_matrix_wrapper(const MatrixType &mat)
120  : mp_matrix(&mat)
121  {}
122 
123  const ActualMatrixType& matrix() const
124  {
125  return *mp_matrix;
126  }
127 
128  void grab(const MatrixType &mat)
129  {
130  mp_matrix = &mat;
131  }
132 
133 protected:
134  const ActualMatrixType *mp_matrix;
135 };
136 
137 }
138 
144 template< typename Derived>
145 class IterativeSolverBase : public SparseSolverBase<Derived>
146 {
147 protected:
149  using Base::m_isInitialized;
150 
151 public:
152  typedef typename internal::traits<Derived>::MatrixType MatrixType;
153  typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
154  typedef typename MatrixType::Scalar Scalar;
155  typedef typename MatrixType::StorageIndex StorageIndex;
156  typedef typename MatrixType::RealScalar RealScalar;
157 
158  enum {
159  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
160  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
161  };
162 
163 public:
164 
165  using Base::derived;
166 
169  {
170  init();
171  }
172 
183  template<typename MatrixDerived>
185  : m_matrixWrapper(A.derived())
186  {
187  init();
188  compute(matrix());
189  }
190 
191 
193 
194  ~IterativeSolverBase() {}
195 
201  template<typename MatrixDerived>
203  {
204  grab(A.derived());
205  m_preconditioner.analyzePattern(matrix());
206  m_isInitialized = true;
207  m_analysisIsOk = true;
208  m_info = m_preconditioner.info();
209  return derived();
210  }
211 
221  template<typename MatrixDerived>
223  {
224  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
225  grab(A.derived());
226  m_preconditioner.factorize(matrix());
227  m_factorizationIsOk = true;
228  m_info = m_preconditioner.info();
229  return derived();
230  }
231 
242  template<typename MatrixDerived>
244  {
245  grab(A.derived());
246  m_preconditioner.compute(matrix());
247  m_isInitialized = true;
248  m_analysisIsOk = true;
249  m_factorizationIsOk = true;
250  m_info = m_preconditioner.info();
251  return derived();
252  }
253 
255  EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return matrix().rows(); }
256 
258  EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return matrix().cols(); }
259 
263  RealScalar tolerance() const { return m_tolerance; }
264 
270  Derived& setTolerance(const RealScalar& tolerance)
271  {
272  m_tolerance = tolerance;
273  return derived();
274  }
275 
277  Preconditioner& preconditioner() { return m_preconditioner; }
278 
280  const Preconditioner& preconditioner() const { return m_preconditioner; }
281 
287  {
288  return (m_maxIterations<0) ? 2*matrix().cols() : m_maxIterations;
289  }
290 
294  Derived& setMaxIterations(Index maxIters)
295  {
296  m_maxIterations = maxIters;
297  return derived();
298  }
299 
302  {
303  eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
304  return m_iterations;
305  }
306 
310  RealScalar error() const
311  {
312  eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
313  return m_error;
314  }
315 
321  template<typename Rhs,typename Guess>
323  solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
324  {
325  eigen_assert(m_isInitialized && "Solver is not initialized.");
326  eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
327  return SolveWithGuess<Derived, Rhs, Guess>(derived(), b.derived(), x0);
328  }
329 
332  {
333  eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
334  return m_info;
335  }
336 
338  template<typename Rhs, typename DestDerived>
339  void _solve_with_guess_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const
340  {
341  eigen_assert(rows()==b.rows());
342 
343  Index rhsCols = b.cols();
344  Index size = b.rows();
345  DestDerived& dest(aDest.derived());
346  typedef typename DestDerived::Scalar DestScalar;
349  // We do not directly fill dest because sparse expressions have to be free of aliasing issue.
350  // For non square least-square problems, b and dest might not have the same size whereas they might alias each-other.
351  typename DestDerived::PlainObject tmp(cols(),rhsCols);
352  ComputationInfo global_info = Success;
353  for(Index k=0; k<rhsCols; ++k)
354  {
355  tb = b.col(k);
356  tx = dest.col(k);
357  derived()._solve_vector_with_guess_impl(tb,tx);
358  tmp.col(k) = tx.sparseView(0);
359 
360  // The call to _solve_vector_with_guess_impl updates m_info, so if it failed for a previous column
361  // we need to restore it to the worst value.
362  if(m_info==NumericalIssue)
363  global_info = NumericalIssue;
364  else if(m_info==NoConvergence)
365  global_info = NoConvergence;
366  }
367  m_info = global_info;
368  dest.swap(tmp);
369  }
370 
371  template<typename Rhs, typename DestDerived>
372  std::enable_if_t<Rhs::ColsAtCompileTime!=1 && DestDerived::ColsAtCompileTime!=1>
373  _solve_with_guess_impl(const Rhs& b, MatrixBase<DestDerived> &aDest) const
374  {
375  eigen_assert(rows()==b.rows());
376 
377  Index rhsCols = b.cols();
378  DestDerived& dest(aDest.derived());
379  ComputationInfo global_info = Success;
380  for(Index k=0; k<rhsCols; ++k)
381  {
382  typename DestDerived::ColXpr xk(dest,k);
383  typename Rhs::ConstColXpr bk(b,k);
384  derived()._solve_vector_with_guess_impl(bk,xk);
385 
386  // The call to _solve_vector_with_guess updates m_info, so if it failed for a previous column
387  // we need to restore it to the worst value.
388  if(m_info==NumericalIssue)
389  global_info = NumericalIssue;
390  else if(m_info==NoConvergence)
391  global_info = NoConvergence;
392  }
393  m_info = global_info;
394  }
395 
396  template<typename Rhs, typename DestDerived>
397  std::enable_if_t<Rhs::ColsAtCompileTime==1 || DestDerived::ColsAtCompileTime==1>
398  _solve_with_guess_impl(const Rhs& b, MatrixBase<DestDerived> &dest) const
399  {
400  derived()._solve_vector_with_guess_impl(b,dest.derived());
401  }
402 
404  template<typename Rhs,typename Dest>
405  void _solve_impl(const Rhs& b, Dest& x) const
406  {
407  x.setZero();
408  derived()._solve_with_guess_impl(b,x);
409  }
410 
411 protected:
412  void init()
413  {
414  m_isInitialized = false;
415  m_analysisIsOk = false;
416  m_factorizationIsOk = false;
417  m_maxIterations = -1;
418  m_tolerance = NumTraits<Scalar>::epsilon();
419  }
420 
421  typedef internal::generic_matrix_wrapper<MatrixType> MatrixWrapper;
422  typedef typename MatrixWrapper::ActualMatrixType ActualMatrixType;
423 
424  const ActualMatrixType& matrix() const
425  {
426  return m_matrixWrapper.matrix();
427  }
428 
429  template<typename InputType>
430  void grab(const InputType &A)
431  {
432  m_matrixWrapper.grab(A);
433  }
434 
435  MatrixWrapper m_matrixWrapper;
436  Preconditioner m_preconditioner;
437 
438  Index m_maxIterations;
439  RealScalar m_tolerance;
440 
441  mutable RealScalar m_error;
442  mutable Index m_iterations;
443  mutable ComputationInfo m_info;
444  mutable bool m_analysisIsOk, m_factorizationIsOk;
445 };
446 
447 } // end namespace Eigen
448 
449 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H
Derived & derived()
Definition: EigenBase.h:48
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:146
IterativeSolverBase()
Definition: IterativeSolverBase.h:168
ComputationInfo info() const
Definition: IterativeSolverBase.h:331
RealScalar error() const
Definition: IterativeSolverBase.h:310
Derived & factorize(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:222
Index maxIterations() const
Definition: IterativeSolverBase.h:286
IterativeSolverBase(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:184
Derived & analyzePattern(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:202
Preconditioner & preconditioner()
Definition: IterativeSolverBase.h:277
const Preconditioner & preconditioner() const
Definition: IterativeSolverBase.h:280
Derived & compute(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:243
Derived & setTolerance(const RealScalar &tolerance)
Definition: IterativeSolverBase.h:270
RealScalar tolerance() const
Definition: IterativeSolverBase.h:263
const SolveWithGuess< Derived, Rhs, Guess > solveWithGuess(const MatrixBase< Rhs > &b, const Guess &x0) const
Definition: IterativeSolverBase.h:323
Index iterations() const
Definition: IterativeSolverBase.h:301
Derived & setMaxIterations(Index maxIters)
Definition: IterativeSolverBase.h:294
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
Pseudo expression representing a solving operation.
Definition: SolveWithGuess.h:17
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:30
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
@ NoConvergence
Definition: Constants.h:448
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48