diff options
-rw-r--r-- | test/jacobisvd.cpp | 413 | ||||
-rw-r--r-- | test/svd_common.h | 454 | ||||
-rw-r--r-- | unsupported/test/bdcsvd.cpp | 233 | ||||
-rw-r--r-- | unsupported/test/jacobisvd.cpp | 198 | ||||
-rw-r--r-- | unsupported/test/svd_common.h | 261 |
5 files changed, 540 insertions, 1019 deletions
diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index cd04db5be..bfcadce95 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> // // This Source Code Form is subject to the terms of the Mozilla @@ -14,273 +14,47 @@ #include "main.h" #include <Eigen/SVD> -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_check_full(const MatrixType& m, const JacobiSVD<MatrixType, QRPreconditioner>& svd) -{ - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - typedef typename MatrixType::Scalar Scalar; - typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; - typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; - - MatrixType sigma = MatrixType::Zero(rows,cols); - sigma.diagonal() = svd.singularValues().template cast<Scalar>(); - MatrixUType u = svd.matrixU(); - MatrixVType v = svd.matrixV(); - - VERIFY_IS_APPROX(m, u * sigma * v.adjoint()); - VERIFY_IS_UNITARY(u); - VERIFY_IS_UNITARY(v); -} - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_compare_to_full(const MatrixType& m, - unsigned int computationOptions, - const JacobiSVD<MatrixType, QRPreconditioner>& referenceSvd) -{ - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - Index diagSize = (std::min)(rows, cols); - - JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions); - - VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues()); - if(computationOptions & ComputeFullU) - VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU()); - if(computationOptions & ComputeThinU) - VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize)); - if(computationOptions & ComputeFullV) - VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV()); - if(computationOptions & ComputeThinV) - VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize)); -} - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions) -{ - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType; - typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType; - - RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols)); - JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions); - - if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8); - else if(internal::is_same<RealScalar,float>::value) svd.setThreshold(1e-4); - - SolutionType x = svd.solve(rhs); - - RealScalar residual = (m*x-rhs).norm(); - // Check that there is no significantly better solution in the neighborhood of x - if(!test_isMuchSmallerThan(residual,rhs.norm())) - { - // If the residual is very small, then we have an exact solution, so we are already good. - for(int k=0;k<x.rows();++k) - { - SolutionType y(x); - y.row(k).array() += 2*NumTraits<RealScalar>::epsilon(); - RealScalar residual_y = (m*y-rhs).norm(); - VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); - - y.row(k) = x.row(k).array() - 2*NumTraits<RealScalar>::epsilon(); - residual_y = (m*y-rhs).norm(); - VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); - } - } - - // evaluate normal equation which works also for least-squares solutions - if(internal::is_same<RealScalar,double>::value) - { - // This test is not stable with single precision. - // This is probably because squaring m signicantly affects the precision. - VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); - } - - // check minimal norm solutions - { - // generate a full-rank m x n problem with m<n - enum { - RankAtCompileTime2 = ColsAtCompileTime==Dynamic ? Dynamic : (ColsAtCompileTime)/2+1, - RowsAtCompileTime3 = ColsAtCompileTime==Dynamic ? Dynamic : ColsAtCompileTime+1 - }; - typedef Matrix<Scalar, RankAtCompileTime2, ColsAtCompileTime> MatrixType2; - typedef Matrix<Scalar, RankAtCompileTime2, 1> RhsType2; - typedef Matrix<Scalar, ColsAtCompileTime, RankAtCompileTime2> MatrixType2T; - Index rank = RankAtCompileTime2==Dynamic ? internal::random<Index>(1,cols) : Index(RankAtCompileTime2); - MatrixType2 m2(rank,cols); - int guard = 0; - do { - m2.setRandom(); - } while(m2.jacobiSvd().setThreshold(test_precision<Scalar>()).rank()!=rank && (++guard)<10); - VERIFY(guard<10); - RhsType2 rhs2 = RhsType2::Random(rank); - // use QR to find a reference minimal norm solution - HouseholderQR<MatrixType2T> qr(m2.adjoint()); - Matrix<Scalar,Dynamic,1> tmp = qr.matrixQR().topLeftCorner(rank,rank).template triangularView<Upper>().adjoint().solve(rhs2); - tmp.conservativeResize(cols); - tmp.tail(cols-rank).setZero(); - SolutionType x21 = qr.householderQ() * tmp; - // now check with SVD - JacobiSVD<MatrixType2, ColPivHouseholderQRPreconditioner> svd2(m2, computationOptions); - SolutionType x22 = svd2.solve(rhs2); - VERIFY_IS_APPROX(m2*x21, rhs2); - VERIFY_IS_APPROX(m2*x22, rhs2); - VERIFY_IS_APPROX(x21, x22); - - // Now check with a rank deficient matrix - typedef Matrix<Scalar, RowsAtCompileTime3, ColsAtCompileTime> MatrixType3; - typedef Matrix<Scalar, RowsAtCompileTime3, 1> RhsType3; - Index rows3 = RowsAtCompileTime3==Dynamic ? internal::random<Index>(rank+1,2*cols) : Index(RowsAtCompileTime3); - Matrix<Scalar,RowsAtCompileTime3,Dynamic> C = Matrix<Scalar,RowsAtCompileTime3,Dynamic>::Random(rows3,rank); - MatrixType3 m3 = C * m2; - RhsType3 rhs3 = C * rhs2; - JacobiSVD<MatrixType3, ColPivHouseholderQRPreconditioner> svd3(m3, computationOptions); - SolutionType x3 = svd3.solve(rhs3); - VERIFY_IS_APPROX(m3*x3, rhs3); - VERIFY_IS_APPROX(m3*x21, rhs3); - VERIFY_IS_APPROX(m2*x3, rhs2); - - VERIFY_IS_APPROX(x21, x3); - } -} - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_test_all_computation_options(const MatrixType& m) -{ - if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) - return; - JacobiSVD<MatrixType, QRPreconditioner> fullSvd(m, ComputeFullU|ComputeFullV); - CALL_SUBTEST(( jacobisvd_check_full(m, fullSvd) )); - CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV) )); - - #if defined __INTEL_COMPILER - // remark #111: statement is unreachable - #pragma warning disable 111 - #endif - if(QRPreconditioner == FullPivHouseholderQRPreconditioner) - return; - - CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU, fullSvd) )); - CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullV, fullSvd) )); - CALL_SUBTEST(( jacobisvd_compare_to_full(m, 0, fullSvd) )); - - if (MatrixType::ColsAtCompileTime == Dynamic) { - // thin U/V are only available with dynamic number of columns - CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) )); - CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinV, fullSvd) )); - CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) )); - CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU , fullSvd) )); - CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) )); - CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV) )); - CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV) )); - CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV) )); - - // test reconstruction - typedef typename MatrixType::Index Index; - Index diagSize = (std::min)(m.rows(), m.cols()); - JacobiSVD<MatrixType, QRPreconditioner> svd(m, ComputeThinU | ComputeThinV); - VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint()); - } -} +#define SVD_DEFAULT(M) JacobiSVD<M> +#define SVD_FOR_MIN_NORM(M) JacobiSVD<M,ColPivHouseholderQRPreconditioner> +#include "svd_common.h" +// Check all variants of JacobiSVD template<typename MatrixType> void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) { MatrixType m = a; if(pickrandom) - { - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; - Index diagSize = (std::min)(a.rows(), a.cols()); - RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4; - s = internal::random<RealScalar>(1,s); - Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize); - for(Index k=0; k<diagSize; ++k) - d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); - m = Matrix<Scalar,Dynamic,Dynamic>::Random(a.rows(),diagSize) * d.asDiagonal() * Matrix<Scalar,Dynamic,Dynamic>::Random(diagSize,a.cols()); - // cancel some coeffs - Index n = internal::random<Index>(0,m.size()-1); - for(Index i=0; i<n; ++i) - m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0); - } + svd_fill_random(m); - CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m) )); - CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m) )); - CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m) )); - CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m) )); + CALL_SUBTEST(( svd_test_all_computation_options<JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> >(m, true) )); // check full only + CALL_SUBTEST(( svd_test_all_computation_options<JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner> >(m, false) )); + CALL_SUBTEST(( svd_test_all_computation_options<JacobiSVD<MatrixType, HouseholderQRPreconditioner> >(m, false) )); + if(m.rows()==m.cols()) + CALL_SUBTEST(( svd_test_all_computation_options<JacobiSVD<MatrixType, NoQRPreconditioner> >(m, false) )); } template<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m) { - typedef typename MatrixType::Scalar Scalar; + svd_verify_assert<JacobiSVD<MatrixType> >(m); typedef typename MatrixType::Index Index; Index rows = m.rows(); Index cols = m.cols(); enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime }; - typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType; - - RhsType rhs(rows); - - JacobiSVD<MatrixType> svd; - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.singularValues()) - VERIFY_RAISES_ASSERT(svd.matrixV()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) MatrixType a = MatrixType::Zero(rows, cols); a.setZero(); - svd.compute(a, 0); - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.matrixV()) - svd.singularValues(); - VERIFY_RAISES_ASSERT(svd.solve(rhs)) if (ColsAtCompileTime == Dynamic) { - svd.compute(a, ComputeThinU); - svd.matrixU(); - VERIFY_RAISES_ASSERT(svd.matrixV()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - - svd.compute(a, ComputeThinV); - svd.matrixV(); - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> svd_fullqr; VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeFullU|ComputeThinV)) VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeThinV)) VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeFullV)) } - else - { - VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU)) - VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV)) - } } template<typename MatrixType> @@ -296,165 +70,17 @@ void jacobisvd_method() VERIFY_IS_APPROX(m.jacobiSvd(ComputeFullU|ComputeFullV).solve(m), m); } -// work around stupid msvc error when constructing at compile time an expression that involves -// a division by zero, even if the numeric type has floating point -template<typename Scalar> -EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); } - -// workaround aggressive optimization in ICC -template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } - -template<typename MatrixType> -void jacobisvd_inf_nan() -{ - // all this function does is verify we don't iterate infinitely on nan/inf values - - JacobiSVD<MatrixType> svd; - typedef typename MatrixType::Scalar Scalar; - Scalar some_inf = Scalar(1) / zero<Scalar>(); - VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf)); - svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); - - Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); - VERIFY(nan != nan); - svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV); - - MatrixType m = MatrixType::Zero(10,10); - m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf; - svd.compute(m, ComputeFullU | ComputeFullV); - - m = MatrixType::Zero(10,10); - m(internal::random<int>(0,9), internal::random<int>(0,9)) = nan; - svd.compute(m, ComputeFullU | ComputeFullV); - - // regression test for bug 791 - m.resize(3,3); - m << 0, 2*NumTraits<Scalar>::epsilon(), 0.5, - 0, -0.5, 0, - nan, 0, 0; - svd.compute(m, ComputeFullU | ComputeFullV); - - m.resize(4,4); - m << 1, 0, 0, 0, - 0, 3, 1, 2e-308, - 1, 0, 1, nan, - 0, nan, nan, 0; - svd.compute(m, ComputeFullU | ComputeFullV); -} - -// Regression test for bug 286: JacobiSVD loops indefinitely with some -// matrices containing denormal numbers. -void jacobisvd_underoverflow() -{ -#if defined __INTEL_COMPILER -// shut up warning #239: floating point underflow -#pragma warning push -#pragma warning disable 239 -#endif - Matrix2d M; - M << -7.90884e-313, -4.94e-324, - 0, 5.60844e-313; - JacobiSVD<Matrix2d> svd; - svd.compute(M,ComputeFullU|ComputeFullV); - jacobisvd_check_full(M,svd); - - VectorXd value_set(9); - value_set << 0, 1, -1, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -4.94e-223, 4.94e-223; - Array4i id(0,0,0,0); - int k = 0; - do - { - M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3)); - svd.compute(M,ComputeFullU|ComputeFullV); - jacobisvd_check_full(M,svd); - - id(k)++; - if(id(k)>=value_set.size()) - { - while(k<3 && id(k)>=value_set.size()) id(++k)++; - id.head(k).setZero(); - k=0; - } - - } while((id<int(value_set.size())).all()); - -#if defined __INTEL_COMPILER -#pragma warning pop -#endif - - // Check for overflow: - Matrix3d M3; - M3 << 4.4331978442502944e+307, -5.8585363752028680e+307, 6.4527017443412964e+307, - 3.7841695601406358e+307, 2.4331702789740617e+306, -3.5235707140272905e+307, - -8.7190887618028355e+307, -7.3453213709232193e+307, -2.4367363684472105e+307; - - JacobiSVD<Matrix3d> svd3; - svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely - jacobisvd_check_full(M3,svd3); -} - -void jacobisvd_preallocate() -{ - Vector3f v(3.f, 2.f, 1.f); - MatrixXf m = v.asDiagonal(); - - internal::set_is_malloc_allowed(false); - VERIFY_RAISES_ASSERT(VectorXf tmp(10);) - JacobiSVD<MatrixXf> svd; - internal::set_is_malloc_allowed(true); - svd.compute(m); - VERIFY_IS_APPROX(svd.singularValues(), v); - - JacobiSVD<MatrixXf> svd2(3,3); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - VERIFY_IS_APPROX(svd2.singularValues(), v); - VERIFY_RAISES_ASSERT(svd2.matrixU()); - VERIFY_RAISES_ASSERT(svd2.matrixV()); - svd2.compute(m, ComputeFullU | ComputeFullV); - VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); - VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - - JacobiSVD<MatrixXf> svd3(3,3,ComputeFullU|ComputeFullV); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - VERIFY_IS_APPROX(svd2.singularValues(), v); - VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); - VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); - internal::set_is_malloc_allowed(false); - svd2.compute(m, ComputeFullU|ComputeFullV); - internal::set_is_malloc_allowed(true); -} - void test_jacobisvd() { CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) )); CALL_SUBTEST_4(( jacobisvd_verify_assert(Matrix4d()) )); CALL_SUBTEST_7(( jacobisvd_verify_assert(MatrixXf(10,12)) )); CALL_SUBTEST_8(( jacobisvd_verify_assert(MatrixXcd(7,5)) )); + + svd_all_trivial_2x2(jacobisvd<Matrix2cd>); + svd_all_trivial_2x2(jacobisvd<Matrix2d>); for(int i = 0; i < g_repeat; i++) { - Matrix2cd m; - m << 0, 1, - 0, 1; - CALL_SUBTEST_1(( jacobisvd(m, false) )); - m << 1, 0, - 1, 0; - CALL_SUBTEST_1(( jacobisvd(m, false) )); - - Matrix2d n; - n << 0, 0, - 0, 0; - CALL_SUBTEST_2(( jacobisvd(n, false) )); - n << 0, 0, - 0, 1; - CALL_SUBTEST_2(( jacobisvd(n, false) )); - CALL_SUBTEST_3(( jacobisvd<Matrix3f>() )); CALL_SUBTEST_4(( jacobisvd<Matrix4d>() )); CALL_SUBTEST_5(( jacobisvd<Matrix<float,3,5> >() )); @@ -473,8 +99,8 @@ void test_jacobisvd() (void) c; // Test on inf/nan matrix - CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() ); - CALL_SUBTEST_10( jacobisvd_inf_nan<MatrixXd>() ); + CALL_SUBTEST_7( (svd_inf_nan<JacobiSVD<MatrixXf>, MatrixXf>()) ); + CALL_SUBTEST_10( (svd_inf_nan<JacobiSVD<MatrixXd>, MatrixXd>()) ); } CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); @@ -488,8 +114,7 @@ void test_jacobisvd() CALL_SUBTEST_7( JacobiSVD<MatrixXf>(10,10) ); // Check that preallocation avoids subsequent mallocs - CALL_SUBTEST_9( jacobisvd_preallocate() ); + CALL_SUBTEST_9( svd_preallocate() ); - // Regression check for bug 286 - CALL_SUBTEST_2( jacobisvd_underoverflow() ); + CALL_SUBTEST_2( svd_underoverflow() ); } diff --git a/test/svd_common.h b/test/svd_common.h new file mode 100644 index 000000000..4631939e5 --- /dev/null +++ b/test/svd_common.h @@ -0,0 +1,454 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef SVD_DEFAULT +#error a macro SVD_DEFAULT(MatrixType) must be defined prior to including svd_common.h +#endif + +#ifndef SVD_FOR_MIN_NORM +#error a macro SVD_FOR_MIN_NORM(MatrixType) must be defined prior to including svd_common.h +#endif + +// Check that the matrix m is properly reconstructed and that the U and V factors are unitary +// The SVD must have already been computed. +template<typename SvdType, typename MatrixType> +void svd_check_full(const MatrixType& m, const SvdType& svd) +{ + typedef typename MatrixType::Index Index; + Index rows = m.rows(); + Index cols = m.cols(); + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime + }; + + typedef typename MatrixType::Scalar Scalar; + typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; + typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; + + MatrixType sigma = MatrixType::Zero(rows,cols); + sigma.diagonal() = svd.singularValues().template cast<Scalar>(); + MatrixUType u = svd.matrixU(); + MatrixVType v = svd.matrixV(); + + VERIFY_IS_APPROX(m, u * sigma * v.adjoint()); + VERIFY_IS_UNITARY(u); + VERIFY_IS_UNITARY(v); +} + +// Compare partial SVD defined by computationOptions to a full SVD referenceSvd +template<typename SvdType, typename MatrixType> +void svd_compare_to_full(const MatrixType& m, + unsigned int computationOptions, + const SvdType& referenceSvd) +{ + typedef typename MatrixType::Index Index; + Index rows = m.rows(); + Index cols = m.cols(); + Index diagSize = (std::min)(rows, cols); + + SvdType svd(m, computationOptions); + + VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues()); + if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU()); + if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize)); + if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV()); + if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize)); +} + +// +template<typename SvdType, typename MatrixType> +void svd_least_square(const MatrixType& m, unsigned int computationOptions) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + Index rows = m.rows(); + Index cols = m.cols(); + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime + }; + + typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType; + typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType; + + RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols)); + SvdType svd(m, computationOptions); + + if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8); + else if(internal::is_same<RealScalar,float>::value) svd.setThreshold(1e-4); + + SolutionType x = svd.solve(rhs); + + RealScalar residual = (m*x-rhs).norm(); + // Check that there is no significantly better solution in the neighborhood of x + if(!test_isMuchSmallerThan(residual,rhs.norm())) + { + // If the residual is very small, then we have an exact solution, so we are already good. + for(int k=0;k<x.rows();++k) + { + SolutionType y(x); + y.row(k).array() += 2*NumTraits<RealScalar>::epsilon(); + RealScalar residual_y = (m*y-rhs).norm(); + VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); + + y.row(k) = x.row(k).array() - 2*NumTraits<RealScalar>::epsilon(); + residual_y = (m*y-rhs).norm(); + VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); + } + } + + // evaluate normal equation which works also for least-squares solutions + if(internal::is_same<RealScalar,double>::value) + { + // This test is not stable with single precision. + // This is probably because squaring m signicantly affects the precision. + VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); + } +} + +// check minimal norm solutions, the inoput matrix m is only used to recover problem size +template<typename MatrixType> +void svd_min_norm(const MatrixType& m, unsigned int computationOptions) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; + Index cols = m.cols(); + + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime + }; + + typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType; + + // generate a full-rank m x n problem with m<n + enum { + RankAtCompileTime2 = ColsAtCompileTime==Dynamic ? Dynamic : (ColsAtCompileTime)/2+1, + RowsAtCompileTime3 = ColsAtCompileTime==Dynamic ? Dynamic : ColsAtCompileTime+1 + }; + typedef Matrix<Scalar, RankAtCompileTime2, ColsAtCompileTime> MatrixType2; + typedef Matrix<Scalar, RankAtCompileTime2, 1> RhsType2; + typedef Matrix<Scalar, ColsAtCompileTime, RankAtCompileTime2> MatrixType2T; + Index rank = RankAtCompileTime2==Dynamic ? internal::random<Index>(1,cols) : Index(RankAtCompileTime2); + MatrixType2 m2(rank,cols); + int guard = 0; + do { + m2.setRandom(); + } while(SVD_FOR_MIN_NORM(MatrixType2)(m2).setThreshold(test_precision<Scalar>()).rank()!=rank && (++guard)<10); + VERIFY(guard<10); + RhsType2 rhs2 = RhsType2::Random(rank); + // use QR to find a reference minimal norm solution + HouseholderQR<MatrixType2T> qr(m2.adjoint()); + Matrix<Scalar,Dynamic,1> tmp = qr.matrixQR().topLeftCorner(rank,rank).template triangularView<Upper>().adjoint().solve(rhs2); + tmp.conservativeResize(cols); + tmp.tail(cols-rank).setZero(); + SolutionType x21 = qr.householderQ() * tmp; + // now check with SVD + SVD_FOR_MIN_NORM(MatrixType2) svd2(m2, computationOptions); + SolutionType x22 = svd2.solve(rhs2); + VERIFY_IS_APPROX(m2*x21, rhs2); + VERIFY_IS_APPROX(m2*x22, rhs2); + VERIFY_IS_APPROX(x21, x22); + + // Now check with a rank deficient matrix + typedef Matrix<Scalar, RowsAtCompileTime3, ColsAtCompileTime> MatrixType3; + typedef Matrix<Scalar, RowsAtCompileTime3, 1> RhsType3; + Index rows3 = RowsAtCompileTime3==Dynamic ? internal::random<Index>(rank+1,2*cols) : Index(RowsAtCompileTime3); + Matrix<Scalar,RowsAtCompileTime3,Dynamic> C = Matrix<Scalar,RowsAtCompileTime3,Dynamic>::Random(rows3,rank); + MatrixType3 m3 = C * m2; + RhsType3 rhs3 = C * rhs2; + SVD_FOR_MIN_NORM(MatrixType3) svd3(m3, computationOptions); + SolutionType x3 = svd3.solve(rhs3); + VERIFY_IS_APPROX(m3*x3, rhs3); + VERIFY_IS_APPROX(m3*x21, rhs3); + VERIFY_IS_APPROX(m2*x3, rhs2); + + VERIFY_IS_APPROX(x21, x3); +} + +// Check full, compare_to_full, least_square, and min_norm for all possible compute-options +template<typename SvdType, typename MatrixType> +void svd_test_all_computation_options(const MatrixType& m, bool full_only) +{ +// if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) +// return; + SvdType fullSvd(m, ComputeFullU|ComputeFullV); + CALL_SUBTEST(( svd_check_full(m, fullSvd) )); + CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeFullU | ComputeFullV) )); + CALL_SUBTEST(( svd_min_norm(m, ComputeFullU | ComputeFullV) )); + + #if defined __INTEL_COMPILER + // remark #111: statement is unreachable + #pragma warning disable 111 + #endif + if(full_only) + return; + + CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullU, fullSvd) )); + CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullV, fullSvd) )); + CALL_SUBTEST(( svd_compare_to_full(m, 0, fullSvd) )); + + if (MatrixType::ColsAtCompileTime == Dynamic) { + // thin U/V are only available with dynamic number of columns + CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) )); + CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinV, fullSvd) )); + CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) )); + CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU , fullSvd) )); + CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) )); + + CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeFullU | ComputeThinV) )); + CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeThinU | ComputeFullV) )); + CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeThinU | ComputeThinV) )); + + CALL_SUBTEST(( svd_min_norm(m, ComputeFullU | ComputeThinV) )); + CALL_SUBTEST(( svd_min_norm(m, ComputeThinU | ComputeFullV) )); + CALL_SUBTEST(( svd_min_norm(m, ComputeThinU | ComputeThinV) )); + + // test reconstruction + typedef typename MatrixType::Index Index; + Index diagSize = (std::min)(m.rows(), m.cols()); + SvdType svd(m, ComputeThinU | ComputeThinV); + VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint()); + } +} + +template<typename MatrixType> +void svd_fill_random(MatrixType &m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + Index diagSize = (std::min)(m.rows(), m.cols()); + RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4; + s = internal::random<RealScalar>(1,s); + Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize); + for(Index k=0; k<diagSize; ++k) + d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); + m = Matrix<Scalar,Dynamic,Dynamic>::Random(m.rows(),diagSize) * d.asDiagonal() * Matrix<Scalar,Dynamic,Dynamic>::Random(diagSize,m.cols()); + // cancel some coeffs + Index n = internal::random<Index>(0,m.size()-1); + for(Index i=0; i<n; ++i) + m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0); +} + + +// work around stupid msvc error when constructing at compile time an expression that involves +// a division by zero, even if the numeric type has floating point +template<typename Scalar> +EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); } + +// workaround aggressive optimization in ICC +template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } + +// all this function does is verify we don't iterate infinitely on nan/inf values +template<typename SvdType, typename MatrixType> +void svd_inf_nan() +{ + SvdType svd; + typedef typename MatrixType::Scalar Scalar; + Scalar some_inf = Scalar(1) / zero<Scalar>(); + VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf)); + svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); + + Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); + VERIFY(nan != nan); + svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV); + + MatrixType m = MatrixType::Zero(10,10); + m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf; + svd.compute(m, ComputeFullU | ComputeFullV); + + m = MatrixType::Zero(10,10); + m(internal::random<int>(0,9), internal::random<int>(0,9)) = nan; + svd.compute(m, ComputeFullU | ComputeFullV); + + // regression test for bug 791 + m.resize(3,3); + m << 0, 2*NumTraits<Scalar>::epsilon(), 0.5, + 0, -0.5, 0, + nan, 0, 0; + svd.compute(m, ComputeFullU | ComputeFullV); + + m.resize(4,4); + m << 1, 0, 0, 0, + 0, 3, 1, 2e-308, + 1, 0, 1, nan, + 0, nan, nan, 0; + svd.compute(m, ComputeFullU | ComputeFullV); +} + +// Regression test for bug 286: JacobiSVD loops indefinitely with some +// matrices containing denormal numbers. +void svd_underoverflow() +{ +#if defined __INTEL_COMPILER +// shut up warning #239: floating point underflow +#pragma warning push +#pragma warning disable 239 +#endif + Matrix2d M; + M << -7.90884e-313, -4.94e-324, + 0, 5.60844e-313; + SVD_DEFAULT(Matrix2d) svd; + svd.compute(M,ComputeFullU|ComputeFullV); + svd_check_full(M,svd); + + // Check all 2x2 matrices made with the following coefficients: + VectorXd value_set(9); + value_set << 0, 1, -1, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -4.94e-223, 4.94e-223; + Array4i id(0,0,0,0); + int k = 0; + do + { + 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); + + id(k)++; + if(id(k)>=value_set.size()) + { + while(k<3 && id(k)>=value_set.size()) id(++k)++; + id.head(k).setZero(); + k=0; + } + + } while((id<int(value_set.size())).all()); + +#if defined __INTEL_COMPILER +#pragma warning pop +#endif + + // Check for overflow: + Matrix3d M3; + M3 << 4.4331978442502944e+307, -5.8585363752028680e+307, 6.4527017443412964e+307, + 3.7841695601406358e+307, 2.4331702789740617e+306, -3.5235707140272905e+307, + -8.7190887618028355e+307, -7.3453213709232193e+307, -2.4367363684472105e+307; + + SVD_DEFAULT(Matrix3d) svd3; + svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely + svd_check_full(M3,svd3); +} + +// void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) + +template<typename MatrixType> +void svd_all_trivial_2x2( void (*cb)(const MatrixType&,bool) ) +{ + MatrixType M; + VectorXd value_set(3); + value_set << 0, 1, -1; + Array4i id(0,0,0,0); + int k = 0; + do + { + M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3)); + + cb(M,false); + + id(k)++; + if(id(k)>=value_set.size()) + { + while(k<3 && id(k)>=value_set.size()) id(++k)++; + id.head(k).setZero(); + k=0; + } + + } while((id<int(value_set.size())).all()); +} + +void svd_preallocate() +{ + Vector3f v(3.f, 2.f, 1.f); + MatrixXf m = v.asDiagonal(); + + internal::set_is_malloc_allowed(false); + VERIFY_RAISES_ASSERT(VectorXf tmp(10);) + SVD_DEFAULT(MatrixXf) svd; + internal::set_is_malloc_allowed(true); + svd.compute(m); + VERIFY_IS_APPROX(svd.singularValues(), v); + + SVD_DEFAULT(MatrixXf) svd2(3,3); + internal::set_is_malloc_allowed(false); + svd2.compute(m); + internal::set_is_malloc_allowed(true); + VERIFY_IS_APPROX(svd2.singularValues(), v); + VERIFY_RAISES_ASSERT(svd2.matrixU()); + VERIFY_RAISES_ASSERT(svd2.matrixV()); + svd2.compute(m, ComputeFullU | ComputeFullV); + VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); + VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); + internal::set_is_malloc_allowed(false); + svd2.compute(m); + internal::set_is_malloc_allowed(true); + + SVD_DEFAULT(MatrixXf) svd3(3,3,ComputeFullU|ComputeFullV); + internal::set_is_malloc_allowed(false); + svd2.compute(m); + internal::set_is_malloc_allowed(true); + VERIFY_IS_APPROX(svd2.singularValues(), v); + VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); + VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); + internal::set_is_malloc_allowed(false); + svd2.compute(m, ComputeFullU|ComputeFullV); + internal::set_is_malloc_allowed(true); +} + +template<typename SvdType,typename MatrixType> +void svd_verify_assert(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; + Index rows = m.rows(); + Index cols = m.cols(); + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime + }; + + typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType; + RhsType rhs(rows); + SvdType svd; + VERIFY_RAISES_ASSERT(svd.matrixU()) + VERIFY_RAISES_ASSERT(svd.singularValues()) + VERIFY_RAISES_ASSERT(svd.matrixV()) + VERIFY_RAISES_ASSERT(svd.solve(rhs)) + MatrixType a = MatrixType::Zero(rows, cols); + a.setZero(); + svd.compute(a, 0); + VERIFY_RAISES_ASSERT(svd.matrixU()) + VERIFY_RAISES_ASSERT(svd.matrixV()) + svd.singularValues(); + VERIFY_RAISES_ASSERT(svd.solve(rhs)) + + if (ColsAtCompileTime == Dynamic) + { + svd.compute(a, ComputeThinU); + svd.matrixU(); + VERIFY_RAISES_ASSERT(svd.matrixV()) + VERIFY_RAISES_ASSERT(svd.solve(rhs)) + svd.compute(a, ComputeThinV); + svd.matrixV(); + VERIFY_RAISES_ASSERT(svd.matrixU()) + VERIFY_RAISES_ASSERT(svd.solve(rhs)) + } + else + { + VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU)) + VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV)) + } +} + +#undef SVD_DEFAULT +#undef SVD_FOR_MIN_NORM diff --git a/unsupported/test/bdcsvd.cpp b/unsupported/test/bdcsvd.cpp index 115a649b0..4ad991522 100644 --- a/unsupported/test/bdcsvd.cpp +++ b/unsupported/test/bdcsvd.cpp @@ -10,204 +10,105 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/ -#include "svd_common.h" +// discard stack allocation as that too bypasses malloc +#define EIGEN_STACK_ALLOCATION_LIMIT 0 +#define EIGEN_RUNTIME_NO_MALLOC + +#include "main.h" +#include <unsupported/Eigen/BDCSVD> #include <iostream> #include <Eigen/LU> -// check if "svd" is the good image of "m" -template<typename MatrixType> -void bdcsvd_check_full(const MatrixType& m, const BDCSVD<MatrixType>& svd) -{ - svd_check_full< MatrixType, BDCSVD< MatrixType > >(m, svd); -} - -// Compare to a reference value -template<typename MatrixType> -void bdcsvd_compare_to_full(const MatrixType& m, - unsigned int computationOptions, - const BDCSVD<MatrixType>& referenceSvd) -{ - svd_compare_to_full< MatrixType, BDCSVD< MatrixType > >(m, computationOptions, referenceSvd); -} // end bdcsvd_compare_to_full +#define SVD_DEFAULT(M) BDCSVD<M> +// #define SVD_FOR_MIN_NORM(M) BDCSVD<M> +#define SVD_FOR_MIN_NORM(M) JacobiSVD<M,ColPivHouseholderQRPreconditioner> +#include "../../test/svd_common.h" -template<typename MatrixType> -void bdcsvd_solve(const MatrixType& m, unsigned int computationOptions) -{ - svd_solve< MatrixType, BDCSVD< MatrixType > >(m, computationOptions); -} // end template bdcsvd_solve - - -// test the computations options -template<typename MatrixType> -void bdcsvd_test_all_computation_options(const MatrixType& m) -{ - BDCSVD<MatrixType> fullSvd(m, ComputeFullU|ComputeFullV); - svd_test_computation_options_1< MatrixType, BDCSVD< MatrixType > >(m, fullSvd); - svd_test_computation_options_2< MatrixType, BDCSVD< MatrixType > >(m, fullSvd); -} // end bdcsvd_test_all_computation_options - - -// Call a test with all the computations options +// Check all variants of JacobiSVD template<typename MatrixType> void bdcsvd(const MatrixType& a = MatrixType(), bool pickrandom = true) { - MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a; - bdcsvd_test_all_computation_options<MatrixType>(m); -} // end template bdcsvd - - -// verify assert -template<typename MatrixType> -void bdcsvd_verify_assert(const MatrixType& m) -{ - svd_verify_assert< MatrixType, BDCSVD< MatrixType > >(m); -}// end template bdcsvd_verify_assert + MatrixType m = a; + if(pickrandom) + svd_fill_random(m); + CALL_SUBTEST(( svd_test_all_computation_options<BDCSVD<MatrixType> >(m, false) )); +} -// test weird values -template<typename MatrixType> -void bdcsvd_inf_nan() -{ - svd_inf_nan< MatrixType, BDCSVD< MatrixType > >(); -}// end template bdcsvd_inf_nan - - - -void bdcsvd_preallocate() -{ - svd_preallocate< BDCSVD< MatrixXf > >(); -} // end bdcsvd_preallocate - +// template<typename MatrixType> +// void bdcsvd_method() +// { +// enum { Size = MatrixType::RowsAtCompileTime }; +// typedef typename MatrixType::RealScalar RealScalar; +// typedef Matrix<RealScalar, Size, 1> RealVecType; +// MatrixType m = MatrixType::Identity(); +// VERIFY_IS_APPROX(m.bdcSvd().singularValues(), RealVecType::Ones()); +// VERIFY_RAISES_ASSERT(m.bdcSvd().matrixU()); +// VERIFY_RAISES_ASSERT(m.bdcSvd().matrixV()); +// VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).solve(m), m); +// } // compare the Singular values returned with Jacobi and Bdc template<typename MatrixType> void compare_bdc_jacobi(const MatrixType& a = MatrixType(), unsigned int computationOptions = 0) { - std::cout << "debut compare" << std::endl; MatrixType m = MatrixType::Random(a.rows(), a.cols()); BDCSVD<MatrixType> bdc_svd(m); JacobiSVD<MatrixType> jacobi_svd(m); VERIFY_IS_APPROX(bdc_svd.singularValues(), jacobi_svd.singularValues()); - if(computationOptions & ComputeFullU) - VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); - if(computationOptions & ComputeThinU) - VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); - if(computationOptions & ComputeFullV) - VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); - if(computationOptions & ComputeThinV) - VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); - std::cout << "fin compare" << std::endl; -} // end template compare_bdc_jacobi - - -// call the tests + if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); + if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); + if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); + if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); +} + void test_bdcsvd() { - // test of Dynamic defined Matrix (42, 42) of float - CALL_SUBTEST_11(( bdcsvd_verify_assert<Matrix<float,Dynamic,Dynamic> > - (Matrix<float,Dynamic,Dynamic>(42,42)) )); - CALL_SUBTEST_11(( compare_bdc_jacobi<Matrix<float,Dynamic,Dynamic> > - (Matrix<float,Dynamic,Dynamic>(42,42), 0) )); - CALL_SUBTEST_11(( bdcsvd<Matrix<float,Dynamic,Dynamic> > - (Matrix<float,Dynamic,Dynamic>(42,42)) )); - - // test of Dynamic defined Matrix (50, 50) of double - CALL_SUBTEST_13(( bdcsvd_verify_assert<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(50,50)) )); - CALL_SUBTEST_13(( compare_bdc_jacobi<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(50,50), 0) )); - CALL_SUBTEST_13(( bdcsvd<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(50, 50)) )); - - // test of Dynamic defined Matrix (22, 22) of complex double - CALL_SUBTEST_14(( bdcsvd_verify_assert<Matrix<std::complex<double>,Dynamic,Dynamic> > - (Matrix<std::complex<double>,Dynamic,Dynamic>(22,22)) )); - CALL_SUBTEST_14(( compare_bdc_jacobi<Matrix<std::complex<double>,Dynamic,Dynamic> > - (Matrix<std::complex<double>, Dynamic, Dynamic> (22,22), 0) )); - CALL_SUBTEST_14(( bdcsvd<Matrix<std::complex<double>,Dynamic,Dynamic> > - (Matrix<std::complex<double>,Dynamic,Dynamic>(22, 22)) )); - - // test of Dynamic defined Matrix (10, 10) of int - //CALL_SUBTEST_15(( bdcsvd_verify_assert<Matrix<int,Dynamic,Dynamic> > - // (Matrix<int,Dynamic,Dynamic>(10,10)) )); - //CALL_SUBTEST_15(( compare_bdc_jacobi<Matrix<int,Dynamic,Dynamic> > - // (Matrix<int,Dynamic,Dynamic>(10,10), 0) )); - //CALL_SUBTEST_15(( bdcsvd<Matrix<int,Dynamic,Dynamic> > - // (Matrix<int,Dynamic,Dynamic>(10, 10)) )); + CALL_SUBTEST_3(( svd_verify_assert<BDCSVD<Matrix3f> >(Matrix3f()) )); + CALL_SUBTEST_4(( svd_verify_assert<BDCSVD<Matrix4d> >(Matrix4d()) )); + 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>); - // test of Dynamic defined Matrix (8, 6) of double - - CALL_SUBTEST_16(( bdcsvd_verify_assert<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(8,6)) )); - CALL_SUBTEST_16(( compare_bdc_jacobi<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(8, 6), 0) )); - CALL_SUBTEST_16(( bdcsvd<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(8, 6)) )); - - - - // test of Dynamic defined Matrix (36, 12) of float - CALL_SUBTEST_17(( compare_bdc_jacobi<Matrix<float,Dynamic,Dynamic> > - (Matrix<float,Dynamic,Dynamic>(36, 12), 0) )); - CALL_SUBTEST_17(( bdcsvd<Matrix<float,Dynamic,Dynamic> > - (Matrix<float,Dynamic,Dynamic>(36, 12)) )); - - // test of Dynamic defined Matrix (5, 8) of double - CALL_SUBTEST_18(( compare_bdc_jacobi<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(5, 8), 0) )); - CALL_SUBTEST_18(( bdcsvd<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(5, 8)) )); - - - // non regression tests - CALL_SUBTEST_3(( bdcsvd_verify_assert(Matrix3f()) )); - CALL_SUBTEST_4(( bdcsvd_verify_assert(Matrix4d()) )); - CALL_SUBTEST_7(( bdcsvd_verify_assert(MatrixXf(10,12)) )); - CALL_SUBTEST_8(( bdcsvd_verify_assert(MatrixXcd(7,5)) )); - - // SUBTESTS 1 and 2 on specifics matrix for(int i = 0; i < g_repeat; i++) { - Matrix2cd m; - m << 0, 1, - 0, 1; - CALL_SUBTEST_1(( bdcsvd(m, false) )); - m << 1, 0, - 1, 0; - CALL_SUBTEST_1(( bdcsvd(m, false) )); - - Matrix2d n; - n << 0, 0, - 0, 0; - CALL_SUBTEST_2(( bdcsvd(n, false) )); - n << 0, 0, - 0, 1; - CALL_SUBTEST_2(( bdcsvd(n, false) )); +// 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); - // Statics matrix don't work with BDSVD yet - // bdc algo on a random 3x3 float matrix - // CALL_SUBTEST_3(( bdcsvd<Matrix3f>() )); - // bdc algo on a random 4x4 double matrix - // CALL_SUBTEST_4(( bdcsvd<Matrix4d>() )); - // bdc algo on a random 3x5 float matrix - // CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() )); - - int r = internal::random<int>(1, 30), - c = internal::random<int>(1, 30); - CALL_SUBTEST_7(( bdcsvd<MatrixXf>(MatrixXf(r,c)) )); - CALL_SUBTEST_8(( bdcsvd<MatrixXcd>(MatrixXcd(r,c)) )); + TEST_SET_BUT_UNUSED_VARIABLE(r) + TEST_SET_BUT_UNUSED_VARIABLE(c) + + CALL_SUBTEST_6(( bdcsvd(Matrix<double,Dynamic,2>(r,2)) )); + CALL_SUBTEST_7(( bdcsvd(MatrixXf(r,c)) )); + CALL_SUBTEST_7(( compare_bdc_jacobi(MatrixXf(r,c)) )); + CALL_SUBTEST_10(( bdcsvd(MatrixXd(r,c)) )); + CALL_SUBTEST_10(( compare_bdc_jacobi(MatrixXd(r,c)) )); + CALL_SUBTEST_8(( bdcsvd(MatrixXcd(r,c)) )); + CALL_SUBTEST_8(( compare_bdc_jacobi(MatrixXcd(r,c)) )); (void) r; (void) c; // Test on inf/nan matrix - CALL_SUBTEST_7( bdcsvd_inf_nan<MatrixXf>() ); + CALL_SUBTEST_7( (svd_inf_nan<BDCSVD<MatrixXf>, MatrixXf>()) ); + CALL_SUBTEST_10( (svd_inf_nan<BDCSVD<MatrixXd>, MatrixXd>()) ); } - CALL_SUBTEST_7(( bdcsvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); - CALL_SUBTEST_8(( bdcsvd<MatrixXcd>(MatrixXcd(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3))) )); + // test matrixbase method +// CALL_SUBTEST_1(( bdcsvd_method<Matrix2cd>() )); +// CALL_SUBTEST_3(( bdcsvd_method<Matrix3f>() )); // Test problem size constructors CALL_SUBTEST_7( BDCSVD<MatrixXf>(10,10) ); -} // end test_bdcsvd + // Check that preallocation avoids subsequent mallocs + CALL_SUBTEST_9( svd_preallocate() ); + + CALL_SUBTEST_2( svd_underoverflow() ); +} + diff --git a/unsupported/test/jacobisvd.cpp b/unsupported/test/jacobisvd.cpp deleted file mode 100644 index b4e884eee..000000000 --- a/unsupported/test/jacobisvd.cpp +++ /dev/null @@ -1,198 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> -// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> -// -// This Source Code Form is subject to the terms of the Mozilla -// Public License v. 2.0. If a copy of the MPL was not distributed -// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include "svd_common.h" - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_check_full(const MatrixType& m, const JacobiSVD<MatrixType, QRPreconditioner>& svd) -{ - svd_check_full<MatrixType, JacobiSVD<MatrixType, QRPreconditioner > >(m, svd); -} - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_compare_to_full(const MatrixType& m, - unsigned int computationOptions, - const JacobiSVD<MatrixType, QRPreconditioner>& referenceSvd) -{ - svd_compare_to_full<MatrixType, JacobiSVD<MatrixType, QRPreconditioner> >(m, computationOptions, referenceSvd); -} - - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions) -{ - svd_solve< MatrixType, JacobiSVD< MatrixType, QRPreconditioner > >(m, computationOptions); -} - - - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_test_all_computation_options(const MatrixType& m) -{ - - if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) - return; - - JacobiSVD< MatrixType, QRPreconditioner > fullSvd(m, ComputeFullU|ComputeFullV); - svd_test_computation_options_1< MatrixType, JacobiSVD< MatrixType, QRPreconditioner > >(m, fullSvd); - - if(QRPreconditioner == FullPivHouseholderQRPreconditioner) - return; - svd_test_computation_options_2< MatrixType, JacobiSVD< MatrixType, QRPreconditioner > >(m, fullSvd); - -} - -template<typename MatrixType> -void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) -{ - MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a; - - jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m); - jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m); - jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m); - jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m); -} - - -template<typename MatrixType> -void jacobisvd_verify_assert(const MatrixType& m) -{ - - svd_verify_assert<MatrixType, JacobiSVD< MatrixType > >(m); - - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - MatrixType a = MatrixType::Zero(rows, cols); - a.setZero(); - - if (ColsAtCompileTime == Dynamic) - { - JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> svd_fullqr; - VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeFullU|ComputeThinV)) - VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeThinV)) - VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeFullV)) - } -} - -template<typename MatrixType> -void jacobisvd_method() -{ - enum { Size = MatrixType::RowsAtCompileTime }; - typedef typename MatrixType::RealScalar RealScalar; - typedef Matrix<RealScalar, Size, 1> RealVecType; - MatrixType m = MatrixType::Identity(); - VERIFY_IS_APPROX(m.jacobiSvd().singularValues(), RealVecType::Ones()); - VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixU()); - VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixV()); - VERIFY_IS_APPROX(m.jacobiSvd(ComputeFullU|ComputeFullV).solve(m), m); -} - - - -template<typename MatrixType> -void jacobisvd_inf_nan() -{ - svd_inf_nan<MatrixType, JacobiSVD< MatrixType > >(); -} - - -// Regression test for bug 286: JacobiSVD loops indefinitely with some -// matrices containing denormal numbers. -void jacobisvd_bug286() -{ -#if defined __INTEL_COMPILER -// shut up warning #239: floating point underflow -#pragma warning push -#pragma warning disable 239 -#endif - Matrix2d M; - M << -7.90884e-313, -4.94e-324, - 0, 5.60844e-313; -#if defined __INTEL_COMPILER -#pragma warning pop -#endif - JacobiSVD<Matrix2d> svd; - svd.compute(M); // just check we don't loop indefinitely -} - - -void jacobisvd_preallocate() -{ - svd_preallocate< JacobiSVD <MatrixXf> >(); -} - -void test_jacobisvd() -{ - CALL_SUBTEST_11(( jacobisvd<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(16, 6)) )); - - CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) )); - CALL_SUBTEST_4(( jacobisvd_verify_assert(Matrix4d()) )); - CALL_SUBTEST_7(( jacobisvd_verify_assert(MatrixXf(10,12)) )); - CALL_SUBTEST_8(( jacobisvd_verify_assert(MatrixXcd(7,5)) )); - - for(int i = 0; i < g_repeat; i++) { - Matrix2cd m; - m << 0, 1, - 0, 1; - CALL_SUBTEST_1(( jacobisvd(m, false) )); - m << 1, 0, - 1, 0; - CALL_SUBTEST_1(( jacobisvd(m, false) )); - - Matrix2d n; - n << 0, 0, - 0, 0; - CALL_SUBTEST_2(( jacobisvd(n, false) )); - n << 0, 0, - 0, 1; - CALL_SUBTEST_2(( jacobisvd(n, false) )); - - CALL_SUBTEST_3(( jacobisvd<Matrix3f>() )); - CALL_SUBTEST_4(( jacobisvd<Matrix4d>() )); - CALL_SUBTEST_5(( jacobisvd<Matrix<float,3,5> >() )); - CALL_SUBTEST_6(( jacobisvd<Matrix<double,Dynamic,2> >(Matrix<double,Dynamic,2>(10,2)) )); - - int r = internal::random<int>(1, 30), - c = internal::random<int>(1, 30); - CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) )); - CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) )); - (void) r; - (void) c; - - // Test on inf/nan matrix - CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() ); - } - - CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); - CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3))) )); - - - // test matrixbase method - CALL_SUBTEST_1(( jacobisvd_method<Matrix2cd>() )); - CALL_SUBTEST_3(( jacobisvd_method<Matrix3f>() )); - - - // Test problem size constructors - CALL_SUBTEST_7( JacobiSVD<MatrixXf>(10,10) ); - - // Check that preallocation avoids subsequent mallocs - CALL_SUBTEST_9( jacobisvd_preallocate() ); - - // Regression check for bug 286 - CALL_SUBTEST_2( jacobisvd_bug286() ); -} diff --git a/unsupported/test/svd_common.h b/unsupported/test/svd_common.h deleted file mode 100644 index 6a96e7c74..000000000 --- a/unsupported/test/svd_common.h +++ /dev/null @@ -1,261 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> -// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> -// -// Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com> -// Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr> -// Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr> -// Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr> -// -// This Source Code Form is subject to the terms of the Mozilla -// Public License v. 2.0. If a copy of the MPL was not distributed -// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. - -// discard stack allocation as that too bypasses malloc -#define EIGEN_STACK_ALLOCATION_LIMIT 0 -#define EIGEN_RUNTIME_NO_MALLOC - -#include "main.h" -#include <unsupported/Eigen/BDCSVD> -#include <Eigen/LU> - - -// check if "svd" is the good image of "m" -template<typename MatrixType, typename SVD> -void svd_check_full(const MatrixType& m, const SVD& svd) -{ - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - typedef typename MatrixType::Scalar Scalar; - typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; - typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; - - - MatrixType sigma = MatrixType::Zero(rows, cols); - sigma.diagonal() = svd.singularValues().template cast<Scalar>(); - MatrixUType u = svd.matrixU(); - MatrixVType v = svd.matrixV(); - VERIFY_IS_APPROX(m, u * sigma * v.adjoint()); - VERIFY_IS_UNITARY(u); - VERIFY_IS_UNITARY(v); -} // end svd_check_full - - - -// Compare to a reference value -template<typename MatrixType, typename SVD> -void svd_compare_to_full(const MatrixType& m, - unsigned int computationOptions, - const SVD& referenceSvd) -{ - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - Index diagSize = (std::min)(rows, cols); - - SVD svd(m, computationOptions); - - VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues()); - if(computationOptions & ComputeFullU) - VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU()); - if(computationOptions & ComputeThinU) - VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize)); - if(computationOptions & ComputeFullV) - VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV()); - if(computationOptions & ComputeThinV) - VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize)); -} // end svd_compare_to_full - - - -template<typename MatrixType, typename SVD> -void svd_solve(const MatrixType& m, unsigned int computationOptions) -{ - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType; - typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType; - - RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols)); - SVD svd(m, computationOptions); - SolutionType x = svd.solve(rhs); - // evaluate normal equation which works also for least-squares solutions - VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); -} // end svd_solve - - -// test computations options -// 2 functions because Jacobisvd can return before the second function -template<typename MatrixType, typename SVD> -void svd_test_computation_options_1(const MatrixType& m, const SVD& fullSvd) -{ - svd_check_full< MatrixType, SVD >(m, fullSvd); - svd_solve< MatrixType, SVD >(m, ComputeFullU | ComputeFullV); -} - - -template<typename MatrixType, typename SVD> -void svd_test_computation_options_2(const MatrixType& m, const SVD& fullSvd) -{ - svd_compare_to_full< MatrixType, SVD >(m, ComputeFullU, fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, ComputeFullV, fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, 0, fullSvd); - - if (MatrixType::ColsAtCompileTime == Dynamic) { - // thin U/V are only available with dynamic number of columns - - svd_compare_to_full< MatrixType, SVD >(m, ComputeFullU|ComputeThinV, fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, ComputeThinV, fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU|ComputeFullV, fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU , fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU|ComputeThinV, fullSvd); - svd_solve<MatrixType, SVD>(m, ComputeFullU | ComputeThinV); - svd_solve<MatrixType, SVD>(m, ComputeThinU | ComputeFullV); - svd_solve<MatrixType, SVD>(m, ComputeThinU | ComputeThinV); - - typedef typename MatrixType::Index Index; - Index diagSize = (std::min)(m.rows(), m.cols()); - SVD svd(m, ComputeThinU | ComputeThinV); - VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint()); - } -} - -template<typename MatrixType, typename SVD> -void svd_verify_assert(const MatrixType& m) -{ - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType; - RhsType rhs(rows); - SVD svd; - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.singularValues()) - VERIFY_RAISES_ASSERT(svd.matrixV()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - MatrixType a = MatrixType::Zero(rows, cols); - a.setZero(); - svd.compute(a, 0); - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.matrixV()) - svd.singularValues(); - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - - if (ColsAtCompileTime == Dynamic) - { - svd.compute(a, ComputeThinU); - svd.matrixU(); - VERIFY_RAISES_ASSERT(svd.matrixV()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - svd.compute(a, ComputeThinV); - svd.matrixV(); - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - } - else - { - VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU)) - VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV)) - } -} - -// work around stupid msvc error when constructing at compile time an expression that involves -// a division by zero, even if the numeric type has floating point -template<typename Scalar> -EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); } - -// workaround aggressive optimization in ICC -template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } - - -template<typename MatrixType, typename SVD> -void svd_inf_nan() -{ - // all this function does is verify we don't iterate infinitely on nan/inf values - - SVD svd; - typedef typename MatrixType::Scalar Scalar; - Scalar some_inf = Scalar(1) / zero<Scalar>(); - VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf)); - svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); - - Scalar some_nan = zero<Scalar> () / zero<Scalar> (); - VERIFY(some_nan != some_nan); - svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV); - - MatrixType m = MatrixType::Zero(10,10); - m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf; - svd.compute(m, ComputeFullU | ComputeFullV); - - m = MatrixType::Zero(10,10); - m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan; - svd.compute(m, ComputeFullU | ComputeFullV); -} - - -template<typename SVD> -void svd_preallocate() -{ - Vector3f v(3.f, 2.f, 1.f); - MatrixXf m = v.asDiagonal(); - - internal::set_is_malloc_allowed(false); - VERIFY_RAISES_ASSERT(VectorXf v(10);) - SVD svd; - internal::set_is_malloc_allowed(true); - svd.compute(m); - VERIFY_IS_APPROX(svd.singularValues(), v); - - SVD svd2(3,3); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - VERIFY_IS_APPROX(svd2.singularValues(), v); - VERIFY_RAISES_ASSERT(svd2.matrixU()); - VERIFY_RAISES_ASSERT(svd2.matrixV()); - svd2.compute(m, ComputeFullU | ComputeFullV); - VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); - VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - - SVD svd3(3,3,ComputeFullU|ComputeFullV); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - VERIFY_IS_APPROX(svd2.singularValues(), v); - VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); - VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); - internal::set_is_malloc_allowed(false); - svd2.compute(m, ComputeFullU|ComputeFullV); - internal::set_is_malloc_allowed(true); -} - - - - - |