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 
13 namespace Eigen {
14 
15 /** \ingroup IterativeLinearSolvers_Module
16   * \brief Base class for linear iterative solvers
17   *
18   * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
19   */
20 template< typename Derived>
21 class IterativeSolverBase : internal::noncopyable
22 {
23 public:
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 
30 public:
31 
derived()32   Derived& derived() { return *static_cast<Derived*>(this); }
derived()33   const Derived& derived() const { return *static_cast<const Derived*>(this); }
34 
35   /** Default constructor. */
IterativeSolverBase()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     */
IterativeSolverBase(const MatrixType & A)52   IterativeSolverBase(const MatrixType& A)
53   {
54     init();
55     compute(A);
56   }
57 
~IterativeSolverBase()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     */
analyzePattern(const MatrixType & A)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     */
factorize(const MatrixType & A)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     */
compute(const MatrixType & A)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 */
rows()115   Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; }
116   /** \internal */
cols()117   Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; }
118 
119   /** \returns the tolerance threshold used by the stopping criteria */
tolerance()120   RealScalar tolerance() const { return m_tolerance; }
121 
122   /** Sets the tolerance threshold used by the stopping criteria */
setTolerance(const RealScalar & tolerance)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. */
preconditioner()130   Preconditioner& preconditioner() { return m_preconditioner; }
131 
132   /** \returns a read-only reference to the preconditioner. */
preconditioner()133   const Preconditioner& preconditioner() const { return m_preconditioner; }
134 
135   /** \returns the max number of iterations */
maxIterations()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 */
setMaxIterations(int maxIters)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 */
iterations()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 */
error()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>
solve(const MatrixBase<Rhs> & b)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>
solve(const SparseMatrixBase<Rhs> & b)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. */
info()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>
_solve_sparse(const Rhs & b,SparseMatrix<DestScalar,DestOptions,DestIndex> & dest)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 
214 protected:
init()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 
235 namespace internal {
236 
237 template<typename Derived, typename Rhs>
238 struct 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