![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_ITERATIVE_SOLVER_BASE_H 00011 #define EIGEN_ITERATIVE_SOLVER_BASE_H 00012 00013 namespace Eigen { 00014 00015 namespace internal { 00016 00017 template<typename MatrixType> 00018 struct is_ref_compatible_impl 00019 { 00020 private: 00021 template <typename T0> 00022 struct any_conversion 00023 { 00024 template <typename T> any_conversion(const volatile T&); 00025 template <typename T> any_conversion(T&); 00026 }; 00027 struct yes {int a[1];}; 00028 struct no {int a[2];}; 00029 00030 template<typename T> 00031 static yes test(const Ref<const T>&, int); 00032 template<typename T> 00033 static no test(any_conversion<T>, ...); 00034 00035 public: 00036 static MatrixType ms_from; 00037 enum { value = sizeof(test<MatrixType>(ms_from, 0))==sizeof(yes) }; 00038 }; 00039 00040 template<typename MatrixType> 00041 struct is_ref_compatible 00042 { 00043 enum { value = is_ref_compatible_impl<typename remove_all<MatrixType>::type>::value }; 00044 }; 00045 00046 template<typename MatrixType, bool MatrixFree = !internal::is_ref_compatible<MatrixType>::value> 00047 class generic_matrix_wrapper; 00048 00049 // We have an explicit matrix at hand, compatible with Ref<> 00050 template<typename MatrixType> 00051 class generic_matrix_wrapper<MatrixType,false> 00052 { 00053 public: 00054 typedef Ref<const MatrixType> ActualMatrixType; 00055 template<int UpLo> struct ConstSelfAdjointViewReturnType { 00056 typedef typename ActualMatrixType::template ConstSelfAdjointViewReturnType<UpLo>::Type Type; 00057 }; 00058 00059 enum { 00060 MatrixFree = false 00061 }; 00062 00063 generic_matrix_wrapper() 00064 : m_dummy(0,0), m_matrix(m_dummy) 00065 {} 00066 00067 template<typename InputType> 00068 generic_matrix_wrapper(const InputType &mat) 00069 : m_matrix(mat) 00070 {} 00071 00072 const ActualMatrixType& matrix() const 00073 { 00074 return m_matrix; 00075 } 00076 00077 template<typename MatrixDerived> 00078 void grab(const EigenBase<MatrixDerived> &mat) 00079 { 00080 m_matrix.~Ref<const MatrixType>(); 00081 ::new (&m_matrix) Ref<const MatrixType>(mat.derived()); 00082 } 00083 00084 void grab(const Ref<const MatrixType> &mat) 00085 { 00086 if(&(mat.derived()) != &m_matrix) 00087 { 00088 m_matrix.~Ref<const MatrixType>(); 00089 ::new (&m_matrix) Ref<const MatrixType>(mat); 00090 } 00091 } 00092 00093 protected: 00094 MatrixType m_dummy; // used to default initialize the Ref<> object 00095 ActualMatrixType m_matrix; 00096 }; 00097 00098 // MatrixType is not compatible with Ref<> -> matrix-free wrapper 00099 template<typename MatrixType> 00100 class generic_matrix_wrapper<MatrixType,true> 00101 { 00102 public: 00103 typedef MatrixType ActualMatrixType; 00104 template<int UpLo> struct ConstSelfAdjointViewReturnType 00105 { 00106 typedef ActualMatrixType Type; 00107 }; 00108 00109 enum { 00110 MatrixFree = true 00111 }; 00112 00113 generic_matrix_wrapper() 00114 : mp_matrix(0) 00115 {} 00116 00117 generic_matrix_wrapper(const MatrixType &mat) 00118 : mp_matrix(&mat) 00119 {} 00120 00121 const ActualMatrixType& matrix() const 00122 { 00123 return *mp_matrix; 00124 } 00125 00126 void grab(const MatrixType &mat) 00127 { 00128 mp_matrix = &mat; 00129 } 00130 00131 protected: 00132 const ActualMatrixType *mp_matrix; 00133 }; 00134 00135 } 00136 00142 template< typename Derived> 00143 class IterativeSolverBase : public SparseSolverBase<Derived> 00144 { 00145 protected: 00146 typedef SparseSolverBase<Derived> Base; 00147 using Base::m_isInitialized; 00148 00149 public: 00150 typedef typename internal::traits<Derived>::MatrixType MatrixType; 00151 typedef typename internal::traits<Derived>::Preconditioner Preconditioner; 00152 typedef typename MatrixType::Scalar Scalar; 00153 typedef typename MatrixType::StorageIndex StorageIndex; 00154 typedef typename MatrixType::RealScalar RealScalar; 00155 00156 enum { 00157 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00158 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00159 }; 00160 00161 public: 00162 00163 using Base::derived; 00164 00166 IterativeSolverBase() 00167 { 00168 init(); 00169 } 00170 00181 template<typename MatrixDerived> 00182 explicit IterativeSolverBase(const EigenBase<MatrixDerived>& A) 00183 : m_matrixWrapper(A.derived()) 00184 { 00185 init(); 00186 compute(matrix()); 00187 } 00188 00189 ~IterativeSolverBase() {} 00190 00196 template<typename MatrixDerived> 00197 Derived& analyzePattern(const EigenBase<MatrixDerived>& A) 00198 { 00199 grab(A.derived()); 00200 m_preconditioner.analyzePattern(matrix()); 00201 m_isInitialized = true; 00202 m_analysisIsOk = true; 00203 m_info = m_preconditioner.info(); 00204 return derived(); 00205 } 00206 00216 template<typename MatrixDerived> 00217 Derived& factorize(const EigenBase<MatrixDerived>& A) 00218 { 00219 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 00220 grab(A.derived()); 00221 m_preconditioner.factorize(matrix()); 00222 m_factorizationIsOk = true; 00223 m_info = m_preconditioner.info(); 00224 return derived(); 00225 } 00226 00237 template<typename MatrixDerived> 00238 Derived& compute(const EigenBase<MatrixDerived>& A) 00239 { 00240 grab(A.derived()); 00241 m_preconditioner.compute(matrix()); 00242 m_isInitialized = true; 00243 m_analysisIsOk = true; 00244 m_factorizationIsOk = true; 00245 m_info = m_preconditioner.info(); 00246 return derived(); 00247 } 00248 00250 Index rows() const { return matrix().rows(); } 00251 00253 Index cols() const { return matrix().cols(); } 00254 00258 RealScalar tolerance() const { return m_tolerance; } 00259 00265 Derived& setTolerance(const RealScalar& tolerance) 00266 { 00267 m_tolerance = tolerance; 00268 return derived(); 00269 } 00270 00272 Preconditioner& preconditioner() { return m_preconditioner; } 00273 00275 const Preconditioner& preconditioner() const { return m_preconditioner; } 00276 00281 Index maxIterations() const 00282 { 00283 return (m_maxIterations<0) ? 2*matrix().cols() : m_maxIterations; 00284 } 00285 00289 Derived& setMaxIterations(Index maxIters) 00290 { 00291 m_maxIterations = maxIters; 00292 return derived(); 00293 } 00294 00296 Index iterations() const 00297 { 00298 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); 00299 return m_iterations; 00300 } 00301 00305 RealScalar error() const 00306 { 00307 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); 00308 return m_error; 00309 } 00310 00316 template<typename Rhs,typename Guess> 00317 inline const SolveWithGuess<Derived, Rhs, Guess> 00318 solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const 00319 { 00320 eigen_assert(m_isInitialized && "Solver is not initialized."); 00321 eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b"); 00322 return SolveWithGuess<Derived, Rhs, Guess>(derived(), b.derived(), x0); 00323 } 00324 00326 ComputationInfo info() const 00327 { 00328 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); 00329 return m_info; 00330 } 00331 00333 template<typename Rhs, typename DestDerived> 00334 void _solve_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const 00335 { 00336 eigen_assert(rows()==b.rows()); 00337 00338 Index rhsCols = b.cols(); 00339 Index size = b.rows(); 00340 DestDerived& dest(aDest.derived()); 00341 typedef typename DestDerived::Scalar DestScalar; 00342 Eigen::Matrix<DestScalar,Dynamic,1> tb(size); 00343 Eigen::Matrix<DestScalar,Dynamic,1> tx(cols()); 00344 // We do not directly fill dest because sparse expressions have to be free of aliasing issue. 00345 // For non square least-square problems, b and dest might not have the same size whereas they might alias each-other. 00346 typename DestDerived::PlainObject tmp(cols(),rhsCols); 00347 for(Index k=0; k<rhsCols; ++k) 00348 { 00349 tb = b.col(k); 00350 tx = derived().solve(tb); 00351 tmp.col(k) = tx.sparseView(0); 00352 } 00353 dest.swap(tmp); 00354 } 00355 00356 protected: 00357 void init() 00358 { 00359 m_isInitialized = false; 00360 m_analysisIsOk = false; 00361 m_factorizationIsOk = false; 00362 m_maxIterations = -1; 00363 m_tolerance = NumTraits<Scalar>::epsilon(); 00364 } 00365 00366 typedef internal::generic_matrix_wrapper<MatrixType> MatrixWrapper; 00367 typedef typename MatrixWrapper::ActualMatrixType ActualMatrixType; 00368 00369 const ActualMatrixType& matrix() const 00370 { 00371 return m_matrixWrapper.matrix(); 00372 } 00373 00374 template<typename InputType> 00375 void grab(const InputType &A) 00376 { 00377 m_matrixWrapper.grab(A); 00378 } 00379 00380 MatrixWrapper m_matrixWrapper; 00381 Preconditioner m_preconditioner; 00382 00383 Index m_maxIterations; 00384 RealScalar m_tolerance; 00385 00386 mutable RealScalar m_error; 00387 mutable Index m_iterations; 00388 mutable ComputationInfo m_info; 00389 mutable bool m_analysisIsOk, m_factorizationIsOk; 00390 }; 00391 00392 } // end namespace Eigen 00393 00394 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H