Eigen  3.3.3
IterativeSolverBase.h
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends