Eigen  3.3.3
BDCSVD.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 // 
00004 // We used the "A Divide-And-Conquer Algorithm for the Bidiagonal SVD"
00005 // research report written by Ming Gu and Stanley C.Eisenstat
00006 // The code variable names correspond to the names they used in their 
00007 // report
00008 //
00009 // Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com>
00010 // Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr>
00011 // Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr>
00012 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr>
00013 // Copyright (C) 2013 Jitse Niesen <jitse@maths.leeds.ac.uk>
00014 // Copyright (C) 2014-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
00015 //
00016 // Source Code Form is subject to the terms of the Mozilla
00017 // Public License v. 2.0. If a copy of the MPL was not distributed
00018 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00019 
00020 #ifndef EIGEN_BDCSVD_H
00021 #define EIGEN_BDCSVD_H
00022 // #define EIGEN_BDCSVD_DEBUG_VERBOSE
00023 // #define EIGEN_BDCSVD_SANITY_CHECKS
00024 
00025 namespace Eigen {
00026 
00027 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00028 IOFormat bdcsvdfmt(8, 0, ", ", "\n", "  [", "]");
00029 #endif
00030   
00031 template<typename _MatrixType> class BDCSVD;
00032 
00033 namespace internal {
00034 
00035 template<typename _MatrixType> 
00036 struct traits<BDCSVD<_MatrixType> >
00037 {
00038   typedef _MatrixType MatrixType;
00039 };  
00040 
00041 } // end namespace internal
00042   
00043   
00066 template<typename _MatrixType> 
00067 class BDCSVD : public SVDBase<BDCSVD<_MatrixType> >
00068 {
00069   typedef SVDBase<BDCSVD> Base;
00070     
00071 public:
00072   using Base::rows;
00073   using Base::cols;
00074   using Base::computeU;
00075   using Base::computeV;
00076   
00077   typedef _MatrixType MatrixType;
00078   typedef typename MatrixType::Scalar Scalar;
00079   typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00080   enum {
00081     RowsAtCompileTime = MatrixType::RowsAtCompileTime, 
00082     ColsAtCompileTime = MatrixType::ColsAtCompileTime, 
00083     DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime), 
00084     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 
00085     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 
00086     MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime, MaxColsAtCompileTime), 
00087     MatrixOptions = MatrixType::Options
00088   };
00089 
00090   typedef typename Base::MatrixUType MatrixUType;
00091   typedef typename Base::MatrixVType MatrixVType;
00092   typedef typename Base::SingularValuesType SingularValuesType;
00093   
00094   typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> MatrixX;
00095   typedef Matrix<RealScalar, Dynamic, Dynamic, ColMajor> MatrixXr;
00096   typedef Matrix<RealScalar, Dynamic, 1> VectorType;
00097   typedef Array<RealScalar, Dynamic, 1> ArrayXr;
00098   typedef Array<Index,1,Dynamic> ArrayXi;
00099   typedef Ref<ArrayXr> ArrayRef;
00100   typedef Ref<ArrayXi> IndicesRef;
00101 
00107   BDCSVD() : m_algoswap(16), m_numIters(0)
00108   {}
00109 
00110 
00117   BDCSVD(Index rows, Index cols, unsigned int computationOptions = 0)
00118     : m_algoswap(16), m_numIters(0)
00119   {
00120     allocate(rows, cols, computationOptions);
00121   }
00122 
00133   BDCSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
00134     : m_algoswap(16), m_numIters(0)
00135   {
00136     compute(matrix, computationOptions);
00137   }
00138 
00139   ~BDCSVD() 
00140   {
00141   }
00142   
00153   BDCSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
00154 
00161   BDCSVD& compute(const MatrixType& matrix)
00162   {
00163     return compute(matrix, this->m_computationOptions);
00164   }
00165 
00166   void setSwitchSize(int s) 
00167   {
00168     eigen_assert(s>3 && "BDCSVD the size of the algo switch has to be greater than 3");
00169     m_algoswap = s;
00170   }
00171  
00172 private:
00173   void allocate(Index rows, Index cols, unsigned int computationOptions);
00174   void divide(Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift);
00175   void computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V);
00176   void computeSingVals(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef& perm, VectorType& singVals, ArrayRef shifts, ArrayRef mus);
00177   void perturbCol0(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef& perm, const VectorType& singVals, const ArrayRef& shifts, const ArrayRef& mus, ArrayRef zhat);
00178   void computeSingVecs(const ArrayRef& zhat, const ArrayRef& diag, const IndicesRef& perm, const VectorType& singVals, const ArrayRef& shifts, const ArrayRef& mus, MatrixXr& U, MatrixXr& V);
00179   void deflation43(Index firstCol, Index shift, Index i, Index size);
00180   void deflation44(Index firstColu , Index firstColm, Index firstRowW, Index firstColW, Index i, Index j, Index size);
00181   void deflation(Index firstCol, Index lastCol, Index k, Index firstRowW, Index firstColW, Index shift);
00182   template<typename HouseholderU, typename HouseholderV, typename NaiveU, typename NaiveV>
00183   void copyUV(const HouseholderU &householderU, const HouseholderV &householderV, const NaiveU &naiveU, const NaiveV &naivev);
00184   void structured_update(Block<MatrixXr,Dynamic,Dynamic> A, const MatrixXr &B, Index n1);
00185   static RealScalar secularEq(RealScalar x, const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const ArrayRef& diagShifted, RealScalar shift);
00186 
00187 protected:
00188   MatrixXr m_naiveU, m_naiveV;
00189   MatrixXr m_computed;
00190   Index m_nRec;
00191   ArrayXr m_workspace;
00192   ArrayXi m_workspaceI;
00193   int m_algoswap;
00194   bool m_isTranspose, m_compU, m_compV;
00195   
00196   using Base::m_singularValues;
00197   using Base::m_diagSize;
00198   using Base::m_computeFullU;
00199   using Base::m_computeFullV;
00200   using Base::m_computeThinU;
00201   using Base::m_computeThinV;
00202   using Base::m_matrixU;
00203   using Base::m_matrixV;
00204   using Base::m_isInitialized;
00205   using Base::m_nonzeroSingularValues;
00206 
00207 public:  
00208   int m_numIters;
00209 }; //end class BDCSVD
00210 
00211 
00212 // Method to allocate and initialize matrix and attributes
00213 template<typename MatrixType>
00214 void BDCSVD<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions)
00215 {
00216   m_isTranspose = (cols > rows);
00217 
00218   if (Base::allocate(rows, cols, computationOptions))
00219     return;
00220   
00221   m_computed = MatrixXr::Zero(m_diagSize + 1, m_diagSize );
00222   m_compU = computeV();
00223   m_compV = computeU();
00224   if (m_isTranspose)
00225     std::swap(m_compU, m_compV);
00226   
00227   if (m_compU) m_naiveU = MatrixXr::Zero(m_diagSize + 1, m_diagSize + 1 );
00228   else         m_naiveU = MatrixXr::Zero(2, m_diagSize + 1 );
00229   
00230   if (m_compV) m_naiveV = MatrixXr::Zero(m_diagSize, m_diagSize);
00231   
00232   m_workspace.resize((m_diagSize+1)*(m_diagSize+1)*3);
00233   m_workspaceI.resize(3*m_diagSize);
00234 }// end allocate
00235 
00236 template<typename MatrixType>
00237 BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOptions) 
00238 {
00239 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00240   std::cout << "\n\n\n======================================================================================================================\n\n\n";
00241 #endif
00242   allocate(matrix.rows(), matrix.cols(), computationOptions);
00243   using std::abs;
00244 
00245   const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
00246   
00247   //**** step -1 - If the problem is too small, directly falls back to JacobiSVD and return
00248   if(matrix.cols() < m_algoswap)
00249   {
00250     // FIXME this line involves temporaries
00251     JacobiSVD<MatrixType> jsvd(matrix,computationOptions);
00252     if(computeU()) m_matrixU = jsvd.matrixU();
00253     if(computeV()) m_matrixV = jsvd.matrixV();
00254     m_singularValues = jsvd.singularValues();
00255     m_nonzeroSingularValues = jsvd.nonzeroSingularValues();
00256     m_isInitialized = true;
00257     return *this;
00258   }
00259   
00260   //**** step 0 - Copy the input matrix and apply scaling to reduce over/under-flows
00261   RealScalar scale = matrix.cwiseAbs().maxCoeff();
00262   if(scale==RealScalar(0)) scale = RealScalar(1);
00263   MatrixX copy;
00264   if (m_isTranspose) copy = matrix.adjoint()/scale;
00265   else               copy = matrix/scale;
00266   
00267   //**** step 1 - Bidiagonalization
00268   // FIXME this line involves temporaries
00269   internal::UpperBidiagonalization<MatrixX> bid(copy);
00270 
00271   //**** step 2 - Divide & Conquer
00272   m_naiveU.setZero();
00273   m_naiveV.setZero();
00274   // FIXME this line involves a temporary matrix
00275   m_computed.topRows(m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose();
00276   m_computed.template bottomRows<1>().setZero();
00277   divide(0, m_diagSize - 1, 0, 0, 0);
00278 
00279   //**** step 3 - Copy singular values and vectors
00280   for (int i=0; i<m_diagSize; i++)
00281   {
00282     RealScalar a = abs(m_computed.coeff(i, i));
00283     m_singularValues.coeffRef(i) = a * scale;
00284     if (a<considerZero)
00285     {
00286       m_nonzeroSingularValues = i;
00287       m_singularValues.tail(m_diagSize - i - 1).setZero();
00288       break;
00289     }
00290     else if (i == m_diagSize - 1)
00291     {
00292       m_nonzeroSingularValues = i + 1;
00293       break;
00294     }
00295   }
00296 
00297 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00298 //   std::cout << "m_naiveU\n" << m_naiveU << "\n\n";
00299 //   std::cout << "m_naiveV\n" << m_naiveV << "\n\n";
00300 #endif
00301   if(m_isTranspose) copyUV(bid.householderV(), bid.householderU(), m_naiveV, m_naiveU);
00302   else              copyUV(bid.householderU(), bid.householderV(), m_naiveU, m_naiveV);
00303 
00304   m_isInitialized = true;
00305   return *this;
00306 }// end compute
00307 
00308 
00309 template<typename MatrixType>
00310 template<typename HouseholderU, typename HouseholderV, typename NaiveU, typename NaiveV>
00311 void BDCSVD<MatrixType>::copyUV(const HouseholderU &householderU, const HouseholderV &householderV, const NaiveU &naiveU, const NaiveV &naiveV)
00312 {
00313   // Note exchange of U and V: m_matrixU is set from m_naiveV and vice versa
00314   if (computeU())
00315   {
00316     Index Ucols = m_computeThinU ? m_diagSize : householderU.cols();
00317     m_matrixU = MatrixX::Identity(householderU.cols(), Ucols);
00318     m_matrixU.topLeftCorner(m_diagSize, m_diagSize) = naiveV.template cast<Scalar>().topLeftCorner(m_diagSize, m_diagSize);
00319     householderU.applyThisOnTheLeft(m_matrixU); // FIXME this line involves a temporary buffer
00320   }
00321   if (computeV())
00322   {
00323     Index Vcols = m_computeThinV ? m_diagSize : householderV.cols();
00324     m_matrixV = MatrixX::Identity(householderV.cols(), Vcols);
00325     m_matrixV.topLeftCorner(m_diagSize, m_diagSize) = naiveU.template cast<Scalar>().topLeftCorner(m_diagSize, m_diagSize);
00326     householderV.applyThisOnTheLeft(m_matrixV); // FIXME this line involves a temporary buffer
00327   }
00328 }
00329 
00338 template<typename MatrixType>
00339 void BDCSVD<MatrixType>::structured_update(Block<MatrixXr,Dynamic,Dynamic> A, const MatrixXr &B, Index n1)
00340 {
00341   Index n = A.rows();
00342   if(n>100)
00343   {
00344     // If the matrices are large enough, let's exploit the sparse structure of A by
00345     // splitting it in half (wrt n1), and packing the non-zero columns.
00346     Index n2 = n - n1;
00347     Map<MatrixXr> A1(m_workspace.data()      , n1, n);
00348     Map<MatrixXr> A2(m_workspace.data()+ n1*n, n2, n);
00349     Map<MatrixXr> B1(m_workspace.data()+  n*n, n,  n);
00350     Map<MatrixXr> B2(m_workspace.data()+2*n*n, n,  n);
00351     Index k1=0, k2=0;
00352     for(Index j=0; j<n; ++j)
00353     {
00354       if( (A.col(j).head(n1).array()!=0).any() )
00355       {
00356         A1.col(k1) = A.col(j).head(n1);
00357         B1.row(k1) = B.row(j);
00358         ++k1;
00359       }
00360       if( (A.col(j).tail(n2).array()!=0).any() )
00361       {
00362         A2.col(k2) = A.col(j).tail(n2);
00363         B2.row(k2) = B.row(j);
00364         ++k2;
00365       }
00366     }
00367   
00368     A.topRows(n1).noalias()    = A1.leftCols(k1) * B1.topRows(k1);
00369     A.bottomRows(n2).noalias() = A2.leftCols(k2) * B2.topRows(k2);
00370   }
00371   else
00372   {
00373     Map<MatrixXr,Aligned> tmp(m_workspace.data(),n,n);
00374     tmp.noalias() = A*B;
00375     A = tmp;
00376   }
00377 }
00378 
00379 // The divide algorithm is done "in place", we are always working on subsets of the same matrix. The divide methods takes as argument the 
00380 // place of the submatrix we are currently working on.
00381 
00382 //@param firstCol : The Index of the first column of the submatrix of m_computed and for m_naiveU;
00383 //@param lastCol : The Index of the last column of the submatrix of m_computed and for m_naiveU; 
00384 // lastCol + 1 - firstCol is the size of the submatrix.
00385 //@param firstRowW : The Index of the first row of the matrix W that we are to change. (see the reference paper section 1 for more information on W)
00386 //@param firstRowW : Same as firstRowW with the column.
00387 //@param shift : Each time one takes the left submatrix, one must add 1 to the shift. Why? Because! We actually want the last column of the U submatrix 
00388 // to become the first column (*coeff) and to shift all the other columns to the right. There are more details on the reference paper.
00389 template<typename MatrixType>
00390 void BDCSVD<MatrixType>::divide (Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift)
00391 {
00392   // requires rows = cols + 1;
00393   using std::pow;
00394   using std::sqrt;
00395   using std::abs;
00396   const Index n = lastCol - firstCol + 1;
00397   const Index k = n/2;
00398   const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
00399   RealScalar alphaK;
00400   RealScalar betaK; 
00401   RealScalar r0; 
00402   RealScalar lambda, phi, c0, s0;
00403   VectorType l, f;
00404   // We use the other algorithm which is more efficient for small 
00405   // matrices.
00406   if (n < m_algoswap)
00407   {
00408     // FIXME this line involves temporaries
00409     JacobiSVD<MatrixXr> b(m_computed.block(firstCol, firstCol, n + 1, n), ComputeFullU | (m_compV ? ComputeFullV : 0));
00410     if (m_compU)
00411       m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() = b.matrixU();
00412     else 
00413     {
00414       m_naiveU.row(0).segment(firstCol, n + 1).real() = b.matrixU().row(0);
00415       m_naiveU.row(1).segment(firstCol, n + 1).real() = b.matrixU().row(n);
00416     }
00417     if (m_compV) m_naiveV.block(firstRowW, firstColW, n, n).real() = b.matrixV();
00418     m_computed.block(firstCol + shift, firstCol + shift, n + 1, n).setZero();
00419     m_computed.diagonal().segment(firstCol + shift, n) = b.singularValues().head(n);
00420     return;
00421   }
00422   // We use the divide and conquer algorithm
00423   alphaK =  m_computed(firstCol + k, firstCol + k);
00424   betaK = m_computed(firstCol + k + 1, firstCol + k);
00425   // The divide must be done in that order in order to have good results. Divide change the data inside the submatrices
00426   // and the divide of the right submatrice reads one column of the left submatrice. That's why we need to treat the 
00427   // right submatrix before the left one. 
00428   divide(k + 1 + firstCol, lastCol, k + 1 + firstRowW, k + 1 + firstColW, shift);
00429   divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1);
00430 
00431   if (m_compU)
00432   {
00433     lambda = m_naiveU(firstCol + k, firstCol + k);
00434     phi = m_naiveU(firstCol + k + 1, lastCol + 1);
00435   } 
00436   else 
00437   {
00438     lambda = m_naiveU(1, firstCol + k);
00439     phi = m_naiveU(0, lastCol + 1);
00440   }
00441   r0 = sqrt((abs(alphaK * lambda) * abs(alphaK * lambda)) + abs(betaK * phi) * abs(betaK * phi));
00442   if (m_compU)
00443   {
00444     l = m_naiveU.row(firstCol + k).segment(firstCol, k);
00445     f = m_naiveU.row(firstCol + k + 1).segment(firstCol + k + 1, n - k - 1);
00446   } 
00447   else 
00448   {
00449     l = m_naiveU.row(1).segment(firstCol, k);
00450     f = m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1);
00451   }
00452   if (m_compV) m_naiveV(firstRowW+k, firstColW) = 1;
00453   if (r0<considerZero)
00454   {
00455     c0 = 1;
00456     s0 = 0;
00457   }
00458   else
00459   {
00460     c0 = alphaK * lambda / r0;
00461     s0 = betaK * phi / r0;
00462   }
00463   
00464 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
00465   assert(m_naiveU.allFinite());
00466   assert(m_naiveV.allFinite());
00467   assert(m_computed.allFinite());
00468 #endif
00469   
00470   if (m_compU)
00471   {
00472     MatrixXr q1 (m_naiveU.col(firstCol + k).segment(firstCol, k + 1));     
00473     // we shiftW Q1 to the right
00474     for (Index i = firstCol + k - 1; i >= firstCol; i--) 
00475       m_naiveU.col(i + 1).segment(firstCol, k + 1) = m_naiveU.col(i).segment(firstCol, k + 1);
00476     // we shift q1 at the left with a factor c0
00477     m_naiveU.col(firstCol).segment( firstCol, k + 1) = (q1 * c0);
00478     // last column = q1 * - s0
00479     m_naiveU.col(lastCol + 1).segment(firstCol, k + 1) = (q1 * ( - s0));
00480     // first column = q2 * s0
00481     m_naiveU.col(firstCol).segment(firstCol + k + 1, n - k) = m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) * s0; 
00482     // q2 *= c0
00483     m_naiveU.col(lastCol + 1).segment(firstCol + k + 1, n - k) *= c0;
00484   } 
00485   else 
00486   {
00487     RealScalar q1 = m_naiveU(0, firstCol + k);
00488     // we shift Q1 to the right
00489     for (Index i = firstCol + k - 1; i >= firstCol; i--) 
00490       m_naiveU(0, i + 1) = m_naiveU(0, i);
00491     // we shift q1 at the left with a factor c0
00492     m_naiveU(0, firstCol) = (q1 * c0);
00493     // last column = q1 * - s0
00494     m_naiveU(0, lastCol + 1) = (q1 * ( - s0));
00495     // first column = q2 * s0
00496     m_naiveU(1, firstCol) = m_naiveU(1, lastCol + 1) *s0; 
00497     // q2 *= c0
00498     m_naiveU(1, lastCol + 1) *= c0;
00499     m_naiveU.row(1).segment(firstCol + 1, k).setZero();
00500     m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1).setZero();
00501   }
00502   
00503 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
00504   assert(m_naiveU.allFinite());
00505   assert(m_naiveV.allFinite());
00506   assert(m_computed.allFinite());
00507 #endif
00508   
00509   m_computed(firstCol + shift, firstCol + shift) = r0;
00510   m_computed.col(firstCol + shift).segment(firstCol + shift + 1, k) = alphaK * l.transpose().real();
00511   m_computed.col(firstCol + shift).segment(firstCol + shift + k + 1, n - k - 1) = betaK * f.transpose().real();
00512 
00513 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00514   ArrayXr tmp1 = (m_computed.block(firstCol+shift, firstCol+shift, n, n)).jacobiSvd().singularValues();
00515 #endif
00516   // Second part: try to deflate singular values in combined matrix
00517   deflation(firstCol, lastCol, k, firstRowW, firstColW, shift);
00518 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00519   ArrayXr tmp2 = (m_computed.block(firstCol+shift, firstCol+shift, n, n)).jacobiSvd().singularValues();
00520   std::cout << "\n\nj1 = " << tmp1.transpose().format(bdcsvdfmt) << "\n";
00521   std::cout << "j2 = " << tmp2.transpose().format(bdcsvdfmt) << "\n\n";
00522   std::cout << "err:      " << ((tmp1-tmp2).abs()>1e-12*tmp2.abs()).transpose() << "\n";
00523   static int count = 0;
00524   std::cout << "# " << ++count << "\n\n";
00525   assert((tmp1-tmp2).matrix().norm() < 1e-14*tmp2.matrix().norm());
00526 //   assert(count<681);
00527 //   assert(((tmp1-tmp2).abs()<1e-13*tmp2.abs()).all());
00528 #endif
00529   
00530   // Third part: compute SVD of combined matrix
00531   MatrixXr UofSVD, VofSVD;
00532   VectorType singVals;
00533   computeSVDofM(firstCol + shift, n, UofSVD, singVals, VofSVD);
00534   
00535 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
00536   assert(UofSVD.allFinite());
00537   assert(VofSVD.allFinite());
00538 #endif
00539   
00540   if (m_compU)
00541     structured_update(m_naiveU.block(firstCol, firstCol, n + 1, n + 1), UofSVD, (n+2)/2);
00542   else
00543   {
00544     Map<Matrix<RealScalar,2,Dynamic>,Aligned> tmp(m_workspace.data(),2,n+1);
00545     tmp.noalias() = m_naiveU.middleCols(firstCol, n+1) * UofSVD;
00546     m_naiveU.middleCols(firstCol, n + 1) = tmp;
00547   }
00548   
00549   if (m_compV)  structured_update(m_naiveV.block(firstRowW, firstColW, n, n), VofSVD, (n+1)/2);
00550   
00551 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
00552   assert(m_naiveU.allFinite());
00553   assert(m_naiveV.allFinite());
00554   assert(m_computed.allFinite());
00555 #endif
00556   
00557   m_computed.block(firstCol + shift, firstCol + shift, n, n).setZero();
00558   m_computed.block(firstCol + shift, firstCol + shift, n, n).diagonal() = singVals;
00559 }// end divide
00560 
00561 // Compute SVD of m_computed.block(firstCol, firstCol, n + 1, n); this block only has non-zeros in
00562 // the first column and on the diagonal and has undergone deflation, so diagonal is in increasing
00563 // order except for possibly the (0,0) entry. The computed SVD is stored U, singVals and V, except
00564 // that if m_compV is false, then V is not computed. Singular values are sorted in decreasing order.
00565 //
00566 // TODO Opportunities for optimization: better root finding algo, better stopping criterion, better
00567 // handling of round-off errors, be consistent in ordering
00568 // For instance, to solve the secular equation using FMM, see http://www.stat.uchicago.edu/~lekheng/courses/302/classics/greengard-rokhlin.pdf
00569 template <typename MatrixType>
00570 void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V)
00571 {
00572   const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
00573   using std::abs;
00574   ArrayRef col0 = m_computed.col(firstCol).segment(firstCol, n);
00575   m_workspace.head(n) =  m_computed.block(firstCol, firstCol, n, n).diagonal();
00576   ArrayRef diag = m_workspace.head(n);
00577   diag(0) = 0;
00578 
00579   // Allocate space for singular values and vectors
00580   singVals.resize(n);
00581   U.resize(n+1, n+1);
00582   if (m_compV) V.resize(n, n);
00583 
00584 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00585   if (col0.hasNaN() || diag.hasNaN())
00586     std::cout << "\n\nHAS NAN\n\n";
00587 #endif
00588   
00589   // Many singular values might have been deflated, the zero ones have been moved to the end,
00590   // but others are interleaved and we must ignore them at this stage.
00591   // To this end, let's compute a permutation skipping them:
00592   Index actual_n = n;
00593   while(actual_n>1 && diag(actual_n-1)==0) --actual_n;
00594   Index m = 0; // size of the deflated problem
00595   for(Index k=0;k<actual_n;++k)
00596     if(abs(col0(k))>considerZero)
00597       m_workspaceI(m++) = k;
00598   Map<ArrayXi> perm(m_workspaceI.data(),m);
00599   
00600   Map<ArrayXr> shifts(m_workspace.data()+1*n, n);
00601   Map<ArrayXr> mus(m_workspace.data()+2*n, n);
00602   Map<ArrayXr> zhat(m_workspace.data()+3*n, n);
00603 
00604 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00605   std::cout << "computeSVDofM using:\n";
00606   std::cout << "  z: " << col0.transpose() << "\n";
00607   std::cout << "  d: " << diag.transpose() << "\n";
00608 #endif
00609   
00610   // Compute singVals, shifts, and mus
00611   computeSingVals(col0, diag, perm, singVals, shifts, mus);
00612   
00613 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00614   std::cout << "  j:        " << (m_computed.block(firstCol, firstCol, n, n)).jacobiSvd().singularValues().transpose().reverse() << "\n\n";
00615   std::cout << "  sing-val: " << singVals.transpose() << "\n";
00616   std::cout << "  mu:       " << mus.transpose() << "\n";
00617   std::cout << "  shift:    " << shifts.transpose() << "\n";
00618   
00619   {
00620     Index actual_n = n;
00621     while(actual_n>1 && abs(col0(actual_n-1))<considerZero) --actual_n;
00622     std::cout << "\n\n    mus:    " << mus.head(actual_n).transpose() << "\n\n";
00623     std::cout << "    check1 (expect0) : " << ((singVals.array()-(shifts+mus)) / singVals.array()).head(actual_n).transpose() << "\n\n";
00624     std::cout << "    check2 (>0)      : " << ((singVals.array()-diag) / singVals.array()).head(actual_n).transpose() << "\n\n";
00625     std::cout << "    check3 (>0)      : " << ((diag.segment(1,actual_n-1)-singVals.head(actual_n-1).array()) / singVals.head(actual_n-1).array()).transpose() << "\n\n\n";
00626     std::cout << "    check4 (>0)      : " << ((singVals.segment(1,actual_n-1)-singVals.head(actual_n-1))).transpose() << "\n\n\n";
00627   }
00628 #endif
00629   
00630 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
00631   assert(singVals.allFinite());
00632   assert(mus.allFinite());
00633   assert(shifts.allFinite());
00634 #endif
00635   
00636   // Compute zhat
00637   perturbCol0(col0, diag, perm, singVals, shifts, mus, zhat);
00638 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
00639   std::cout << "  zhat: " << zhat.transpose() << "\n";
00640 #endif
00641   
00642 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
00643   assert(zhat.allFinite());
00644 #endif
00645   
00646   computeSingVecs(zhat, diag, perm, singVals, shifts, mus, U, V);
00647   
00648 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
00649   std::cout << "U^T U: " << (U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() << "\n";
00650   std::cout << "V^T V: " << (V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() << "\n";
00651 #endif
00652   
00653 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
00654   assert(U.allFinite());
00655   assert(V.allFinite());
00656   assert((U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() < 1e-14 * n);
00657   assert((V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() < 1e-14 * n);
00658   assert(m_naiveU.allFinite());
00659   assert(m_naiveV.allFinite());
00660   assert(m_computed.allFinite());
00661 #endif
00662   
00663   // Because of deflation, the singular values might not be completely sorted.
00664   // Fortunately, reordering them is a O(n) problem
00665   for(Index i=0; i<actual_n-1; ++i)
00666   {
00667     if(singVals(i)>singVals(i+1))
00668     {
00669       using std::swap;
00670       swap(singVals(i),singVals(i+1));
00671       U.col(i).swap(U.col(i+1));
00672       if(m_compV) V.col(i).swap(V.col(i+1));
00673     }
00674   }
00675   
00676   // Reverse order so that singular values in increased order
00677   // Because of deflation, the zeros singular-values are already at the end
00678   singVals.head(actual_n).reverseInPlace();
00679   U.leftCols(actual_n).rowwise().reverseInPlace();
00680   if (m_compV) V.leftCols(actual_n).rowwise().reverseInPlace();
00681   
00682 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00683   JacobiSVD<MatrixXr> jsvd(m_computed.block(firstCol, firstCol, n, n) );
00684   std::cout << "  * j:        " << jsvd.singularValues().transpose() << "\n\n";
00685   std::cout << "  * sing-val: " << singVals.transpose() << "\n";
00686 //   std::cout << "  * err:      " << ((jsvd.singularValues()-singVals)>1e-13*singVals.norm()).transpose() << "\n";
00687 #endif
00688 }
00689 
00690 template <typename MatrixType>
00691 typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::secularEq(RealScalar mu, const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const ArrayRef& diagShifted, RealScalar shift)
00692 {
00693   Index m = perm.size();
00694   RealScalar res = 1;
00695   for(Index i=0; i<m; ++i)
00696   {
00697     Index j = perm(i);
00698     res += numext::abs2(col0(j)) / ((diagShifted(j) - mu) * (diag(j) + shift + mu));
00699   }
00700   return res;
00701 
00702 }
00703 
00704 template <typename MatrixType>
00705 void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm,
00706                                          VectorType& singVals, ArrayRef shifts, ArrayRef mus)
00707 {
00708   using std::abs;
00709   using std::swap;
00710 
00711   Index n = col0.size();
00712   Index actual_n = n;
00713   while(actual_n>1 && col0(actual_n-1)==0) --actual_n;
00714 
00715   for (Index k = 0; k < n; ++k)
00716   {
00717     if (col0(k) == 0 || actual_n==1)
00718     {
00719       // if col0(k) == 0, then entry is deflated, so singular value is on diagonal
00720       // if actual_n==1, then the deflated problem is already diagonalized
00721       singVals(k) = k==0 ? col0(0) : diag(k);
00722       mus(k) = 0;
00723       shifts(k) = k==0 ? col0(0) : diag(k);
00724       continue;
00725     } 
00726 
00727     // otherwise, use secular equation to find singular value
00728     RealScalar left = diag(k);
00729     RealScalar right; // was: = (k != actual_n-1) ? diag(k+1) : (diag(actual_n-1) + col0.matrix().norm());
00730     if(k==actual_n-1)
00731       right = (diag(actual_n-1) + col0.matrix().norm());
00732     else
00733     {
00734       // Skip deflated singular values
00735       Index l = k+1;
00736       while(col0(l)==0) { ++l; eigen_internal_assert(l<actual_n); }
00737       right = diag(l);
00738     }
00739 
00740     // first decide whether it's closer to the left end or the right end
00741     RealScalar mid = left + (right-left) / 2;
00742     RealScalar fMid = secularEq(mid, col0, diag, perm, diag, 0);
00743 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00744     std::cout << right-left << "\n";
00745     std::cout << "fMid = " << fMid << " " << secularEq(mid-left, col0, diag, perm, diag-left, left) << " " << secularEq(mid-right, col0, diag, perm, diag-right, right)   << "\n";
00746     std::cout << "     = " << secularEq(0.1*(left+right), col0, diag, perm, diag, 0)
00747               << " "       << secularEq(0.2*(left+right), col0, diag, perm, diag, 0)
00748               << " "       << secularEq(0.3*(left+right), col0, diag, perm, diag, 0)
00749               << " "       << secularEq(0.4*(left+right), col0, diag, perm, diag, 0)
00750               << " "       << secularEq(0.49*(left+right), col0, diag, perm, diag, 0)
00751               << " "       << secularEq(0.5*(left+right), col0, diag, perm, diag, 0)
00752               << " "       << secularEq(0.51*(left+right), col0, diag, perm, diag, 0)
00753               << " "       << secularEq(0.6*(left+right), col0, diag, perm, diag, 0)
00754               << " "       << secularEq(0.7*(left+right), col0, diag, perm, diag, 0)
00755               << " "       << secularEq(0.8*(left+right), col0, diag, perm, diag, 0)
00756               << " "       << secularEq(0.9*(left+right), col0, diag, perm, diag, 0) << "\n";
00757 #endif
00758     RealScalar shift = (k == actual_n-1 || fMid > 0) ? left : right;
00759     
00760     // measure everything relative to shift
00761     Map<ArrayXr> diagShifted(m_workspace.data()+4*n, n);
00762     diagShifted = diag - shift;
00763     
00764     // initial guess
00765     RealScalar muPrev, muCur;
00766     if (shift == left)
00767     {
00768       muPrev = (right - left) * RealScalar(0.1);
00769       if (k == actual_n-1) muCur = right - left;
00770       else                 muCur = (right - left) * RealScalar(0.5);
00771     }
00772     else
00773     {
00774       muPrev = -(right - left) * RealScalar(0.1);
00775       muCur = -(right - left) * RealScalar(0.5);
00776     }
00777 
00778     RealScalar fPrev = secularEq(muPrev, col0, diag, perm, diagShifted, shift);
00779     RealScalar fCur = secularEq(muCur, col0, diag, perm, diagShifted, shift);
00780     if (abs(fPrev) < abs(fCur))
00781     {
00782       swap(fPrev, fCur);
00783       swap(muPrev, muCur);
00784     }
00785 
00786     // rational interpolation: fit a function of the form a / mu + b through the two previous
00787     // iterates and use its zero to compute the next iterate
00788     bool useBisection = fPrev*fCur>0;
00789     while (fCur!=0 && abs(muCur - muPrev) > 8 * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection)
00790     {
00791       ++m_numIters;
00792 
00793       // Find a and b such that the function f(mu) = a / mu + b matches the current and previous samples.
00794       RealScalar a = (fCur - fPrev) / (1/muCur - 1/muPrev);
00795       RealScalar b = fCur - a / muCur;
00796       // And find mu such that f(mu)==0:
00797       RealScalar muZero = -a/b;
00798       RealScalar fZero = secularEq(muZero, col0, diag, perm, diagShifted, shift);
00799       
00800       muPrev = muCur;
00801       fPrev = fCur;
00802       muCur = muZero;
00803       fCur = fZero;
00804       
00805       
00806       if (shift == left  && (muCur < 0 || muCur > right - left)) useBisection = true;
00807       if (shift == right && (muCur < -(right - left) || muCur > 0)) useBisection = true;
00808       if (abs(fCur)>abs(fPrev)) useBisection = true;
00809     }
00810 
00811     // fall back on bisection method if rational interpolation did not work
00812     if (useBisection)
00813     {
00814 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
00815       std::cout << "useBisection for k = " << k << ", actual_n = " << actual_n << "\n";
00816 #endif
00817       RealScalar leftShifted, rightShifted;
00818       if (shift == left)
00819       {
00820         leftShifted = (std::numeric_limits<RealScalar>::min)();
00821         // I don't understand why the case k==0 would be special there:
00822         // if (k == 0) rightShifted = right - left; else 
00823         rightShifted = (k==actual_n-1) ? right : ((right - left) * RealScalar(0.6)); // theoretically we can take 0.5, but let's be safe
00824       }
00825       else
00826       {
00827         leftShifted = -(right - left) * RealScalar(0.6);
00828         rightShifted = -(std::numeric_limits<RealScalar>::min)();
00829       }
00830       
00831       RealScalar fLeft = secularEq(leftShifted, col0, diag, perm, diagShifted, shift);
00832 
00833 #if defined EIGEN_INTERNAL_DEBUGGING || defined EIGEN_BDCSVD_DEBUG_VERBOSE
00834       RealScalar fRight = secularEq(rightShifted, col0, diag, perm, diagShifted, shift);
00835 #endif
00836 
00837 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
00838       if(!(fLeft * fRight<0))
00839       {
00840         std::cout << "fLeft: " << leftShifted << " - " << diagShifted.head(10).transpose()  << "\n ; " << bool(left==shift) << " " << (left-shift) << "\n";
00841         std::cout << k << " : " <<  fLeft << " * " << fRight << " == " << fLeft * fRight << "  ;  " << left << " - " << right << " -> " <<  leftShifted << " " << rightShifted << "   shift=" << shift << "\n";
00842       }
00843 #endif
00844       eigen_internal_assert(fLeft * fRight < 0);
00845       
00846       while (rightShifted - leftShifted > 2 * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(leftShifted), abs(rightShifted)))
00847       {
00848         RealScalar midShifted = (leftShifted + rightShifted) / 2;
00849         fMid = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
00850         if (fLeft * fMid < 0)
00851         {
00852           rightShifted = midShifted;
00853         }
00854         else
00855         {
00856           leftShifted = midShifted;
00857           fLeft = fMid;
00858         }
00859       }
00860 
00861       muCur = (leftShifted + rightShifted) / 2;
00862     }
00863       
00864     singVals[k] = shift + muCur;
00865     shifts[k] = shift;
00866     mus[k] = muCur;
00867 
00868     // perturb singular value slightly if it equals diagonal entry to avoid division by zero later
00869     // (deflation is supposed to avoid this from happening)
00870     // - this does no seem to be necessary anymore -
00871 //     if (singVals[k] == left) singVals[k] *= 1 + NumTraits<RealScalar>::epsilon();
00872 //     if (singVals[k] == right) singVals[k] *= 1 - NumTraits<RealScalar>::epsilon();
00873   }
00874 }
00875 
00876 
00877 // zhat is perturbation of col0 for which singular vectors can be computed stably (see Section 3.1)
00878 template <typename MatrixType>
00879 void BDCSVD<MatrixType>::perturbCol0
00880    (const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const VectorType& singVals,
00881     const ArrayRef& shifts, const ArrayRef& mus, ArrayRef zhat)
00882 {
00883   using std::sqrt;
00884   Index n = col0.size();
00885   Index m = perm.size();
00886   if(m==0)
00887   {
00888     zhat.setZero();
00889     return;
00890   }
00891   Index last = perm(m-1);
00892   // The offset permits to skip deflated entries while computing zhat
00893   for (Index k = 0; k < n; ++k)
00894   {
00895     if (col0(k) == 0) // deflated
00896       zhat(k) = 0;
00897     else
00898     {
00899       // see equation (3.6)
00900       RealScalar dk = diag(k);
00901       RealScalar prod = (singVals(last) + dk) * (mus(last) + (shifts(last) - dk));
00902 
00903       for(Index l = 0; l<m; ++l)
00904       {
00905         Index i = perm(l);
00906         if(i!=k)
00907         {
00908           Index j = i<k ? i : perm(l-1);
00909           prod *= ((singVals(j)+dk) / ((diag(i)+dk))) * ((mus(j)+(shifts(j)-dk)) / ((diag(i)-dk)));
00910 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00911           if(i!=k && std::abs(((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) - 1) > 0.9 )
00912             std::cout << "     " << ((singVals(j)+dk)*(mus(j)+(shifts(j)-dk)))/((diag(i)+dk)*(diag(i)-dk)) << " == (" << (singVals(j)+dk) << " * " << (mus(j)+(shifts(j)-dk))
00913                        << ") / (" << (diag(i)+dk) << " * " << (diag(i)-dk) << ")\n";
00914 #endif
00915         }
00916       }
00917 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
00918       std::cout << "zhat(" << k << ") =  sqrt( " << prod << ")  ;  " << (singVals(last) + dk) << " * " << mus(last) + shifts(last) << " - " << dk << "\n";
00919 #endif
00920       RealScalar tmp = sqrt(prod);
00921       zhat(k) = col0(k) > 0 ? tmp : -tmp;
00922     }
00923   }
00924 }
00925 
00926 // compute singular vectors
00927 template <typename MatrixType>
00928 void BDCSVD<MatrixType>::computeSingVecs
00929    (const ArrayRef& zhat, const ArrayRef& diag, const IndicesRef &perm, const VectorType& singVals,
00930     const ArrayRef& shifts, const ArrayRef& mus, MatrixXr& U, MatrixXr& V)
00931 {
00932   Index n = zhat.size();
00933   Index m = perm.size();
00934   
00935   for (Index k = 0; k < n; ++k)
00936   {
00937     if (zhat(k) == 0)
00938     {
00939       U.col(k) = VectorType::Unit(n+1, k);
00940       if (m_compV) V.col(k) = VectorType::Unit(n, k);
00941     }
00942     else
00943     {
00944       U.col(k).setZero();
00945       for(Index l=0;l<m;++l)
00946       {
00947         Index i = perm(l);
00948         U(i,k) = zhat(i)/(((diag(i) - shifts(k)) - mus(k)) )/( (diag(i) + singVals[k]));
00949       }
00950       U(n,k) = 0;      
00951       U.col(k).normalize();
00952     
00953       if (m_compV)
00954       {
00955         V.col(k).setZero();
00956         for(Index l=1;l<m;++l)
00957         {
00958           Index i = perm(l);
00959           V(i,k) = diag(i) * zhat(i) / (((diag(i) - shifts(k)) - mus(k)) )/( (diag(i) + singVals[k]));
00960         }
00961         V(0,k) = -1;
00962         V.col(k).normalize();
00963       }
00964     }
00965   }
00966   U.col(n) = VectorType::Unit(n+1, n);
00967 }
00968 
00969 
00970 // page 12_13
00971 // i >= 1, di almost null and zi non null.
00972 // We use a rotation to zero out zi applied to the left of M
00973 template <typename MatrixType>
00974 void BDCSVD<MatrixType>::deflation43(Index firstCol, Index shift, Index i, Index size)
00975 {
00976   using std::abs;
00977   using std::sqrt;
00978   using std::pow;
00979   Index start = firstCol + shift;
00980   RealScalar c = m_computed(start, start);
00981   RealScalar s = m_computed(start+i, start);
00982   RealScalar r = sqrt(numext::abs2(c) + numext::abs2(s));
00983   if (r == 0)
00984   {
00985     m_computed(start+i, start+i) = 0;
00986     return;
00987   }
00988   m_computed(start,start) = r;  
00989   m_computed(start+i, start) = 0;
00990   m_computed(start+i, start+i) = 0;
00991   
00992   JacobiRotation<RealScalar> J(c/r,-s/r);
00993   if (m_compU)  m_naiveU.middleRows(firstCol, size+1).applyOnTheRight(firstCol, firstCol+i, J);
00994   else          m_naiveU.applyOnTheRight(firstCol, firstCol+i, J);
00995 }// end deflation 43
00996 
00997 
00998 // page 13
00999 // i,j >= 1, i!=j and |di - dj| < epsilon * norm2(M)
01000 // We apply two rotations to have zj = 0;
01001 // TODO deflation44 is still broken and not properly tested
01002 template <typename MatrixType>
01003 void BDCSVD<MatrixType>::deflation44(Index firstColu , Index firstColm, Index firstRowW, Index firstColW, Index i, Index j, Index size)
01004 {
01005   using std::abs;
01006   using std::sqrt;
01007   using std::conj;
01008   using std::pow;
01009   RealScalar c = m_computed(firstColm+i, firstColm);
01010   RealScalar s = m_computed(firstColm+j, firstColm);
01011   RealScalar r = sqrt(numext::abs2(c) + numext::abs2(s));
01012 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
01013   std::cout << "deflation 4.4: " << i << "," << j << " -> " << c << " " << s << " " << r << " ; "
01014     << m_computed(firstColm + i-1, firstColm)  << " "
01015     << m_computed(firstColm + i, firstColm)  << " "
01016     << m_computed(firstColm + i+1, firstColm) << " "
01017     << m_computed(firstColm + i+2, firstColm) << "\n";
01018   std::cout << m_computed(firstColm + i-1, firstColm + i-1)  << " "
01019     << m_computed(firstColm + i, firstColm+i)  << " "
01020     << m_computed(firstColm + i+1, firstColm+i+1) << " "
01021     << m_computed(firstColm + i+2, firstColm+i+2) << "\n";
01022 #endif
01023   if (r==0)
01024   {
01025     m_computed(firstColm + i, firstColm + i) = m_computed(firstColm + j, firstColm + j);
01026     return;
01027   }
01028   c/=r;
01029   s/=r;
01030   m_computed(firstColm + i, firstColm) = r;  
01031   m_computed(firstColm + j, firstColm + j) = m_computed(firstColm + i, firstColm + i);
01032   m_computed(firstColm + j, firstColm) = 0;
01033 
01034   JacobiRotation<RealScalar> J(c,-s);
01035   if (m_compU)  m_naiveU.middleRows(firstColu, size+1).applyOnTheRight(firstColu + i, firstColu + j, J);
01036   else          m_naiveU.applyOnTheRight(firstColu+i, firstColu+j, J);
01037   if (m_compV)  m_naiveV.middleRows(firstRowW, size).applyOnTheRight(firstColW + i, firstColW + j, J);
01038 }// end deflation 44
01039 
01040 
01041 // acts on block from (firstCol+shift, firstCol+shift) to (lastCol+shift, lastCol+shift) [inclusive]
01042 template <typename MatrixType>
01043 void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index firstRowW, Index firstColW, Index shift)
01044 {
01045   using std::sqrt;
01046   using std::abs;
01047   const Index length = lastCol + 1 - firstCol;
01048   
01049   Block<MatrixXr,Dynamic,1> col0(m_computed, firstCol+shift, firstCol+shift, length, 1);
01050   Diagonal<MatrixXr> fulldiag(m_computed);
01051   VectorBlock<Diagonal<MatrixXr>,Dynamic> diag(fulldiag, firstCol+shift, length);
01052   
01053   const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
01054   RealScalar maxDiag = diag.tail((std::max)(Index(1),length-1)).cwiseAbs().maxCoeff();
01055   RealScalar epsilon_strict = numext::maxi<RealScalar>(considerZero,NumTraits<RealScalar>::epsilon() * maxDiag);
01056   RealScalar epsilon_coarse = 8 * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(col0.cwiseAbs().maxCoeff(), maxDiag);
01057   
01058 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
01059   assert(m_naiveU.allFinite());
01060   assert(m_naiveV.allFinite());
01061   assert(m_computed.allFinite());
01062 #endif
01063 
01064 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE  
01065   std::cout << "\ndeflate:" << diag.head(k+1).transpose() << "  |  " << diag.segment(k+1,length-k-1).transpose() << "\n";
01066 #endif
01067   
01068   //condition 4.1
01069   if (diag(0) < epsilon_coarse)
01070   { 
01071 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
01072     std::cout << "deflation 4.1, because " << diag(0) << " < " << epsilon_coarse << "\n";
01073 #endif
01074     diag(0) = epsilon_coarse;
01075   }
01076 
01077   //condition 4.2
01078   for (Index i=1;i<length;++i)
01079     if (abs(col0(i)) < epsilon_strict)
01080     {
01081 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
01082       std::cout << "deflation 4.2, set z(" << i << ") to zero because " << abs(col0(i)) << " < " << epsilon_strict << "  (diag(" << i << ")=" << diag(i) << ")\n";
01083 #endif
01084       col0(i) = 0;
01085     }
01086 
01087   //condition 4.3
01088   for (Index i=1;i<length; i++)
01089     if (diag(i) < epsilon_coarse)
01090     {
01091 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
01092       std::cout << "deflation 4.3, cancel z(" << i << ")=" << col0(i) << " because diag(" << i << ")=" << diag(i) << " < " << epsilon_coarse << "\n";
01093 #endif
01094       deflation43(firstCol, shift, i, length);
01095     }
01096 
01097 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
01098   assert(m_naiveU.allFinite());
01099   assert(m_naiveV.allFinite());
01100   assert(m_computed.allFinite());
01101 #endif
01102 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
01103   std::cout << "to be sorted: " << diag.transpose() << "\n\n";
01104 #endif
01105   {
01106     // Check for total deflation
01107     // If we have a total deflation, then we have to consider col0(0)==diag(0) as a singular value during sorting
01108     bool total_deflation = (col0.tail(length-1).array()<considerZero).all();
01109     
01110     // Sort the diagonal entries, since diag(1:k-1) and diag(k:length) are already sorted, let's do a sorted merge.
01111     // First, compute the respective permutation.
01112     Index *permutation = m_workspaceI.data();
01113     {
01114       permutation[0] = 0;
01115       Index p = 1;
01116       
01117       // Move deflated diagonal entries at the end.
01118       for(Index i=1; i<length; ++i)
01119         if(abs(diag(i))<considerZero)
01120           permutation[p++] = i;
01121         
01122       Index i=1, j=k+1;
01123       for( ; p < length; ++p)
01124       {
01125              if (i > k)             permutation[p] = j++;
01126         else if (j >= length)       permutation[p] = i++;
01127         else if (diag(i) < diag(j)) permutation[p] = j++;
01128         else                        permutation[p] = i++;
01129       }
01130     }
01131     
01132     // If we have a total deflation, then we have to insert diag(0) at the right place
01133     if(total_deflation)
01134     {
01135       for(Index i=1; i<length; ++i)
01136       {
01137         Index pi = permutation[i];
01138         if(abs(diag(pi))<considerZero || diag(0)<diag(pi))
01139           permutation[i-1] = permutation[i];
01140         else
01141         {
01142           permutation[i-1] = 0;
01143           break;
01144         }
01145       }
01146     }
01147     
01148     // Current index of each col, and current column of each index
01149     Index *realInd = m_workspaceI.data()+length;
01150     Index *realCol = m_workspaceI.data()+2*length;
01151     
01152     for(int pos = 0; pos< length; pos++)
01153     {
01154       realCol[pos] = pos;
01155       realInd[pos] = pos;
01156     }
01157     
01158     for(Index i = total_deflation?0:1; i < length; i++)
01159     {
01160       const Index pi = permutation[length - (total_deflation ? i+1 : i)];
01161       const Index J = realCol[pi];
01162       
01163       using std::swap;
01164       // swap diagonal and first column entries:
01165       swap(diag(i), diag(J));
01166       if(i!=0 && J!=0) swap(col0(i), col0(J));
01167 
01168       // change columns
01169       if (m_compU) m_naiveU.col(firstCol+i).segment(firstCol, length + 1).swap(m_naiveU.col(firstCol+J).segment(firstCol, length + 1));
01170       else         m_naiveU.col(firstCol+i).segment(0, 2)                .swap(m_naiveU.col(firstCol+J).segment(0, 2));
01171       if (m_compV) m_naiveV.col(firstColW + i).segment(firstRowW, length).swap(m_naiveV.col(firstColW + J).segment(firstRowW, length));
01172 
01173       //update real pos
01174       const Index realI = realInd[i];
01175       realCol[realI] = J;
01176       realCol[pi] = i;
01177       realInd[J] = realI;
01178       realInd[i] = pi;
01179     }
01180   }
01181 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
01182   std::cout << "sorted: " << diag.transpose().format(bdcsvdfmt) << "\n";
01183   std::cout << "      : " << col0.transpose() << "\n\n";
01184 #endif
01185     
01186   //condition 4.4
01187   {
01188     Index i = length-1;
01189     while(i>0 && (abs(diag(i))<considerZero || abs(col0(i))<considerZero)) --i;
01190     for(; i>1;--i)
01191        if( (diag(i) - diag(i-1)) < NumTraits<RealScalar>::epsilon()*maxDiag )
01192       {
01193 #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
01194         std::cout << "deflation 4.4 with i = " << i << " because " << (diag(i) - diag(i-1)) << " < " << NumTraits<RealScalar>::epsilon()*diag(i) << "\n";
01195 #endif
01196         eigen_internal_assert(abs(diag(i) - diag(i-1))<epsilon_coarse && " diagonal entries are not properly sorted");
01197         deflation44(firstCol, firstCol + shift, firstRowW, firstColW, i-1, i, length);
01198       }
01199   }
01200   
01201 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
01202   for(Index j=2;j<length;++j)
01203     assert(diag(j-1)<=diag(j) || abs(diag(j))<considerZero);
01204 #endif
01205   
01206 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
01207   assert(m_naiveU.allFinite());
01208   assert(m_naiveV.allFinite());
01209   assert(m_computed.allFinite());
01210 #endif
01211 }//end deflation
01212 
01213 #ifndef __CUDACC__
01214 
01220 template<typename Derived>
01221 BDCSVD<typename MatrixBase<Derived>::PlainObject>
01222 MatrixBase<Derived>::bdcSvd(unsigned int computationOptions) const
01223 {
01224   return BDCSVD<PlainObject>(*this, computationOptions);
01225 }
01226 #endif
01227 
01228 } // end namespace Eigen
01229 
01230 #endif
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends