10 #ifndef EIGEN_JACOBISVD_H 11 #define EIGEN_JACOBISVD_H 18 template<
typename MatrixType,
int QRPreconditioner,
19 bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
20 struct svd_precondition_2x2_block_to_be_real {};
29 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
31 template<
typename MatrixType,
int QRPreconditioner,
int Case>
32 struct qr_preconditioner_should_do_anything
34 enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
35 MatrixType::ColsAtCompileTime != Dynamic &&
36 MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
37 b = MatrixType::RowsAtCompileTime != Dynamic &&
38 MatrixType::ColsAtCompileTime != Dynamic &&
39 MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
41 (Case == PreconditionIfMoreColsThanRows &&
bool(a)) ||
42 (Case == PreconditionIfMoreRowsThanCols &&
bool(b)) )
46 template<
typename MatrixType,
int QRPreconditioner,
int Case,
47 bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
48 >
struct qr_preconditioner_impl {};
50 template<
typename MatrixType,
int QRPreconditioner,
int Case>
51 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
54 typedef typename MatrixType::Index Index;
55 void allocate(
const JacobiSVD<MatrixType, QRPreconditioner>&) {}
56 bool run(JacobiSVD<MatrixType, QRPreconditioner>&,
const MatrixType&)
64 template<
typename MatrixType>
68 typedef typename MatrixType::Index Index;
69 typedef typename MatrixType::Scalar Scalar;
72 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
73 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
75 typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
77 void allocate(
const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
79 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
82 ::new (&m_qr) QRType(svd.rows(), svd.cols());
84 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
89 if(matrix.rows() > matrix.cols())
92 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).
template triangularView<Upper>();
93 if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
94 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
100 typedef FullPivHouseholderQR<MatrixType> QRType;
102 WorkspaceType m_workspace;
105 template<
typename MatrixType>
109 typedef typename MatrixType::Index Index;
110 typedef typename MatrixType::Scalar Scalar;
113 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
114 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
115 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
116 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
117 Options = MatrixType::Options
119 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
120 TransposeTypeWithSameStorageOrder;
122 void allocate(
const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
124 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
127 ::new (&m_qr) QRType(svd.cols(), svd.rows());
129 m_adjoint.resize(svd.cols(), svd.rows());
130 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
135 if(matrix.cols() > matrix.rows())
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();
147 typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
149 TransposeTypeWithSameStorageOrder m_adjoint;
150 typename internal::plain_row_type<MatrixType>::type m_workspace;
155 template<
typename MatrixType>
159 typedef typename MatrixType::Index Index;
161 void allocate(
const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
163 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
166 ::new (&m_qr) QRType(svd.rows(), svd.cols());
168 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
169 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
174 if(matrix.rows() > matrix.cols())
176 m_qr.compute(matrix);
177 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).
template triangularView<Upper>();
178 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
179 else if(svd.m_computeThinU)
181 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
182 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
184 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
191 typedef ColPivHouseholderQR<MatrixType> QRType;
193 typename internal::plain_col_type<MatrixType>::type m_workspace;
196 template<
typename MatrixType>
200 typedef typename MatrixType::Index Index;
201 typedef typename MatrixType::Scalar Scalar;
204 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
205 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
206 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
207 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
208 Options = MatrixType::Options
211 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
212 TransposeTypeWithSameStorageOrder;
214 void allocate(
const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
216 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
219 ::new (&m_qr) QRType(svd.cols(), svd.rows());
221 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
222 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
223 m_adjoint.resize(svd.cols(), svd.rows());
228 if(matrix.cols() > matrix.rows())
230 m_adjoint = matrix.adjoint();
231 m_qr.compute(m_adjoint);
233 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).
template triangularView<Upper>().adjoint();
234 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
235 else if(svd.m_computeThinV)
237 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
238 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
240 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
247 typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
249 TransposeTypeWithSameStorageOrder m_adjoint;
250 typename internal::plain_row_type<MatrixType>::type m_workspace;
255 template<
typename MatrixType>
259 typedef typename MatrixType::Index Index;
261 void allocate(
const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
263 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
266 ::new (&m_qr) QRType(svd.rows(), svd.cols());
268 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
269 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
274 if(matrix.rows() > matrix.cols())
276 m_qr.compute(matrix);
277 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).
template triangularView<Upper>();
278 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
279 else if(svd.m_computeThinU)
281 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
282 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
284 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
290 typedef HouseholderQR<MatrixType> QRType;
292 typename internal::plain_col_type<MatrixType>::type m_workspace;
295 template<
typename MatrixType>
299 typedef typename MatrixType::Index Index;
300 typedef typename MatrixType::Scalar Scalar;
303 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
304 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
305 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
306 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
307 Options = MatrixType::Options
310 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
311 TransposeTypeWithSameStorageOrder;
313 void allocate(
const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
315 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
318 ::new (&m_qr) QRType(svd.cols(), svd.rows());
320 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
321 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
322 m_adjoint.resize(svd.cols(), svd.rows());
327 if(matrix.cols() > matrix.rows())
329 m_adjoint = matrix.adjoint();
330 m_qr.compute(m_adjoint);
332 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).
template triangularView<Upper>().adjoint();
333 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
334 else if(svd.m_computeThinV)
336 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
337 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
339 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
346 typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
348 TransposeTypeWithSameStorageOrder m_adjoint;
349 typename internal::plain_row_type<MatrixType>::type m_workspace;
357 template<
typename MatrixType,
int QRPreconditioner>
358 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
360 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
361 typedef typename SVD::Index Index;
362 typedef typename MatrixType::RealScalar RealScalar;
363 static bool run(
typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) {
return true; }
366 template<
typename MatrixType,
int QRPreconditioner>
367 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
369 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
370 typedef typename SVD::Index Index;
371 typedef typename MatrixType::Scalar Scalar;
372 typedef typename MatrixType::RealScalar RealScalar;
373 static bool run(
typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry)
379 JacobiRotation<Scalar> rot;
380 RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
382 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
383 const RealScalar precision = NumTraits<Scalar>::epsilon();
388 work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0);
390 if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
393 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
394 work_matrix.row(p) *= z;
395 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
397 if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
399 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
400 work_matrix.row(q) *= z;
401 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
407 rot.c() = conj(work_matrix.coeff(p,p)) / n;
408 rot.s() = work_matrix.coeff(q,p) / n;
409 work_matrix.applyOnTheLeft(p,q,rot);
410 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
411 if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
413 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
414 work_matrix.col(q) *= z;
415 if(svd.computeV()) svd.m_matrixV.col(q) *= z;
417 if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
419 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
420 work_matrix.row(q) *= z;
421 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
426 maxDiagEntry = max EIGEN_EMPTY (maxDiagEntry,max EIGEN_EMPTY (abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q))));
428 RealScalar threshold = max EIGEN_EMPTY (considerAsZero, precision * maxDiagEntry);
429 return abs(work_matrix.coeff(p,q))>threshold || abs(work_matrix.coeff(q,p)) > threshold;
433 template<
typename MatrixType,
typename RealScalar,
typename Index>
434 void real_2x2_jacobi_svd(
const MatrixType& matrix, Index p, Index q,
435 JacobiRotation<RealScalar> *j_left,
436 JacobiRotation<RealScalar> *j_right)
440 Matrix<RealScalar,2,2> m;
441 m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)),
442 numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q));
443 JacobiRotation<RealScalar> rot1;
444 RealScalar t = m.coeff(0,0) + m.coeff(1,1);
445 RealScalar d = m.coeff(1,0) - m.coeff(0,1);
446 if(d == RealScalar(0))
448 rot1.s() = RealScalar(0);
449 rot1.c() = RealScalar(1);
455 RealScalar u = t / d;
456 RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u));
457 rot1.s() = RealScalar(1) / tmp;
460 m.applyOnTheLeft(0,1,rot1);
461 j_right->makeJacobi(m,0,1);
462 *j_left = rot1 * j_right->transpose();
520 template<
typename _MatrixType,
int QRPreconditioner>
class JacobiSVD
524 typedef _MatrixType MatrixType;
525 typedef typename MatrixType::Scalar Scalar;
526 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
527 typedef typename MatrixType::Index Index;
529 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
530 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
531 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
532 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
533 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
534 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
535 MatrixOptions = MatrixType::Options
538 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
539 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>
541 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
542 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>
544 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
545 typedef typename internal::plain_row_type<MatrixType>::type RowType;
546 typedef typename internal::plain_col_type<MatrixType>::type ColType;
547 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
548 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
557 : m_isInitialized(false),
558 m_isAllocated(false),
559 m_usePrescribedThreshold(false),
560 m_computationOptions(0),
561 m_rows(-1), m_cols(-1), m_diagSize(0)
571 JacobiSVD(Index rows, Index cols,
unsigned int computationOptions = 0)
572 : m_isInitialized(false),
573 m_isAllocated(false),
574 m_usePrescribedThreshold(false),
575 m_computationOptions(0),
576 m_rows(-1), m_cols(-1)
578 allocate(rows, cols, computationOptions);
591 JacobiSVD(
const MatrixType& matrix,
unsigned int computationOptions = 0)
592 : m_isInitialized(false),
593 m_isAllocated(false),
594 m_usePrescribedThreshold(false),
595 m_computationOptions(0),
596 m_rows(-1), m_cols(-1)
598 compute(matrix, computationOptions);
611 JacobiSVD& compute(
const MatrixType& matrix,
unsigned int computationOptions);
621 return compute(matrix, m_computationOptions);
635 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
636 eigen_assert(computeU() &&
"This JacobiSVD decomposition didn't compute U. Did you ask for it?");
651 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
652 eigen_assert(computeV() &&
"This JacobiSVD decomposition didn't compute V. Did you ask for it?");
663 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
664 return m_singularValues;
668 inline bool computeU()
const {
return m_computeFullU || m_computeThinU; }
670 inline bool computeV()
const {
return m_computeFullV || m_computeThinV; }
681 template<
typename Rhs>
682 inline const internal::solve_retval<JacobiSVD, Rhs>
685 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
686 eigen_assert(computeU() && computeV() &&
"JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
687 return internal::solve_retval<JacobiSVD, Rhs>(*
this, b.derived());
693 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
694 return m_nonzeroSingularValues;
706 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
707 if(m_singularValues.size()==0)
return 0;
708 RealScalar premultiplied_threshold = m_singularValues.coeff(0) * threshold();
709 Index i = m_nonzeroSingularValues-1;
710 while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i;
730 m_usePrescribedThreshold =
true;
731 m_prescribedThreshold = threshold;
745 m_usePrescribedThreshold =
false;
755 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
756 return m_usePrescribedThreshold ? m_prescribedThreshold
760 inline Index rows()
const {
return m_rows; }
761 inline Index cols()
const {
return m_cols; }
764 void allocate(Index rows, Index cols,
unsigned int computationOptions);
766 static void check_template_parameters()
768 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
772 MatrixUType m_matrixU;
773 MatrixVType m_matrixV;
774 SingularValuesType m_singularValues;
775 WorkMatrixType m_workMatrix;
776 bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold;
777 bool m_computeFullU, m_computeThinU;
778 bool m_computeFullV, m_computeThinV;
779 unsigned int m_computationOptions;
780 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
781 RealScalar m_prescribedThreshold;
783 template<
typename __MatrixType,
int _QRPreconditioner,
bool _IsComplex>
784 friend struct internal::svd_precondition_2x2_block_to_be_real;
785 template<
typename __MatrixType,
int _QRPreconditioner,
int _Case,
bool _DoAnything>
786 friend struct internal::qr_preconditioner_impl;
788 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
789 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
790 MatrixType m_scaledMatrix;
793 template<
typename MatrixType,
int QRPreconditioner>
796 eigen_assert(rows >= 0 && cols >= 0);
801 computationOptions == m_computationOptions)
808 m_isInitialized =
false;
809 m_isAllocated =
true;
810 m_computationOptions = computationOptions;
811 m_computeFullU = (computationOptions &
ComputeFullU) != 0;
812 m_computeThinU = (computationOptions &
ComputeThinU) != 0;
813 m_computeFullV = (computationOptions &
ComputeFullV) != 0;
814 m_computeThinV = (computationOptions &
ComputeThinV) != 0;
815 eigen_assert(!(m_computeFullU && m_computeThinU) &&
"JacobiSVD: you can't ask for both full and thin U");
816 eigen_assert(!(m_computeFullV && m_computeThinV) &&
"JacobiSVD: you can't ask for both full and thin V");
817 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
818 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
821 eigen_assert(!(m_computeThinU || m_computeThinV) &&
822 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " 823 "Use the ColPivHouseholderQR preconditioner instead.");
825 m_diagSize = (std::min)(m_rows, m_cols);
826 m_singularValues.resize(m_diagSize);
827 if(RowsAtCompileTime==Dynamic)
828 m_matrixU.resize(m_rows, m_computeFullU ? m_rows
829 : m_computeThinU ? m_diagSize
831 if(ColsAtCompileTime==Dynamic)
832 m_matrixV.resize(m_cols, m_computeFullV ? m_cols
833 : m_computeThinV ? m_diagSize
835 m_workMatrix.resize(m_diagSize, m_diagSize);
837 if(m_cols>m_rows) m_qr_precond_morecols.allocate(*
this);
838 if(m_rows>m_cols) m_qr_precond_morerows.allocate(*
this);
839 if(m_rows!=m_cols) m_scaledMatrix.resize(rows,cols);
842 template<
typename MatrixType,
int QRPreconditioner>
846 check_template_parameters();
850 allocate(matrix.rows(), matrix.cols(), computationOptions);
857 const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
860 RealScalar scale = matrix.cwiseAbs().maxCoeff();
861 if(scale==RealScalar(0)) scale = RealScalar(1);
867 m_scaledMatrix = matrix / scale;
868 m_qr_precond_morecols.run(*
this, m_scaledMatrix);
869 m_qr_precond_morerows.run(*
this, m_scaledMatrix);
873 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
874 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
875 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
876 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
877 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
881 RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff();
883 bool finished =
false;
890 for(Index p = 1; p < m_diagSize; ++p)
892 for(Index q = 0; q < p; ++q)
897 RealScalar threshold = max EIGEN_EMPTY (considerAsZero, precision * maxDiagEntry);
898 if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
903 if(internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *
this, p, q, maxDiagEntry))
906 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
909 m_workMatrix.applyOnTheLeft(p,q,j_left);
910 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.
transpose());
912 m_workMatrix.applyOnTheRight(p,q,j_right);
913 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
916 maxDiagEntry = max EIGEN_EMPTY (maxDiagEntry,max EIGEN_EMPTY (abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q))));
925 for(Index i = 0; i < m_diagSize; ++i)
927 RealScalar a = abs(m_workMatrix.coeff(i,i));
928 m_singularValues.coeffRef(i) = a;
929 if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
934 m_nonzeroSingularValues = m_diagSize;
935 for(Index i = 0; i < m_diagSize; i++)
938 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
939 if(maxRemainingSingularValue == RealScalar(0))
941 m_nonzeroSingularValues = i;
947 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
948 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
949 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
953 m_singularValues *= scale;
955 m_isInitialized =
true;
960 template<
typename _MatrixType,
int QRPreconditioner,
typename Rhs>
961 struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
962 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
965 EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
967 template<
typename Dest>
void evalTo(Dest& dst)
const 969 eigen_assert(rhs().rows() == dec().rows());
975 Index rank = dec().rank();
977 tmp.
noalias() = dec().matrixU().leftCols(rank).adjoint() * rhs();
978 tmp = dec().singularValues().
head(rank).asDiagonal().inverse() * tmp;
979 dst = dec().matrixV().
leftCols(rank) * tmp;
991 template<
typename Derived>
1000 #endif // EIGEN_JACOBISVD_H Index rank() const
Definition: JacobiSVD.h:703
Definition: Constants.h:365
JacobiSVD< PlainObject > jacobiSvd(unsigned int computationOptions=0) const
Definition: JacobiSVD.h:993
Definition: Constants.h:359
bool computeV() const
Definition: JacobiSVD.h:670
Definition: Constants.h:329
const internal::solve_retval< JacobiSVD, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: JacobiSVD.h:683
RealScalar threshold() const
Definition: JacobiSVD.h:753
NoAlias< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >, Eigen::MatrixBase > noalias()
Rotation given by a cosine-sine pair.
Definition: ForwardDeclarations.h:228
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Definition: Constants.h:331
const SingularValuesType & singularValues() const
Definition: JacobiSVD.h:661
ColsBlockXpr leftCols(Index n)
Definition: DenseBase.h:527
SegmentReturnType head(Index n)
Definition: DenseBase.h:806
JacobiRotation transpose() const
Definition: Jacobi.h:59
JacobiSVD()
Default Constructor.
Definition: JacobiSVD.h:556
const MatrixUType & matrixU() const
Definition: JacobiSVD.h:633
JacobiSVD(Index rows, Index cols, unsigned int computationOptions=0)
Default Constructor with memory preallocation.
Definition: JacobiSVD.h:571
JacobiSVD & compute(const MatrixType &matrix, unsigned int computationOptions)
Method performing the decomposition of given matrix using custom options.
Definition: JacobiSVD.h:844
Definition: Constants.h:333
Index nonzeroSingularValues() const
Definition: JacobiSVD.h:691
Definition: Eigen_Colamd.h:50
const MatrixVType & matrixV() const
Definition: JacobiSVD.h:649
bool computeU() const
Definition: JacobiSVD.h:668
Two-sided Jacobi SVD decomposition of a rectangular matrix.
Definition: ForwardDeclarations.h:224
JacobiSVD & setThreshold(Default_t)
Definition: JacobiSVD.h:743
JacobiSVD & setThreshold(const RealScalar &threshold)
Definition: JacobiSVD.h:728
JacobiSVD(const MatrixType &matrix, unsigned int computationOptions=0)
Constructor performing the decomposition of given matrix.
Definition: JacobiSVD.h:591
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
Definition: Constants.h:327
Definition: Constants.h:361
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Definition: Constants.h:363
JacobiSVD & compute(const MatrixType &matrix)
Method performing the decomposition of given matrix using current options.
Definition: JacobiSVD.h:619