![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> 00005 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // 00007 // Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com> 00008 // Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr> 00009 // Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr> 00010 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr> 00011 // 00012 // This Source Code Form is subject to the terms of the Mozilla 00013 // Public License v. 2.0. If a copy of the MPL was not distributed 00014 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00015 00016 #ifndef EIGEN_SVDBASE_H 00017 #define EIGEN_SVDBASE_H 00018 00019 namespace Eigen { 00047 template<typename Derived> 00048 class SVDBase 00049 { 00050 00051 public: 00052 typedef typename internal::traits<Derived>::MatrixType MatrixType; 00053 typedef typename MatrixType::Scalar Scalar; 00054 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 00055 typedef typename MatrixType::StorageIndex StorageIndex; 00056 typedef Eigen::Index Index; 00057 enum { 00058 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00059 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00060 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime), 00061 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00062 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00063 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime), 00064 MatrixOptions = MatrixType::Options 00065 }; 00066 00067 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixUType; 00068 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> MatrixVType; 00069 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType; 00070 00071 Derived& derived() { return *static_cast<Derived*>(this); } 00072 const Derived& derived() const { return *static_cast<const Derived*>(this); } 00073 00083 const MatrixUType& matrixU() const 00084 { 00085 eigen_assert(m_isInitialized && "SVD is not initialized."); 00086 eigen_assert(computeU() && "This SVD decomposition didn't compute U. Did you ask for it?"); 00087 return m_matrixU; 00088 } 00089 00099 const MatrixVType& matrixV() const 00100 { 00101 eigen_assert(m_isInitialized && "SVD is not initialized."); 00102 eigen_assert(computeV() && "This SVD decomposition didn't compute V. Did you ask for it?"); 00103 return m_matrixV; 00104 } 00105 00111 const SingularValuesType& singularValues() const 00112 { 00113 eigen_assert(m_isInitialized && "SVD is not initialized."); 00114 return m_singularValues; 00115 } 00116 00118 Index nonzeroSingularValues() const 00119 { 00120 eigen_assert(m_isInitialized && "SVD is not initialized."); 00121 return m_nonzeroSingularValues; 00122 } 00123 00130 inline Index rank() const 00131 { 00132 using std::abs; 00133 eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 00134 if(m_singularValues.size()==0) return 0; 00135 RealScalar premultiplied_threshold = numext::maxi<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)()); 00136 Index i = m_nonzeroSingularValues-1; 00137 while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i; 00138 return i+1; 00139 } 00140 00155 Derived& setThreshold(const RealScalar& threshold) 00156 { 00157 m_usePrescribedThreshold = true; 00158 m_prescribedThreshold = threshold; 00159 return derived(); 00160 } 00161 00170 Derived& setThreshold(Default_t) 00171 { 00172 m_usePrescribedThreshold = false; 00173 return derived(); 00174 } 00175 00180 RealScalar threshold() const 00181 { 00182 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 00183 return m_usePrescribedThreshold ? m_prescribedThreshold 00184 : (std::max<Index>)(1,m_diagSize)*NumTraits<Scalar>::epsilon(); 00185 } 00186 00188 inline bool computeU() const { return m_computeFullU || m_computeThinU; } 00190 inline bool computeV() const { return m_computeFullV || m_computeThinV; } 00191 00192 inline Index rows() const { return m_rows; } 00193 inline Index cols() const { return m_cols; } 00194 00204 template<typename Rhs> 00205 inline const Solve<Derived, Rhs> 00206 solve(const MatrixBase<Rhs>& b) const 00207 { 00208 eigen_assert(m_isInitialized && "SVD is not initialized."); 00209 eigen_assert(computeU() && computeV() && "SVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); 00210 return Solve<Derived, Rhs>(derived(), b.derived()); 00211 } 00212 00213 #ifndef EIGEN_PARSED_BY_DOXYGEN 00214 template<typename RhsType, typename DstType> 00215 EIGEN_DEVICE_FUNC 00216 void _solve_impl(const RhsType &rhs, DstType &dst) const; 00217 #endif 00218 00219 protected: 00220 00221 static void check_template_parameters() 00222 { 00223 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00224 } 00225 00226 // return true if already allocated 00227 bool allocate(Index rows, Index cols, unsigned int computationOptions) ; 00228 00229 MatrixUType m_matrixU; 00230 MatrixVType m_matrixV; 00231 SingularValuesType m_singularValues; 00232 bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold; 00233 bool m_computeFullU, m_computeThinU; 00234 bool m_computeFullV, m_computeThinV; 00235 unsigned int m_computationOptions; 00236 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize; 00237 RealScalar m_prescribedThreshold; 00238 00243 SVDBase() 00244 : m_isInitialized(false), 00245 m_isAllocated(false), 00246 m_usePrescribedThreshold(false), 00247 m_computationOptions(0), 00248 m_rows(-1), m_cols(-1), m_diagSize(0) 00249 { 00250 check_template_parameters(); 00251 } 00252 00253 00254 }; 00255 00256 #ifndef EIGEN_PARSED_BY_DOXYGEN 00257 template<typename Derived> 00258 template<typename RhsType, typename DstType> 00259 void SVDBase<Derived>::_solve_impl(const RhsType &rhs, DstType &dst) const 00260 { 00261 eigen_assert(rhs.rows() == rows()); 00262 00263 // A = U S V^* 00264 // So A^{-1} = V S^{-1} U^* 00265 00266 Matrix<Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp; 00267 Index l_rank = rank(); 00268 tmp.noalias() = m_matrixU.leftCols(l_rank).adjoint() * rhs; 00269 tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp; 00270 dst = m_matrixV.leftCols(l_rank) * tmp; 00271 } 00272 #endif 00273 00274 template<typename MatrixType> 00275 bool SVDBase<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions) 00276 { 00277 eigen_assert(rows >= 0 && cols >= 0); 00278 00279 if (m_isAllocated && 00280 rows == m_rows && 00281 cols == m_cols && 00282 computationOptions == m_computationOptions) 00283 { 00284 return true; 00285 } 00286 00287 m_rows = rows; 00288 m_cols = cols; 00289 m_isInitialized = false; 00290 m_isAllocated = true; 00291 m_computationOptions = computationOptions; 00292 m_computeFullU = (computationOptions & ComputeFullU) != 0; 00293 m_computeThinU = (computationOptions & ComputeThinU) != 0; 00294 m_computeFullV = (computationOptions & ComputeFullV) != 0; 00295 m_computeThinV = (computationOptions & ComputeThinV) != 0; 00296 eigen_assert(!(m_computeFullU && m_computeThinU) && "SVDBase: you can't ask for both full and thin U"); 00297 eigen_assert(!(m_computeFullV && m_computeThinV) && "SVDBase: you can't ask for both full and thin V"); 00298 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && 00299 "SVDBase: thin U and V are only available when your matrix has a dynamic number of columns."); 00300 00301 m_diagSize = (std::min)(m_rows, m_cols); 00302 m_singularValues.resize(m_diagSize); 00303 if(RowsAtCompileTime==Dynamic) 00304 m_matrixU.resize(m_rows, m_computeFullU ? m_rows : m_computeThinU ? m_diagSize : 0); 00305 if(ColsAtCompileTime==Dynamic) 00306 m_matrixV.resize(m_cols, m_computeFullV ? m_cols : m_computeThinV ? m_diagSize : 0); 00307 00308 return false; 00309 } 00310 00311 }// end namespace 00312 00313 #endif // EIGEN_SVDBASE_H