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 |
|
---|
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 | template<typename InputDerived>
|
---|
53 | IterativeSolverBase(const EigenBase<InputDerived>& A)
|
---|
54 | {
|
---|
55 | init();
|
---|
56 | compute(A.derived());
|
---|
57 | }
|
---|
58 |
|
---|
59 | ~IterativeSolverBase() {}
|
---|
60 |
|
---|
61 | /** Initializes the iterative solver for the sparcity pattern of the matrix \a A for further solving \c Ax=b problems.
|
---|
62 | *
|
---|
63 | * Currently, this function mostly call analyzePattern on the preconditioner. In the future
|
---|
64 | * we might, for instance, implement column reodering for faster matrix vector products.
|
---|
65 | */
|
---|
66 | template<typename InputDerived>
|
---|
67 | Derived& analyzePattern(const EigenBase<InputDerived>& A)
|
---|
68 | {
|
---|
69 | grabInput(A.derived());
|
---|
70 | m_preconditioner.analyzePattern(*mp_matrix);
|
---|
71 | m_isInitialized = true;
|
---|
72 | m_analysisIsOk = true;
|
---|
73 | m_info = Success;
|
---|
74 | return derived();
|
---|
75 | }
|
---|
76 |
|
---|
77 | /** Initializes the iterative solver with the numerical values of the matrix \a A for further solving \c Ax=b problems.
|
---|
78 | *
|
---|
79 | * Currently, this function mostly call factorize on the preconditioner.
|
---|
80 | *
|
---|
81 | * \warning this class stores a reference to the matrix A as well as some
|
---|
82 | * precomputed values that depend on it. Therefore, if \a A is changed
|
---|
83 | * this class becomes invalid. Call compute() to update it with the new
|
---|
84 | * matrix A, or modify a copy of A.
|
---|
85 | */
|
---|
86 | template<typename InputDerived>
|
---|
87 | Derived& factorize(const EigenBase<InputDerived>& A)
|
---|
88 | {
|
---|
89 | grabInput(A.derived());
|
---|
90 | eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
|
---|
91 | m_preconditioner.factorize(*mp_matrix);
|
---|
92 | m_factorizationIsOk = true;
|
---|
93 | m_info = Success;
|
---|
94 | return derived();
|
---|
95 | }
|
---|
96 |
|
---|
97 | /** Initializes the iterative solver with the matrix \a A for further solving \c Ax=b problems.
|
---|
98 | *
|
---|
99 | * Currently, this function mostly initialized/compute the preconditioner. In the future
|
---|
100 | * we might, for instance, implement column reodering for faster matrix vector products.
|
---|
101 | *
|
---|
102 | * \warning this class stores a reference to the matrix A as well as some
|
---|
103 | * precomputed values that depend on it. Therefore, if \a A is changed
|
---|
104 | * this class becomes invalid. Call compute() to update it with the new
|
---|
105 | * matrix A, or modify a copy of A.
|
---|
106 | */
|
---|
107 | template<typename InputDerived>
|
---|
108 | Derived& compute(const EigenBase<InputDerived>& A)
|
---|
109 | {
|
---|
110 | grabInput(A.derived());
|
---|
111 | m_preconditioner.compute(*mp_matrix);
|
---|
112 | m_isInitialized = true;
|
---|
113 | m_analysisIsOk = true;
|
---|
114 | m_factorizationIsOk = true;
|
---|
115 | m_info = Success;
|
---|
116 | return derived();
|
---|
117 | }
|
---|
118 |
|
---|
119 | /** \internal */
|
---|
120 | Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; }
|
---|
121 | /** \internal */
|
---|
122 | Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; }
|
---|
123 |
|
---|
124 | /** \returns the tolerance threshold used by the stopping criteria */
|
---|
125 | RealScalar tolerance() const { return m_tolerance; }
|
---|
126 |
|
---|
127 | /** Sets the tolerance threshold used by the stopping criteria */
|
---|
128 | Derived& setTolerance(const RealScalar& tolerance)
|
---|
129 | {
|
---|
130 | m_tolerance = tolerance;
|
---|
131 | return derived();
|
---|
132 | }
|
---|
133 |
|
---|
134 | /** \returns a read-write reference to the preconditioner for custom configuration. */
|
---|
135 | Preconditioner& preconditioner() { return m_preconditioner; }
|
---|
136 |
|
---|
137 | /** \returns a read-only reference to the preconditioner. */
|
---|
138 | const Preconditioner& preconditioner() const { return m_preconditioner; }
|
---|
139 |
|
---|
140 | /** \returns the max number of iterations */
|
---|
141 | int maxIterations() const
|
---|
142 | {
|
---|
143 | return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations;
|
---|
144 | }
|
---|
145 |
|
---|
146 | /** Sets the max number of iterations */
|
---|
147 | Derived& setMaxIterations(int maxIters)
|
---|
148 | {
|
---|
149 | m_maxIterations = maxIters;
|
---|
150 | return derived();
|
---|
151 | }
|
---|
152 |
|
---|
153 | /** \returns the number of iterations performed during the last solve */
|
---|
154 | int iterations() const
|
---|
155 | {
|
---|
156 | eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
|
---|
157 | return m_iterations;
|
---|
158 | }
|
---|
159 |
|
---|
160 | /** \returns the tolerance error reached during the last solve */
|
---|
161 | RealScalar error() const
|
---|
162 | {
|
---|
163 | eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
|
---|
164 | return m_error;
|
---|
165 | }
|
---|
166 |
|
---|
167 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
|
---|
168 | *
|
---|
169 | * \sa compute()
|
---|
170 | */
|
---|
171 | template<typename Rhs> inline const internal::solve_retval<Derived, Rhs>
|
---|
172 | solve(const MatrixBase<Rhs>& b) const
|
---|
173 | {
|
---|
174 | eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
|
---|
175 | eigen_assert(rows()==b.rows()
|
---|
176 | && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
|
---|
177 | return internal::solve_retval<Derived, Rhs>(derived(), b.derived());
|
---|
178 | }
|
---|
179 |
|
---|
180 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
|
---|
181 | *
|
---|
182 | * \sa compute()
|
---|
183 | */
|
---|
184 | template<typename Rhs>
|
---|
185 | inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs>
|
---|
186 | solve(const SparseMatrixBase<Rhs>& b) const
|
---|
187 | {
|
---|
188 | eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
|
---|
189 | eigen_assert(rows()==b.rows()
|
---|
190 | && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
|
---|
191 | return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived());
|
---|
192 | }
|
---|
193 |
|
---|
194 | /** \returns Success if the iterations converged, and NoConvergence otherwise. */
|
---|
195 | ComputationInfo info() const
|
---|
196 | {
|
---|
197 | eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
|
---|
198 | return m_info;
|
---|
199 | }
|
---|
200 |
|
---|
201 | /** \internal */
|
---|
202 | template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
|
---|
203 | void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
|
---|
204 | {
|
---|
205 | eigen_assert(rows()==b.rows());
|
---|
206 |
|
---|
207 | int rhsCols = b.cols();
|
---|
208 | int size = b.rows();
|
---|
209 | Eigen::Matrix<DestScalar,Dynamic,1> tb(size);
|
---|
210 | Eigen::Matrix<DestScalar,Dynamic,1> tx(size);
|
---|
211 | for(int k=0; k<rhsCols; ++k)
|
---|
212 | {
|
---|
213 | tb = b.col(k);
|
---|
214 | tx = derived().solve(tb);
|
---|
215 | dest.col(k) = tx.sparseView(0);
|
---|
216 | }
|
---|
217 | }
|
---|
218 |
|
---|
219 | protected:
|
---|
220 |
|
---|
221 | template<typename InputDerived>
|
---|
222 | void grabInput(const EigenBase<InputDerived>& A)
|
---|
223 | {
|
---|
224 | // we const cast to prevent the creation of a MatrixType temporary by the compiler.
|
---|
225 | grabInput_impl(A.const_cast_derived());
|
---|
226 | }
|
---|
227 |
|
---|
228 | template<typename InputDerived>
|
---|
229 | void grabInput_impl(const EigenBase<InputDerived>& A)
|
---|
230 | {
|
---|
231 | m_copyMatrix = A;
|
---|
232 | mp_matrix = &m_copyMatrix;
|
---|
233 | }
|
---|
234 |
|
---|
235 | void grabInput_impl(MatrixType& A)
|
---|
236 | {
|
---|
237 | if(MatrixType::RowsAtCompileTime==Dynamic && MatrixType::ColsAtCompileTime==Dynamic)
|
---|
238 | m_copyMatrix.resize(0,0);
|
---|
239 | mp_matrix = &A;
|
---|
240 | }
|
---|
241 |
|
---|
242 | void init()
|
---|
243 | {
|
---|
244 | m_isInitialized = false;
|
---|
245 | m_analysisIsOk = false;
|
---|
246 | m_factorizationIsOk = false;
|
---|
247 | m_maxIterations = -1;
|
---|
248 | m_tolerance = NumTraits<Scalar>::epsilon();
|
---|
249 | }
|
---|
250 | MatrixType m_copyMatrix;
|
---|
251 | const MatrixType* mp_matrix;
|
---|
252 | Preconditioner m_preconditioner;
|
---|
253 |
|
---|
254 | int m_maxIterations;
|
---|
255 | RealScalar m_tolerance;
|
---|
256 |
|
---|
257 | mutable RealScalar m_error;
|
---|
258 | mutable int m_iterations;
|
---|
259 | mutable ComputationInfo m_info;
|
---|
260 | mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk;
|
---|
261 | };
|
---|
262 |
|
---|
263 | namespace internal {
|
---|
264 |
|
---|
265 | template<typename Derived, typename Rhs>
|
---|
266 | struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs>
|
---|
267 | : sparse_solve_retval_base<IterativeSolverBase<Derived>, Rhs>
|
---|
268 | {
|
---|
269 | typedef IterativeSolverBase<Derived> Dec;
|
---|
270 | EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
|
---|
271 |
|
---|
272 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
273 | {
|
---|
274 | dec().derived()._solve_sparse(rhs(),dst);
|
---|
275 | }
|
---|
276 | };
|
---|
277 |
|
---|
278 | } // end namespace internal
|
---|
279 |
|
---|
280 | } // end namespace Eigen
|
---|
281 |
|
---|
282 | #endif // EIGEN_ITERATIVE_SOLVER_BASE_H
|
---|