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 // discard stack allocation as that too bypasses malloc
     12 #define EIGEN_STACK_ALLOCATION_LIMIT 0
     13 #define EIGEN_RUNTIME_NO_MALLOC
     14 #include "main.h"
     15 #include <Eigen/SVD>
     16 
     17 template<typename MatrixType, int QRPreconditioner>
     18 void jacobisvd_check_full(const MatrixType& m, const JacobiSVD<MatrixType, QRPreconditioner>& svd)
     19 {
     20   typedef typename MatrixType::Index Index;
     21   Index rows = m.rows();
     22   Index cols = m.cols();
     23 
     24   enum {
     25     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
     26     ColsAtCompileTime = MatrixType::ColsAtCompileTime
     27   };
     28 
     29   typedef typename MatrixType::Scalar Scalar;
     30   typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType;
     31   typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType;
     32 
     33   MatrixType sigma = MatrixType::Zero(rows,cols);
     34   sigma.diagonal() = svd.singularValues().template cast<Scalar>();
     35   MatrixUType u = svd.matrixU();
     36   MatrixVType v = svd.matrixV();
     37 
     38   VERIFY_IS_APPROX(m, u * sigma * v.adjoint());
     39   VERIFY_IS_UNITARY(u);
     40   VERIFY_IS_UNITARY(v);
     41 }
     42 
     43 template<typename MatrixType, int QRPreconditioner>
     44 void jacobisvd_compare_to_full(const MatrixType& m,
     45                                unsigned int computationOptions,
     46                                const JacobiSVD<MatrixType, QRPreconditioner>& referenceSvd)
     47 {
     48   typedef typename MatrixType::Index Index;
     49   Index rows = m.rows();
     50   Index cols = m.cols();
     51   Index diagSize = (std::min)(rows, cols);
     52 
     53   JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions);
     54 
     55   VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues());
     56   if(computationOptions & ComputeFullU)
     57     VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU());
     58   if(computationOptions & ComputeThinU)
     59     VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize));
     60   if(computationOptions & ComputeFullV)
     61     VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV());
     62   if(computationOptions & ComputeThinV)
     63     VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize));
     64 }
     65 
     66 template<typename MatrixType, int QRPreconditioner>
     67 void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions)
     68 {
     69   typedef typename MatrixType::Scalar Scalar;
     70   typedef typename MatrixType::RealScalar RealScalar;
     71   typedef typename MatrixType::Index Index;
     72   Index rows = m.rows();
     73   Index cols = m.cols();
     74 
     75   enum {
     76     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
     77     ColsAtCompileTime = MatrixType::ColsAtCompileTime
     78   };
     79 
     80   typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType;
     81   typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType;
     82 
     83   RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols));
     84   JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions);
     85 
     86        if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8);
     87   else if(internal::is_same<RealScalar,float>::value)  svd.setThreshold(1e-4);
     88 
     89   SolutionType x = svd.solve(rhs);
     90 
     91   RealScalar residual = (m*x-rhs).norm();
     92   // Check that there is no significantly better solution in the neighborhood of x
     93   if(!test_isMuchSmallerThan(residual,rhs.norm()))
     94   {
     95     // If the residual is very small, then we have an exact solution, so we are already good.
     96     for(int k=0;k<x.rows();++k)
     97     {
     98       SolutionType y(x);
     99       y.row(k).array() += 2*NumTraits<RealScalar>::epsilon();
    100       RealScalar residual_y = (m*y-rhs).norm();
    101       VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
    102 
    103       y.row(k) = x.row(k).array() - 2*NumTraits<RealScalar>::epsilon();
    104       residual_y = (m*y-rhs).norm();
    105       VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
    106     }
    107   }
    108 
    109   // evaluate normal equation which works also for least-squares solutions
    110   if(internal::is_same<RealScalar,double>::value)
    111   {
    112     // This test is not stable with single precision.
    113     // This is probably because squaring m signicantly affects the precision.
    114     VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs);
    115   }
    116 
    117   // check minimal norm solutions
    118   {
    119     // generate a full-rank m x n problem with m<n
    120     enum {
    121       RankAtCompileTime2 = ColsAtCompileTime==Dynamic ? Dynamic : (ColsAtCompileTime)/2+1,
    122       RowsAtCompileTime3 = ColsAtCompileTime==Dynamic ? Dynamic : ColsAtCompileTime+1
    123     };
    124     typedef Matrix<Scalar, RankAtCompileTime2, ColsAtCompileTime> MatrixType2;
    125     typedef Matrix<Scalar, RankAtCompileTime2, 1> RhsType2;
    126     typedef Matrix<Scalar, ColsAtCompileTime, RankAtCompileTime2> MatrixType2T;
    127     Index rank = RankAtCompileTime2==Dynamic ? internal::random<Index>(1,cols) : Index(RankAtCompileTime2);
    128     MatrixType2 m2(rank,cols);
    129     int guard = 0;
    130     do {
    131       m2.setRandom();
    132     } while(m2.jacobiSvd().setThreshold(test_precision<Scalar>()).rank()!=rank && (++guard)<10);
    133     VERIFY(guard<10);
    134     RhsType2 rhs2 = RhsType2::Random(rank);
    135     // use QR to find a reference minimal norm solution
    136     HouseholderQR<MatrixType2T> qr(m2.adjoint());
    137     Matrix<Scalar,Dynamic,1> tmp = qr.matrixQR().topLeftCorner(rank,rank).template triangularView<Upper>().adjoint().solve(rhs2);
    138     tmp.conservativeResize(cols);
    139     tmp.tail(cols-rank).setZero();
    140     SolutionType x21 = qr.householderQ() * tmp;
    141     // now check with SVD
    142     JacobiSVD<MatrixType2, ColPivHouseholderQRPreconditioner> svd2(m2, computationOptions);
    143     SolutionType x22 = svd2.solve(rhs2);
    144     VERIFY_IS_APPROX(m2*x21, rhs2);
    145     VERIFY_IS_APPROX(m2*x22, rhs2);
    146     VERIFY_IS_APPROX(x21, x22);
    147 
    148     // Now check with a rank deficient matrix
    149     typedef Matrix<Scalar, RowsAtCompileTime3, ColsAtCompileTime> MatrixType3;
    150     typedef Matrix<Scalar, RowsAtCompileTime3, 1> RhsType3;
    151     Index rows3 = RowsAtCompileTime3==Dynamic ? internal::random<Index>(rank+1,2*cols) : Index(RowsAtCompileTime3);
    152     Matrix<Scalar,RowsAtCompileTime3,Dynamic> C = Matrix<Scalar,RowsAtCompileTime3,Dynamic>::Random(rows3,rank);
    153     MatrixType3 m3 = C * m2;
    154     RhsType3 rhs3 = C * rhs2;
    155     JacobiSVD<MatrixType3, ColPivHouseholderQRPreconditioner> svd3(m3, computationOptions);
    156     SolutionType x3 = svd3.solve(rhs3);
    157     if(svd3.rank()!=rank) {
    158       std::cout << m3 << "\n\n";
    159       std::cout << svd3.singularValues().transpose() << "\n";
    160     std::cout << svd3.rank() << " == " << rank << "\n";
    161     std::cout << x21.norm() << " == " << x3.norm() << "\n";
    162     }
    163 //     VERIFY_IS_APPROX(m3*x3, rhs3);
    164     VERIFY_IS_APPROX(m3*x21, rhs3);
    165     VERIFY_IS_APPROX(m2*x3, rhs2);
    166 
    167     VERIFY_IS_APPROX(x21, x3);
    168   }
    169 }
    170 
    171 template<typename MatrixType, int QRPreconditioner>
    172 void jacobisvd_test_all_computation_options(const MatrixType& m)
    173 {
    174   if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols())
    175     return;
    176   JacobiSVD<MatrixType, QRPreconditioner> fullSvd(m, ComputeFullU|ComputeFullV);
    177   CALL_SUBTEST(( jacobisvd_check_full(m, fullSvd) ));
    178   CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV) ));
    179 
    180   #if defined __INTEL_COMPILER
    181   // remark #111: statement is unreachable
    182   #pragma warning disable 111
    183   #endif
    184   if(QRPreconditioner == FullPivHouseholderQRPreconditioner)
    185     return;
    186 
    187   CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU, fullSvd) ));
    188   CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullV, fullSvd) ));
    189   CALL_SUBTEST(( jacobisvd_compare_to_full(m, 0, fullSvd) ));
    190 
    191   if (MatrixType::ColsAtCompileTime == Dynamic) {
    192     // thin U/V are only available with dynamic number of columns
    193     CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) ));
    194     CALL_SUBTEST(( jacobisvd_compare_to_full(m,              ComputeThinV, fullSvd) ));
    195     CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) ));
    196     CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU             , fullSvd) ));
    197     CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) ));
    198     CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV) ));
    199     CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV) ));
    200     CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV) ));
    201 
    202     // test reconstruction
    203     typedef typename MatrixType::Index Index;
    204     Index diagSize = (std::min)(m.rows(), m.cols());
    205     JacobiSVD<MatrixType, QRPreconditioner> svd(m, ComputeThinU | ComputeThinV);
    206     VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint());
    207   }
    208 }
    209 
    210 template<typename MatrixType>
    211 void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true)
    212 {
    213   MatrixType m = a;
    214   if(pickrandom)
    215   {
    216     typedef typename MatrixType::Scalar Scalar;
    217     typedef typename MatrixType::RealScalar RealScalar;
    218     typedef typename MatrixType::Index Index;
    219     Index diagSize = (std::min)(a.rows(), a.cols());
    220     RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4;
    221     s = internal::random<RealScalar>(1,s);
    222     Matrix<RealScalar,Dynamic,1> d =  Matrix<RealScalar,Dynamic,1>::Random(diagSize);
    223     for(Index k=0; k<diagSize; ++k)
    224       d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s));
    225     m = Matrix<Scalar,Dynamic,Dynamic>::Random(a.rows(),diagSize) * d.asDiagonal() * Matrix<Scalar,Dynamic,Dynamic>::Random(diagSize,a.cols());
    226     // cancel some coeffs
    227     Index n  = internal::random<Index>(0,m.size()-1);
    228     for(Index i=0; i<n; ++i)
    229       m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0);
    230   }
    231 
    232   CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m) ));
    233   CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m) ));
    234   CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m) ));
    235   CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m) ));
    236 }
    237 
    238 template<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m)
    239 {
    240   typedef typename MatrixType::Scalar Scalar;
    241   typedef typename MatrixType::Index Index;
    242   Index rows = m.rows();
    243   Index cols = m.cols();
    244 
    245   enum {
    246     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
    247     ColsAtCompileTime = MatrixType::ColsAtCompileTime
    248   };
    249 
    250   typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType;
    251 
    252   RhsType rhs(rows);
    253 
    254   JacobiSVD<MatrixType> svd;
    255   VERIFY_RAISES_ASSERT(svd.matrixU())
    256   VERIFY_RAISES_ASSERT(svd.singularValues())
    257   VERIFY_RAISES_ASSERT(svd.matrixV())
    258   VERIFY_RAISES_ASSERT(svd.solve(rhs))
    259 
    260   MatrixType a = MatrixType::Zero(rows, cols);
    261   a.setZero();
    262   svd.compute(a, 0);
    263   VERIFY_RAISES_ASSERT(svd.matrixU())
    264   VERIFY_RAISES_ASSERT(svd.matrixV())
    265   svd.singularValues();
    266   VERIFY_RAISES_ASSERT(svd.solve(rhs))
    267 
    268   if (ColsAtCompileTime == Dynamic)
    269   {
    270     svd.compute(a, ComputeThinU);
    271     svd.matrixU();
    272     VERIFY_RAISES_ASSERT(svd.matrixV())
    273     VERIFY_RAISES_ASSERT(svd.solve(rhs))
    274 
    275     svd.compute(a, ComputeThinV);
    276     svd.matrixV();
    277     VERIFY_RAISES_ASSERT(svd.matrixU())
    278     VERIFY_RAISES_ASSERT(svd.solve(rhs))
    279 
    280     JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> svd_fullqr;
    281     VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeFullU|ComputeThinV))
    282     VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeThinV))
    283     VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeFullV))
    284   }
    285   else
    286   {
    287     VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU))
    288     VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV))
    289   }
    290 }
    291 
    292 template<typename MatrixType>
    293 void jacobisvd_method()
    294 {
    295   enum { Size = MatrixType::RowsAtCompileTime };
    296   typedef typename MatrixType::RealScalar RealScalar;
    297   typedef Matrix<RealScalar, Size, 1> RealVecType;
    298   MatrixType m = MatrixType::Identity();
    299   VERIFY_IS_APPROX(m.jacobiSvd().singularValues(), RealVecType::Ones());
    300   VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixU());
    301   VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixV());
    302   VERIFY_IS_APPROX(m.jacobiSvd(ComputeFullU|ComputeFullV).solve(m), m);
    303 }
    304 
    305 // work around stupid msvc error when constructing at compile time an expression that involves
    306 // a division by zero, even if the numeric type has floating point
    307 template<typename Scalar>
    308 EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); }
    309 
    310 // workaround aggressive optimization in ICC
    311 template<typename T> EIGEN_DONT_INLINE  T sub(T a, T b) { return a - b; }
    312 
    313 template<typename MatrixType>
    314 void jacobisvd_inf_nan()
    315 {
    316   // all this function does is verify we don't iterate infinitely on nan/inf values
    317 
    318   JacobiSVD<MatrixType> svd;
    319   typedef typename MatrixType::Scalar Scalar;
    320   Scalar some_inf = Scalar(1) / zero<Scalar>();
    321   VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf));
    322   svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV);
    323 
    324   Scalar some_nan = zero<Scalar>() / zero<Scalar>();
    325   VERIFY(some_nan != some_nan);
    326   svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV);
    327 
    328   MatrixType m = MatrixType::Zero(10,10);
    329   m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf;
    330   svd.compute(m, ComputeFullU | ComputeFullV);
    331 
    332   m = MatrixType::Zero(10,10);
    333   m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan;
    334   svd.compute(m, ComputeFullU | ComputeFullV);
    335 }
    336 
    337 // Regression test for bug 286: JacobiSVD loops indefinitely with some
    338 // matrices containing denormal numbers.
    339 void jacobisvd_bug286()
    340 {
    341 #if defined __INTEL_COMPILER
    342 // shut up warning #239: floating point underflow
    343 #pragma warning push
    344 #pragma warning disable 239
    345 #endif
    346   Matrix2d M;
    347   M << -7.90884e-313, -4.94e-324,
    348                  0, 5.60844e-313;
    349 #if defined __INTEL_COMPILER
    350 #pragma warning pop
    351 #endif
    352   JacobiSVD<Matrix2d> svd;
    353   svd.compute(M); // just check we don't loop indefinitely
    354 }
    355 
    356 void jacobisvd_preallocate()
    357 {
    358   Vector3f v(3.f, 2.f, 1.f);
    359   MatrixXf m = v.asDiagonal();
    360 
    361   internal::set_is_malloc_allowed(false);
    362   VERIFY_RAISES_ASSERT(VectorXf tmp(10);)
    363   JacobiSVD<MatrixXf> svd;
    364   internal::set_is_malloc_allowed(true);
    365   svd.compute(m);
    366   VERIFY_IS_APPROX(svd.singularValues(), v);
    367 
    368   JacobiSVD<MatrixXf> svd2(3,3);
    369   internal::set_is_malloc_allowed(false);
    370   svd2.compute(m);
    371   internal::set_is_malloc_allowed(true);
    372   VERIFY_IS_APPROX(svd2.singularValues(), v);
    373   VERIFY_RAISES_ASSERT(svd2.matrixU());
    374   VERIFY_RAISES_ASSERT(svd2.matrixV());
    375   svd2.compute(m, ComputeFullU | ComputeFullV);
    376   VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
    377   VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
    378   internal::set_is_malloc_allowed(false);
    379   svd2.compute(m);
    380   internal::set_is_malloc_allowed(true);
    381 
    382   JacobiSVD<MatrixXf> svd3(3,3,ComputeFullU|ComputeFullV);
    383   internal::set_is_malloc_allowed(false);
    384   svd2.compute(m);
    385   internal::set_is_malloc_allowed(true);
    386   VERIFY_IS_APPROX(svd2.singularValues(), v);
    387   VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
    388   VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
    389   internal::set_is_malloc_allowed(false);
    390   svd2.compute(m, ComputeFullU|ComputeFullV);
    391   internal::set_is_malloc_allowed(true);
    392 }
    393 
    394 void test_jacobisvd()
    395 {
    396   CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) ));
    397   CALL_SUBTEST_4(( jacobisvd_verify_assert(Matrix4d()) ));
    398   CALL_SUBTEST_7(( jacobisvd_verify_assert(MatrixXf(10,12)) ));
    399   CALL_SUBTEST_8(( jacobisvd_verify_assert(MatrixXcd(7,5)) ));
    400 
    401   for(int i = 0; i < g_repeat; i++) {
    402     Matrix2cd m;
    403     m << 0, 1,
    404          0, 1;
    405     CALL_SUBTEST_1(( jacobisvd(m, false) ));
    406     m << 1, 0,
    407          1, 0;
    408     CALL_SUBTEST_1(( jacobisvd(m, false) ));
    409 
    410     Matrix2d n;
    411     n << 0, 0,
    412          0, 0;
    413     CALL_SUBTEST_2(( jacobisvd(n, false) ));
    414     n << 0, 0,
    415          0, 1;
    416     CALL_SUBTEST_2(( jacobisvd(n, false) ));
    417 
    418     CALL_SUBTEST_3(( jacobisvd<Matrix3f>() ));
    419     CALL_SUBTEST_4(( jacobisvd<Matrix4d>() ));
    420     CALL_SUBTEST_5(( jacobisvd<Matrix<float,3,5> >() ));
    421     CALL_SUBTEST_6(( jacobisvd<Matrix<double,Dynamic,2> >(Matrix<double,Dynamic,2>(10,2)) ));
    422 
    423     int r = internal::random<int>(1, 30),
    424         c = internal::random<int>(1, 30);
    425 
    426     TEST_SET_BUT_UNUSED_VARIABLE(r)
    427     TEST_SET_BUT_UNUSED_VARIABLE(c)
    428 
    429     CALL_SUBTEST_10(( jacobisvd<MatrixXd>(MatrixXd(r,c)) ));
    430     CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) ));
    431     CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) ));
    432     (void) r;
    433     (void) c;
    434 
    435     // Test on inf/nan matrix
    436     CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() );
    437   }
    438 
    439   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))) ));
    440   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))) ));
    441 
    442   // test matrixbase method
    443   CALL_SUBTEST_1(( jacobisvd_method<Matrix2cd>() ));
    444   CALL_SUBTEST_3(( jacobisvd_method<Matrix3f>() ));
    445 
    446   // Test problem size constructors
    447   CALL_SUBTEST_7( JacobiSVD<MatrixXf>(10,10) );
    448 
    449   // Check that preallocation avoids subsequent mallocs
    450   CALL_SUBTEST_9( jacobisvd_preallocate() );
    451 
    452   // Regression check for bug 286
    453   CALL_SUBTEST_2( jacobisvd_bug286() );
    454 }
    455