Eigen  3.3.3
JacobiSVD.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 // Copyright (C) 2013-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_JACOBISVD_H
00012 #define EIGEN_JACOBISVD_H
00013 
00014 namespace Eigen { 
00015 
00016 namespace internal {
00017 // forward declaration (needed by ICC)
00018 // the empty body is required by MSVC
00019 template<typename MatrixType, int QRPreconditioner,
00020          bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
00021 struct svd_precondition_2x2_block_to_be_real {};
00022 
00023 /*** QR preconditioners (R-SVD)
00024  ***
00025  *** Their role is to reduce the problem of computing the SVD to the case of a square matrix.
00026  *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for
00027  *** JacobiSVD which by itself is only able to work on square matrices.
00028  ***/
00029 
00030 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
00031 
00032 template<typename MatrixType, int QRPreconditioner, int Case>
00033 struct qr_preconditioner_should_do_anything
00034 {
00035   enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
00036              MatrixType::ColsAtCompileTime != Dynamic &&
00037              MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
00038          b = MatrixType::RowsAtCompileTime != Dynamic &&
00039              MatrixType::ColsAtCompileTime != Dynamic &&
00040              MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
00041          ret = !( (QRPreconditioner == NoQRPreconditioner) ||
00042                   (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
00043                   (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
00044   };
00045 };
00046 
00047 template<typename MatrixType, int QRPreconditioner, int Case,
00048          bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
00049 > struct qr_preconditioner_impl {};
00050 
00051 template<typename MatrixType, int QRPreconditioner, int Case>
00052 class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
00053 {
00054 public:
00055   void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {}
00056   bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
00057   {
00058     return false;
00059   }
00060 };
00061 
00062 /*** preconditioner using FullPivHouseholderQR ***/
00063 
00064 template<typename MatrixType>
00065 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
00066 {
00067 public:
00068   typedef typename MatrixType::Scalar Scalar;
00069   enum
00070   {
00071     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00072     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
00073   };
00074   typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
00075 
00076   void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
00077   {
00078     if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
00079     {
00080       m_qr.~QRType();
00081       ::new (&m_qr) QRType(svd.rows(), svd.cols());
00082     }
00083     if (svd.m_computeFullU) m_workspace.resize(svd.rows());
00084   }
00085 
00086   bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00087   {
00088     if(matrix.rows() > matrix.cols())
00089     {
00090       m_qr.compute(matrix);
00091       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
00092       if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace);
00093       if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
00094       return true;
00095     }
00096     return false;
00097   }
00098 private:
00099   typedef FullPivHouseholderQR<MatrixType> QRType;
00100   QRType m_qr;
00101   WorkspaceType m_workspace;
00102 };
00103 
00104 template<typename MatrixType>
00105 class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
00106 {
00107 public:
00108   typedef typename MatrixType::Scalar Scalar;
00109   enum
00110   {
00111     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00112     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00113     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00114     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00115     TrOptions = RowsAtCompileTime==1 ? (MatrixType::Options & ~(RowMajor))
00116               : ColsAtCompileTime==1 ? (MatrixType::Options |   RowMajor)
00117               : MatrixType::Options
00118   };
00119   typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, TrOptions, MaxColsAtCompileTime, MaxRowsAtCompileTime>
00120           TransposeTypeWithSameStorageOrder;
00121 
00122   void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
00123   {
00124     if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
00125     {
00126       m_qr.~QRType();
00127       ::new (&m_qr) QRType(svd.cols(), svd.rows());
00128     }
00129     m_adjoint.resize(svd.cols(), svd.rows());
00130     if (svd.m_computeFullV) m_workspace.resize(svd.cols());
00131   }
00132 
00133   bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00134   {
00135     if(matrix.cols() > matrix.rows())
00136     {
00137       m_adjoint = matrix.adjoint();
00138       m_qr.compute(m_adjoint);
00139       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
00140       if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace);
00141       if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
00142       return true;
00143     }
00144     else return false;
00145   }
00146 private:
00147   typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
00148   QRType m_qr;
00149   TransposeTypeWithSameStorageOrder m_adjoint;
00150   typename internal::plain_row_type<MatrixType>::type m_workspace;
00151 };
00152 
00153 /*** preconditioner using ColPivHouseholderQR ***/
00154 
00155 template<typename MatrixType>
00156 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
00157 {
00158 public:
00159   void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
00160   {
00161     if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
00162     {
00163       m_qr.~QRType();
00164       ::new (&m_qr) QRType(svd.rows(), svd.cols());
00165     }
00166     if (svd.m_computeFullU) m_workspace.resize(svd.rows());
00167     else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
00168   }
00169 
00170   bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00171   {
00172     if(matrix.rows() > matrix.cols())
00173     {
00174       m_qr.compute(matrix);
00175       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
00176       if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
00177       else if(svd.m_computeThinU)
00178       {
00179         svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
00180         m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
00181       }
00182       if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation();
00183       return true;
00184     }
00185     return false;
00186   }
00187 
00188 private:
00189   typedef ColPivHouseholderQR<MatrixType> QRType;
00190   QRType m_qr;
00191   typename internal::plain_col_type<MatrixType>::type m_workspace;
00192 };
00193 
00194 template<typename MatrixType>
00195 class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
00196 {
00197 public:
00198   typedef typename MatrixType::Scalar Scalar;
00199   enum
00200   {
00201     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00202     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00203     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00204     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00205     TrOptions = RowsAtCompileTime==1 ? (MatrixType::Options & ~(RowMajor))
00206               : ColsAtCompileTime==1 ? (MatrixType::Options |   RowMajor)
00207               : MatrixType::Options
00208   };
00209 
00210   typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, TrOptions, MaxColsAtCompileTime, MaxRowsAtCompileTime>
00211           TransposeTypeWithSameStorageOrder;
00212 
00213   void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
00214   {
00215     if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
00216     {
00217       m_qr.~QRType();
00218       ::new (&m_qr) QRType(svd.cols(), svd.rows());
00219     }
00220     if (svd.m_computeFullV) m_workspace.resize(svd.cols());
00221     else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
00222     m_adjoint.resize(svd.cols(), svd.rows());
00223   }
00224 
00225   bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00226   {
00227     if(matrix.cols() > matrix.rows())
00228     {
00229       m_adjoint = matrix.adjoint();
00230       m_qr.compute(m_adjoint);
00231 
00232       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
00233       if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
00234       else if(svd.m_computeThinV)
00235       {
00236         svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
00237         m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
00238       }
00239       if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation();
00240       return true;
00241     }
00242     else return false;
00243   }
00244 
00245 private:
00246   typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
00247   QRType m_qr;
00248   TransposeTypeWithSameStorageOrder m_adjoint;
00249   typename internal::plain_row_type<MatrixType>::type m_workspace;
00250 };
00251 
00252 /*** preconditioner using HouseholderQR ***/
00253 
00254 template<typename MatrixType>
00255 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
00256 {
00257 public:
00258   void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
00259   {
00260     if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols())
00261     {
00262       m_qr.~QRType();
00263       ::new (&m_qr) QRType(svd.rows(), svd.cols());
00264     }
00265     if (svd.m_computeFullU) m_workspace.resize(svd.rows());
00266     else if (svd.m_computeThinU) m_workspace.resize(svd.cols());
00267   }
00268 
00269   bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00270   {
00271     if(matrix.rows() > matrix.cols())
00272     {
00273       m_qr.compute(matrix);
00274       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>();
00275       if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace);
00276       else if(svd.m_computeThinU)
00277       {
00278         svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols());
00279         m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace);
00280       }
00281       if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols());
00282       return true;
00283     }
00284     return false;
00285   }
00286 private:
00287   typedef HouseholderQR<MatrixType> QRType;
00288   QRType m_qr;
00289   typename internal::plain_col_type<MatrixType>::type m_workspace;
00290 };
00291 
00292 template<typename MatrixType>
00293 class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
00294 {
00295 public:
00296   typedef typename MatrixType::Scalar Scalar;
00297   enum
00298   {
00299     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00300     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00301     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00302     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00303     Options = MatrixType::Options
00304   };
00305 
00306   typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime>
00307           TransposeTypeWithSameStorageOrder;
00308 
00309   void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
00310   {
00311     if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols())
00312     {
00313       m_qr.~QRType();
00314       ::new (&m_qr) QRType(svd.cols(), svd.rows());
00315     }
00316     if (svd.m_computeFullV) m_workspace.resize(svd.cols());
00317     else if (svd.m_computeThinV) m_workspace.resize(svd.rows());
00318     m_adjoint.resize(svd.cols(), svd.rows());
00319   }
00320 
00321   bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix)
00322   {
00323     if(matrix.cols() > matrix.rows())
00324     {
00325       m_adjoint = matrix.adjoint();
00326       m_qr.compute(m_adjoint);
00327 
00328       svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint();
00329       if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace);
00330       else if(svd.m_computeThinV)
00331       {
00332         svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows());
00333         m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace);
00334       }
00335       if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows());
00336       return true;
00337     }
00338     else return false;
00339   }
00340 
00341 private:
00342   typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType;
00343   QRType m_qr;
00344   TransposeTypeWithSameStorageOrder m_adjoint;
00345   typename internal::plain_row_type<MatrixType>::type m_workspace;
00346 };
00347 
00348 /*** 2x2 SVD implementation
00349  ***
00350  *** JacobiSVD consists in performing a series of 2x2 SVD subproblems
00351  ***/
00352 
00353 template<typename MatrixType, int QRPreconditioner>
00354 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
00355 {
00356   typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
00357   typedef typename MatrixType::RealScalar RealScalar;
00358   static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; }
00359 };
00360 
00361 template<typename MatrixType, int QRPreconditioner>
00362 struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
00363 {
00364   typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
00365   typedef typename MatrixType::Scalar Scalar;
00366   typedef typename MatrixType::RealScalar RealScalar;
00367   static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry)
00368   {
00369     using std::sqrt;
00370     using std::abs;
00371     Scalar z;
00372     JacobiRotation<Scalar> rot;
00373     RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
00374 
00375     const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
00376     const RealScalar precision = NumTraits<Scalar>::epsilon();
00377 
00378     if(n==0)
00379     {
00380       // make sure first column is zero
00381       work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0);
00382 
00383       if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
00384       {
00385         // work_matrix.coeff(p,q) can be zero if work_matrix.coeff(q,p) is not zero but small enough to underflow when computing n
00386         z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
00387         work_matrix.row(p) *= z;
00388         if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
00389       }
00390       if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
00391       {
00392         z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
00393         work_matrix.row(q) *= z;
00394         if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
00395       }
00396       // otherwise the second row is already zero, so we have nothing to do.
00397     }
00398     else
00399     {
00400       rot.c() = conj(work_matrix.coeff(p,p)) / n;
00401       rot.s() = work_matrix.coeff(q,p) / n;
00402       work_matrix.applyOnTheLeft(p,q,rot);
00403       if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint());
00404       if(abs(numext::imag(work_matrix.coeff(p,q)))>considerAsZero)
00405       {
00406         z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
00407         work_matrix.col(q) *= z;
00408         if(svd.computeV()) svd.m_matrixV.col(q) *= z;
00409       }
00410       if(abs(numext::imag(work_matrix.coeff(q,q)))>considerAsZero)
00411       {
00412         z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
00413         work_matrix.row(q) *= z;
00414         if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z);
00415       }
00416     }
00417 
00418     // update largest diagonal entry
00419     maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q))));
00420     // and check whether the 2x2 block is already diagonal
00421     RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
00422     return abs(work_matrix.coeff(p,q))>threshold || abs(work_matrix.coeff(q,p)) > threshold;
00423   }
00424 };
00425 
00426 template<typename _MatrixType, int QRPreconditioner> 
00427 struct traits<JacobiSVD<_MatrixType,QRPreconditioner> >
00428 {
00429   typedef _MatrixType MatrixType;
00430 };
00431 
00432 } // end namespace internal
00433 
00487 template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
00488  : public SVDBase<JacobiSVD<_MatrixType,QRPreconditioner> >
00489 {
00490     typedef SVDBase<JacobiSVD> Base;
00491   public:
00492 
00493     typedef _MatrixType MatrixType;
00494     typedef typename MatrixType::Scalar Scalar;
00495     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00496     enum {
00497       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00498       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00499       DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime),
00500       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00501       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
00502       MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime),
00503       MatrixOptions = MatrixType::Options
00504     };
00505 
00506     typedef typename Base::MatrixUType MatrixUType;
00507     typedef typename Base::MatrixVType MatrixVType;
00508     typedef typename Base::SingularValuesType SingularValuesType;
00509     
00510     typedef typename internal::plain_row_type<MatrixType>::type RowType;
00511     typedef typename internal::plain_col_type<MatrixType>::type ColType;
00512     typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
00513                    MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
00514             WorkMatrixType;
00515 
00521     JacobiSVD()
00522     {}
00523 
00524 
00531     JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
00532     {
00533       allocate(rows, cols, computationOptions);
00534     }
00535 
00546     explicit JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
00547     {
00548       compute(matrix, computationOptions);
00549     }
00550 
00561     JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
00562 
00569     JacobiSVD& compute(const MatrixType& matrix)
00570     {
00571       return compute(matrix, m_computationOptions);
00572     }
00573 
00574     using Base::computeU;
00575     using Base::computeV;
00576     using Base::rows;
00577     using Base::cols;
00578     using Base::rank;
00579 
00580   private:
00581     void allocate(Index rows, Index cols, unsigned int computationOptions);
00582 
00583   protected:
00584     using Base::m_matrixU;
00585     using Base::m_matrixV;
00586     using Base::m_singularValues;
00587     using Base::m_isInitialized;
00588     using Base::m_isAllocated;
00589     using Base::m_usePrescribedThreshold;
00590     using Base::m_computeFullU;
00591     using Base::m_computeThinU;
00592     using Base::m_computeFullV;
00593     using Base::m_computeThinV;
00594     using Base::m_computationOptions;
00595     using Base::m_nonzeroSingularValues;
00596     using Base::m_rows;
00597     using Base::m_cols;
00598     using Base::m_diagSize;
00599     using Base::m_prescribedThreshold;
00600     WorkMatrixType m_workMatrix;
00601 
00602     template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
00603     friend struct internal::svd_precondition_2x2_block_to_be_real;
00604     template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything>
00605     friend struct internal::qr_preconditioner_impl;
00606 
00607     internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
00608     internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
00609     MatrixType m_scaledMatrix;
00610 };
00611 
00612 template<typename MatrixType, int QRPreconditioner>
00613 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
00614 {
00615   eigen_assert(rows >= 0 && cols >= 0);
00616 
00617   if (m_isAllocated &&
00618       rows == m_rows &&
00619       cols == m_cols &&
00620       computationOptions == m_computationOptions)
00621   {
00622     return;
00623   }
00624 
00625   m_rows = rows;
00626   m_cols = cols;
00627   m_isInitialized = false;
00628   m_isAllocated = true;
00629   m_computationOptions = computationOptions;
00630   m_computeFullU = (computationOptions & ComputeFullU) != 0;
00631   m_computeThinU = (computationOptions & ComputeThinU) != 0;
00632   m_computeFullV = (computationOptions & ComputeFullV) != 0;
00633   m_computeThinV = (computationOptions & ComputeThinV) != 0;
00634   eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
00635   eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
00636   eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
00637               "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
00638   if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
00639   {
00640       eigen_assert(!(m_computeThinU || m_computeThinV) &&
00641               "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
00642               "Use the ColPivHouseholderQR preconditioner instead.");
00643   }
00644   m_diagSize = (std::min)(m_rows, m_cols);
00645   m_singularValues.resize(m_diagSize);
00646   if(RowsAtCompileTime==Dynamic)
00647     m_matrixU.resize(m_rows, m_computeFullU ? m_rows
00648                             : m_computeThinU ? m_diagSize
00649                             : 0);
00650   if(ColsAtCompileTime==Dynamic)
00651     m_matrixV.resize(m_cols, m_computeFullV ? m_cols
00652                             : m_computeThinV ? m_diagSize
00653                             : 0);
00654   m_workMatrix.resize(m_diagSize, m_diagSize);
00655   
00656   if(m_cols>m_rows)   m_qr_precond_morecols.allocate(*this);
00657   if(m_rows>m_cols)   m_qr_precond_morerows.allocate(*this);
00658   if(m_rows!=m_cols)  m_scaledMatrix.resize(rows,cols);
00659 }
00660 
00661 template<typename MatrixType, int QRPreconditioner>
00662 JacobiSVD<MatrixType, QRPreconditioner>&
00663 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
00664 {
00665   using std::abs;
00666   allocate(matrix.rows(), matrix.cols(), computationOptions);
00667 
00668   // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
00669   // only worsening the precision of U and V as we accumulate more rotations
00670   const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();
00671 
00672   // limit for denormal numbers to be considered zero in order to avoid infinite loops (see bug 286)
00673   const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
00674 
00675   // Scaling factor to reduce over/under-flows
00676   RealScalar scale = matrix.cwiseAbs().maxCoeff();
00677   if(scale==RealScalar(0)) scale = RealScalar(1);
00678   
00679   /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
00680 
00681   if(m_rows!=m_cols)
00682   {
00683     m_scaledMatrix = matrix / scale;
00684     m_qr_precond_morecols.run(*this, m_scaledMatrix);
00685     m_qr_precond_morerows.run(*this, m_scaledMatrix);
00686   }
00687   else
00688   {
00689     m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
00690     if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
00691     if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
00692     if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
00693     if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
00694   }
00695 
00696   /*** step 2. The main Jacobi SVD iteration. ***/
00697   RealScalar maxDiagEntry = m_workMatrix.cwiseAbs().diagonal().maxCoeff();
00698 
00699   bool finished = false;
00700   while(!finished)
00701   {
00702     finished = true;
00703 
00704     // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
00705 
00706     for(Index p = 1; p < m_diagSize; ++p)
00707     {
00708       for(Index q = 0; q < p; ++q)
00709       {
00710         // if this 2x2 sub-matrix is not diagonal already...
00711         // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
00712         // keep us iterating forever. Similarly, small denormal numbers are considered zero.
00713         RealScalar threshold = numext::maxi<RealScalar>(considerAsZero, precision * maxDiagEntry);
00714         if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
00715         {
00716           finished = false;
00717           // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
00718           // the complex to real operation returns true if the updated 2x2 block is not already diagonal
00719           if(internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q, maxDiagEntry))
00720           {
00721             JacobiRotation<RealScalar> j_left, j_right;
00722             internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
00723 
00724             // accumulate resulting Jacobi rotations
00725             m_workMatrix.applyOnTheLeft(p,q,j_left);
00726             if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
00727 
00728             m_workMatrix.applyOnTheRight(p,q,j_right);
00729             if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right);
00730 
00731             // keep track of the largest diagonal coefficient
00732             maxDiagEntry = numext::maxi<RealScalar>(maxDiagEntry,numext::maxi<RealScalar>(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q))));
00733           }
00734         }
00735       }
00736     }
00737   }
00738 
00739   /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
00740 
00741   for(Index i = 0; i < m_diagSize; ++i)
00742   {
00743     // For a complex matrix, some diagonal coefficients might note have been
00744     // treated by svd_precondition_2x2_block_to_be_real, and the imaginary part
00745     // of some diagonal entry might not be null.
00746     if(NumTraits<Scalar>::IsComplex && abs(numext::imag(m_workMatrix.coeff(i,i)))>considerAsZero)
00747     {
00748       RealScalar a = abs(m_workMatrix.coeff(i,i));
00749       m_singularValues.coeffRef(i) = abs(a);
00750       if(computeU()) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
00751     }
00752     else
00753     {
00754       // m_workMatrix.coeff(i,i) is already real, no difficulty:
00755       RealScalar a = numext::real(m_workMatrix.coeff(i,i));
00756       m_singularValues.coeffRef(i) = abs(a);
00757       if(computeU() && (a<RealScalar(0))) m_matrixU.col(i) = -m_matrixU.col(i);
00758     }
00759   }
00760   
00761   m_singularValues *= scale;
00762 
00763   /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
00764 
00765   m_nonzeroSingularValues = m_diagSize;
00766   for(Index i = 0; i < m_diagSize; i++)
00767   {
00768     Index pos;
00769     RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
00770     if(maxRemainingSingularValue == RealScalar(0))
00771     {
00772       m_nonzeroSingularValues = i;
00773       break;
00774     }
00775     if(pos)
00776     {
00777       pos += i;
00778       std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
00779       if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i));
00780       if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
00781     }
00782   }
00783 
00784   m_isInitialized = true;
00785   return *this;
00786 }
00787 
00795 template<typename Derived>
00796 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
00797 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
00798 {
00799   return JacobiSVD<PlainObject>(*this, computationOptions);
00800 }
00801 
00802 } // end namespace Eigen
00803 
00804 #endif // EIGEN_JACOBISVD_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends