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 static void run(
typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
365 template<
typename MatrixType,
int QRPreconditioner>
366 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
368 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
369 typedef typename MatrixType::Scalar Scalar;
370 typedef typename MatrixType::RealScalar RealScalar;
371 typedef typename SVD::Index Index;
372 static void run(
typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q)
376 JacobiRotation<Scalar> rot;
377 RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
380 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
381 work_matrix.row(p) *= z;
382 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
383 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
384 work_matrix.row(q) *= z;
385 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
389 rot.c() = conj(work_matrix.coeff(p,p)) / n;
390 rot.s() = work_matrix.coeff(q,p) / n;
391 work_matrix.applyOnTheLeft(p,q,rot);
392 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
393 if(work_matrix.coeff(p,q) != Scalar(0))
395 Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
396 work_matrix.col(q) *= z;
397 if(svd.computeV()) svd.m_matrixV.col(q) *= z;
399 if(work_matrix.coeff(q,q) != Scalar(0))
401 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
402 work_matrix.row(q) *= z;
403 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
409 template<
typename MatrixType,
typename RealScalar,
typename Index>
410 void real_2x2_jacobi_svd(
const MatrixType& matrix, Index p, Index q,
411 JacobiRotation<RealScalar> *j_left,
412 JacobiRotation<RealScalar> *j_right)
415 Matrix<RealScalar,2,2> m;
416 m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)),
417 numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q));
418 JacobiRotation<RealScalar> rot1;
419 RealScalar t = m.coeff(0,0) + m.coeff(1,1);
420 RealScalar d = m.coeff(1,0) - m.coeff(0,1);
421 if(t == RealScalar(0))
423 rot1.c() = RealScalar(0);
424 rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
428 RealScalar u = d / t;
429 rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u));
430 rot1.s() = rot1.c() * u;
432 m.applyOnTheLeft(0,1,rot1);
433 j_right->makeJacobi(m,0,1);
434 *j_left = rot1 * j_right->transpose();
492 template<
typename _MatrixType,
int QRPreconditioner>
class JacobiSVD
496 typedef _MatrixType MatrixType;
497 typedef typename MatrixType::Scalar Scalar;
498 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
499 typedef typename MatrixType::Index Index;
501 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
502 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
503 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
504 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
505 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
506 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
507 MatrixOptions = MatrixType::Options
510 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
511 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>
513 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
514 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>
516 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
517 typedef typename internal::plain_row_type<MatrixType>::type RowType;
518 typedef typename internal::plain_col_type<MatrixType>::type ColType;
519 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
520 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
529 : m_isInitialized(false),
530 m_isAllocated(false),
531 m_computationOptions(0),
532 m_rows(-1), m_cols(-1)
542 JacobiSVD(Index rows, Index cols,
unsigned int computationOptions = 0)
543 : m_isInitialized(false),
544 m_isAllocated(false),
545 m_computationOptions(0),
546 m_rows(-1), m_cols(-1)
548 allocate(rows, cols, computationOptions);
561 JacobiSVD(
const MatrixType& matrix,
unsigned int computationOptions = 0)
562 : m_isInitialized(false),
563 m_isAllocated(false),
564 m_computationOptions(0),
565 m_rows(-1), m_cols(-1)
567 compute(matrix, computationOptions);
580 JacobiSVD& compute(
const MatrixType& matrix,
unsigned int computationOptions);
590 return compute(matrix, m_computationOptions);
604 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
605 eigen_assert(computeU() &&
"This JacobiSVD decomposition didn't compute U. Did you ask for it?");
620 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
621 eigen_assert(computeV() &&
"This JacobiSVD decomposition didn't compute V. Did you ask for it?");
632 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
633 return m_singularValues;
637 inline bool computeU()
const {
return m_computeFullU || m_computeThinU; }
639 inline bool computeV()
const {
return m_computeFullV || m_computeThinV; }
650 template<
typename Rhs>
651 inline const internal::solve_retval<JacobiSVD, Rhs>
654 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
655 eigen_assert(computeU() && computeV() &&
"JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
656 return internal::solve_retval<JacobiSVD, Rhs>(*
this, b.derived());
662 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
663 return m_nonzeroSingularValues;
666 inline Index rows()
const {
return m_rows; }
667 inline Index cols()
const {
return m_cols; }
670 void allocate(Index rows, Index cols,
unsigned int computationOptions);
673 MatrixUType m_matrixU;
674 MatrixVType m_matrixV;
675 SingularValuesType m_singularValues;
676 WorkMatrixType m_workMatrix;
677 bool m_isInitialized, m_isAllocated;
678 bool m_computeFullU, m_computeThinU;
679 bool m_computeFullV, m_computeThinV;
680 unsigned int m_computationOptions;
681 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
683 template<
typename __MatrixType,
int _QRPreconditioner,
bool _IsComplex>
684 friend struct internal::svd_precondition_2x2_block_to_be_real;
685 template<
typename __MatrixType,
int _QRPreconditioner,
int _Case,
bool _DoAnything>
686 friend struct internal::qr_preconditioner_impl;
688 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
689 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
692 template<
typename MatrixType,
int QRPreconditioner>
693 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols,
unsigned int computationOptions)
695 eigen_assert(rows >= 0 && cols >= 0);
700 computationOptions == m_computationOptions)
707 m_isInitialized =
false;
708 m_isAllocated =
true;
709 m_computationOptions = computationOptions;
710 m_computeFullU = (computationOptions &
ComputeFullU) != 0;
711 m_computeThinU = (computationOptions &
ComputeThinU) != 0;
712 m_computeFullV = (computationOptions &
ComputeFullV) != 0;
713 m_computeThinV = (computationOptions &
ComputeThinV) != 0;
714 eigen_assert(!(m_computeFullU && m_computeThinU) &&
"JacobiSVD: you can't ask for both full and thin U");
715 eigen_assert(!(m_computeFullV && m_computeThinV) &&
"JacobiSVD: you can't ask for both full and thin V");
716 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==
Dynamic) &&
717 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
720 eigen_assert(!(m_computeThinU || m_computeThinV) &&
721 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
722 "Use the ColPivHouseholderQR preconditioner instead.");
724 m_diagSize = (std::min)(m_rows, m_cols);
725 m_singularValues.resize(m_diagSize);
727 m_matrixU.resize(m_rows, m_computeFullU ? m_rows
728 : m_computeThinU ? m_diagSize
731 m_matrixV.resize(m_cols, m_computeFullV ? m_cols
732 : m_computeThinV ? m_diagSize
734 m_workMatrix.resize(m_diagSize, m_diagSize);
736 if(m_cols>m_rows) m_qr_precond_morecols.allocate(*
this);
737 if(m_rows>m_cols) m_qr_precond_morerows.allocate(*
this);
740 template<
typename MatrixType,
int QRPreconditioner>
741 JacobiSVD<MatrixType, QRPreconditioner>&
745 allocate(matrix.rows(), matrix.cols(), computationOptions);
752 const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
756 if(!m_qr_precond_morecols.run(*
this, matrix) && !m_qr_precond_morerows.run(*
this, matrix))
758 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
759 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
760 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
761 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
762 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
767 bool finished =
false;
774 for(Index p = 1; p < m_diagSize; ++p)
776 for(Index q = 0; q < p; ++q)
782 RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)),
783 abs(m_workMatrix.coeff(q,q))));
784 if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold)
789 internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *
this, p, q);
791 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
794 m_workMatrix.applyOnTheLeft(p,q,j_left);
795 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.
transpose());
797 m_workMatrix.applyOnTheRight(p,q,j_right);
798 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
806 for(Index i = 0; i < m_diagSize; ++i)
808 RealScalar a = abs(m_workMatrix.coeff(i,i));
809 m_singularValues.coeffRef(i) = a;
810 if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
815 m_nonzeroSingularValues = m_diagSize;
816 for(Index i = 0; i < m_diagSize; i++)
819 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
820 if(maxRemainingSingularValue == RealScalar(0))
822 m_nonzeroSingularValues = i;
828 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
829 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
830 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
834 m_isInitialized =
true;
839 template<
typename _MatrixType,
int QRPreconditioner,
typename Rhs>
840 struct solve_retval<
JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
841 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
844 EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
846 template<typename Dest>
void evalTo(Dest& dst)
const
848 eigen_assert(rhs().rows() == dec().rows());
854 Index nonzeroSingVals = dec().nonzeroSingularValues();
856 tmp.
noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs();
857 tmp = dec().singularValues().
head(nonzeroSingVals).asDiagonal().inverse() * tmp;
858 dst = dec().matrixV().
leftCols(nonzeroSingVals) * tmp;
870 template<
typename Derived>
871 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
879 #endif // EIGEN_JACOBISVD_H