BasicPreconditioners.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_BASIC_PRECONDITIONERS_H
00026 #define EIGEN_BASIC_PRECONDITIONERS_H
00027 
00028 namespace Eigen { 
00029 
00047 template <typename _Scalar>
00048 class DiagonalPreconditioner
00049 {
00050     typedef _Scalar Scalar;
00051     typedef Matrix<Scalar,Dynamic,1> Vector;
00052     typedef typename Vector::Index Index;
00053 
00054   public:
00055     typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
00056 
00057     DiagonalPreconditioner() : m_isInitialized(false) {}
00058 
00059     template<typename MatrixType>
00060     DiagonalPreconditioner(const MatrixType& mat) : m_invdiag(mat.cols())
00061     {
00062       compute(mat);
00063     }
00064 
00065     Index rows() const { return m_invdiag.size(); }
00066     Index cols() const { return m_invdiag.size(); }
00067     
00068     template<typename MatrixType>
00069     DiagonalPreconditioner& analyzePattern(const MatrixType& )
00070     {
00071       return *this;
00072     }
00073     
00074     template<typename MatrixType>
00075     DiagonalPreconditioner& factorize(const MatrixType& mat)
00076     {
00077       m_invdiag.resize(mat.cols());
00078       for(int j=0; j<mat.outerSize(); ++j)
00079       {
00080         typename MatrixType::InnerIterator it(mat,j);
00081         while(it && it.index()!=j) ++it;
00082         if(it && it.index()==j)
00083           m_invdiag(j) = Scalar(1)/it.value();
00084         else
00085           m_invdiag(j) = 0;
00086       }
00087       m_isInitialized = true;
00088       return *this;
00089     }
00090     
00091     template<typename MatrixType>
00092     DiagonalPreconditioner& compute(const MatrixType& mat)
00093     {
00094       return factorize(mat);
00095     }
00096 
00097     template<typename Rhs, typename Dest>
00098     void _solve(const Rhs& b, Dest& x) const
00099     {
00100       x = m_invdiag.array() * b.array() ;
00101     }
00102 
00103     template<typename Rhs> inline const internal::solve_retval<DiagonalPreconditioner, Rhs>
00104     solve(const MatrixBase<Rhs>& b) const
00105     {
00106       eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
00107       eigen_assert(m_invdiag.size()==b.rows()
00108                 && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
00109       return internal::solve_retval<DiagonalPreconditioner, Rhs>(*this, b.derived());
00110     }
00111 
00112   protected:
00113     Vector m_invdiag;
00114     bool m_isInitialized;
00115 };
00116 
00117 namespace internal {
00118 
00119 template<typename _MatrixType, typename Rhs>
00120 struct solve_retval<DiagonalPreconditioner<_MatrixType>, Rhs>
00121   : solve_retval_base<DiagonalPreconditioner<_MatrixType>, Rhs>
00122 {
00123   typedef DiagonalPreconditioner<_MatrixType> Dec;
00124   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00125 
00126   template<typename Dest> void evalTo(Dest& dst) const
00127   {
00128     dec()._solve(rhs(),dst);
00129   }
00130 };
00131 
00132 }
00133 
00139 class IdentityPreconditioner
00140 {
00141   public:
00142 
00143     IdentityPreconditioner() {}
00144 
00145     template<typename MatrixType>
00146     IdentityPreconditioner(const MatrixType& ) {}
00147     
00148     template<typename MatrixType>
00149     IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
00150     
00151     template<typename MatrixType>
00152     IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
00153 
00154     template<typename MatrixType>
00155     IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
00156     
00157     template<typename Rhs>
00158     inline const Rhs& solve(const Rhs& b) const { return b; }
00159 };
00160 
00161 } // end namespace Eigen
00162 
00163 #endif // EIGEN_BASIC_PRECONDITIONERS_H