Path Tracer
SVDBase.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
5 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com>
8 // Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr>
9 // Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr>
10 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr>
11 //
12 // This Source Code Form is subject to the terms of the Mozilla
13 // Public License v. 2.0. If a copy of the MPL was not distributed
14 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
15 
16 #ifndef EIGEN_SVDBASE_H
17 #define EIGEN_SVDBASE_H
18 
19 namespace Eigen {
20 
21 namespace internal {
22 template<typename Derived> struct traits<SVDBase<Derived> >
23  : traits<Derived>
24 {
25  typedef MatrixXpr XprKind;
26  typedef SolverStorage StorageKind;
27  typedef int StorageIndex;
28  enum { Flags = 0 };
29 };
30 }
31 
59 template<typename Derived> class SVDBase
60  : public SolverBase<SVDBase<Derived> >
61 {
62 public:
63 
64  template<typename Derived_>
65  friend struct internal::solve_assertion;
66 
67  typedef typename internal::traits<Derived>::MatrixType MatrixType;
68  typedef typename MatrixType::Scalar Scalar;
69  typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
70  typedef typename Eigen::internal::traits<SVDBase>::StorageIndex StorageIndex;
71  typedef Eigen::Index Index;
72  enum {
73  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
74  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
75  DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
76  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
77  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
78  MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
79  MatrixOptions = MatrixType::Options
80  };
81 
82  typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixUType;
83  typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> MatrixVType;
84  typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
85 
86  Derived& derived() { return *static_cast<Derived*>(this); }
87  const Derived& derived() const { return *static_cast<const Derived*>(this); }
88 
98  const MatrixUType& matrixU() const
99  {
100  eigen_assert(m_isInitialized && "SVD is not initialized.");
101  eigen_assert(computeU() && "This SVD decomposition didn't compute U. Did you ask for it?");
102  return m_matrixU;
103  }
104 
114  const MatrixVType& matrixV() const
115  {
116  eigen_assert(m_isInitialized && "SVD is not initialized.");
117  eigen_assert(computeV() && "This SVD decomposition didn't compute V. Did you ask for it?");
118  return m_matrixV;
119  }
120 
126  const SingularValuesType& singularValues() const
127  {
128  eigen_assert(m_isInitialized && "SVD is not initialized.");
129  return m_singularValues;
130  }
131 
133  Index nonzeroSingularValues() const
134  {
135  eigen_assert(m_isInitialized && "SVD is not initialized.");
136  return m_nonzeroSingularValues;
137  }
138 
145  inline Index rank() const
146  {
147  using std::abs;
148  eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
149  if(m_singularValues.size()==0) return 0;
150  RealScalar premultiplied_threshold = numext::maxi<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)());
151  Index i = m_nonzeroSingularValues-1;
152  while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i;
153  return i+1;
154  }
155 
170  Derived& setThreshold(const RealScalar& threshold)
171  {
172  m_usePrescribedThreshold = true;
173  m_prescribedThreshold = threshold;
174  return derived();
175  }
176 
185  Derived& setThreshold(Default_t)
186  {
187  m_usePrescribedThreshold = false;
188  return derived();
189  }
190 
195  RealScalar threshold() const
196  {
197  eigen_assert(m_isInitialized || m_usePrescribedThreshold);
198  // this temporary is needed to workaround a MSVC issue
199  Index diagSize = (std::max<Index>)(1,m_diagSize);
200  return m_usePrescribedThreshold ? m_prescribedThreshold
201  : RealScalar(diagSize)*NumTraits<Scalar>::epsilon();
202  }
203 
205  inline bool computeU() const { return m_computeFullU || m_computeThinU; }
207  inline bool computeV() const { return m_computeFullV || m_computeThinV; }
208 
209  inline Index rows() const { return m_rows; }
210  inline Index cols() const { return m_cols; }
211 
212  #ifdef EIGEN_PARSED_BY_DOXYGEN
213 
222  template<typename Rhs>
223  inline const Solve<Derived, Rhs>
224  solve(const MatrixBase<Rhs>& b) const;
225  #endif
226 
227  #ifndef EIGEN_PARSED_BY_DOXYGEN
228  template<typename RhsType, typename DstType>
229  void _solve_impl(const RhsType &rhs, DstType &dst) const;
230 
231  template<bool Conjugate, typename RhsType, typename DstType>
232  void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
233  #endif
234 
235 protected:
236 
237  static void check_template_parameters()
238  {
239  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
240  }
241 
242  template<bool Transpose_, typename Rhs>
243  void _check_solve_assertion(const Rhs& b) const {
244  EIGEN_ONLY_USED_FOR_DEBUG(b);
245  eigen_assert(m_isInitialized && "SVD is not initialized.");
246  eigen_assert(computeU() && computeV() && "SVDBase::solve(): Both unitaries U and V are required to be computed (thin unitaries suffice).");
247  eigen_assert((Transpose_?cols():rows())==b.rows() && "SVDBase::solve(): invalid number of rows of the right hand side matrix b");
248  }
249 
250  // return true if already allocated
251  bool allocate(Index rows, Index cols, unsigned int computationOptions) ;
252 
253  MatrixUType m_matrixU;
254  MatrixVType m_matrixV;
255  SingularValuesType m_singularValues;
256  bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold;
257  bool m_computeFullU, m_computeThinU;
258  bool m_computeFullV, m_computeThinV;
259  unsigned int m_computationOptions;
260  Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
261  RealScalar m_prescribedThreshold;
262 
267  SVDBase()
268  : m_isInitialized(false),
269  m_isAllocated(false),
270  m_usePrescribedThreshold(false),
271  m_computeFullU(false),
272  m_computeThinU(false),
273  m_computeFullV(false),
274  m_computeThinV(false),
275  m_computationOptions(0),
276  m_rows(-1), m_cols(-1), m_diagSize(0)
277  {
278  check_template_parameters();
279  }
280 
281 
282 };
283 
284 #ifndef EIGEN_PARSED_BY_DOXYGEN
285 template<typename Derived>
286 template<typename RhsType, typename DstType>
287 void SVDBase<Derived>::_solve_impl(const RhsType &rhs, DstType &dst) const
288 {
289  // A = U S V^*
290  // So A^{-1} = V S^{-1} U^*
291 
292  Matrix<typename RhsType::Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp;
293  Index l_rank = rank();
294  tmp.noalias() = m_matrixU.leftCols(l_rank).adjoint() * rhs;
295  tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp;
296  dst = m_matrixV.leftCols(l_rank) * tmp;
297 }
298 
299 template<typename Derived>
300 template<bool Conjugate, typename RhsType, typename DstType>
301 void SVDBase<Derived>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
302 {
303  // A = U S V^*
304  // So A^{-*} = U S^{-1} V^*
305  // And A^{-T} = U_conj S^{-1} V^T
306  Matrix<typename RhsType::Scalar, Dynamic, RhsType::ColsAtCompileTime, 0, MatrixType::MaxRowsAtCompileTime, RhsType::MaxColsAtCompileTime> tmp;
307  Index l_rank = rank();
308 
309  tmp.noalias() = m_matrixV.leftCols(l_rank).transpose().template conjugateIf<Conjugate>() * rhs;
310  tmp = m_singularValues.head(l_rank).asDiagonal().inverse() * tmp;
311  dst = m_matrixU.template conjugateIf<!Conjugate>().leftCols(l_rank) * tmp;
312 }
313 #endif
314 
315 template<typename MatrixType>
316 bool SVDBase<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions)
317 {
318  eigen_assert(rows >= 0 && cols >= 0);
319 
320  if (m_isAllocated &&
321  rows == m_rows &&
322  cols == m_cols &&
323  computationOptions == m_computationOptions)
324  {
325  return true;
326  }
327 
328  m_rows = rows;
329  m_cols = cols;
330  m_isInitialized = false;
331  m_isAllocated = true;
332  m_computationOptions = computationOptions;
333  m_computeFullU = (computationOptions & ComputeFullU) != 0;
334  m_computeThinU = (computationOptions & ComputeThinU) != 0;
335  m_computeFullV = (computationOptions & ComputeFullV) != 0;
336  m_computeThinV = (computationOptions & ComputeThinV) != 0;
337  eigen_assert(!(m_computeFullU && m_computeThinU) && "SVDBase: you can't ask for both full and thin U");
338  eigen_assert(!(m_computeFullV && m_computeThinV) && "SVDBase: you can't ask for both full and thin V");
339  eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
340  "SVDBase: thin U and V are only available when your matrix has a dynamic number of columns.");
341 
342  m_diagSize = (std::min)(m_rows, m_cols);
343  m_singularValues.resize(m_diagSize);
344  if(RowsAtCompileTime==Dynamic)
345  m_matrixU.resize(m_rows, m_computeFullU ? m_rows : m_computeThinU ? m_diagSize : 0);
346  if(ColsAtCompileTime==Dynamic)
347  m_matrixV.resize(m_cols, m_computeFullV ? m_cols : m_computeThinV ? m_diagSize : 0);
348 
349  return false;
350 }
351 
352 }// end namespace
353 
354 #endif // EIGEN_SVDBASE_H
Eigen::MatrixXpr
Definition: Constants.h:518
Eigen
Namespace containing all symbols from the Eigen library.
Definition: LDLT.h:16
Eigen::ComputeFullV
@ ComputeFullV
Definition: Constants.h:396
Eigen::SVDBase
Base class of SVD algorithms.
Definition: SVDBase.h:61
Eigen::SVDBase::singularValues
const SingularValuesType & singularValues() const
Definition: SVDBase.h:126
Eigen::ComputeFullU
@ ComputeFullU
Definition: Constants.h:392
Eigen::internal::solve_assertion
Definition: SolverBase.h:18
Eigen::SVDBase::threshold
RealScalar threshold() const
Definition: SVDBase.h:195
Eigen::SVDBase::rank
Index rank() const
Definition: SVDBase.h:145
Eigen::SVDBase::setThreshold
Derived & setThreshold(const RealScalar &threshold)
Definition: SVDBase.h:170
Eigen::SVDBase::setThreshold
Derived & setThreshold(Default_t)
Definition: SVDBase.h:185
Eigen::SVDBase::nonzeroSingularValues
Index nonzeroSingularValues() const
Definition: SVDBase.h:133
Eigen::internal::true_type
Definition: Meta.h:63
Eigen::SVDBase::matrixU
const MatrixUType & matrixU() const
Definition: SVDBase.h:98
Eigen::SolverStorage
Definition: Constants.h:509
Eigen::ComputeThinU
@ ComputeThinU
Definition: Constants.h:394
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:21
Eigen::SVDBase::computeU
bool computeU() const
Definition: SVDBase.h:205
Eigen::SVDBase::matrixV
const MatrixVType & matrixV() const
Definition: SVDBase.h:114
Eigen::SVDBase::computeV
bool computeV() const
Definition: SVDBase.h:207
Eigen::ComputeThinV
@ ComputeThinV
Definition: Constants.h:398
Eigen::internal::traits
Definition: ForwardDeclarations.h:17
Eigen::SVDBase::SVDBase
SVDBase()
Default Constructor.
Definition: SVDBase.h:267
Eigen::Matrix< Scalar, RowsAtCompileTime, RowsAtCompileTime, MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime >
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:213
Eigen::SVDBase::Index
Eigen::Index Index
Definition: SVDBase.h:71
Eigen::SolverBase
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:69
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:42
Generated by   doxygen 1.8.19