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 internal {
00029
00030
00031 template<typename MatrixType, int QRPreconditioner,
00032 bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
00033 struct svd_precondition_2x2_block_to_be_real {};
00034
00035
00036
00037
00038
00039
00040
00041
00042 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
00043
00044 template<typename MatrixType, int QRPreconditioner, int Case>
00045 struct qr_preconditioner_should_do_anything
00046 {
00047 enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
00048 MatrixType::ColsAtCompileTime != Dynamic &&
00049 MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
00050 b = MatrixType::RowsAtCompileTime != Dynamic &&
00051 MatrixType::ColsAtCompileTime != Dynamic &&
00052 MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
00053 ret = !( (QRPreconditioner == NoQRPreconditioner) ||
00054 (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
00055 (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
00056 };
00057 };
00058
00059 template<typename MatrixType, int QRPreconditioner, int Case,
00060 bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
00061 > struct qr_preconditioner_impl {};
00062
00063 template<typename MatrixType, int QRPreconditioner, int Case>
00064 struct qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
00065 {
00066 static bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
00067 {
00068 return false;
00069 }
00070 };
00071
00072
00073
00074 template<typename MatrixType>
00075 struct qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
00076 {
00077 static bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00078 {
00079 if(matrix.rows() > matrix.cols())
00080 {
00081 FullPivHouseholderQR<MatrixType> qr(matrix);
00082 svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
00083 if(svd.m_computeFullU) svd.m_matrixU = qr.matrixQ();
00084 if(svd.computeV()) svd.m_matrixV = qr.colsPermutation();
00085 return true;
00086 }
00087 return false;
00088 }
00089 };
00090
00091 template<typename MatrixType>
00092 struct qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
00093 {
00094 static bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00095 {
00096 if(matrix.cols() > matrix.rows())
00097 {
00098 typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
00099 MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
00100 TransposeTypeWithSameStorageOrder;
00101 FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
00102 svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
00103 if(svd.m_computeFullV) svd.m_matrixV = qr.matrixQ();
00104 if(svd.computeU()) svd.m_matrixU = qr.colsPermutation();
00105 return true;
00106 }
00107 else return false;
00108 }
00109 };
00110
00111
00112
00113 template<typename MatrixType>
00114 struct qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
00115 {
00116 static bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00117 {
00118 if(matrix.rows() > matrix.cols())
00119 {
00120 ColPivHouseholderQR<MatrixType> qr(matrix);
00121 svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
00122 if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ();
00123 else if(svd.m_computeThinU) {
00124 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
00125 qr.householderQ().applyThisOnTheLeft(svd.m_matrixU);
00126 }
00127 if(svd.computeV()) svd.m_matrixV = qr.colsPermutation();
00128 return true;
00129 }
00130 return false;
00131 }
00132 };
00133
00134 template<typename MatrixType>
00135 struct qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
00136 {
00137 static bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00138 {
00139 if(matrix.cols() > matrix.rows())
00140 {
00141 typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
00142 MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
00143 TransposeTypeWithSameStorageOrder;
00144 ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
00145 svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
00146 if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ();
00147 else if(svd.m_computeThinV) {
00148 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
00149 qr.householderQ().applyThisOnTheLeft(svd.m_matrixV);
00150 }
00151 if(svd.computeU()) svd.m_matrixU = qr.colsPermutation();
00152 return true;
00153 }
00154 else return false;
00155 }
00156 };
00157
00158
00159
00160 template<typename MatrixType>
00161 struct qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
00162 {
00163 static bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00164 {
00165 if(matrix.rows() > matrix.cols())
00166 {
00167 HouseholderQR<MatrixType> qr(matrix);
00168 svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
00169 if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ();
00170 else if(svd.m_computeThinU) {
00171 svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
00172 qr.householderQ().applyThisOnTheLeft(svd.m_matrixU);
00173 }
00174 if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
00175 return true;
00176 }
00177 return false;
00178 }
00179 };
00180
00181 template<typename MatrixType>
00182 struct qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
00183 {
00184 static bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00185 {
00186 if(matrix.cols() > matrix.rows())
00187 {
00188 typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime,
00189 MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime>
00190 TransposeTypeWithSameStorageOrder;
00191 HouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint());
00192 svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
00193 if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ();
00194 else if(svd.m_computeThinV) {
00195 svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
00196 qr.householderQ().applyThisOnTheLeft(svd.m_matrixV);
00197 }
00198 if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
00199 return true;
00200 }
00201 else return false;
00202 }
00203 };
00204
00205
00206
00207
00208
00209
00210 template<typename MatrixType, int QRPreconditioner>
00211 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
00212 {
00213 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
00214 typedef typename SVD::Index Index;
00215 static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {}
00216 };
00217
00218 template<typename MatrixType, int QRPreconditioner>
00219 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
00220 {
00221 typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
00222 typedef typename MatrixType::Scalar Scalar;
00223 typedef typename MatrixType::RealScalar RealScalar;
00224 typedef typename SVD::Index Index;
00225 static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q)
00226 {
00227 Scalar z;
00228 JacobiRotation<Scalar> rot;
00229 RealScalar n = sqrt(abs2(work_matrix.coeff(p,p)) + abs2(work_matrix.coeff(q,p)));
00230 if(n==0)
00231 {
00232 z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
00233 work_matrix.row(p) *= z;
00234 if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
00235 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
00236 work_matrix.row(q) *= z;
00237 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
00238 }
00239 else
00240 {
00241 rot.c() = conj(work_matrix.coeff(p,p)) / n;
00242 rot.s() = work_matrix.coeff(q,p) / n;
00243 work_matrix.applyOnTheLeft(p,q,rot);
00244 if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
00245 if(work_matrix.coeff(p,q) != Scalar(0))
00246 {
00247 Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
00248 work_matrix.col(q) *= z;
00249 if(svd.computeV()) svd.m_matrixV.col(q) *= z;
00250 }
00251 if(work_matrix.coeff(q,q) != Scalar(0))
00252 {
00253 z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
00254 work_matrix.row(q) *= z;
00255 if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
00256 }
00257 }
00258 }
00259 };
00260
00261 template<typename MatrixType, typename RealScalar, typename Index>
00262 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
00263 JacobiRotation<RealScalar> *j_left,
00264 JacobiRotation<RealScalar> *j_right)
00265 {
00266 Matrix<RealScalar,2,2> m;
00267 m << real(matrix.coeff(p,p)), real(matrix.coeff(p,q)),
00268 real(matrix.coeff(q,p)), real(matrix.coeff(q,q));
00269 JacobiRotation<RealScalar> rot1;
00270 RealScalar t = m.coeff(0,0) + m.coeff(1,1);
00271 RealScalar d = m.coeff(1,0) - m.coeff(0,1);
00272 if(t == RealScalar(0))
00273 {
00274 rot1.c() = 0;
00275 rot1.s() = d > 0 ? 1 : -1;
00276 }
00277 else
00278 {
00279 RealScalar u = d / t;
00280 rot1.c() = RealScalar(1) / sqrt(1 + abs2(u));
00281 rot1.s() = rot1.c() * u;
00282 }
00283 m.applyOnTheLeft(0,1,rot1);
00284 j_right->makeJacobi(m,0,1);
00285 *j_left = rot1 * j_right->transpose();
00286 }
00287
00288 }
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
00344 {
00345 public:
00346
00347 typedef _MatrixType MatrixType;
00348 typedef typename MatrixType::Scalar Scalar;
00349 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00350 typedef typename MatrixType::Index Index;
00351 enum {
00352 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00353 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00354 DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
00355 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00356 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00357 MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
00358 MatrixOptions = MatrixType::Options
00359 };
00360
00361 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime,
00362 MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>
00363 MatrixUType;
00364 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime,
00365 MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>
00366 MatrixVType;
00367 typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
00368 typedef typename internal::plain_row_type<MatrixType>::type RowType;
00369 typedef typename internal::plain_col_type<MatrixType>::type ColType;
00370 typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
00371 MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
00372 WorkMatrixType;
00373
00374
00375
00376
00377
00378
00379 JacobiSVD() : m_isInitialized(false) {}
00380
00381
00382
00383
00384
00385
00386
00387
00388 JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
00389 {
00390 allocate(rows, cols, computationOptions);
00391 }
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403 JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
00404 {
00405 compute(matrix, computationOptions);
00406 }
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418 JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions = 0);
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429 const MatrixUType& matrixU() const
00430 {
00431 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00432 eigen_assert(computeU() && "This JacobiSVD decomposition didn't compute U. Did you ask for it?");
00433 return m_matrixU;
00434 }
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445 const MatrixVType& matrixV() const
00446 {
00447 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00448 eigen_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?");
00449 return m_matrixV;
00450 }
00451
00452
00453
00454
00455
00456
00457 const SingularValuesType& singularValues() const
00458 {
00459 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00460 return m_singularValues;
00461 }
00462
00463
00464 inline bool computeU() const { return m_computeFullU || m_computeThinU; }
00465
00466 inline bool computeV() const { return m_computeFullV || m_computeThinV; }
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 template<typename Rhs>
00478 inline const internal::solve_retval<JacobiSVD, Rhs>
00479 solve(const MatrixBase<Rhs>& b) const
00480 {
00481 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00482 eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice).");
00483 return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived());
00484 }
00485
00486
00487 Index nonzeroSingularValues() const
00488 {
00489 eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
00490 return m_nonzeroSingularValues;
00491 }
00492
00493 inline Index rows() const { return m_rows; }
00494 inline Index cols() const { return m_cols; }
00495
00496 private:
00497 void allocate(Index rows, Index cols, unsigned int computationOptions = 0);
00498
00499 protected:
00500 MatrixUType m_matrixU;
00501 MatrixVType m_matrixV;
00502 SingularValuesType m_singularValues;
00503 WorkMatrixType m_workMatrix;
00504 bool m_isInitialized;
00505 bool m_computeFullU, m_computeThinU;
00506 bool m_computeFullV, m_computeThinV;
00507 Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
00508
00509 template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
00510 friend struct internal::svd_precondition_2x2_block_to_be_real;
00511 template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
00512 friend struct internal::qr_preconditioner_impl;
00513 };
00514
00515 template<typename MatrixType, int QRPreconditioner>
00516 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
00517 {
00518 m_rows = rows;
00519 m_cols = cols;
00520 m_isInitialized = false;
00521 m_computeFullU = (computationOptions & ComputeFullU) != 0;
00522 m_computeThinU = (computationOptions & ComputeThinU) != 0;
00523 m_computeFullV = (computationOptions & ComputeFullV) != 0;
00524 m_computeThinV = (computationOptions & ComputeThinV) != 0;
00525 eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
00526 eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
00527 eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
00528 "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
00529 if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
00530 {
00531 eigen_assert(!(m_computeThinU || m_computeThinV) &&
00532 "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
00533 "Use the ColPivHouseholderQR preconditioner instead.");
00534 }
00535 m_diagSize = std::min(m_rows, m_cols);
00536 m_singularValues.resize(m_diagSize);
00537 m_matrixU.resize(m_rows, m_computeFullU ? m_rows
00538 : m_computeThinU ? m_diagSize
00539 : 0);
00540 m_matrixV.resize(m_cols, m_computeFullV ? m_cols
00541 : m_computeThinV ? m_diagSize
00542 : 0);
00543 m_workMatrix.resize(m_diagSize, m_diagSize);
00544 }
00545
00546 template<typename MatrixType, int QRPreconditioner>
00547 JacobiSVD<MatrixType, QRPreconditioner>&
00548 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
00549 {
00550 allocate(matrix.rows(), matrix.cols(), computationOptions);
00551
00552
00553
00554 const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
00555
00556
00557
00558 if(!internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows>::run(*this, matrix)
00559 && !internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols>::run(*this, matrix))
00560 {
00561 m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
00562 if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
00563 if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
00564 if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
00565 if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
00566 }
00567
00568
00569
00570 bool finished = false;
00571 while(!finished)
00572 {
00573 finished = true;
00574
00575
00576
00577 for(Index p = 1; p < m_diagSize; ++p)
00578 {
00579 for(Index q = 0; q < p; ++q)
00580 {
00581
00582
00583
00584 if(std::max(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p)))
00585 > std::max(internal::abs(m_workMatrix.coeff(p,p)),internal::abs(m_workMatrix.coeff(q,q)))*precision)
00586 {
00587 finished = false;
00588
00589
00590 internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q);
00591 JacobiRotation<RealScalar> j_left, j_right;
00592 internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
00593
00594
00595 m_workMatrix.applyOnTheLeft(p,q,j_left);
00596 if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
00597
00598 m_workMatrix.applyOnTheRight(p,q,j_right);
00599 if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
00600 }
00601 }
00602 }
00603 }
00604
00605
00606
00607 for(Index i = 0; i < m_diagSize; ++i)
00608 {
00609 RealScalar a = internal::abs(m_workMatrix.coeff(i,i));
00610 m_singularValues.coeffRef(i) = a;
00611 if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
00612 }
00613
00614
00615
00616 m_nonzeroSingularValues = m_diagSize;
00617 for(Index i = 0; i < m_diagSize; i++)
00618 {
00619 Index pos;
00620 RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
00621 if(maxRemainingSingularValue == RealScalar(0))
00622 {
00623 m_nonzeroSingularValues = i;
00624 break;
00625 }
00626 if(pos)
00627 {
00628 pos += i;
00629 std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
00630 if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
00631 if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
00632 }
00633 }
00634
00635 m_isInitialized = true;
00636 return *this;
00637 }
00638
00639 namespace internal {
00640 template<typename _MatrixType, int QRPreconditioner, typename Rhs>
00641 struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
00642 : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
00643 {
00644 typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType;
00645 EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs)
00646
00647 template<typename Dest> void evalTo(Dest& dst) const
00648 {
00649 eigen_assert(rhs().rows() == dec().rows());
00650
00651
00652
00653
00654 Index diagSize = std::min(dec().rows(), dec().cols());
00655 typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize);
00656
00657 Index nonzeroSingVals = dec().nonzeroSingularValues();
00658 invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse();
00659 invertedSingVals.tail(diagSize - nonzeroSingVals).setZero();
00660
00661 dst = dec().matrixV().leftCols(diagSize)
00662 * invertedSingVals.asDiagonal()
00663 * dec().matrixU().leftCols(diagSize).adjoint()
00664 * rhs();
00665 }
00666 };
00667 }
00668
00669 template<typename Derived>
00670 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
00671 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
00672 {
00673 return JacobiSVD<PlainObject>(*this, computationOptions);
00674 }
00675
00676
00677
00678 #endif // EIGEN_JACOBISVD_H