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 (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 namespace internal { 16 17 template<typename MatrixType> 18 struct is_ref_compatible_impl 19 { 20 private: 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 35 public: 36 static MatrixType ms_from; 37 enum { value = sizeof(test<MatrixType>(ms_from, 0))==sizeof(yes) }; 38 }; 39 40 template<typename MatrixType> 41 struct is_ref_compatible 42 { 43 enum { value = is_ref_compatible_impl<typename remove_all<MatrixType>::type>::value }; 44 }; 45 46 template<typename MatrixType, bool MatrixFree = !internal::is_ref_compatible<MatrixType>::value> 47 class generic_matrix_wrapper; 48 49 // We have an explicit matrix at hand, compatible with Ref<> 50 template<typename MatrixType> 51 class generic_matrix_wrapper<MatrixType,false> 52 { 53 public: 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 93 protected: 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 99 template<typename MatrixType> 100 class generic_matrix_wrapper<MatrixType,true> 101 { 102 public: 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 131 protected: 132 const ActualMatrixType *mp_matrix; 133 }; 134 135 } 136 137 /** \ingroup IterativeLinearSolvers_Module 138 * \brief Base class for linear iterative solvers 139 * 140 * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner 141 */ 142 template< typename Derived> 143 class IterativeSolverBase : public SparseSolverBase<Derived> 144 { 145 protected: 146 typedef SparseSolverBase<Derived> Base; 147 using Base::m_isInitialized; 148 149 public: 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 161 public: 162 163 using Base::derived; 164 165 /** Default constructor. */ 166 IterativeSolverBase() 167 { 168 init(); 169 } 170 171 /** Initialize the solver with matrix \a A for further \c Ax=b solving. 172 * 173 * This constructor is a shortcut for the default constructor followed 174 * by a call to compute(). 175 * 176 * \warning this class stores a reference to the matrix A as well as some 177 * precomputed values that depend on it. Therefore, if \a A is changed 178 * this class becomes invalid. Call compute() to update it with the new 179 * matrix A, or modify a copy of A. 180 */ 181 template<typename MatrixDerived> 182 explicit IterativeSolverBase(const EigenBase<MatrixDerived>& A) 183 : m_matrixWrapper(A.derived()) 184 { 185 init(); 186 compute(matrix()); 187 } 188 189 ~IterativeSolverBase() {} 190 191 /** Initializes the iterative solver for the sparsity pattern of the matrix \a A for further solving \c Ax=b problems. 192 * 193 * Currently, this function mostly calls analyzePattern on the preconditioner. In the future 194 * we might, for instance, implement column reordering for faster matrix vector products. 195 */ 196 template<typename MatrixDerived> 197 Derived& analyzePattern(const EigenBase<MatrixDerived>& A) 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 207 /** Initializes the iterative solver with the numerical values of the matrix \a A for further solving \c Ax=b problems. 208 * 209 * Currently, this function mostly calls factorize on the preconditioner. 210 * 211 * \warning this class stores a reference to the matrix A as well as some 212 * precomputed values that depend on it. Therefore, if \a A is changed 213 * this class becomes invalid. Call compute() to update it with the new 214 * matrix A, or modify a copy of A. 215 */ 216 template<typename MatrixDerived> 217 Derived& factorize(const EigenBase<MatrixDerived>& A) 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 227 /** Initializes the iterative solver with the matrix \a A for further solving \c Ax=b problems. 228 * 229 * Currently, this function mostly initializes/computes the preconditioner. In the future 230 * we might, for instance, implement column reordering for faster matrix vector products. 231 * 232 * \warning this class stores a reference to the matrix A as well as some 233 * precomputed values that depend on it. Therefore, if \a A is changed 234 * this class becomes invalid. Call compute() to update it with the new 235 * matrix A, or modify a copy of A. 236 */ 237 template<typename MatrixDerived> 238 Derived& compute(const EigenBase<MatrixDerived>& A) 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 249 /** \internal */ 250 Index rows() const { return matrix().rows(); } 251 252 /** \internal */ 253 Index cols() const { return matrix().cols(); } 254 255 /** \returns the tolerance threshold used by the stopping criteria. 256 * \sa setTolerance() 257 */ 258 RealScalar tolerance() const { return m_tolerance; } 259 260 /** Sets the tolerance threshold used by the stopping criteria. 261 * 262 * This value is used as an upper bound to the relative residual error: |Ax-b|/|b|. 263 * The default value is the machine precision given by NumTraits<Scalar>::epsilon() 264 */ 265 Derived& setTolerance(const RealScalar& tolerance) 266 { 267 m_tolerance = tolerance; 268 return derived(); 269 } 270 271 /** \returns a read-write reference to the preconditioner for custom configuration. */ 272 Preconditioner& preconditioner() { return m_preconditioner; } 273 274 /** \returns a read-only reference to the preconditioner. */ 275 const Preconditioner& preconditioner() const { return m_preconditioner; } 276 277 /** \returns the max number of iterations. 278 * It is either the value setted by setMaxIterations or, by default, 279 * twice the number of columns of the matrix. 280 */ 281 Index maxIterations() const 282 { 283 return (m_maxIterations<0) ? 2*matrix().cols() : m_maxIterations; 284 } 285 286 /** Sets the max number of iterations. 287 * Default is twice the number of columns of the matrix. 288 */ 289 Derived& setMaxIterations(Index maxIters) 290 { 291 m_maxIterations = maxIters; 292 return derived(); 293 } 294 295 /** \returns the number of iterations performed during the last solve */ 296 Index iterations() const 297 { 298 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); 299 return m_iterations; 300 } 301 302 /** \returns the tolerance error reached during the last solve. 303 * It is a close approximation of the true relative residual error |Ax-b|/|b|. 304 */ 305 RealScalar error() const 306 { 307 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); 308 return m_error; 309 } 310 311 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A 312 * and \a x0 as an initial solution. 313 * 314 * \sa solve(), compute() 315 */ 316 template<typename Rhs,typename Guess> 317 inline const SolveWithGuess<Derived, Rhs, 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 325 /** \returns Success if the iterations converged, and NoConvergence otherwise. */ 326 ComputationInfo info() const 327 { 328 eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); 329 return m_info; 330 } 331 332 /** \internal */ 333 template<typename Rhs, typename DestDerived> 334 void _solve_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; 342 Eigen::Matrix<DestScalar,Dynamic,1> tb(size); 343 Eigen::Matrix<DestScalar,Dynamic,1> tx(cols()); 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 for(Index k=0; k<rhsCols; ++k) 348 { 349 tb = b.col(k); 350 tx = derived().solve(tb); 351 tmp.col(k) = tx.sparseView(0); 352 } 353 dest.swap(tmp); 354 } 355 356 protected: 357 void init() 358 { 359 m_isInitialized = false; 360 m_analysisIsOk = false; 361 m_factorizationIsOk = false; 362 m_maxIterations = -1; 363 m_tolerance = NumTraits<Scalar>::epsilon(); 364 } 365 366 typedef internal::generic_matrix_wrapper<MatrixType> MatrixWrapper; 367 typedef typename MatrixWrapper::ActualMatrixType ActualMatrixType; 368 369 const ActualMatrixType& matrix() const 370 { 371 return m_matrixWrapper.matrix(); 372 } 373 374 template<typename InputType> 375 void grab(const InputType &A) 376 { 377 m_matrixWrapper.grab(A); 378 } 379 380 MatrixWrapper m_matrixWrapper; 381 Preconditioner m_preconditioner; 382 383 Index m_maxIterations; 384 RealScalar m_tolerance; 385 386 mutable RealScalar m_error; 387 mutable Index m_iterations; 388 mutable ComputationInfo m_info; 389 mutable bool m_analysisIsOk, m_factorizationIsOk; 390 }; 391 392 } // end namespace Eigen 393 394 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H 395
