diff options
| -rw-r--r-- | test/svd_common.h | 6 | ||||
| -rw-r--r-- | unsupported/Eigen/src/BDCSVD/BDCSVD.h | 35 | ||||
| -rw-r--r-- | unsupported/test/bdcsvd.cpp | 10 |
3 files changed, 28 insertions, 23 deletions
diff --git a/test/svd_common.h b/test/svd_common.h index 4631939e5..b880efce5 100644 --- a/test/svd_common.h +++ b/test/svd_common.h @@ -301,7 +301,7 @@ void svd_underoverflow() 0, 5.60844e-313; SVD_DEFAULT(Matrix2d) svd; svd.compute(M,ComputeFullU|ComputeFullV); - svd_check_full(M,svd); + CALL_SUBTEST( svd_check_full(M,svd) ); // Check all 2x2 matrices made with the following coefficients: VectorXd value_set(9); @@ -312,7 +312,7 @@ void svd_underoverflow() { M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3)); svd.compute(M,ComputeFullU|ComputeFullV); - svd_check_full(M,svd); + CALL_SUBTEST( svd_check_full(M,svd) ); id(k)++; if(id(k)>=value_set.size()) @@ -336,7 +336,7 @@ void svd_underoverflow() SVD_DEFAULT(Matrix3d) svd3; svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely - svd_check_full(M3,svd3); + CALL_SUBTEST( svd_check_full(M3,svd3) ); } // void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) diff --git a/unsupported/Eigen/src/BDCSVD/BDCSVD.h b/unsupported/Eigen/src/BDCSVD/BDCSVD.h index b13ee415e..5bf9b0ae2 100644 --- a/unsupported/Eigen/src/BDCSVD/BDCSVD.h +++ b/unsupported/Eigen/src/BDCSVD/BDCSVD.h @@ -236,26 +236,29 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign { allocate(matrix.rows(), matrix.cols(), computationOptions); using std::abs; - - //**** step 1 Bidiagonalization - MatrixType copy; - if (m_isTranspose) copy = matrix.adjoint(); - else copy = matrix; + //**** step 0 - Copy the input matrix and apply scaling to reduce over/under-flows + RealScalar scale = matrix.cwiseAbs().maxCoeff(); + if(scale==RealScalar(0)) scale = RealScalar(1); + MatrixX copy; + if (m_isTranspose) copy = matrix.adjoint()/scale; + else copy = matrix/scale; + + //**** step 1 - Bidiagonalization internal::UpperBidiagonalization<MatrixX> bid(copy); - //**** step 2 Divide + //**** step 2 - Divide & Conquer m_naiveU.setZero(); m_naiveV.setZero(); m_computed.topRows(m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose(); m_computed.template bottomRows<1>().setZero(); divide(0, m_diagSize - 1, 0, 0, 0); - //**** step 3 copy + //**** step 3 - Copy singular values and vectors for (int i=0; i<m_diagSize; i++) { RealScalar a = abs(m_computed.coeff(i, i)); - m_singularValues.coeffRef(i) = a; + m_singularValues.coeffRef(i) = a * scale; if (a == 0) { m_nonzeroSingularValues = i; @@ -268,11 +271,13 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign break; } } -// std::cout << "m_naiveU\n" << m_naiveU << "\n\n"; -// std::cout << "m_naiveV\n" << m_naiveV << "\n\n"; +#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE + std::cout << "m_naiveU\n" << m_naiveU << "\n\n"; + std::cout << "m_naiveV\n" << m_naiveV << "\n\n"; +#endif if(m_isTranspose) copyUV(bid.householderV(), bid.householderU(), m_naiveV, m_naiveU); else copyUV(bid.householderU(), bid.householderV(), m_naiveU, m_naiveV); -// std::cout << "DONE\n"; + m_isInitialized = true; return *this; }// end compute @@ -569,13 +574,13 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec computeSingVecs(zhat, diag, singVals, shifts, mus, U, V); #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE - std::cout << "U^T U: " << (U.transpose() * U - MatrixType(MatrixType::Identity(U.cols(),U.cols()))).norm() << "\n"; - std::cout << "V^T V: " << (V.transpose() * V - MatrixType(MatrixType::Identity(V.cols(),V.cols()))).norm() << "\n"; + std::cout << "U^T U: " << (U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() << "\n"; + std::cout << "V^T V: " << (V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() << "\n"; #endif #ifdef EIGEN_BDCSVD_SANITY_CHECKS - assert((U.transpose() * U - MatrixType(MatrixType::Identity(U.cols(),U.cols()))).norm() < 1e-14 * n); - assert((V.transpose() * V - MatrixType(MatrixType::Identity(V.cols(),V.cols()))).norm() < 1e-14 * n); + 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()); assert(m_naiveV.allFinite()); assert(m_computed.allFinite()); diff --git a/unsupported/test/bdcsvd.cpp b/unsupported/test/bdcsvd.cpp index 4ad991522..95cdb6a2e 100644 --- a/unsupported/test/bdcsvd.cpp +++ b/unsupported/test/bdcsvd.cpp @@ -70,13 +70,13 @@ void test_bdcsvd() CALL_SUBTEST_7(( svd_verify_assert<BDCSVD<MatrixXf> >(MatrixXf(10,12)) )); CALL_SUBTEST_8(( svd_verify_assert<BDCSVD<MatrixXcd> >(MatrixXcd(7,5)) )); -// svd_all_trivial_2x2(bdcsvd<Matrix2cd>); -// svd_all_trivial_2x2(bdcsvd<Matrix2d>); + CALL_SUBTEST_1(( svd_all_trivial_2x2(bdcsvd<Matrix2cd>) )); + CALL_SUBTEST_1(( svd_all_trivial_2x2(bdcsvd<Matrix2d>) )); for(int i = 0; i < g_repeat; i++) { -// CALL_SUBTEST_3(( bdcsvd<Matrix3f>() )); -// CALL_SUBTEST_4(( bdcsvd<Matrix4d>() )); -// CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() )); + CALL_SUBTEST_3(( bdcsvd<Matrix3f>() )); + CALL_SUBTEST_4(( bdcsvd<Matrix4d>() )); + CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() )); int r = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2), c = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2); |