Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
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 }
00162
00163 #endif // EIGEN_BASIC_PRECONDITIONERS_H