$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2011 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 00020 template< typename Derived> 00021 class IterativeSolverBase : internal::noncopyable 00022 { 00023 public: 00024 typedef typename internal::traits<Derived>::MatrixType MatrixType; 00025 typedef typename internal::traits<Derived>::Preconditioner Preconditioner; 00026 typedef typename MatrixType::Scalar Scalar; 00027 typedef typename MatrixType::Index Index; 00028 typedef typename MatrixType::RealScalar RealScalar; 00029 00030 public: 00031 00032 Derived& derived() { return *static_cast<Derived*>(this); } 00033 const Derived& derived() const { return *static_cast<const Derived*>(this); } 00034 00036 IterativeSolverBase() 00037 : mp_matrix(0) 00038 { 00039 init(); 00040 } 00041 00052 IterativeSolverBase(const MatrixType& A) 00053 { 00054 init(); 00055 compute(A); 00056 } 00057 00058 ~IterativeSolverBase() {} 00059 00065 Derived& analyzePattern(const MatrixType& A) 00066 { 00067 m_preconditioner.analyzePattern(A); 00068 m_isInitialized = true; 00069 m_analysisIsOk = true; 00070 m_info = Success; 00071 return derived(); 00072 } 00073 00083 Derived& factorize(const MatrixType& A) 00084 { 00085 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 00086 mp_matrix = &A; 00087 m_preconditioner.factorize(A); 00088 m_factorizationIsOk = true; 00089 m_info = Success; 00090 return derived(); 00091 } 00092 00103 Derived& compute(const MatrixType& A) 00104 { 00105 mp_matrix = &A; 00106 m_preconditioner.compute(A); 00107 m_isInitialized = true; 00108 m_analysisIsOk = true; 00109 m_factorizationIsOk = true; 00110 m_info = Success; 00111 return derived(); 00112 } 00113 00115 Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; } 00117 Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; } 00118 00120 RealScalar tolerance() const { return m_tolerance; } 00121 00123 Derived& setTolerance(const RealScalar& tolerance) 00124 { 00125 m_tolerance = tolerance; 00126 return derived(); 00127 } 00128 00130 Preconditioner& preconditioner() { return m_preconditioner; } 00131 00133 const Preconditioner& preconditioner() const { return m_preconditioner; } 00134 00136 int maxIterations() const 00137 { 00138 return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations; 00139 } 00140 00142 Derived& setMaxIterations(int maxIters) 00143 { 00144 m_maxIterations = maxIters; 00145 return derived(); 00146 } 00147 00149 int iterations() const 00150 { 00151 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); 00152 return m_iterations; 00153 } 00154 00156 RealScalar error() const 00157 { 00158 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); 00159 return m_error; 00160 } 00161 00166 template<typename Rhs> inline const internal::solve_retval<Derived, Rhs> 00167 solve(const MatrixBase<Rhs>& b) const 00168 { 00169 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); 00170 eigen_assert(rows()==b.rows() 00171 && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b"); 00172 return internal::solve_retval<Derived, Rhs>(derived(), b.derived()); 00173 } 00174 00179 template<typename Rhs> 00180 inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs> 00181 solve(const SparseMatrixBase<Rhs>& b) const 00182 { 00183 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); 00184 eigen_assert(rows()==b.rows() 00185 && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b"); 00186 return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived()); 00187 } 00188 00190 ComputationInfo info() const 00191 { 00192 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); 00193 return m_info; 00194 } 00195 00197 template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> 00198 void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const 00199 { 00200 eigen_assert(rows()==b.rows()); 00201 00202 int rhsCols = b.cols(); 00203 int size = b.rows(); 00204 Eigen::Matrix<DestScalar,Dynamic,1> tb(size); 00205 Eigen::Matrix<DestScalar,Dynamic,1> tx(size); 00206 for(int k=0; k<rhsCols; ++k) 00207 { 00208 tb = b.col(k); 00209 tx = derived().solve(tb); 00210 dest.col(k) = tx.sparseView(0); 00211 } 00212 } 00213 00214 protected: 00215 void init() 00216 { 00217 m_isInitialized = false; 00218 m_analysisIsOk = false; 00219 m_factorizationIsOk = false; 00220 m_maxIterations = -1; 00221 m_tolerance = NumTraits<Scalar>::epsilon(); 00222 } 00223 const MatrixType* mp_matrix; 00224 Preconditioner m_preconditioner; 00225 00226 int m_maxIterations; 00227 RealScalar m_tolerance; 00228 00229 mutable RealScalar m_error; 00230 mutable int m_iterations; 00231 mutable ComputationInfo m_info; 00232 mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk; 00233 }; 00234 00235 namespace internal { 00236 00237 template<typename Derived, typename Rhs> 00238 struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs> 00239 : sparse_solve_retval_base<IterativeSolverBase<Derived>, Rhs> 00240 { 00241 typedef IterativeSolverBase<Derived> Dec; 00242 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 00243 00244 template<typename Dest> void evalTo(Dest& dst) const 00245 { 00246 dec().derived()._solve_sparse(rhs(),dst); 00247 } 00248 }; 00249 00250 } // end namespace internal 00251 00252 } // end namespace Eigen 00253 00254 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H