IterativeSolverBase.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_ITERATIVE_SOLVER_BASE_H
00026 #define EIGEN_ITERATIVE_SOLVER_BASE_H
00027 
00028 namespace Eigen { 
00029 
00035 template< typename Derived>
00036 class IterativeSolverBase : internal::noncopyable
00037 {
00038 public:
00039   typedef typename internal::traits<Derived>::MatrixType MatrixType;
00040   typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
00041   typedef typename MatrixType::Scalar Scalar;
00042   typedef typename MatrixType::Index Index;
00043   typedef typename MatrixType::RealScalar RealScalar;
00044 
00045 public:
00046 
00047   Derived& derived() { return *static_cast<Derived*>(this); }
00048   const Derived& derived() const { return *static_cast<const Derived*>(this); }
00049 
00051   IterativeSolverBase()
00052     : mp_matrix(0)
00053   {
00054     init();
00055   }
00056 
00067   IterativeSolverBase(const MatrixType& A)
00068   {
00069     init();
00070     compute(A);
00071   }
00072 
00073   ~IterativeSolverBase() {}
00074   
00080   Derived& analyzePattern(const MatrixType& A)
00081   {
00082     m_preconditioner.analyzePattern(A);
00083     m_isInitialized = true;
00084     m_analysisIsOk = true;
00085     m_info = Success;
00086     return derived();
00087   }
00088   
00098   Derived& factorize(const MatrixType& A)
00099   {
00100     eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 
00101     mp_matrix = &A;
00102     m_preconditioner.factorize(A);
00103     m_factorizationIsOk = true;
00104     m_info = Success;
00105     return derived();
00106   }
00107 
00118   Derived& compute(const MatrixType& A)
00119   {
00120     mp_matrix = &A;
00121     m_preconditioner.compute(A);
00122     m_isInitialized = true;
00123     m_analysisIsOk = true;
00124     m_factorizationIsOk = true;
00125     m_info = Success;
00126     return derived();
00127   }
00128 
00130   Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; }
00132   Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; }
00133 
00135   RealScalar tolerance() const { return m_tolerance; }
00136   
00138   Derived& setTolerance(RealScalar tolerance)
00139   {
00140     m_tolerance = tolerance;
00141     return derived();
00142   }
00143 
00145   Preconditioner& preconditioner() { return m_preconditioner; }
00146   
00148   const Preconditioner& preconditioner() const { return m_preconditioner; }
00149 
00151   int maxIterations() const
00152   {
00153     return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations;
00154   }
00155   
00157   Derived& setMaxIterations(int maxIters)
00158   {
00159     m_maxIterations = maxIters;
00160     return derived();
00161   }
00162 
00164   int iterations() const
00165   {
00166     eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
00167     return m_iterations;
00168   }
00169 
00171   RealScalar error() const
00172   {
00173     eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
00174     return m_error;
00175   }
00176 
00181   template<typename Rhs> inline const internal::solve_retval<Derived, Rhs>
00182   solve(const MatrixBase<Rhs>& b) const
00183   {
00184     eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
00185     eigen_assert(rows()==b.rows()
00186               && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
00187     return internal::solve_retval<Derived, Rhs>(derived(), b.derived());
00188   }
00189   
00194   template<typename Rhs>
00195   inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs>
00196   solve(const SparseMatrixBase<Rhs>& b) const
00197   {
00198     eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
00199     eigen_assert(rows()==b.rows()
00200               && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
00201     return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived());
00202   }
00203 
00205   ComputationInfo info() const
00206   {
00207     eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
00208     return m_info;
00209   }
00210   
00212   template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
00213   void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00214   {
00215     eigen_assert(rows()==b.rows());
00216     
00217     int rhsCols = b.cols();
00218     int size = b.rows();
00219     Eigen::Matrix<DestScalar,Dynamic,1> tb(size);
00220     Eigen::Matrix<DestScalar,Dynamic,1> tx(size);
00221     for(int k=0; k<rhsCols; ++k)
00222     {
00223       tb = b.col(k);
00224       tx = derived().solve(tb);
00225       dest.col(k) = tx.sparseView(0);
00226     }
00227   }
00228 
00229 protected:
00230   void init()
00231   {
00232     m_isInitialized = false;
00233     m_analysisIsOk = false;
00234     m_factorizationIsOk = false;
00235     m_maxIterations = -1;
00236     m_tolerance = NumTraits<Scalar>::epsilon();
00237   }
00238   const MatrixType* mp_matrix;
00239   Preconditioner m_preconditioner;
00240 
00241   int m_maxIterations;
00242   RealScalar m_tolerance;
00243   
00244   mutable RealScalar m_error;
00245   mutable int m_iterations;
00246   mutable ComputationInfo m_info;
00247   mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk;
00248 };
00249 
00250 namespace internal {
00251  
00252 template<typename Derived, typename Rhs>
00253 struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs>
00254   : sparse_solve_retval_base<IterativeSolverBase<Derived>, Rhs>
00255 {
00256   typedef IterativeSolverBase<Derived> Dec;
00257   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00258 
00259   template<typename Dest> void evalTo(Dest& dst) const
00260   {
00261     dec().derived()._solve_sparse(rhs(),dst);
00262   }
00263 };
00264 
00265 } // end namespace internal
00266 
00267 } // end namespace Eigen
00268 
00269 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H