16#ifndef EIGEN_SVDBASE_H
17#define EIGEN_SVDBASE_H
19#include "./InternalHeaderCheck.h"
24template<
typename Derived>
struct traits<SVDBase<Derived> >
27 typedef MatrixXpr XprKind;
28 typedef SolverStorage StorageKind;
29 typedef int StorageIndex;
69 template<
typename Derived_>
70 friend struct internal::solve_assertion;
72 typedef typename internal::traits<Derived>::MatrixType MatrixType;
73 typedef typename MatrixType::Scalar Scalar;
74 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
75 typedef typename Eigen::internal::traits<SVDBase>::StorageIndex StorageIndex;
78 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
79 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
80 DiagSizeAtCompileTime = internal::min_size_prefer_dynamic(RowsAtCompileTime,ColsAtCompileTime),
81 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
82 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
83 MaxDiagSizeAtCompileTime = internal::min_size_prefer_fixed(MaxRowsAtCompileTime,MaxColsAtCompileTime),
84 MatrixOptions = MatrixType::Options
89 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
91 Derived& derived() {
return *
static_cast<Derived*
>(
this); }
92 const Derived& derived()
const {
return *
static_cast<const Derived*
>(
this); }
105 _check_compute_assertions();
106 eigen_assert(
computeU() &&
"This SVD decomposition didn't compute U. Did you ask for it?");
121 _check_compute_assertions();
122 eigen_assert(
computeV() &&
"This SVD decomposition didn't compute V. Did you ask for it?");
133 _check_compute_assertions();
134 return m_singularValues;
140 _check_compute_assertions();
141 return m_nonzeroSingularValues;
153 _check_compute_assertions();
154 if(m_singularValues.size()==0)
return 0;
155 RealScalar premultiplied_threshold = numext::maxi<RealScalar>(m_singularValues.coeff(0) *
threshold(), (std::numeric_limits<RealScalar>::min)());
156 Index i = m_nonzeroSingularValues-1;
157 while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i;
177 m_usePrescribedThreshold =
true;
192 m_usePrescribedThreshold =
false;
202 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
204 Index diagSize = (std::max<Index>)(1,m_diagSize);
205 return m_usePrescribedThreshold ? m_prescribedThreshold
210 inline bool computeU()
const {
return m_computeFullU || m_computeThinU; }
212 inline bool computeV()
const {
return m_computeFullV || m_computeThinV; }
214 inline Index rows()
const {
return m_rows; }
215 inline Index cols()
const {
return m_cols; }
217 #ifdef EIGEN_PARSED_BY_DOXYGEN
227 template<
typename Rhs>
228 inline const Solve<Derived, Rhs>
240 eigen_assert(m_isInitialized &&
"SVD is not initialized.");
244 #ifndef EIGEN_PARSED_BY_DOXYGEN
245 template<
typename RhsType,
typename DstType>
246 void _solve_impl(
const RhsType &rhs, DstType &dst)
const;
248 template<
bool Conjugate,
typename RhsType,
typename DstType>
249 void _solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const;
254 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
256 void _check_compute_assertions()
const {
257 eigen_assert(m_isInitialized &&
"SVD is not initialized.");
260 template<
bool Transpose_,
typename Rhs>
261 void _check_solve_assertion(
const Rhs& b)
const {
262 EIGEN_ONLY_USED_FOR_DEBUG(b);
263 _check_compute_assertions();
264 eigen_assert(
computeU() &&
computeV() &&
"SVDBase::solve(): Both unitaries U and V are required to be computed (thin unitaries suffice).");
265 eigen_assert((Transpose_?cols():rows())==b.rows() &&
"SVDBase::solve(): invalid number of rows of the right hand side matrix b");
269 bool allocate(
Index rows,
Index cols,
unsigned int computationOptions) ;
271 MatrixUType m_matrixU;
272 MatrixVType m_matrixV;
273 SingularValuesType m_singularValues;
275 bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold;
276 bool m_computeFullU, m_computeThinU;
277 bool m_computeFullV, m_computeThinV;
278 unsigned int m_computationOptions;
279 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
280 RealScalar m_prescribedThreshold;
288 m_isInitialized(false),
289 m_isAllocated(false),
290 m_usePrescribedThreshold(false),
291 m_computeFullU(false),
292 m_computeThinU(false),
293 m_computeFullV(false),
294 m_computeThinV(false),
295 m_computationOptions(0),
296 m_rows(-1), m_cols(-1), m_diagSize(0)
302#ifndef EIGEN_PARSED_BY_DOXYGEN
303template<
typename Derived>
304template<
typename RhsType,
typename DstType>
305void SVDBase<Derived>::_solve_impl(
const RhsType &rhs, DstType &dst)
const
310 Matrix<typename RhsType::Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp;
311 Index l_rank = rank();
312 tmp.noalias() = m_matrixU.leftCols(l_rank).adjoint() * rhs;
313 tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp;
314 dst = m_matrixV.leftCols(l_rank) * tmp;
317template<
typename Derived>
318template<
bool Conjugate,
typename RhsType,
typename DstType>
319void SVDBase<Derived>::_solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const
324 Matrix<typename RhsType::Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp;
325 Index l_rank = rank();
327 tmp.noalias() = m_matrixV.leftCols(l_rank).transpose().template conjugateIf<Conjugate>() * rhs;
328 tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp;
329 dst = m_matrixU.template conjugateIf<!Conjugate>().leftCols(l_rank) * tmp;
333template<
typename MatrixType>
334bool SVDBase<MatrixType>::allocate(
Index rows,
Index cols,
unsigned int computationOptions)
336 eigen_assert(rows >= 0 && cols >= 0);
341 computationOptions == m_computationOptions)
349 m_isInitialized =
false;
350 m_isAllocated =
true;
351 m_computationOptions = computationOptions;
352 m_computeFullU = (computationOptions &
ComputeFullU) != 0;
353 m_computeThinU = (computationOptions &
ComputeThinU) != 0;
354 m_computeFullV = (computationOptions &
ComputeFullV) != 0;
355 m_computeThinV = (computationOptions &
ComputeThinV) != 0;
356 eigen_assert(!(m_computeFullU && m_computeThinU) &&
"SVDBase: you can't ask for both full and thin U");
357 eigen_assert(!(m_computeFullV && m_computeThinV) &&
"SVDBase: you can't ask for both full and thin V");
358 eigen_assert(internal::check_implication(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==
Dynamic) &&
359 "SVDBase: thin U and V are only available when your matrix has a dynamic number of columns.");
361 m_diagSize = (std::min)(m_rows, m_cols);
362 m_singularValues.resize(m_diagSize);
364 m_matrixU.resize(m_rows, m_computeFullU ? m_rows : m_computeThinU ? m_diagSize : 0);
366 m_matrixV.resize(m_cols, m_computeFullV ? m_cols : m_computeThinV ? m_diagSize : 0);
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
Base class of SVD algorithms.
Definition: SVDBase.h:66
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SVDBase.h:238
Derived & setThreshold(const RealScalar &threshold)
Definition: SVDBase.h:175
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Index rank() const
Definition: SVDBase.h:150
bool computeV() const
Definition: SVDBase.h:212
Eigen::Index Index
Definition: SVDBase.h:76
bool computeU() const
Definition: SVDBase.h:210
Derived & setThreshold(Default_t)
Definition: SVDBase.h:190
RealScalar threshold() const
Definition: SVDBase.h:200
SVDBase()
Default Constructor.
Definition: SVDBase.h:286
const SingularValuesType & singularValues() const
Definition: SVDBase.h:131
const MatrixUType & matrixU() const
Definition: SVDBase.h:103
const MatrixVType & matrixV() const
Definition: SVDBase.h:119
Index nonzeroSingularValues() const
Definition: SVDBase.h:138
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
ComputationInfo
Definition: Constants.h:442
@ Success
Definition: Constants.h:444
@ ComputeFullV
Definition: Constants.h:399
@ ComputeThinV
Definition: Constants.h:401
@ ComputeFullU
Definition: Constants.h:395
@ ComputeThinU
Definition: Constants.h:397
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:235