Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 1 | // This file is part of Eigen, a lightweight C++ template library |
| 2 | // for linear algebra. |
| 3 | // |
| 4 | // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> |
| 5 | // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> |
| 6 | // |
| 7 | // This Source Code Form is subject to the terms of the Mozilla |
| 8 | // Public License v. 2.0. If a copy of the MPL was not distributed |
| 9 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| 10 | |
| 11 | #include "svd_common.h" |
| 12 | |
| 13 | template<typename MatrixType, int QRPreconditioner> |
| 14 | void jacobisvd_check_full(const MatrixType& m, const JacobiSVD<MatrixType, QRPreconditioner>& svd) |
| 15 | { |
| 16 | svd_check_full<MatrixType, JacobiSVD<MatrixType, QRPreconditioner > >(m, svd); |
| 17 | } |
| 18 | |
| 19 | template<typename MatrixType, int QRPreconditioner> |
| 20 | void jacobisvd_compare_to_full(const MatrixType& m, |
| 21 | unsigned int computationOptions, |
| 22 | const JacobiSVD<MatrixType, QRPreconditioner>& referenceSvd) |
| 23 | { |
| 24 | svd_compare_to_full<MatrixType, JacobiSVD<MatrixType, QRPreconditioner> >(m, computationOptions, referenceSvd); |
| 25 | } |
| 26 | |
| 27 | |
| 28 | template<typename MatrixType, int QRPreconditioner> |
| 29 | void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions) |
| 30 | { |
| 31 | svd_solve< MatrixType, JacobiSVD< MatrixType, QRPreconditioner > >(m, computationOptions); |
| 32 | } |
| 33 | |
| 34 | |
| 35 | |
| 36 | template<typename MatrixType, int QRPreconditioner> |
| 37 | void jacobisvd_test_all_computation_options(const MatrixType& m) |
| 38 | { |
| 39 | |
| 40 | if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) |
| 41 | return; |
| 42 | |
| 43 | JacobiSVD< MatrixType, QRPreconditioner > fullSvd(m, ComputeFullU|ComputeFullV); |
| 44 | svd_test_computation_options_1< MatrixType, JacobiSVD< MatrixType, QRPreconditioner > >(m, fullSvd); |
| 45 | |
| 46 | if(QRPreconditioner == FullPivHouseholderQRPreconditioner) |
| 47 | return; |
| 48 | svd_test_computation_options_2< MatrixType, JacobiSVD< MatrixType, QRPreconditioner > >(m, fullSvd); |
| 49 | |
| 50 | } |
| 51 | |
| 52 | template<typename MatrixType> |
| 53 | void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) |
| 54 | { |
| 55 | MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a; |
| 56 | |
| 57 | jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m); |
| 58 | jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m); |
| 59 | jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m); |
| 60 | jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m); |
| 61 | } |
| 62 | |
| 63 | |
| 64 | template<typename MatrixType> |
| 65 | void jacobisvd_verify_assert(const MatrixType& m) |
| 66 | { |
| 67 | |
| 68 | svd_verify_assert<MatrixType, JacobiSVD< MatrixType > >(m); |
| 69 | |
| 70 | typedef typename MatrixType::Index Index; |
| 71 | Index rows = m.rows(); |
| 72 | Index cols = m.cols(); |
| 73 | |
| 74 | enum { |
| 75 | RowsAtCompileTime = MatrixType::RowsAtCompileTime, |
| 76 | ColsAtCompileTime = MatrixType::ColsAtCompileTime |
| 77 | }; |
| 78 | |
| 79 | MatrixType a = MatrixType::Zero(rows, cols); |
| 80 | a.setZero(); |
| 81 | |
| 82 | if (ColsAtCompileTime == Dynamic) |
| 83 | { |
| 84 | JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> svd_fullqr; |
| 85 | VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeFullU|ComputeThinV)) |
| 86 | VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeThinV)) |
| 87 | VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeFullV)) |
| 88 | } |
| 89 | } |
| 90 | |
| 91 | template<typename MatrixType> |
| 92 | void jacobisvd_method() |
| 93 | { |
| 94 | enum { Size = MatrixType::RowsAtCompileTime }; |
| 95 | typedef typename MatrixType::RealScalar RealScalar; |
| 96 | typedef Matrix<RealScalar, Size, 1> RealVecType; |
| 97 | MatrixType m = MatrixType::Identity(); |
| 98 | VERIFY_IS_APPROX(m.jacobiSvd().singularValues(), RealVecType::Ones()); |
| 99 | VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixU()); |
| 100 | VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixV()); |
| 101 | VERIFY_IS_APPROX(m.jacobiSvd(ComputeFullU|ComputeFullV).solve(m), m); |
| 102 | } |
| 103 | |
| 104 | |
| 105 | |
| 106 | template<typename MatrixType> |
| 107 | void jacobisvd_inf_nan() |
| 108 | { |
| 109 | svd_inf_nan<MatrixType, JacobiSVD< MatrixType > >(); |
| 110 | } |
| 111 | |
| 112 | |
| 113 | // Regression test for bug 286: JacobiSVD loops indefinitely with some |
| 114 | // matrices containing denormal numbers. |
| 115 | void jacobisvd_bug286() |
| 116 | { |
| 117 | #if defined __INTEL_COMPILER |
| 118 | // shut up warning #239: floating point underflow |
| 119 | #pragma warning push |
| 120 | #pragma warning disable 239 |
| 121 | #endif |
| 122 | Matrix2d M; |
| 123 | M << -7.90884e-313, -4.94e-324, |
| 124 | 0, 5.60844e-313; |
| 125 | #if defined __INTEL_COMPILER |
| 126 | #pragma warning pop |
| 127 | #endif |
| 128 | JacobiSVD<Matrix2d> svd; |
| 129 | svd.compute(M); // just check we don't loop indefinitely |
| 130 | } |
| 131 | |
| 132 | |
| 133 | void jacobisvd_preallocate() |
| 134 | { |
| 135 | svd_preallocate< JacobiSVD <MatrixXf> >(); |
| 136 | } |
| 137 | |
| 138 | void test_jacobisvd() |
| 139 | { |
| 140 | CALL_SUBTEST_11(( jacobisvd<Matrix<double,Dynamic,Dynamic> > |
| 141 | (Matrix<double,Dynamic,Dynamic>(16, 6)) )); |
| 142 | |
| 143 | CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) )); |
| 144 | CALL_SUBTEST_4(( jacobisvd_verify_assert(Matrix4d()) )); |
| 145 | CALL_SUBTEST_7(( jacobisvd_verify_assert(MatrixXf(10,12)) )); |
| 146 | CALL_SUBTEST_8(( jacobisvd_verify_assert(MatrixXcd(7,5)) )); |
| 147 | |
| 148 | for(int i = 0; i < g_repeat; i++) { |
| 149 | Matrix2cd m; |
| 150 | m << 0, 1, |
| 151 | 0, 1; |
| 152 | CALL_SUBTEST_1(( jacobisvd(m, false) )); |
| 153 | m << 1, 0, |
| 154 | 1, 0; |
| 155 | CALL_SUBTEST_1(( jacobisvd(m, false) )); |
| 156 | |
| 157 | Matrix2d n; |
| 158 | n << 0, 0, |
| 159 | 0, 0; |
| 160 | CALL_SUBTEST_2(( jacobisvd(n, false) )); |
| 161 | n << 0, 0, |
| 162 | 0, 1; |
| 163 | CALL_SUBTEST_2(( jacobisvd(n, false) )); |
| 164 | |
| 165 | CALL_SUBTEST_3(( jacobisvd<Matrix3f>() )); |
| 166 | CALL_SUBTEST_4(( jacobisvd<Matrix4d>() )); |
| 167 | CALL_SUBTEST_5(( jacobisvd<Matrix<float,3,5> >() )); |
| 168 | CALL_SUBTEST_6(( jacobisvd<Matrix<double,Dynamic,2> >(Matrix<double,Dynamic,2>(10,2)) )); |
| 169 | |
| 170 | int r = internal::random<int>(1, 30), |
| 171 | c = internal::random<int>(1, 30); |
| 172 | CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) )); |
| 173 | CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) )); |
| 174 | (void) r; |
| 175 | (void) c; |
| 176 | |
| 177 | // Test on inf/nan matrix |
| 178 | CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() ); |
| 179 | } |
| 180 | |
| 181 | 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))) )); |
| 182 | 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))) )); |
| 183 | |
| 184 | |
| 185 | // test matrixbase method |
| 186 | CALL_SUBTEST_1(( jacobisvd_method<Matrix2cd>() )); |
| 187 | CALL_SUBTEST_3(( jacobisvd_method<Matrix3f>() )); |
| 188 | |
| 189 | |
| 190 | // Test problem size constructors |
| 191 | CALL_SUBTEST_7( JacobiSVD<MatrixXf>(10,10) ); |
| 192 | |
| 193 | // Check that preallocation avoids subsequent mallocs |
| 194 | CALL_SUBTEST_9( jacobisvd_preallocate() ); |
| 195 | |
| 196 | // Regression check for bug 286 |
| 197 | CALL_SUBTEST_2( jacobisvd_bug286() ); |
| 198 | } |