From 8ba8d90063ab6b297b1ba912f806eaa71cbdf003 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Fri, 8 Oct 2010 10:42:40 -0400 Subject: add option to compute thin U/V. By default nothing is computed. You have to ask explicitly for thin/full U/V if you want them. --- test/jacobisvd.cpp | 105 +++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 74 insertions(+), 31 deletions(-) (limited to 'test/jacobisvd.cpp') diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 56b91c983..32b25f9d2 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -28,7 +28,7 @@ #include template -void svd_with_qr_preconditioner(const MatrixType& m = MatrixType(), bool pickrandom = true) +void jacobisvd_check_full(const MatrixType& m, const JacobiSVD& svd) { typedef typename MatrixType::Index Index; Index rows = m.rows(); @@ -46,33 +46,76 @@ void svd_with_qr_preconditioner(const MatrixType& m = MatrixType(), bool pickran typedef Matrix ColVectorType; typedef Matrix InputVectorType; - MatrixType a; - if(pickrandom) a = MatrixType::Random(rows,cols); - else a = m; - - JacobiSVD svd(a, ComputeU|ComputeV); MatrixType sigma = MatrixType::Zero(rows,cols); sigma.diagonal() = svd.singularValues().template cast(); MatrixUType u = svd.matrixU(); MatrixVType v = svd.matrixV(); - //std::cout << "a\n" << a << std::endl; - //std::cout << "b\n" << u * sigma * v.adjoint() << std::endl; - - VERIFY_IS_APPROX(a, u * sigma * v.adjoint()); + VERIFY_IS_APPROX(m, u * sigma * v.adjoint()); VERIFY_IS_UNITARY(u); VERIFY_IS_UNITARY(v); } +template +void jacobisvd_compare_to_full(const MatrixType& m, + unsigned int computationOptions, + const JacobiSVD& referenceSvd) +{ + typedef typename MatrixType::Index Index; + Index rows = m.rows(); + Index cols = m.cols(); + Index diagSize = std::min(rows, cols); + + JacobiSVD svd(m, computationOptions); + + VERIFY_IS_EQUAL(svd.singularValues(), referenceSvd.singularValues()); + if(computationOptions & ComputeFullU) + VERIFY_IS_EQUAL(svd.matrixU(), referenceSvd.matrixU()); + if(computationOptions & ComputeThinU) + VERIFY_IS_EQUAL(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize)); + if(computationOptions & ComputeFullV) + VERIFY_IS_EQUAL(svd.matrixV(), referenceSvd.matrixV()); + if(computationOptions & ComputeThinV) + VERIFY_IS_EQUAL(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize)); +} + +template +void jacobisvd_test_all_computation_options(const MatrixType& m) +{ + if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) + return; + JacobiSVD fullSvd(m, ComputeFullU|ComputeFullV); + + jacobisvd_check_full(m, fullSvd); + + if(QRPreconditioner == FullPivHouseholderQRPreconditioner) + return; + + jacobisvd_compare_to_full(m, ComputeFullU, fullSvd); + jacobisvd_compare_to_full(m, ComputeFullV, fullSvd); + jacobisvd_compare_to_full(m, 0, fullSvd); + + if (MatrixType::ColsAtCompileTime == Dynamic) { + // thin U/V are only available with dynamic number of columns + jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd); + jacobisvd_compare_to_full(m, ComputeThinV, fullSvd); + jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd); + jacobisvd_compare_to_full(m, ComputeThinU , fullSvd); + jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd); + } +} + template -void svd(const MatrixType& m = MatrixType(), bool pickrandom = true) +void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) { - svd_with_qr_preconditioner(m, pickrandom); - svd_with_qr_preconditioner(m, pickrandom); - svd_with_qr_preconditioner(m, pickrandom); + MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a; + jacobisvd_test_all_computation_options(m); + jacobisvd_test_all_computation_options(m); + jacobisvd_test_all_computation_options(m); + jacobisvd_test_all_computation_options(m); } -template void svd_verify_assert() +template void jacobisvd_verify_assert() { MatrixType tmp; @@ -93,29 +136,29 @@ void test_jacobisvd() Matrix2cd m; m << 0, 1, 0, 1; - CALL_SUBTEST_1(( svd(m, false) )); + CALL_SUBTEST_1(( jacobisvd(m, false) )); m << 1, 0, 1, 0; - CALL_SUBTEST_1(( svd(m, false) )); + CALL_SUBTEST_1(( jacobisvd(m, false) )); Matrix2d n; n << 1, 1, 1, -1; - CALL_SUBTEST_2(( svd(n, false) )); - CALL_SUBTEST_3(( svd() )); - CALL_SUBTEST_4(( svd() )); - CALL_SUBTEST_5(( svd >() )); - CALL_SUBTEST_6(( svd >(Matrix(10,2)) )); - - CALL_SUBTEST_7(( svd(MatrixXf(50,50)) )); - CALL_SUBTEST_8(( svd(MatrixXcd(14,7)) )); + CALL_SUBTEST_2(( jacobisvd(n, false) )); + CALL_SUBTEST_3(( jacobisvd() )); + CALL_SUBTEST_4(( jacobisvd() )); + CALL_SUBTEST_5(( jacobisvd >() )); + CALL_SUBTEST_6(( jacobisvd >(Matrix(10,2)) )); + + CALL_SUBTEST_7(( jacobisvd(MatrixXf(50,50)) )); + CALL_SUBTEST_8(( jacobisvd(MatrixXcd(14,7)) )); } - CALL_SUBTEST_9(( svd(MatrixXf(300,200)) )); - CALL_SUBTEST_10(( svd(MatrixXcd(100,150)) )); + CALL_SUBTEST_9(( jacobisvd(MatrixXf(300,200)) )); + CALL_SUBTEST_10(( jacobisvd(MatrixXcd(100,150)) )); - CALL_SUBTEST_3(( svd_verify_assert() )); - CALL_SUBTEST_3(( svd_verify_assert() )); - CALL_SUBTEST_9(( svd_verify_assert() )); - CALL_SUBTEST_11(( svd_verify_assert() )); + CALL_SUBTEST_3(( jacobisvd_verify_assert() )); + CALL_SUBTEST_3(( jacobisvd_verify_assert() )); + CALL_SUBTEST_9(( jacobisvd_verify_assert() )); + CALL_SUBTEST_11(( jacobisvd_verify_assert() )); // Test problem size constructors CALL_SUBTEST_12( JacobiSVD(10, 20) ); -- cgit v1.2.3