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