Home | History | Annotate | Download | only in test
      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 (at) inria.fr>
      5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1 (at) 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 }
    199