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 if(work_matrix.coeff(q,q)!=Scalar(0))
384 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
387 work_matrix.row(q) *= z;
388 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
392 rot.c() = conj(work_matrix.coeff(p,p)) / n;
393 rot.s() = work_matrix.coeff(q,p) / n;
394 work_matrix.applyOnTheLeft(p,q,rot);
395 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
396 if(work_matrix.coeff(p,q) != Scalar(0))
398 Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
399 work_matrix.col(q) *= z;
400 if(svd.computeV()) svd.m_matrixV.col(q) *= z;
402 if(work_matrix.coeff(q,q) != Scalar(0))
404 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
405 work_matrix.row(q) *= z;
406 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
412 template<
typename MatrixType,
typename RealScalar,
typename Index>
413 void real_2x2_jacobi_svd(
const MatrixType& matrix, Index p, Index q,
414 JacobiRotation<RealScalar> *j_left,
415 JacobiRotation<RealScalar> *j_right)
418 Matrix<RealScalar,2,2> m;
419 m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)),
420 numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q));
421 JacobiRotation<RealScalar> rot1;
422 RealScalar t = m.coeff(0,0) + m.coeff(1,1);
423 RealScalar d = m.coeff(1,0) - m.coeff(0,1);
424 if(t == RealScalar(0))
426 rot1.c() = RealScalar(0);
427 rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
431 RealScalar u = d / t;
432 rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u));
433 rot1.s() = rot1.c() * u;
435 m.applyOnTheLeft(0,1,rot1);
436 j_right->makeJacobi(m,0,1);
437 *j_left = rot1 * j_right->transpose();
495 template<
typename _MatrixType,
int QRPreconditioner>
class JacobiSVD
499 typedef _MatrixType MatrixType;
500 typedef typename MatrixType::Scalar Scalar;
501 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
502 typedef typename MatrixType::Index Index;
504 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
505 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
506 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
507 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
508 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
509 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
510 MatrixOptions = MatrixType::Options
513 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
514 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>
516 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
517 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>
519 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
520 typedef typename internal::plain_row_type<MatrixType>::type RowType;
521 typedef typename internal::plain_col_type<MatrixType>::type ColType;
522 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
523 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
532 : m_isInitialized(false),
533 m_isAllocated(false),
534 m_computationOptions(0),
535 m_rows(-1), m_cols(-1)
545 JacobiSVD(Index rows, Index cols,
unsigned int computationOptions = 0)
546 : m_isInitialized(false),
547 m_isAllocated(false),
548 m_computationOptions(0),
549 m_rows(-1), m_cols(-1)
551 allocate(rows, cols, computationOptions);
564 JacobiSVD(
const MatrixType& matrix,
unsigned int computationOptions = 0)
565 : m_isInitialized(false),
566 m_isAllocated(false),
567 m_computationOptions(0),
568 m_rows(-1), m_cols(-1)
570 compute(matrix, computationOptions);
583 JacobiSVD& compute(
const MatrixType& matrix,
unsigned int computationOptions);
593 return compute(matrix, m_computationOptions);
607 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
608 eigen_assert(computeU() &&
"This JacobiSVD decomposition didn't compute U. Did you ask for it?");
623 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
624 eigen_assert(computeV() &&
"This JacobiSVD decomposition didn't compute V. Did you ask for it?");
635 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
636 return m_singularValues;
640 inline bool computeU()
const {
return m_computeFullU || m_computeThinU; }
642 inline bool computeV()
const {
return m_computeFullV || m_computeThinV; }
653 template<
typename Rhs>
654 inline const internal::solve_retval<JacobiSVD, Rhs>
657 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
658 eigen_assert(computeU() && computeV() &&
"JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
659 return internal::solve_retval<JacobiSVD, Rhs>(*
this, b.derived());
665 eigen_assert(m_isInitialized &&
"JacobiSVD is not initialized.");
666 return m_nonzeroSingularValues;
669 inline Index rows()
const {
return m_rows; }
670 inline Index cols()
const {
return m_cols; }
673 void allocate(Index rows, Index cols,
unsigned int computationOptions);
676 MatrixUType m_matrixU;
677 MatrixVType m_matrixV;
678 SingularValuesType m_singularValues;
679 WorkMatrixType m_workMatrix;
680 bool m_isInitialized, m_isAllocated;
681 bool m_computeFullU, m_computeThinU;
682 bool m_computeFullV, m_computeThinV;
683 unsigned int m_computationOptions;
684 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
686 template<
typename __MatrixType,
int _QRPreconditioner,
bool _IsComplex>
687 friend struct internal::svd_precondition_2x2_block_to_be_real;
688 template<
typename __MatrixType,
int _QRPreconditioner,
int _Case,
bool _DoAnything>
689 friend struct internal::qr_preconditioner_impl;
691 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
692 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
695 template<
typename MatrixType,
int QRPreconditioner>
696 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols,
unsigned int computationOptions)
698 eigen_assert(rows >= 0 && cols >= 0);
703 computationOptions == m_computationOptions)
710 m_isInitialized =
false;
711 m_isAllocated =
true;
712 m_computationOptions = computationOptions;
713 m_computeFullU = (computationOptions &
ComputeFullU) != 0;
714 m_computeThinU = (computationOptions &
ComputeThinU) != 0;
715 m_computeFullV = (computationOptions &
ComputeFullV) != 0;
716 m_computeThinV = (computationOptions &
ComputeThinV) != 0;
717 eigen_assert(!(m_computeFullU && m_computeThinU) &&
"JacobiSVD: you can't ask for both full and thin U");
718 eigen_assert(!(m_computeFullV && m_computeThinV) &&
"JacobiSVD: you can't ask for both full and thin V");
719 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==
Dynamic) &&
720 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
723 eigen_assert(!(m_computeThinU || m_computeThinV) &&
724 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
725 "Use the ColPivHouseholderQR preconditioner instead.");
727 m_diagSize = (std::min)(m_rows, m_cols);
728 m_singularValues.resize(m_diagSize);
730 m_matrixU.resize(m_rows, m_computeFullU ? m_rows
731 : m_computeThinU ? m_diagSize
734 m_matrixV.resize(m_cols, m_computeFullV ? m_cols
735 : m_computeThinV ? m_diagSize
737 m_workMatrix.resize(m_diagSize, m_diagSize);
739 if(m_cols>m_rows) m_qr_precond_morecols.allocate(*
this);
740 if(m_rows>m_cols) m_qr_precond_morerows.allocate(*
this);
743 template<
typename MatrixType,
int QRPreconditioner>
744 JacobiSVD<MatrixType, QRPreconditioner>&
748 allocate(matrix.rows(), matrix.cols(), computationOptions);
755 const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
759 if(!m_qr_precond_morecols.run(*
this, matrix) && !m_qr_precond_morerows.run(*
this, matrix))
761 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
762 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
763 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
764 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
765 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
770 bool finished =
false;
777 for(Index p = 1; p < m_diagSize; ++p)
779 for(Index q = 0; q < p; ++q)
785 RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)),
786 abs(m_workMatrix.coeff(q,q))));
787 if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold)
792 internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *
this, p, q);
794 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
797 m_workMatrix.applyOnTheLeft(p,q,j_left);
798 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.
transpose());
800 m_workMatrix.applyOnTheRight(p,q,j_right);
801 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
809 for(Index i = 0; i < m_diagSize; ++i)
811 RealScalar a = abs(m_workMatrix.coeff(i,i));
812 m_singularValues.coeffRef(i) = a;
813 if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
818 m_nonzeroSingularValues = m_diagSize;
819 for(Index i = 0; i < m_diagSize; i++)
822 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
823 if(maxRemainingSingularValue == RealScalar(0))
825 m_nonzeroSingularValues = i;
831 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
832 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
833 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
837 m_isInitialized =
true;
842 template<
typename _MatrixType,
int QRPreconditioner,
typename Rhs>
843 struct solve_retval<
JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
844 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
847 EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
849 template<typename Dest>
void evalTo(Dest& dst)
const
851 eigen_assert(rhs().rows() == dec().rows());
857 Index nonzeroSingVals = dec().nonzeroSingularValues();
859 tmp.
noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs();
860 tmp = dec().singularValues().
head(nonzeroSingVals).asDiagonal().inverse() * tmp;
861 dst = dec().matrixV().
leftCols(nonzeroSingVals) * tmp;
873 template<
typename Derived>
874 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
882 #endif // EIGEN_JACOBISVD_H
Definition: Constants.h:365
JacobiSVD< PlainObject > jacobiSvd(unsigned int computationOptions=0) const
Definition: JacobiSVD.h:875
Definition: Constants.h:359
bool computeV() const
Definition: JacobiSVD.h:642
Definition: Constants.h:329
const internal::solve_retval< JacobiSVD, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: JacobiSVD.h:655
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
const int Dynamic
Definition: Constants.h:21
Definition: Constants.h:331
const SingularValuesType & singularValues() const
Definition: JacobiSVD.h:633
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:531
const MatrixUType & matrixU() const
Definition: JacobiSVD.h:605
JacobiSVD(Index rows, Index cols, unsigned int computationOptions=0)
Default Constructor with memory preallocation.
Definition: JacobiSVD.h:545
JacobiSVD & compute(const MatrixType &matrix, unsigned int computationOptions)
Method performing the decomposition of given matrix using custom options.
Definition: JacobiSVD.h:745
Definition: Constants.h:333
Index nonzeroSingularValues() const
Definition: JacobiSVD.h:663
const MatrixVType & matrixV() const
Definition: JacobiSVD.h:621
bool computeU() const
Definition: JacobiSVD.h:640
Two-sided Jacobi SVD decomposition of a rectangular matrix.
Definition: ForwardDeclarations.h:224
JacobiSVD(const MatrixType &matrix, unsigned int computationOptions=0)
Constructor performing the decomposition of given matrix.
Definition: JacobiSVD.h:564
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:591