Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
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
13namespace Eigen {
14
15namespace internal {
16
17template<typename MatrixType>
18struct is_ref_compatible_impl
19{
20private:
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
35public:
36 static MatrixType ms_from;
37 enum { value = sizeof(test<MatrixType>(ms_from, 0))==sizeof(yes) };
38};
39
40template<typename MatrixType>
41struct is_ref_compatible
42{
43 enum { value = is_ref_compatible_impl<typename remove_all<MatrixType>::type>::value };
44};
45
46template<typename MatrixType, bool MatrixFree = !internal::is_ref_compatible<MatrixType>::value>
47class generic_matrix_wrapper;
48
49// We have an explicit matrix at hand, compatible with Ref<>
50template<typename MatrixType>
51class generic_matrix_wrapper<MatrixType,false>
52{
53public:
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
93protected:
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
99template<typename MatrixType>
100class generic_matrix_wrapper<MatrixType,true>
101{
102public:
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
131protected:
132 const ActualMatrixType *mp_matrix;
133};
134
135}
136
142template< typename Derived>
144{
145protected:
147 using Base::m_isInitialized;
148
149public:
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
161public:
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
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 EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return matrix().rows(); }
251
253 EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { 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
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
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 DestDerived>
334 void _solve_with_guess_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const
335 {
336 eigen_assert(rows()==b.rows());
337
338 Index rhsCols = b.cols();
339 Index size = b.rows();
340 DestDerived& dest(aDest.derived());
341 typedef typename DestDerived::Scalar DestScalar;
344 // We do not directly fill dest because sparse expressions have to be free of aliasing issue.
345 // For non square least-square problems, b and dest might not have the same size whereas they might alias each-other.
346 typename DestDerived::PlainObject tmp(cols(),rhsCols);
347 ComputationInfo global_info = Success;
348 for(Index k=0; k<rhsCols; ++k)
349 {
350 tb = b.col(k);
351 tx = dest.col(k);
352 derived()._solve_vector_with_guess_impl(tb,tx);
353 tmp.col(k) = tx.sparseView(0);
354
355 // The call to _solve_vector_with_guess_impl updates m_info, so if it failed for a previous column
356 // we need to restore it to the worst value.
357 if(m_info==NumericalIssue)
358 global_info = NumericalIssue;
359 else if(m_info==NoConvergence)
360 global_info = NoConvergence;
361 }
362 m_info = global_info;
363 dest.swap(tmp);
364 }
365
366 template<typename Rhs, typename DestDerived>
367 typename internal::enable_if<Rhs::ColsAtCompileTime!=1 && DestDerived::ColsAtCompileTime!=1>::type
368 _solve_with_guess_impl(const Rhs& b, MatrixBase<DestDerived> &aDest) const
369 {
370 eigen_assert(rows()==b.rows());
371
372 Index rhsCols = b.cols();
373 DestDerived& dest(aDest.derived());
374 ComputationInfo global_info = Success;
375 for(Index k=0; k<rhsCols; ++k)
376 {
377 typename DestDerived::ColXpr xk(dest,k);
378 typename Rhs::ConstColXpr bk(b,k);
379 derived()._solve_vector_with_guess_impl(bk,xk);
380
381 // The call to _solve_vector_with_guess updates m_info, so if it failed for a previous column
382 // we need to restore it to the worst value.
383 if(m_info==NumericalIssue)
384 global_info = NumericalIssue;
385 else if(m_info==NoConvergence)
386 global_info = NoConvergence;
387 }
388 m_info = global_info;
389 }
390
391 template<typename Rhs, typename DestDerived>
392 typename internal::enable_if<Rhs::ColsAtCompileTime==1 || DestDerived::ColsAtCompileTime==1>::type
393 _solve_with_guess_impl(const Rhs& b, MatrixBase<DestDerived> &dest) const
394 {
395 derived()._solve_vector_with_guess_impl(b,dest.derived());
396 }
397
399 template<typename Rhs,typename Dest>
400 void _solve_impl(const Rhs& b, Dest& x) const
401 {
402 x.setZero();
403 derived()._solve_with_guess_impl(b,x);
404 }
405
406protected:
407 void init()
408 {
409 m_isInitialized = false;
410 m_analysisIsOk = false;
411 m_factorizationIsOk = false;
412 m_maxIterations = -1;
413 m_tolerance = NumTraits<Scalar>::epsilon();
414 }
415
416 typedef internal::generic_matrix_wrapper<MatrixType> MatrixWrapper;
417 typedef typename MatrixWrapper::ActualMatrixType ActualMatrixType;
418
419 const ActualMatrixType& matrix() const
420 {
421 return m_matrixWrapper.matrix();
422 }
423
424 template<typename InputType>
425 void grab(const InputType &A)
426 {
427 m_matrixWrapper.grab(A);
428 }
429
430 MatrixWrapper m_matrixWrapper;
431 Preconditioner m_preconditioner;
432
433 Index m_maxIterations;
434 RealScalar m_tolerance;
435
436 mutable RealScalar m_error;
437 mutable Index m_iterations;
438 mutable ComputationInfo m_info;
439 mutable bool m_analysisIsOk, m_factorizationIsOk;
440};
441
442} // end namespace Eigen
443
444#endif // EIGEN_ITERATIVE_SOLVER_BASE_H
Derived & derived()
Definition: EigenBase.h:46
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:144
IterativeSolverBase()
Definition: IterativeSolverBase.h:166
ComputationInfo info() const
Definition: IterativeSolverBase.h:326
RealScalar error() const
Definition: IterativeSolverBase.h:305
Index maxIterations() const
Definition: IterativeSolverBase.h:281
Derived & setMaxIterations(Index maxIters)
Definition: IterativeSolverBase.h:289
Derived & compute(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:238
IterativeSolverBase(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:182
const Preconditioner & preconditioner() const
Definition: IterativeSolverBase.h:275
Derived & analyzePattern(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:197
Derived & factorize(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:217
Preconditioner & preconditioner()
Definition: IterativeSolverBase.h:272
Derived & setTolerance(const RealScalar &tolerance)
Definition: IterativeSolverBase.h:265
RealScalar tolerance() const
Definition: IterativeSolverBase.h:258
const SolveWithGuess< Derived, Rhs, Guess > solveWithGuess(const MatrixBase< Rhs > &b, const Guess &x0) const
Definition: IterativeSolverBase.h:318
Index iterations() const
Definition: IterativeSolverBase.h:296
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Pseudo expression representing a solving operation.
Definition: SolveWithGuess.h:42
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:28
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
ComputationInfo
Definition: Constants.h:440
@ NumericalIssue
Definition: Constants.h:444
@ Success
Definition: Constants.h:442
@ NoConvergence
Definition: Constants.h:446
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: EigenBase.h:30
Derived & derived()
Definition: EigenBase.h:46