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())
81 m_qr = FullPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
83 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
86 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
88 if(matrix.rows() > matrix.cols())
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();
99 FullPivHouseholderQR<MatrixType> m_qr;
100 WorkspaceType m_workspace;
103 template<
typename MatrixType>
107 typedef typename MatrixType::Index Index;
108 typedef typename MatrixType::Scalar Scalar;
111 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
112 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
113 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
114 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
115 Options = MatrixType::Options
117 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
118 TransposeTypeWithSameStorageOrder;
120 void allocate(
const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
122 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
124 m_qr = FullPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
126 m_adjoint.resize(svd.cols(), svd.rows());
127 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
130 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
132 if(matrix.cols() > matrix.rows())
134 m_adjoint = matrix.adjoint();
135 m_qr.compute(m_adjoint);
136 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).
template triangularView<Upper>().adjoint();
137 if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
138 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
144 FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
145 TransposeTypeWithSameStorageOrder m_adjoint;
146 typename internal::plain_row_type<MatrixType>::type m_workspace;
151 template<
typename MatrixType>
155 typedef typename MatrixType::Index Index;
157 void allocate(
const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
159 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
161 m_qr = ColPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
163 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
164 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
167 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
169 if(matrix.rows() > matrix.cols())
171 m_qr.compute(matrix);
172 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).
template triangularView<Upper>();
173 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
174 else if(svd.m_computeThinU)
176 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
177 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
179 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
186 ColPivHouseholderQR<MatrixType> m_qr;
187 typename internal::plain_col_type<MatrixType>::type m_workspace;
190 template<
typename MatrixType>
194 typedef typename MatrixType::Index Index;
195 typedef typename MatrixType::Scalar Scalar;
198 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
199 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
200 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
201 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
202 Options = MatrixType::Options
205 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
206 TransposeTypeWithSameStorageOrder;
208 void allocate(
const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
210 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
212 m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
214 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
215 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
216 m_adjoint.resize(svd.cols(), svd.rows());
219 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
221 if(matrix.cols() > matrix.rows())
223 m_adjoint = matrix.adjoint();
224 m_qr.compute(m_adjoint);
226 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).
template triangularView<Upper>().adjoint();
227 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
228 else if(svd.m_computeThinV)
230 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
231 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
233 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
240 ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
241 TransposeTypeWithSameStorageOrder m_adjoint;
242 typename internal::plain_row_type<MatrixType>::type m_workspace;
247 template<
typename MatrixType>
251 typedef typename MatrixType::Index Index;
253 void allocate(
const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
255 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
257 m_qr = HouseholderQR<MatrixType>(svd.rows(), svd.cols());
259 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
260 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
263 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
265 if(matrix.rows() > matrix.cols())
267 m_qr.compute(matrix);
268 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).
template triangularView<Upper>();
269 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
270 else if(svd.m_computeThinU)
272 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
273 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
275 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
281 HouseholderQR<MatrixType> m_qr;
282 typename internal::plain_col_type<MatrixType>::type m_workspace;
285 template<
typename MatrixType>
289 typedef typename MatrixType::Index Index;
290 typedef typename MatrixType::Scalar Scalar;
293 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
294 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
295 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
296 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
297 Options = MatrixType::Options
300 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
301 TransposeTypeWithSameStorageOrder;
303 void allocate(
const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
305 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
307 m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
309 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
310 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
311 m_adjoint.resize(svd.cols(), svd.rows());
314 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd,
const MatrixType& matrix)
316 if(matrix.cols() > matrix.rows())
318 m_adjoint = matrix.adjoint();
319 m_qr.compute(m_adjoint);
321 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).
template triangularView<Upper>().adjoint();
322 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
323 else if(svd.m_computeThinV)
325 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
326 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
328 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
335 HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
336 TransposeTypeWithSameStorageOrder m_adjoint;
337 typename internal::plain_row_type<MatrixType>::type m_workspace;
345 template<
typename MatrixType,
int QRPreconditioner>
346 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
348 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
349 typedef typename SVD::Index Index;
350 static void run(
typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
353 template<
typename MatrixType,
int QRPreconditioner>
354 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
356 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
357 typedef typename MatrixType::Scalar Scalar;
358 typedef typename MatrixType::RealScalar RealScalar;
359 typedef typename SVD::Index Index;
360 static void run(
typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q)
363 JacobiRotation<Scalar> rot;
364 RealScalar n = sqrt(abs2(work_matrix.coeff(p,p)) + abs2(work_matrix.coeff(q,p)));
367 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
368 work_matrix.row(p) *= z;
369 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
370 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
371 work_matrix.row(q) *= z;
372 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
376 rot.c() = conj(work_matrix.coeff(p,p)) / n;
377 rot.s() = work_matrix.coeff(q,p) / n;
378 work_matrix.applyOnTheLeft(p,q,rot);
379 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
380 if(work_matrix.coeff(p,q) != Scalar(0))
382 Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
383 work_matrix.col(q) *= z;
384 if(svd.computeV()) svd.m_matrixV.col(q) *= z;
386 if(work_matrix.coeff(q,q) != Scalar(0))
388 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
389 work_matrix.row(q) *= z;
390 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
396 template<
typename MatrixType,
typename RealScalar,
typename Index>
397 void real_2x2_jacobi_svd(
const MatrixType& matrix, Index p, Index q,
398 JacobiRotation<RealScalar> *j_left,
399 JacobiRotation<RealScalar> *j_right)
401 Matrix<RealScalar,2,2> m;
402 m << real(matrix.coeff(p,p)), real(matrix.coeff(p,q)),
403 real(matrix.coeff(q,p)), real(matrix.coeff(q,q));
404 JacobiRotation<RealScalar> rot1;
405 RealScalar t = m.coeff(0,0) + m.coeff(1,1);
406 RealScalar d = m.coeff(1,0) - m.coeff(0,1);
407 if(t == RealScalar(0))
409 rot1.c() = RealScalar(0);
410 rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
414 RealScalar u = d / t;
415 rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + abs2(u));
416 rot1.s() = rot1.c() * u;
418 m.applyOnTheLeft(0,1,rot1);
419 j_right->makeJacobi(m,0,1);
420 *j_left = rot1 * j_right->transpose();
478 template<
typename _MatrixType,
int QRPreconditioner>
class JacobiSVD
482 typedef _MatrixType MatrixType;
483 typedef typename MatrixType::Scalar Scalar;
484 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
485 typedef typename MatrixType::Index Index;
487 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
488 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
489 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
490 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
491 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
492 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
493 MatrixOptions = MatrixType::Options
496 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
497 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>
499 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
500 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>
502 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
503 typedef typename internal::plain_row_type<MatrixType>::type RowType;
504 typedef typename internal::plain_col_type<MatrixType>::type ColType;
505 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
506 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
515 : m_isInitialized(false),
516 m_isAllocated(false),
517 m_computationOptions(0),
518 m_rows(-1), m_cols(-1)
528 JacobiSVD(Index rows, Index cols,
unsigned int computationOptions = 0)
529 : m_isInitialized(false),
530 m_isAllocated(false),
531 m_computationOptions(0),
532 m_rows(-1), m_cols(-1)
534 allocate(rows, cols, computationOptions);
547 JacobiSVD(
const MatrixType& matrix,
unsigned int computationOptions = 0)
548 : m_isInitialized(false),
549 m_isAllocated(false),
550 m_computationOptions(0),
551 m_rows(-1), m_cols(-1)
553 compute(matrix, computationOptions);
566 JacobiSVD&
compute(
const MatrixType& matrix,
unsigned int computationOptions);
576 return compute(matrix, m_computationOptions);
590 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
591 eigen_assert(
computeU() &&
"This JacobiSVD decomposition didn't compute U. Did you ask for it?");
606 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
607 eigen_assert(
computeV() &&
"This JacobiSVD decomposition didn't compute V. Did you ask for it?");
618 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
619 return m_singularValues;
623 inline bool computeU()
const {
return m_computeFullU || m_computeThinU; }
625 inline bool computeV()
const {
return m_computeFullV || m_computeThinV; }
636 template<
typename Rhs>
637 inline const internal::solve_retval<JacobiSVD, Rhs>
640 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
641 eigen_assert(
computeU() &&
computeV() &&
"JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
642 return internal::solve_retval<JacobiSVD, Rhs>(*
this, b.derived());
648 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
649 return m_nonzeroSingularValues;
652 inline Index rows()
const {
return m_rows; }
653 inline Index cols()
const {
return m_cols; }
656 void allocate(Index rows, Index cols,
unsigned int computationOptions);
659 MatrixUType m_matrixU;
660 MatrixVType m_matrixV;
661 SingularValuesType m_singularValues;
662 WorkMatrixType m_workMatrix;
663 bool m_isInitialized, m_isAllocated;
664 bool m_computeFullU, m_computeThinU;
665 bool m_computeFullV, m_computeThinV;
666 unsigned int m_computationOptions;
667 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
669 template<
typename __MatrixType,
int _QRPreconditioner,
bool _IsComplex>
670 friend struct internal::svd_precondition_2x2_block_to_be_real;
671 template<
typename __MatrixType,
int _QRPreconditioner,
int _Case,
bool _DoAnything>
672 friend struct internal::qr_preconditioner_impl;
674 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
675 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
678 template<
typename MatrixType,
int QRPreconditioner>
679 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols,
unsigned int computationOptions)
681 eigen_assert(rows >= 0 && cols >= 0);
686 computationOptions == m_computationOptions)
693 m_isInitialized =
false;
694 m_isAllocated =
true;
695 m_computationOptions = computationOptions;
696 m_computeFullU = (computationOptions &
ComputeFullU) != 0;
697 m_computeThinU = (computationOptions &
ComputeThinU) != 0;
698 m_computeFullV = (computationOptions &
ComputeFullV) != 0;
699 m_computeThinV = (computationOptions &
ComputeThinV) != 0;
700 eigen_assert(!(m_computeFullU && m_computeThinU) &&
"JacobiSVD: you can't ask for both full and thin U");
701 eigen_assert(!(m_computeFullV && m_computeThinV) &&
"JacobiSVD: you can't ask for both full and thin V");
702 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==
Dynamic) &&
703 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
706 eigen_assert(!(m_computeThinU || m_computeThinV) &&
707 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
708 "Use the ColPivHouseholderQR preconditioner instead.");
710 m_diagSize = (std::min)(m_rows, m_cols);
711 m_singularValues.resize(m_diagSize);
712 m_matrixU.resize(m_rows, m_computeFullU ? m_rows
713 : m_computeThinU ? m_diagSize
715 m_matrixV.resize(m_cols, m_computeFullV ? m_cols
716 : m_computeThinV ? m_diagSize
718 m_workMatrix.resize(m_diagSize, m_diagSize);
720 if(m_cols>m_rows) m_qr_precond_morecols.allocate(*
this);
721 if(m_rows>m_cols) m_qr_precond_morerows.allocate(*
this);
724 template<
typename MatrixType,
int QRPreconditioner>
725 JacobiSVD<MatrixType, QRPreconditioner>&
728 allocate(matrix.rows(), matrix.cols(), computationOptions);
735 const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
739 if(!m_qr_precond_morecols.run(*
this, matrix) && !m_qr_precond_morerows.run(*
this, matrix))
741 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
742 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
743 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
744 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
745 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
750 bool finished =
false;
757 for(Index p = 1; p < m_diagSize; ++p)
759 for(Index q = 0; q < p; ++q)
765 RealScalar threshold = (max)(considerAsZero, precision * (max)(internal::abs(m_workMatrix.coeff(p,p)),
766 internal::abs(m_workMatrix.coeff(q,q))));
767 if((max)(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p))) > threshold)
772 internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *
this, p, q);
774 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
777 m_workMatrix.applyOnTheLeft(p,q,j_left);
778 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.
transpose());
780 m_workMatrix.applyOnTheRight(p,q,j_right);
781 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
789 for(Index i = 0; i < m_diagSize; ++i)
791 RealScalar a = internal::abs(m_workMatrix.coeff(i,i));
792 m_singularValues.coeffRef(i) = a;
793 if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
798 m_nonzeroSingularValues = m_diagSize;
799 for(Index i = 0; i < m_diagSize; i++)
802 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
803 if(maxRemainingSingularValue == RealScalar(0))
805 m_nonzeroSingularValues = i;
811 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
812 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
813 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
817 m_isInitialized =
true;
822 template<
typename _MatrixType,
int QRPreconditioner,
typename Rhs>
823 struct solve_retval<
JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
824 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
827 EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
829 template<typename Dest>
void evalTo(Dest& dst)
const
831 eigen_assert(rhs().rows() == dec().rows());
837 Index diagSize = (std::min)(dec().rows(), dec().cols());
838 Index nonzeroSingVals = dec().nonzeroSingularValues();
840 tmp.
noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs();
841 tmp = dec().singularValues().
head(nonzeroSingVals).asDiagonal().inverse() * tmp;
842 dst = dec().matrixV().
leftCols(nonzeroSingVals) * tmp;
854 template<
typename Derived>
855 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
863 #endif // EIGEN_JACOBISVD_H