JacobiSVD.h
Go to the documentation of this file.
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) 2013-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_JACOBISVD_H
12 #define EIGEN_JACOBISVD_H
13 
14 namespace Eigen {
15 
16 namespace internal {
17 // forward declaration (needed by ICC)
18 // the empty body is required by MSVC
19 template<typename MatrixType, int QRPreconditioner,
22 
23 /*** QR preconditioners (R-SVD)
24  ***
25  *** Their role is to reduce the problem of computing the SVD to the case of a square matrix.
26  *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for
27  *** JacobiSVD which by itself is only able to work on square matrices.
28  ***/
29 
31 
32 template<typename MatrixType, int QRPreconditioner, int Case>
34 {
35  enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
36  MatrixType::ColsAtCompileTime != Dynamic &&
37  MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
38  b = MatrixType::RowsAtCompileTime != Dynamic &&
39  MatrixType::ColsAtCompileTime != Dynamic &&
40  MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
41  ret = !( (QRPreconditioner == NoQRPreconditioner) ||
42  (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
43  (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
44  };
45 };
46 
47 template<typename MatrixType, int QRPreconditioner, int Case,
50 
51 template<typename MatrixType, int QRPreconditioner, int Case>
52 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
53 {
54 public:
57  {
58  return false;
59  }
60 };
61 
62 /*** preconditioner using FullPivHouseholderQR ***/
63 
64 template<typename MatrixType>
65 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
66 {
67 public:
68  typedef typename MatrixType::Scalar Scalar;
69  enum
70  {
71  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
72  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
73  };
75 
77  {
78  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
79  {
80  m_qr.~QRType();
81  ::new (&m_qr) QRType(svd.rows(), svd.cols());
82  }
83  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
84  }
85 
87  {
88  if(matrix.rows() > matrix.cols())
89  {
90  m_qr.compute(matrix);
91  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
92  if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
93  if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
94  return true;
95  }
96  return false;
97  }
98 private:
100  QRType m_qr;
101  WorkspaceType m_workspace;
102 };
103 
104 template<typename MatrixType>
106 {
107 public:
108  typedef typename MatrixType::Scalar Scalar;
109  enum
110  {
111  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
112  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
113  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
114  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
115  TrOptions = RowsAtCompileTime==1 ? (MatrixType::Options & ~(RowMajor))
116  : ColsAtCompileTime==1 ? (MatrixType::Options | RowMajor)
117  : MatrixType::Options
118  };
121 
123  {
124  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
125  {
126  m_qr.~QRType();
127  ::new (&m_qr) QRType(svd.cols(), svd.rows());
128  }
129  m_adjoint.resize(svd.cols(), svd.rows());
130  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
131  }
132 
134  {
135  if(matrix.cols() > matrix.rows())
136  {
137  m_adjoint = matrix.adjoint();
138  m_qr.compute(m_adjoint);
139  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
140  if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
141  if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
142  return true;
143  }
144  else return false;
145  }
146 private:
148  QRType m_qr;
151 };
152 
153 /*** preconditioner using ColPivHouseholderQR ***/
154 
155 template<typename MatrixType>
156 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
157 {
158 public:
160  {
161  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
162  {
163  m_qr.~QRType();
164  ::new (&m_qr) QRType(svd.rows(), svd.cols());
165  }
166  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
167  else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
168  }
169 
171  {
172  if(matrix.rows() > matrix.cols())
173  {
174  m_qr.compute(matrix);
175  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
176  if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
177  else if(svd.m_computeThinU)
178  {
179  svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
180  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
181  }
182  if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
183  return true;
184  }
185  return false;
186  }
187 
188 private:
190  QRType m_qr;
192 };
193 
194 template<typename MatrixType>
196 {
197 public:
198  typedef typename MatrixType::Scalar Scalar;
199  enum
200  {
201  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
202  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
203  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
204  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
205  TrOptions = RowsAtCompileTime==1 ? (MatrixType::Options & ~(RowMajor))
206  : ColsAtCompileTime==1 ? (MatrixType::Options | RowMajor)
207  : MatrixType::Options
208  };
209 
212 
214  {
215  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
216  {
217  m_qr.~QRType();
218  ::new (&m_qr) QRType(svd.cols(), svd.rows());
219  }
220  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
221  else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
222  m_adjoint.resize(svd.cols(), svd.rows());
223  }
224 
226  {
227  if(matrix.cols() > matrix.rows())
228  {
229  m_adjoint = matrix.adjoint();
230  m_qr.compute(m_adjoint);
231 
232  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
233  if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
234  else if(svd.m_computeThinV)
235  {
236  svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
237  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
238  }
239  if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
240  return true;
241  }
242  else return false;
243  }
244 
245 private:
247  QRType m_qr;
250 };
251 
252 /*** preconditioner using HouseholderQR ***/
253 
254 template<typename MatrixType>
255 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
256 {
257 public:
259  {
260  if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
261  {
262  m_qr.~QRType();
263  ::new (&m_qr) QRType(svd.rows(), svd.cols());
264  }
265  if (svd.m_computeFullU) m_workspace.resize(svd.rows());
266  else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
267  }
268 
270  {
271  if(matrix.rows() > matrix.cols())
272  {
273  m_qr.compute(matrix);
274  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
275  if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
276  else if(svd.m_computeThinU)
277  {
278  svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
279  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
280  }
281  if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
282  return true;
283  }
284  return false;
285  }
286 private:
288  QRType m_qr;
290 };
291 
292 template<typename MatrixType>
294 {
295 public:
296  typedef typename MatrixType::Scalar Scalar;
297  enum
298  {
299  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
300  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
301  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
302  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
303  Options = MatrixType::Options
304  };
305 
308 
310  {
311  if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
312  {
313  m_qr.~QRType();
314  ::new (&m_qr) QRType(svd.cols(), svd.rows());
315  }
316  if (svd.m_computeFullV) m_workspace.resize(svd.cols());
317  else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
318  m_adjoint.resize(svd.cols(), svd.rows());
319  }
320 
322  {
323  if(matrix.cols() > matrix.rows())
324  {
325  m_adjoint = matrix.adjoint();
326  m_qr.compute(m_adjoint);
327 
328  svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
329  if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
330  else if(svd.m_computeThinV)
331  {
332  svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
333  m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
334  }
335  if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
336  return true;
337  }
338  else return false;
339  }
340 
341 private:
343  QRType m_qr;
346 };
347 
348 /*** 2x2 SVD implementation
349  ***
350  *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
351  ***/
352 
353 template<typename MatrixType, int QRPreconditioner>
354 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
355 {
358  static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; }
359 };
360 
361 template<typename MatrixType, int QRPreconditioner>
362 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
363 {
365  typedef typename MatrixType::Scalar Scalar;
367  static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry)
368  {
369  using std::sqrt;
370  using std::abs;
371  Scalar z;
373  RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
374 
375  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
376  const RealScalar precision = NumTraits<Scalar>::epsilon();
377 
378  if(n==0)
379  {
380  // make sure first column is zero
381  work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0);
382 
383  if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
384  {
385  // work_matrix.coeff(p,q) can be zero if work_matrix.coeff(q,p) is not zero but small enough to underflow when computing n
386  z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
387  work_matrix.row(p) *= z;
388  if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
389  }
390  if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
391  {
392  z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
393  work_matrix.row(q) *= z;
394  if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
395  }
396  // otherwise the second row is already zero, so we have nothing to do.
397  }
398  else
399  {
400  rot.c() = conj(work_matrix.coeff(p,p)) / n;
401  rot.s() = work_matrix.coeff(q,p) / n;
402  work_matrix.applyOnTheLeft(p,q,rot);
403  if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
404  if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
405  {
406  z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
407  work_matrix.col(q) *= z;
408  if(svd.computeV()) svd.m_matrixV.col(q) *= z;
409  }
410  if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
411  {
412  z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
413  work_matrix.row(q) *= z;
414  if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
415  }
416  }
417 
418  // update largest diagonal entry
419  maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q))));
420  // and check whether the 2x2 block is already diagonal
421  RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
422  return abs(work_matrix.coeff(p,q))>threshold || abs(work_matrix.coeff(q,p)) > threshold;
423  }
424 };
425 
426 template<typename _MatrixType, int QRPreconditioner>
427 struct traits<JacobiSVD<_MatrixType,QRPreconditioner> >
428 {
429  typedef _MatrixType MatrixType;
430 };
431 
432 } // end namespace internal
433 
487 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
488  : public SVDBase<JacobiSVD<_MatrixType,QRPreconditioner> >
489 {
491  public:
492 
493  typedef _MatrixType MatrixType;
494  typedef typename MatrixType::Scalar Scalar;
496  enum {
497  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
498  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
499  DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
500  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
501  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
502  MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
503  MatrixOptions = MatrixType::Options
504  };
505 
506  typedef typename Base::MatrixUType MatrixUType;
507  typedef typename Base::MatrixVType MatrixVType;
509 
512  typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
513  MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
515 
522  {}
523 
524 
531  JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
532  {
533  allocate(rows, cols, computationOptions);
534  }
535 
546  explicit JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
547  {
548  compute(matrix, computationOptions);
549  }
550 
561  JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
562 
569  JacobiSVD& compute(const MatrixType& matrix)
570  {
571  return compute(matrix, m_computationOptions);
572  }
573 
574  using Base::computeU;
575  using Base::computeV;
576  using Base::rows;
577  using Base::cols;
578  using Base::rank;
579 
580  private:
581  void allocate(Index rows, Index cols, unsigned int computationOptions);
582 
583  protected:
584  using Base::m_matrixU;
585  using Base::m_matrixV;
586  using Base::m_singularValues;
587  using Base::m_isInitialized;
588  using Base::m_isAllocated;
589  using Base::m_usePrescribedThreshold;
590  using Base::m_computeFullU;
591  using Base::m_computeThinU;
592  using Base::m_computeFullV;
593  using Base::m_computeThinV;
594  using Base::m_computationOptions;
595  using Base::m_nonzeroSingularValues;
596  using Base::m_rows;
597  using Base::m_cols;
598  using Base::m_diagSize;
599  using Base::m_prescribedThreshold;
601 
602  template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
604  template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
606 
609  MatrixType m_scaledMatrix;
610 };
611 
612 template<typename MatrixType, int QRPreconditioner>
614 {
615  eigen_assert(rows >= 0 && cols >= 0);
616 
617  if (m_isAllocated &&
618  rows == m_rows &&
619  cols == m_cols &&
620  computationOptions == m_computationOptions)
621  {
622  return;
623  }
624 
625  m_rows = rows;
626  m_cols = cols;
627  m_isInitialized = false;
628  m_isAllocated = true;
629  m_computationOptions = computationOptions;
630  m_computeFullU = (computationOptions & ComputeFullU) != 0;
631  m_computeThinU = (computationOptions & ComputeThinU) != 0;
632  m_computeFullV = (computationOptions & ComputeFullV) != 0;
633  m_computeThinV = (computationOptions & ComputeThinV) != 0;
634  eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
635  eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
636  eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
637  "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
638  if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
639  {
640  eigen_assert(!(m_computeThinU || m_computeThinV) &&
641  "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
642  "Use the ColPivHouseholderQR preconditioner instead.");
643  }
644  m_diagSize = (std::min)(m_rows, m_cols);
645  m_singularValues.resize(m_diagSize);
646  if(RowsAtCompileTime==Dynamic)
647  m_matrixU.resize(m_rows, m_computeFullU ? m_rows
648  : m_computeThinU ? m_diagSize
649  : 0);
650  if(ColsAtCompileTime==Dynamic)
651  m_matrixV.resize(m_cols, m_computeFullV ? m_cols
652  : m_computeThinV ? m_diagSize
653  : 0);
654  m_workMatrix.resize(m_diagSize, m_diagSize);
655 
656  if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
657  if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
658  if(m_rows!=m_cols) m_scaledMatrix.resize(rows,cols);
659 }
660 
661 template<typename MatrixType, int QRPreconditioner>
664 {
665  using std::abs;
666  allocate(matrix.rows(), matrix.cols(), computationOptions);
667 
668  // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
669  // only worsening the precision of U and V as we accumulate more rotations
671 
672  // limit for denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
673  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
674 
675  // Scaling factor to reduce over/under-flows
676  RealScalar scale = matrix.cwiseAbs().maxCoeff();
677  if(scale==RealScalar(0)) scale = RealScalar(1);
678 
679  /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
680 
681  if(m_rows!=m_cols)
682  {
683  m_scaledMatrix = matrix / scale;
684  m_qr_precond_morecols.run(*this, m_scaledMatrix);
685  m_qr_precond_morerows.run(*this, m_scaledMatrix);
686  }
687  else
688  {
689  m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
690  if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
691  if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
692  if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
693  if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
694  }
695 
696  /*** step 2. The main Jacobi SVD iteration. ***/
697  RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff();
698 
699  bool finished = false;
700  while(!finished)
701  {
702  finished = true;
703 
704  // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
705 
706  for(Index p = 1; p < m_diagSize; ++p)
707  {
708  for(Index q = 0; q < p; ++q)
709  {
710  // if this 2x2 sub-matrix is not diagonal already...
711  // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
712  // keep us iterating forever. Similarly, small denormal numbers are considered zero.
713  RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
714  if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
715  {
716  finished = false;
717  // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
718  // the complex to real operation returns true if the updated 2x2 block is not already diagonal
720  {
721  JacobiRotation<RealScalar> j_left, j_right;
722  internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
723 
724  // accumulate resulting Jacobi rotations
725  m_workMatrix.applyOnTheLeft(p,q,j_left);
726  if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
727 
728  m_workMatrix.applyOnTheRight(p,q,j_right);
729  if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
730 
731  // keep track of the largest diagonal coefficient
732  maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q))));
733  }
734  }
735  }
736  }
737  }
738 
739  /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
740 
741  for(Index i = 0; i < m_diagSize; ++i)
742  {
743  // For a complex matrix, some diagonal coefficients might note have been
744  // treated by svd_precondition_2x2_block_to_be_real, and the imaginary part
745  // of some diagonal entry might not be null.
746  if(NumTraits<Scalar>::IsComplex && abs(numext::imag(m_workMatrix.coeff(i,i)))>considerAsZero)
747  {
748  RealScalar a = abs(m_workMatrix.coeff(i,i));
749  m_singularValues.coeffRef(i) = abs(a);
750  if(computeU()) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
751  }
752  else
753  {
754  // m_workMatrix.coeff(i,i) is already real, no difficulty:
755  RealScalar a = numext::real(m_workMatrix.coeff(i,i));
756  m_singularValues.coeffRef(i) = abs(a);
757  if(computeU() && (a<RealScalar(0))) m_matrixU.col(i) = -m_matrixU.col(i);
758  }
759  }
760 
761  m_singularValues *= scale;
762 
763  /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
764 
765  m_nonzeroSingularValues = m_diagSize;
766  for(Index i = 0; i < m_diagSize; i++)
767  {
768  Index pos;
769  RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
770  if(maxRemainingSingularValue == RealScalar(0))
771  {
772  m_nonzeroSingularValues = i;
773  break;
774  }
775  if(pos)
776  {
777  pos += i;
778  std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
779  if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
780  if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
781  }
782  }
783 
784  m_isInitialized = true;
785  return *this;
786 }
787 
795 template<typename Derived>
797 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
798 {
799  return JacobiSVD<PlainObject>(*this, computationOptions);
800 }
801 
802 } // end namespace Eigen
803 
804 #endif // EIGEN_JACOBISVD_H
Matrix< Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime > WorkMatrixType
Definition: JacobiSVD.h:514
bool run(JacobiSVD< MatrixType, QRPreconditioner > &, const MatrixType &)
Definition: JacobiSVD.h:56
Scalar & c()
Definition: Jacobi.h:45
int EIGEN_BLAS_FUNC() rot(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps)
SCALAR Scalar
Definition: bench_gemm.cpp:33
cout<< "Here is the matrix m:"<< endl<< m<< endl;JacobiSVD< MatrixXf > svd(m, ComputeThinU|ComputeThinV)
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
internal::plain_col_type< MatrixType >::type ColType
Definition: JacobiSVD.h:511
static bool run(typename SVD::WorkMatrixType &, SVD &, Index, Index, RealScalar &)
Definition: JacobiSVD.h:358
float real
Definition: datatypes.h:10
Matrix< Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime > TransposeTypeWithSameStorageOrder
Definition: JacobiSVD.h:307
Scalar * b
Definition: benchVecAdd.cpp:17
void allocate(const JacobiSVD< MatrixType, QRPreconditioner > &)
Definition: JacobiSVD.h:55
#define min(a, b)
Definition: datatypes.h:19
int n
SVDBase< JacobiSVD > Base
Definition: JacobiSVD.h:490
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Rotation given by a cosine-sine pair.
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
static bool run(typename SVD::WorkMatrixType &work_matrix, SVD &svd, Index p, Index q, RealScalar &maxDiagEntry)
Definition: JacobiSVD.h:367
JacobiRotation transpose() const
Definition: Jacobi.h:59
#define EIGEN_IMPLIES(a, b)
Definition: Macros.h:902
#define EIGEN_SIZE_MIN_PREFER_FIXED(a, b)
Definition: Macros.h:889
static double epsilon
Definition: testRot3.cpp:39
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
Array33i a
MatrixType::Scalar Scalar
Definition: JacobiSVD.h:494
Matrix< Scalar, ColsAtCompileTime, RowsAtCompileTime, TrOptions, MaxColsAtCompileTime, MaxRowsAtCompileTime > TransposeTypeWithSameStorageOrder
Definition: JacobiSVD.h:211
bool run(JacobiSVD< MatrixType, ColPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:170
void real_2x2_jacobi_svd(const MatrixType &matrix, Index p, Index q, JacobiRotation< RealScalar > *j_left, JacobiRotation< RealScalar > *j_right)
Definition: RealSvd2x2.h:19
Base class of SVD algorithms.
Definition: SVDBase.h:48
Matrix< Scalar, ColsAtCompileTime, RowsAtCompileTime, TrOptions, MaxColsAtCompileTime, MaxRowsAtCompileTime > TransposeTypeWithSameStorageOrder
Definition: JacobiSVD.h:120
JacobiSVD()
Default Constructor.
Definition: JacobiSVD.h:521
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
#define eigen_assert(x)
Definition: Macros.h:579
internal::qr_preconditioner_impl< MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows > m_qr_precond_morecols
Definition: JacobiSVD.h:607
_MatrixType MatrixType
Definition: JacobiSVD.h:493
JacobiRotation adjoint() const
Definition: Jacobi.h:62
EIGEN_DEVICE_FUNC const Scalar & q
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Abs2ReturnType abs2() const
Base::MatrixVType MatrixVType
Definition: JacobiSVD.h:507
bool run(JacobiSVD< MatrixType, HouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:269
bool run(JacobiSVD< MatrixType, FullPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:86
bool run(JacobiSVD< MatrixType, HouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:321
MatrixType m_scaledMatrix
Definition: JacobiSVD.h:609
DenseIndex ret
Definition: level1_impl.h:59
JacobiSVD & compute(const MatrixType &matrix)
Method performing the decomposition of given matrix using current options.
Definition: JacobiSVD.h:569
cout precision(2)
void allocate(Index rows, Index cols, unsigned int computationOptions)
Definition: JacobiSVD.h:613
Two-sided Jacobi SVD decomposition of a rectangular matrix.
float * p
Scalar & s()
Definition: Jacobi.h:47
internal::plain_row_type< MatrixType >::type RowType
Definition: JacobiSVD.h:510
NumTraits< typename MatrixType::Scalar >::Real RealScalar
Definition: JacobiSVD.h:495
JacobiSVD & compute(const MatrixType &matrix, unsigned int computationOptions)
Method performing the decomposition of given matrix using custom options.
Definition: JacobiSVD.h:663
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics scale
bool run(JacobiSVD< MatrixType, FullPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:133
EIGEN_DEVICE_FUNC const ImagReturnType imag() const
Base::SingularValuesType SingularValuesType
Definition: JacobiSVD.h:508
internal::qr_preconditioner_impl< MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols > m_qr_precond_morerows
Definition: JacobiSVD.h:608
const int Dynamic
Definition: Constants.h:21
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
#define EIGEN_SIZE_MIN_PREFER_DYNAMIC(a, b)
Definition: Macros.h:881
#define abs(x)
Definition: datatypes.h:17
JacobiSVD(Index rows, Index cols, unsigned int computationOptions=0)
Default Constructor with memory preallocation.
Definition: JacobiSVD.h:531
bool run(JacobiSVD< MatrixType, ColPivHouseholderQRPreconditioner > &svd, const MatrixType &matrix)
Definition: JacobiSVD.h:225
WorkMatrixType m_workMatrix
Definition: JacobiSVD.h:600
void swap(mpfr::mpreal &x, mpfr::mpreal &y)
Definition: mpreal.h:2986
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
JacobiSVD(const MatrixType &matrix, unsigned int computationOptions=0)
Constructor performing the decomposition of given matrix.
Definition: JacobiSVD.h:546
JacobiSVD< PlainObject > jacobiSvd(unsigned int computationOptions=0) const
Definition: JacobiSVD.h:797
Base::MatrixUType MatrixUType
Definition: JacobiSVD.h:506


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:26