IterativeSolverBase.h revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba 
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2011 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
13namespace Eigen {
14
15/** \ingroup IterativeLinearSolvers_Module
16  * \brief Base class for linear iterative solvers
17  *
18  * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
19  */
20template< typename Derived>
21class IterativeSolverBase : internal::noncopyable
22{
23public:
24  typedef typename internal::traits<Derived>::MatrixType MatrixType;
25  typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
26  typedef typename MatrixType::Scalar Scalar;
27  typedef typename MatrixType::Index Index;
28  typedef typename MatrixType::RealScalar RealScalar;
29
30public:
31
32  Derived& derived() { return *static_cast<Derived*>(this); }
33  const Derived& derived() const { return *static_cast<const Derived*>(this); }
34
35  /** Default constructor. */
36  IterativeSolverBase()
37    : mp_matrix(0)
38  {
39    init();
40  }
41
42  /** Initialize the solver with matrix \a A for further \c Ax=b solving.
43    *
44    * This constructor is a shortcut for the default constructor followed
45    * by a call to compute().
46    *
47    * \warning this class stores a reference to the matrix A as well as some
48    * precomputed values that depend on it. Therefore, if \a A is changed
49    * this class becomes invalid. Call compute() to update it with the new
50    * matrix A, or modify a copy of A.
51    */
52  IterativeSolverBase(const MatrixType& A)
53  {
54    init();
55    compute(A);
56  }
57
58  ~IterativeSolverBase() {}
59
60  /** Initializes the iterative solver for the sparcity pattern of the matrix \a A for further solving \c Ax=b problems.
61    *
62    * Currently, this function mostly call analyzePattern on the preconditioner. In the future
63    * we might, for instance, implement column reodering for faster matrix vector products.
64    */
65  Derived& analyzePattern(const MatrixType& A)
66  {
67    m_preconditioner.analyzePattern(A);
68    m_isInitialized = true;
69    m_analysisIsOk = true;
70    m_info = Success;
71    return derived();
72  }
73
74  /** Initializes the iterative solver with the numerical values of the matrix \a A for further solving \c Ax=b problems.
75    *
76    * Currently, this function mostly call factorize on the preconditioner.
77    *
78    * \warning this class stores a reference to the matrix A as well as some
79    * precomputed values that depend on it. Therefore, if \a A is changed
80    * this class becomes invalid. Call compute() to update it with the new
81    * matrix A, or modify a copy of A.
82    */
83  Derived& factorize(const MatrixType& A)
84  {
85    eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
86    mp_matrix = &A;
87    m_preconditioner.factorize(A);
88    m_factorizationIsOk = true;
89    m_info = Success;
90    return derived();
91  }
92
93  /** Initializes the iterative solver with the matrix \a A for further solving \c Ax=b problems.
94    *
95    * Currently, this function mostly initialized/compute the preconditioner. In the future
96    * we might, for instance, implement column reodering for faster matrix vector products.
97    *
98    * \warning this class stores a reference to the matrix A as well as some
99    * precomputed values that depend on it. Therefore, if \a A is changed
100    * this class becomes invalid. Call compute() to update it with the new
101    * matrix A, or modify a copy of A.
102    */
103  Derived& compute(const MatrixType& A)
104  {
105    mp_matrix = &A;
106    m_preconditioner.compute(A);
107    m_isInitialized = true;
108    m_analysisIsOk = true;
109    m_factorizationIsOk = true;
110    m_info = Success;
111    return derived();
112  }
113
114  /** \internal */
115  Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; }
116  /** \internal */
117  Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; }
118
119  /** \returns the tolerance threshold used by the stopping criteria */
120  RealScalar tolerance() const { return m_tolerance; }
121
122  /** Sets the tolerance threshold used by the stopping criteria */
123  Derived& setTolerance(const RealScalar& tolerance)
124  {
125    m_tolerance = tolerance;
126    return derived();
127  }
128
129  /** \returns a read-write reference to the preconditioner for custom configuration. */
130  Preconditioner& preconditioner() { return m_preconditioner; }
131
132  /** \returns a read-only reference to the preconditioner. */
133  const Preconditioner& preconditioner() const { return m_preconditioner; }
134
135  /** \returns the max number of iterations */
136  int maxIterations() const
137  {
138    return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations;
139  }
140
141  /** Sets the max number of iterations */
142  Derived& setMaxIterations(int maxIters)
143  {
144    m_maxIterations = maxIters;
145    return derived();
146  }
147
148  /** \returns the number of iterations performed during the last solve */
149  int iterations() const
150  {
151    eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
152    return m_iterations;
153  }
154
155  /** \returns the tolerance error reached during the last solve */
156  RealScalar error() const
157  {
158    eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
159    return m_error;
160  }
161
162  /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
163    *
164    * \sa compute()
165    */
166  template<typename Rhs> inline const internal::solve_retval<Derived, Rhs>
167  solve(const MatrixBase<Rhs>& b) const
168  {
169    eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
170    eigen_assert(rows()==b.rows()
171              && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
172    return internal::solve_retval<Derived, Rhs>(derived(), b.derived());
173  }
174
175  /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
176    *
177    * \sa compute()
178    */
179  template<typename Rhs>
180  inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs>
181  solve(const SparseMatrixBase<Rhs>& b) const
182  {
183    eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
184    eigen_assert(rows()==b.rows()
185              && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
186    return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived());
187  }
188
189  /** \returns Success if the iterations converged, and NoConvergence otherwise. */
190  ComputationInfo info() const
191  {
192    eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
193    return m_info;
194  }
195
196  /** \internal */
197  template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
198  void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
199  {
200    eigen_assert(rows()==b.rows());
201
202    int rhsCols = b.cols();
203    int size = b.rows();
204    Eigen::Matrix<DestScalar,Dynamic,1> tb(size);
205    Eigen::Matrix<DestScalar,Dynamic,1> tx(size);
206    for(int k=0; k<rhsCols; ++k)
207    {
208      tb = b.col(k);
209      tx = derived().solve(tb);
210      dest.col(k) = tx.sparseView(0);
211    }
212  }
213
214protected:
215  void init()
216  {
217    m_isInitialized = false;
218    m_analysisIsOk = false;
219    m_factorizationIsOk = false;
220    m_maxIterations = -1;
221    m_tolerance = NumTraits<Scalar>::epsilon();
222  }
223  const MatrixType* mp_matrix;
224  Preconditioner m_preconditioner;
225
226  int m_maxIterations;
227  RealScalar m_tolerance;
228
229  mutable RealScalar m_error;
230  mutable int m_iterations;
231  mutable ComputationInfo m_info;
232  mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk;
233};
234
235namespace internal {
236
237template<typename Derived, typename Rhs>
238struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs>
239  : sparse_solve_retval_base<IterativeSolverBase<Derived>, Rhs>
240{
241  typedef IterativeSolverBase<Derived> Dec;
242  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
243
244  template<typename Dest> void evalTo(Dest& dst) const
245  {
246    dec().derived()._solve_sparse(rhs(),dst);
247  }
248};
249
250} // end namespace internal
251
252} // end namespace Eigen
253
254#endif // EIGEN_ITERATIVE_SOLVER_BASE_H
255