From 4b7c3abbea97c604d56f732cd5ba5f8ff8356814 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 17 Oct 2014 15:32:55 +0200 Subject: Fix D&C SVD wrt zero matrices --- unsupported/Eigen/src/BDCSVD/BDCSVD.h | 30 +++++++++--------------------- 1 file changed, 9 insertions(+), 21 deletions(-) (limited to 'unsupported') diff --git a/unsupported/Eigen/src/BDCSVD/BDCSVD.h b/unsupported/Eigen/src/BDCSVD/BDCSVD.h index e8e795589..e2551cf88 100644 --- a/unsupported/Eigen/src/BDCSVD/BDCSVD.h +++ b/unsupported/Eigen/src/BDCSVD/BDCSVD.h @@ -194,7 +194,7 @@ public: }; //end class BDCSVD -// Methode to allocate ans initialize matrix and attributes +// Method to allocate and initialize matrix and attributes template void BDCSVD::allocate(Index rows, Index cols, unsigned int computationOptions) { @@ -215,24 +215,6 @@ void BDCSVD::allocate(Index rows, Index cols, unsigned int computati if (m_compV) m_naiveV = MatrixXr::Zero(m_diagSize, m_diagSize); }// end allocate -// Methode which compute the BDCSVD for the int -template<> -BDCSVD >& BDCSVD >::compute(const MatrixType& matrix, unsigned int computationOptions) -{ - allocate(matrix.rows(), matrix.cols(), computationOptions); - m_nonzeroSingularValues = 0; - m_computed = Matrix::Zero(rows(), cols()); - - m_singularValues.head(m_diagSize).setZero(); - - if (m_computeFullU) m_matrixU.setZero(rows(), rows()); - if (m_computeFullV) m_matrixV.setZero(cols(), cols()); - m_isInitialized = true; - return *this; -} - - -// Methode which compute the BDCSVD template BDCSVD& BDCSVD::compute(const MatrixType& matrix, unsigned int computationOptions) { @@ -612,6 +594,8 @@ void BDCSVD::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec #endif #ifdef EIGEN_BDCSVD_SANITY_CHECKS + assert(U.allFinite()); + assert(V.allFinite()); assert((U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() < 1e-14 * n); assert((V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() < 1e-14 * n); assert(m_naiveU.allFinite()); @@ -836,8 +820,12 @@ void BDCSVD::perturbCol0 using std::sqrt; Index n = col0.size(); Index m = perm.size(); - Index last = perm(m-1); - + if(m==0) + { + zhat.setZero(); + return; + } + Index last = last = perm(m-1); // The offset permits to skip deflated entries while computing zhat for (Index k = 0; k < n; ++k) { -- cgit v1.2.3