Eigen  3.3.3
BasicPreconditioners.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_BASIC_PRECONDITIONERS_H
00011 #define EIGEN_BASIC_PRECONDITIONERS_H
00012 
00013 namespace Eigen { 
00014 
00035 template <typename _Scalar>
00036 class DiagonalPreconditioner
00037 {
00038     typedef _Scalar Scalar;
00039     typedef Matrix<Scalar,Dynamic,1> Vector;
00040   public:
00041     typedef typename Vector::StorageIndex StorageIndex;
00042     enum {
00043       ColsAtCompileTime = Dynamic,
00044       MaxColsAtCompileTime = Dynamic
00045     };
00046 
00047     DiagonalPreconditioner() : m_isInitialized(false) {}
00048 
00049     template<typename MatType>
00050     explicit DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols())
00051     {
00052       compute(mat);
00053     }
00054 
00055     Index rows() const { return m_invdiag.size(); }
00056     Index cols() const { return m_invdiag.size(); }
00057     
00058     template<typename MatType>
00059     DiagonalPreconditioner& analyzePattern(const MatType& )
00060     {
00061       return *this;
00062     }
00063     
00064     template<typename MatType>
00065     DiagonalPreconditioner& factorize(const MatType& mat)
00066     {
00067       m_invdiag.resize(mat.cols());
00068       for(int j=0; j<mat.outerSize(); ++j)
00069       {
00070         typename MatType::InnerIterator it(mat,j);
00071         while(it && it.index()!=j) ++it;
00072         if(it && it.index()==j && it.value()!=Scalar(0))
00073           m_invdiag(j) = Scalar(1)/it.value();
00074         else
00075           m_invdiag(j) = Scalar(1);
00076       }
00077       m_isInitialized = true;
00078       return *this;
00079     }
00080     
00081     template<typename MatType>
00082     DiagonalPreconditioner& compute(const MatType& mat)
00083     {
00084       return factorize(mat);
00085     }
00086 
00088     template<typename Rhs, typename Dest>
00089     void _solve_impl(const Rhs& b, Dest& x) const
00090     {
00091       x = m_invdiag.array() * b.array() ;
00092     }
00093 
00094     template<typename Rhs> inline const Solve<DiagonalPreconditioner, Rhs>
00095     solve(const MatrixBase<Rhs>& b) const
00096     {
00097       eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
00098       eigen_assert(m_invdiag.size()==b.rows()
00099                 && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
00100       return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived());
00101     }
00102     
00103     ComputationInfo info() { return Success; }
00104 
00105   protected:
00106     Vector m_invdiag;
00107     bool m_isInitialized;
00108 };
00109 
00127 template <typename _Scalar>
00128 class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar>
00129 {
00130     typedef _Scalar Scalar;
00131     typedef typename NumTraits<Scalar>::Real RealScalar;
00132     typedef DiagonalPreconditioner<_Scalar> Base;
00133     using Base::m_invdiag;
00134   public:
00135 
00136     LeastSquareDiagonalPreconditioner() : Base() {}
00137 
00138     template<typename MatType>
00139     explicit LeastSquareDiagonalPreconditioner(const MatType& mat) : Base()
00140     {
00141       compute(mat);
00142     }
00143 
00144     template<typename MatType>
00145     LeastSquareDiagonalPreconditioner& analyzePattern(const MatType& )
00146     {
00147       return *this;
00148     }
00149     
00150     template<typename MatType>
00151     LeastSquareDiagonalPreconditioner& factorize(const MatType& mat)
00152     {
00153       // Compute the inverse squared-norm of each column of mat
00154       m_invdiag.resize(mat.cols());
00155       for(Index j=0; j<mat.outerSize(); ++j)
00156       {
00157         RealScalar sum = mat.innerVector(j).squaredNorm();
00158         if(sum>0)
00159           m_invdiag(j) = RealScalar(1)/sum;
00160         else
00161           m_invdiag(j) = RealScalar(1);
00162       }
00163       Base::m_isInitialized = true;
00164       return *this;
00165     }
00166     
00167     template<typename MatType>
00168     LeastSquareDiagonalPreconditioner& compute(const MatType& mat)
00169     {
00170       return factorize(mat);
00171     }
00172     
00173     ComputationInfo info() { return Success; }
00174 
00175   protected:
00176 };
00177 
00185 class IdentityPreconditioner
00186 {
00187   public:
00188 
00189     IdentityPreconditioner() {}
00190 
00191     template<typename MatrixType>
00192     explicit IdentityPreconditioner(const MatrixType& ) {}
00193     
00194     template<typename MatrixType>
00195     IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
00196     
00197     template<typename MatrixType>
00198     IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
00199 
00200     template<typename MatrixType>
00201     IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
00202     
00203     template<typename Rhs>
00204     inline const Rhs& solve(const Rhs& b) const { return b; }
00205     
00206     ComputationInfo info() { return Success; }
00207 };
00208 
00209 } // end namespace Eigen
00210 
00211 #endif // EIGEN_BASIC_PRECONDITIONERS_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends