$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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