00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_JACOBISVD_H
00026 #define EIGEN_JACOBISVD_H
00027
00028 namespace Eigen {
00029
00030 namespace internal {
00031
00032
00033 template<typename MatrixType, int QRPreconditioner,
00034 bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
00035 struct svd_precondition_2x2_block_to_be_real {};
00036
00037
00038
00039
00040
00041
00042
00043
00044 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
00045
00046 template<typename MatrixType, int QRPreconditioner, int Case>
00047 struct qr_preconditioner_should_do_anything
00048 {
00049 enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
00050 MatrixType::ColsAtCompileTime != Dynamic &&
00051 MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
00052 b = MatrixType::RowsAtCompileTime != Dynamic &&
00053 MatrixType::ColsAtCompileTime != Dynamic &&
00054 MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
00055 ret = !( (QRPreconditioner == NoQRPreconditioner) ||
00056 (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
00057 (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
00058 };
00059 };
00060
00061 template<typename MatrixType, int QRPreconditioner, int Case,
00062 bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
00063 > struct qr_preconditioner_impl {};
00064
00065 template<typename MatrixType, int QRPreconditioner, int Case>
00066 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
00067 {
00068 public:
00069 typedef typename MatrixType::Index Index;
00070 void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {}
00071 bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
00072 {
00073 return false;
00074 }
00075 };
00076
00077
00078
00079 template<typename MatrixType>
00080 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
00081 {
00082 public:
00083 typedef typename MatrixType::Index Index;
00084 typedef typename MatrixType::Scalar Scalar;
00085 enum
00086 {
00087 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00088 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
00089 };
00090 typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
00091
00092 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
00093 {
00094 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
00095 {
00096 m_qr = FullPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
00097 }
00098 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
00099 }
00100
00101 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00102 {
00103 if(matrix.rows() > matrix.cols())
00104 {
00105 m_qr.compute(matrix);
00106 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
00107 if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
00108 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
00109 return true;
00110 }
00111 return false;
00112 }
00113 private:
00114 FullPivHouseholderQR<MatrixType> m_qr;
00115 WorkspaceType m_workspace;
00116 };
00117
00118 template<typename MatrixType>
00119 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
00120 {
00121 public:
00122 typedef typename MatrixType::Index Index;
00123 typedef typename MatrixType::Scalar Scalar;
00124 enum
00125 {
00126 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00127 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00128 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00129 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00130 Options = MatrixType::Options
00131 };
00132 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
00133 TransposeTypeWithSameStorageOrder;
00134
00135 void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
00136 {
00137 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
00138 {
00139 m_qr = FullPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
00140 }
00141 m_adjoint.resize(svd.cols(), svd.rows());
00142 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
00143 }
00144
00145 bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00146 {
00147 if(matrix.cols() > matrix.rows())
00148 {
00149 m_adjoint = matrix.adjoint();
00150 m_qr.compute(m_adjoint);
00151 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
00152 if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
00153 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
00154 return true;
00155 }
00156 else return false;
00157 }
00158 private:
00159 FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
00160 TransposeTypeWithSameStorageOrder m_adjoint;
00161 typename internal::plain_row_type<MatrixType>::type m_workspace;
00162 };
00163
00164
00165
00166 template<typename MatrixType>
00167 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
00168 {
00169 public:
00170 typedef typename MatrixType::Index Index;
00171
00172 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
00173 {
00174 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
00175 {
00176 m_qr = ColPivHouseholderQR<MatrixType>(svd.rows(), svd.cols());
00177 }
00178 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
00179 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
00180 }
00181
00182 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00183 {
00184 if(matrix.rows() > matrix.cols())
00185 {
00186 m_qr.compute(matrix);
00187 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
00188 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
00189 else if(svd.m_computeThinU)
00190 {
00191 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
00192 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
00193 }
00194 if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
00195 return true;
00196 }
00197 return false;
00198 }
00199
00200 private:
00201 ColPivHouseholderQR<MatrixType> m_qr;
00202 typename internal::plain_col_type<MatrixType>::type m_workspace;
00203 };
00204
00205 template<typename MatrixType>
00206 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
00207 {
00208 public:
00209 typedef typename MatrixType::Index Index;
00210 typedef typename MatrixType::Scalar Scalar;
00211 enum
00212 {
00213 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00214 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00215 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00216 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00217 Options = MatrixType::Options
00218 };
00219
00220 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
00221 TransposeTypeWithSameStorageOrder;
00222
00223 void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
00224 {
00225 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
00226 {
00227 m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
00228 }
00229 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
00230 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
00231 m_adjoint.resize(svd.cols(), svd.rows());
00232 }
00233
00234 bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00235 {
00236 if(matrix.cols() > matrix.rows())
00237 {
00238 m_adjoint = matrix.adjoint();
00239 m_qr.compute(m_adjoint);
00240
00241 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
00242 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
00243 else if(svd.m_computeThinV)
00244 {
00245 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
00246 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
00247 }
00248 if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
00249 return true;
00250 }
00251 else return false;
00252 }
00253
00254 private:
00255 ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
00256 TransposeTypeWithSameStorageOrder m_adjoint;
00257 typename internal::plain_row_type<MatrixType>::type m_workspace;
00258 };
00259
00260
00261
00262 template<typename MatrixType>
00263 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
00264 {
00265 public:
00266 typedef typename MatrixType::Index Index;
00267
00268 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
00269 {
00270 if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
00271 {
00272 m_qr = HouseholderQR<MatrixType>(svd.rows(), svd.cols());
00273 }
00274 if (svd.m_computeFullU) m_workspace.resize(svd.rows());
00275 else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
00276 }
00277
00278 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00279 {
00280 if(matrix.rows() > matrix.cols())
00281 {
00282 m_qr.compute(matrix);
00283 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
00284 if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
00285 else if(svd.m_computeThinU)
00286 {
00287 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
00288 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
00289 }
00290 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
00291 return true;
00292 }
00293 return false;
00294 }
00295 private:
00296 HouseholderQR<MatrixType> m_qr;
00297 typename internal::plain_col_type<MatrixType>::type m_workspace;
00298 };
00299
00300 template<typename MatrixType>
00301 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
00302 {
00303 public:
00304 typedef typename MatrixType::Index Index;
00305 typedef typename MatrixType::Scalar Scalar;
00306 enum
00307 {
00308 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00309 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00310 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00311 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00312 Options = MatrixType::Options
00313 };
00314
00315 typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
00316 TransposeTypeWithSameStorageOrder;
00317
00318 void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
00319 {
00320 if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
00321 {
00322 m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows());
00323 }
00324 if (svd.m_computeFullV) m_workspace.resize(svd.cols());
00325 else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
00326 m_adjoint.resize(svd.cols(), svd.rows());
00327 }
00328
00329 bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00330 {
00331 if(matrix.cols() > matrix.rows())
00332 {
00333 m_adjoint = matrix.adjoint();
00334 m_qr.compute(m_adjoint);
00335
00336 svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
00337 if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
00338 else if(svd.m_computeThinV)
00339 {
00340 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
00341 m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
00342 }
00343 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
00344 return true;
00345 }
00346 else return false;
00347 }
00348
00349 private:
00350 HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr;
00351 TransposeTypeWithSameStorageOrder m_adjoint;
00352 typename internal::plain_row_type<MatrixType>::type m_workspace;
00353 };
00354
00355
00356
00357
00358
00359
00360 template<typename MatrixType, int QRPreconditioner>
00361 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
00362 {
00363 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
00364 typedef typename SVD::Index Index;
00365 static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
00366 };
00367
00368 template<typename MatrixType, int QRPreconditioner>
00369 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
00370 {
00371 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
00372 typedef typename MatrixType::Scalar Scalar;
00373 typedef typename MatrixType::RealScalar RealScalar;
00374 typedef typename SVD::Index Index;
00375 static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q)
00376 {
00377 Scalar z;
00378 JacobiRotation<Scalar> rot;
00379 RealScalar n = sqrt(abs2(work_matrix.coeff(p,p)) + abs2(work_matrix.coeff(q,p)));
00380 if(n==0)
00381 {
00382 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
00383 work_matrix.row(p) *= z;
00384 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
00385 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
00386 work_matrix.row(q) *= z;
00387 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
00388 }
00389 else
00390 {
00391 rot.c() = conj(work_matrix.coeff(p,p)) / n;
00392 rot.s() = work_matrix.coeff(q,p) / n;
00393 work_matrix.applyOnTheLeft(p,q,rot);
00394 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
00395 if(work_matrix.coeff(p,q) != Scalar(0))
00396 {
00397 Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
00398 work_matrix.col(q) *= z;
00399 if(svd.computeV()) svd.m_matrixV.col(q) *= z;
00400 }
00401 if(work_matrix.coeff(q,q) != Scalar(0))
00402 {
00403 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
00404 work_matrix.row(q) *= z;
00405 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
00406 }
00407 }
00408 }
00409 };
00410
00411 template<typename MatrixType, typename RealScalar, typename Index>
00412 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
00413 JacobiRotation<RealScalar> *j_left,
00414 JacobiRotation<RealScalar> *j_right)
00415 {
00416 Matrix<RealScalar,2,2> m;
00417 m << real(matrix.coeff(p,p)), real(matrix.coeff(p,q)),
00418 real(matrix.coeff(q,p)), real(matrix.coeff(q,q));
00419 JacobiRotation<RealScalar> rot1;
00420 RealScalar t = m.coeff(0,0) + m.coeff(1,1);
00421 RealScalar d = m.coeff(1,0) - m.coeff(0,1);
00422 if(t == RealScalar(0))
00423 {
00424 rot1.c() = RealScalar(0);
00425 rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
00426 }
00427 else
00428 {
00429 RealScalar u = d / t;
00430 rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + abs2(u));
00431 rot1.s() = rot1.c() * u;
00432 }
00433 m.applyOnTheLeft(0,1,rot1);
00434 j_right->makeJacobi(m,0,1);
00435 *j_left = rot1 * j_right->transpose();
00436 }
00437
00438 }
00439
00493 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
00494 {
00495 public:
00496
00497 typedef _MatrixType MatrixType;
00498 typedef typename MatrixType::Scalar Scalar;
00499 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00500 typedef typename MatrixType::Index Index;
00501 enum {
00502 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00503 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00504 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
00505 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00506 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00507 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
00508 MatrixOptions = MatrixType::Options
00509 };
00510
00511 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
00512 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>
00513 MatrixUType;
00514 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
00515 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>
00516 MatrixVType;
00517 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
00518 typedef typename internal::plain_row_type<MatrixType>::type RowType;
00519 typedef typename internal::plain_col_type<MatrixType>::type ColType;
00520 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
00521 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
00522 WorkMatrixType;
00523
00529 JacobiSVD()
00530 : m_isInitialized(false),
00531 m_isAllocated(false),
00532 m_computationOptions(0),
00533 m_rows(-1), m_cols(-1)
00534 {}
00535
00536
00543 JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
00544 : m_isInitialized(false),
00545 m_isAllocated(false),
00546 m_computationOptions(0),
00547 m_rows(-1), m_cols(-1)
00548 {
00549 allocate(rows, cols, computationOptions);
00550 }
00551
00562 JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
00563 : m_isInitialized(false),
00564 m_isAllocated(false),
00565 m_computationOptions(0),
00566 m_rows(-1), m_cols(-1)
00567 {
00568 compute(matrix, computationOptions);
00569 }
00570
00581 JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
00582
00589 JacobiSVD& compute(const MatrixType& matrix)
00590 {
00591 return compute(matrix, m_computationOptions);
00592 }
00593
00603 const MatrixUType& matrixU() const
00604 {
00605 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00606 eigen_assert(computeU() && "This JacobiSVD decomposition didn't compute U. Did you ask for it?");
00607 return m_matrixU;
00608 }
00609
00619 const MatrixVType& matrixV() const
00620 {
00621 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00622 eigen_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?");
00623 return m_matrixV;
00624 }
00625
00631 const SingularValuesType& singularValues() const
00632 {
00633 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00634 return m_singularValues;
00635 }
00636
00638 inline bool computeU() const { return m_computeFullU || m_computeThinU; }
00640 inline bool computeV() const { return m_computeFullV || m_computeThinV; }
00641
00651 template<typename Rhs>
00652 inline const internal::solve_retval<JacobiSVD, Rhs>
00653 solve(const MatrixBase<Rhs>& b) const
00654 {
00655 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00656 eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
00657 return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived());
00658 }
00659
00661 Index nonzeroSingularValues() const
00662 {
00663 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00664 return m_nonzeroSingularValues;
00665 }
00666
00667 inline Index rows() const { return m_rows; }
00668 inline Index cols() const { return m_cols; }
00669
00670 private:
00671 void allocate(Index rows, Index cols, unsigned int computationOptions);
00672
00673 protected:
00674 MatrixUType m_matrixU;
00675 MatrixVType m_matrixV;
00676 SingularValuesType m_singularValues;
00677 WorkMatrixType m_workMatrix;
00678 bool m_isInitialized, m_isAllocated;
00679 bool m_computeFullU, m_computeThinU;
00680 bool m_computeFullV, m_computeThinV;
00681 unsigned int m_computationOptions;
00682 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
00683
00684 template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
00685 friend struct internal::svd_precondition_2x2_block_to_be_real;
00686 template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
00687 friend struct internal::qr_preconditioner_impl;
00688
00689 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
00690 internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
00691 };
00692
00693 template<typename MatrixType, int QRPreconditioner>
00694 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
00695 {
00696 eigen_assert(rows >= 0 && cols >= 0);
00697
00698 if (m_isAllocated &&
00699 rows == m_rows &&
00700 cols == m_cols &&
00701 computationOptions == m_computationOptions)
00702 {
00703 return;
00704 }
00705
00706 m_rows = rows;
00707 m_cols = cols;
00708 m_isInitialized = false;
00709 m_isAllocated = true;
00710 m_computationOptions = computationOptions;
00711 m_computeFullU = (computationOptions & ComputeFullU) != 0;
00712 m_computeThinU = (computationOptions & ComputeThinU) != 0;
00713 m_computeFullV = (computationOptions & ComputeFullV) != 0;
00714 m_computeThinV = (computationOptions & ComputeThinV) != 0;
00715 eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
00716 eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
00717 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
00718 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
00719 if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
00720 {
00721 eigen_assert(!(m_computeThinU || m_computeThinV) &&
00722 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
00723 "Use the ColPivHouseholderQR preconditioner instead.");
00724 }
00725 m_diagSize = (std::min)(m_rows, m_cols);
00726 m_singularValues.resize(m_diagSize);
00727 m_matrixU.resize(m_rows, m_computeFullU ? m_rows
00728 : m_computeThinU ? m_diagSize
00729 : 0);
00730 m_matrixV.resize(m_cols, m_computeFullV ? m_cols
00731 : m_computeThinV ? m_diagSize
00732 : 0);
00733 m_workMatrix.resize(m_diagSize, m_diagSize);
00734
00735 m_qr_precond_morecols.allocate(*this);
00736 m_qr_precond_morerows.allocate(*this);
00737 }
00738
00739 template<typename MatrixType, int QRPreconditioner>
00740 JacobiSVD<MatrixType, QRPreconditioner>&
00741 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
00742 {
00743 allocate(matrix.rows(), matrix.cols(), computationOptions);
00744
00745
00746
00747 const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
00748
00749
00750 const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min();
00751
00752
00753
00754 if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix))
00755 {
00756 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
00757 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
00758 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
00759 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
00760 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
00761 }
00762
00763
00764
00765 bool finished = false;
00766 while(!finished)
00767 {
00768 finished = true;
00769
00770
00771
00772 for(Index p = 1; p < m_diagSize; ++p)
00773 {
00774 for(Index q = 0; q < p; ++q)
00775 {
00776
00777
00778
00779 using std::max;
00780 RealScalar threshold = (max)(considerAsZero, precision * (max)(internal::abs(m_workMatrix.coeff(p,p)),
00781 internal::abs(m_workMatrix.coeff(q,q))));
00782 if((max)(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p))) > threshold)
00783 {
00784 finished = false;
00785
00786
00787 internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q);
00788 JacobiRotation<RealScalar> j_left, j_right;
00789 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
00790
00791
00792 m_workMatrix.applyOnTheLeft(p,q,j_left);
00793 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
00794
00795 m_workMatrix.applyOnTheRight(p,q,j_right);
00796 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
00797 }
00798 }
00799 }
00800 }
00801
00802
00803
00804 for(Index i = 0; i < m_diagSize; ++i)
00805 {
00806 RealScalar a = internal::abs(m_workMatrix.coeff(i,i));
00807 m_singularValues.coeffRef(i) = a;
00808 if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
00809 }
00810
00811
00812
00813 m_nonzeroSingularValues = m_diagSize;
00814 for(Index i = 0; i < m_diagSize; i++)
00815 {
00816 Index pos;
00817 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
00818 if(maxRemainingSingularValue == RealScalar(0))
00819 {
00820 m_nonzeroSingularValues = i;
00821 break;
00822 }
00823 if(pos)
00824 {
00825 pos += i;
00826 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
00827 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
00828 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
00829 }
00830 }
00831
00832 m_isInitialized = true;
00833 return *this;
00834 }
00835
00836 namespace internal {
00837 template<typename _MatrixType, int QRPreconditioner, typename Rhs>
00838 struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
00839 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
00840 {
00841 typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType;
00842 EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
00843
00844 template<typename Dest> void evalTo(Dest& dst) const
00845 {
00846 eigen_assert(rhs().rows() == dec().rows());
00847
00848
00849
00850
00851 Index diagSize = (std::min)(dec().rows(), dec().cols());
00852 typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize);
00853
00854 Index nonzeroSingVals = dec().nonzeroSingularValues();
00855 invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse();
00856 invertedSingVals.tail(diagSize - nonzeroSingVals).setZero();
00857
00858 dst = dec().matrixV().leftCols(diagSize)
00859 * invertedSingVals.asDiagonal()
00860 * dec().matrixU().leftCols(diagSize).adjoint()
00861 * rhs();
00862 }
00863 };
00864 }
00865
00873 template<typename Derived>
00874 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
00875 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
00876 {
00877 return JacobiSVD<PlainObject>(*this, computationOptions);
00878 }
00879
00880 }
00881
00882 #endif // EIGEN_JACOBISVD_H