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 // 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_BASIC_PRECONDITIONERS_H
00011 #define EIGEN_BASIC_PRECONDITIONERS_H
00012 
00013 namespace Eigen { 
00014 
00032 template <typename _Scalar>
00033 class DiagonalPreconditioner
00034 {
00035     typedef _Scalar Scalar;
00036     typedef Matrix<Scalar,Dynamic,1> Vector;
00037     typedef typename Vector::Index Index;
00038 
00039   public:
00040     // this typedef is only to export the scalar type and compile-time dimensions to solve_retval
00041     typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
00042 
00043     DiagonalPreconditioner() : m_isInitialized(false) {}
00044 
00045     template<typename MatType>
00046     DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols())
00047     {
00048       compute(mat);
00049     }
00050 
00051     Index rows() const { return m_invdiag.size(); }
00052     Index cols() const { return m_invdiag.size(); }
00053     
00054     template<typename MatType>
00055     DiagonalPreconditioner& analyzePattern(const MatType& )
00056     {
00057       return *this;
00058     }
00059     
00060     template<typename MatType>
00061     DiagonalPreconditioner& factorize(const MatType& mat)
00062     {
00063       m_invdiag.resize(mat.cols());
00064       for(int j=0; j<mat.outerSize(); ++j)
00065       {
00066         typename MatType::InnerIterator it(mat,j);
00067         while(it && it.index()!=j) ++it;
00068         if(it && it.index()==j)
00069           m_invdiag(j) = Scalar(1)/it.value();
00070         else
00071           m_invdiag(j) = 0;
00072       }
00073       m_isInitialized = true;
00074       return *this;
00075     }
00076     
00077     template<typename MatType>
00078     DiagonalPreconditioner& compute(const MatType& mat)
00079     {
00080       return factorize(mat);
00081     }
00082 
00083     template<typename Rhs, typename Dest>
00084     void _solve(const Rhs& b, Dest& x) const
00085     {
00086       x = m_invdiag.array() * b.array() ;
00087     }
00088 
00089     template<typename Rhs> inline const internal::solve_retval<DiagonalPreconditioner, Rhs>
00090     solve(const MatrixBase<Rhs>& b) const
00091     {
00092       eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
00093       eigen_assert(m_invdiag.size()==b.rows()
00094                 && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
00095       return internal::solve_retval<DiagonalPreconditioner, Rhs>(*this, b.derived());
00096     }
00097 
00098   protected:
00099     Vector m_invdiag;
00100     bool m_isInitialized;
00101 };
00102 
00103 namespace internal {
00104 
00105 template<typename _MatrixType, typename Rhs>
00106 struct solve_retval<DiagonalPreconditioner<_MatrixType>, Rhs>
00107   : solve_retval_base<DiagonalPreconditioner<_MatrixType>, Rhs>
00108 {
00109   typedef DiagonalPreconditioner<_MatrixType> Dec;
00110   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00111 
00112   template<typename Dest> void evalTo(Dest& dst) const
00113   {
00114     dec()._solve(rhs(),dst);
00115   }
00116 };
00117 
00118 }
00119 
00125 class IdentityPreconditioner
00126 {
00127   public:
00128 
00129     IdentityPreconditioner() {}
00130 
00131     template<typename MatrixType>
00132     IdentityPreconditioner(const MatrixType& ) {}
00133     
00134     template<typename MatrixType>
00135     IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
00136     
00137     template<typename MatrixType>
00138     IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
00139 
00140     template<typename MatrixType>
00141     IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
00142     
00143     template<typename Rhs>
00144     inline const Rhs& solve(const Rhs& b) const { return b; }
00145 };
00146 
00147 } // end namespace Eigen
00148 
00149 #endif // EIGEN_BASIC_PRECONDITIONERS_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Wed Sep 16 2015 07:10:19