Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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
15namespace Eigen {
16
17namespace internal {
18
19template<typename MatrixType>
20struct is_ref_compatible_impl
21{
22private:
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
37public:
38 static MatrixType ms_from;
39 enum { value = sizeof(test<MatrixType>(ms_from, 0))==sizeof(yes) };
40};
41
42template<typename MatrixType>
43struct is_ref_compatible
44{
45 enum { value = is_ref_compatible_impl<typename remove_all<MatrixType>::type>::value };
46};
47
48template<typename MatrixType, bool MatrixFree = !internal::is_ref_compatible<MatrixType>::value>
49class generic_matrix_wrapper;
50
51// We have an explicit matrix at hand, compatible with Ref<>
52template<typename MatrixType>
53class generic_matrix_wrapper<MatrixType,false>
54{
55public:
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 m_matrix.~Ref<const MatrixType>();
83 ::new (&m_matrix) Ref<const MatrixType>(mat.derived());
84 }
85
86 void grab(const Ref<const MatrixType> &mat)
87 {
88 if(&(mat.derived()) != &m_matrix)
89 {
90 m_matrix.~Ref<const MatrixType>();
91 ::new (&m_matrix) Ref<const MatrixType>(mat);
92 }
93 }
94
95protected:
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
101template<typename MatrixType>
102class generic_matrix_wrapper<MatrixType,true>
103{
104public:
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
133protected:
134 const ActualMatrixType *mp_matrix;
135};
136
137}
138
144template< typename Derived>
146{
147protected:
149 using Base::m_isInitialized;
150
151public:
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
163public:
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
192
198 template<typename MatrixDerived>
200 {
201 grab(A.derived());
202 m_preconditioner.analyzePattern(matrix());
203 m_isInitialized = true;
204 m_analysisIsOk = true;
205 m_info = m_preconditioner.info();
206 return derived();
207 }
208
218 template<typename MatrixDerived>
220 {
221 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
222 grab(A.derived());
223 m_preconditioner.factorize(matrix());
224 m_factorizationIsOk = true;
225 m_info = m_preconditioner.info();
226 return derived();
227 }
228
239 template<typename MatrixDerived>
241 {
242 grab(A.derived());
243 m_preconditioner.compute(matrix());
244 m_isInitialized = true;
245 m_analysisIsOk = true;
246 m_factorizationIsOk = true;
247 m_info = m_preconditioner.info();
248 return derived();
249 }
250
252 EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return matrix().rows(); }
253
255 EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return matrix().cols(); }
256
260 RealScalar tolerance() const { return m_tolerance; }
261
267 Derived& setTolerance(const RealScalar& tolerance)
268 {
269 m_tolerance = tolerance;
270 return derived();
271 }
272
274 Preconditioner& preconditioner() { return m_preconditioner; }
275
277 const Preconditioner& preconditioner() const { return m_preconditioner; }
278
284 {
285 return (m_maxIterations<0) ? 2*matrix().cols() : m_maxIterations;
286 }
287
291 Derived& setMaxIterations(Index maxIters)
292 {
293 m_maxIterations = maxIters;
294 return derived();
295 }
296
299 {
300 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
301 return m_iterations;
302 }
303
307 RealScalar error() const
308 {
309 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
310 return m_error;
311 }
312
318 template<typename Rhs,typename Guess>
320 solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
321 {
322 eigen_assert(m_isInitialized && "Solver is not initialized.");
323 eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
324 return SolveWithGuess<Derived, Rhs, Guess>(derived(), b.derived(), x0);
325 }
326
329 {
330 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
331 return m_info;
332 }
333
335 template<typename Rhs, typename DestDerived>
336 void _solve_with_guess_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const
337 {
338 eigen_assert(rows()==b.rows());
339
340 Index rhsCols = b.cols();
341 Index size = b.rows();
342 DestDerived& dest(aDest.derived());
343 typedef typename DestDerived::Scalar DestScalar;
346 // We do not directly fill dest because sparse expressions have to be free of aliasing issue.
347 // For non square least-square problems, b and dest might not have the same size whereas they might alias each-other.
348 typename DestDerived::PlainObject tmp(cols(),rhsCols);
349 ComputationInfo global_info = Success;
350 for(Index k=0; k<rhsCols; ++k)
351 {
352 tb = b.col(k);
353 tx = dest.col(k);
354 derived()._solve_vector_with_guess_impl(tb,tx);
355 tmp.col(k) = tx.sparseView(0);
356
357 // The call to _solve_vector_with_guess_impl updates m_info, so if it failed for a previous column
358 // we need to restore it to the worst value.
359 if(m_info==NumericalIssue)
360 global_info = NumericalIssue;
361 else if(m_info==NoConvergence)
362 global_info = NoConvergence;
363 }
364 m_info = global_info;
365 dest.swap(tmp);
366 }
367
368 template<typename Rhs, typename DestDerived>
369 typename internal::enable_if<Rhs::ColsAtCompileTime!=1 && DestDerived::ColsAtCompileTime!=1>::type
370 _solve_with_guess_impl(const Rhs& b, MatrixBase<DestDerived> &aDest) const
371 {
372 eigen_assert(rows()==b.rows());
373
374 Index rhsCols = b.cols();
375 DestDerived& dest(aDest.derived());
376 ComputationInfo global_info = Success;
377 for(Index k=0; k<rhsCols; ++k)
378 {
379 typename DestDerived::ColXpr xk(dest,k);
380 typename Rhs::ConstColXpr bk(b,k);
381 derived()._solve_vector_with_guess_impl(bk,xk);
382
383 // The call to _solve_vector_with_guess updates m_info, so if it failed for a previous column
384 // we need to restore it to the worst value.
385 if(m_info==NumericalIssue)
386 global_info = NumericalIssue;
387 else if(m_info==NoConvergence)
388 global_info = NoConvergence;
389 }
390 m_info = global_info;
391 }
392
393 template<typename Rhs, typename DestDerived>
394 typename internal::enable_if<Rhs::ColsAtCompileTime==1 || DestDerived::ColsAtCompileTime==1>::type
395 _solve_with_guess_impl(const Rhs& b, MatrixBase<DestDerived> &dest) const
396 {
397 derived()._solve_vector_with_guess_impl(b,dest.derived());
398 }
399
401 template<typename Rhs,typename Dest>
402 void _solve_impl(const Rhs& b, Dest& x) const
403 {
404 x.setZero();
405 derived()._solve_with_guess_impl(b,x);
406 }
407
408protected:
409 void init()
410 {
411 m_isInitialized = false;
412 m_analysisIsOk = false;
413 m_factorizationIsOk = false;
414 m_maxIterations = -1;
415 m_tolerance = NumTraits<Scalar>::epsilon();
416 }
417
418 typedef internal::generic_matrix_wrapper<MatrixType> MatrixWrapper;
419 typedef typename MatrixWrapper::ActualMatrixType ActualMatrixType;
420
421 const ActualMatrixType& matrix() const
422 {
423 return m_matrixWrapper.matrix();
424 }
425
426 template<typename InputType>
427 void grab(const InputType &A)
428 {
429 m_matrixWrapper.grab(A);
430 }
431
432 MatrixWrapper m_matrixWrapper;
433 Preconditioner m_preconditioner;
434
435 Index m_maxIterations;
436 RealScalar m_tolerance;
437
438 mutable RealScalar m_error;
439 mutable Index m_iterations;
440 mutable ComputationInfo m_info;
441 mutable bool m_analysisIsOk, m_factorizationIsOk;
442};
443
444} // end namespace Eigen
445
446#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:328
RealScalar error() const
Definition: IterativeSolverBase.h:307
Index maxIterations() const
Definition: IterativeSolverBase.h:283
Derived & setMaxIterations(Index maxIters)
Definition: IterativeSolverBase.h:291
Derived & compute(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:240
IterativeSolverBase(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:184
const Preconditioner & preconditioner() const
Definition: IterativeSolverBase.h:277
Derived & analyzePattern(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:199
Derived & factorize(const EigenBase< MatrixDerived > &A)
Definition: IterativeSolverBase.h:219
Preconditioner & preconditioner()
Definition: IterativeSolverBase.h:274
Derived & setTolerance(const RealScalar &tolerance)
Definition: IterativeSolverBase.h:267
RealScalar tolerance() const
Definition: IterativeSolverBase.h:260
const SolveWithGuess< Derived, Rhs, Guess > solveWithGuess(const MatrixBase< Rhs > &b, const Guess &x0) const
Definition: IterativeSolverBase.h:320
Index iterations() const
Definition: IterativeSolverBase.h:298
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: B01_Experimental.dox:1
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