![]() |
Eigen
3.3.3
|
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