1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com>
5 // Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr>
6 // Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr>
7 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr>
8 //
9 // This Source Code Form is subject to the terms of the Mozilla
10 // Public License v. 2.0. If a copy of the MPL was not distributed
11 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/
12 
13 // Bench to compare the efficiency of SVD algorithms
14 
15 #include <iostream>
16 #include <bench/BenchTimer.h>
17 #include <unsupported/Eigen/SVD>
18 
19 
20 using namespace Eigen;
21 using namespace std;
22 
23 // number of computations of each algorithm before the print of the time
24 #ifndef REPEAT
25 #define REPEAT 10
26 #endif
27 
28 // number of tests of the same type
29 #ifndef NUMBER_SAMPLE
30 #define NUMBER_SAMPLE 2
31 #endif
32 
33 template<typename MatrixType>
34 void bench_svd(const MatrixType& a = MatrixType())
35 {
36   MatrixType m = MatrixType::Random(a.rows(), a.cols());
37   BenchTimer timerJacobi;
38   BenchTimer timerBDC;
39   timerJacobi.reset();
40   timerBDC.reset();
41 
42   cout << " Only compute Singular Values" <<endl;
43   for (int k=1; k<=NUMBER_SAMPLE; ++k)
44   {
45     timerBDC.start();
46     for (int i=0; i<REPEAT; ++i)
47     {
48       BDCSVD<MatrixType> bdc_matrix(m);
49     }
50     timerBDC.stop();
51 
52     timerJacobi.start();
53     for (int i=0; i<REPEAT; ++i)
54     {
55       JacobiSVD<MatrixType> jacobi_matrix(m);
56     }
57     timerJacobi.stop();
58 
59 
60     cout << "Sample " << k << " : " << REPEAT << " computations :  Jacobi : " << fixed << timerJacobi.value() << "s ";
61     cout << " || " << " BDC : " << timerBDC.value() << "s " <<endl <<endl;
62 
63     if (timerBDC.value() >= timerJacobi.value())
64       cout << "KO : BDC is " <<  timerJacobi.value() / timerBDC.value() << "  times faster than Jacobi" <<endl;
65     else
66       cout << "OK : BDC is " << timerJacobi.value() / timerBDC.value() << "  times faster than Jacobi"  <<endl;
67 
68   }
69   cout << "       =================" <<endl;
70   std::cout<< std::endl;
71   timerJacobi.reset();
72   timerBDC.reset();
73   cout << " Computes rotaion matrix" <<endl;
74   for (int k=1; k<=NUMBER_SAMPLE; ++k)
75   {
76     timerBDC.start();
77     for (int i=0; i<REPEAT; ++i)
78     {
79       BDCSVD<MatrixType> bdc_matrix(m, ComputeFullU|ComputeFullV);
80     }
81     timerBDC.stop();
82 
83     timerJacobi.start();
84     for (int i=0; i<REPEAT; ++i)
85     {
86       JacobiSVD<MatrixType> jacobi_matrix(m, ComputeFullU|ComputeFullV);
87     }
88     timerJacobi.stop();
89 
90 
91     cout << "Sample " << k << " : " << REPEAT << " computations :  Jacobi : " << fixed << timerJacobi.value() << "s ";
92     cout << " || " << " BDC : " << timerBDC.value() << "s " <<endl <<endl;
93 
94     if (timerBDC.value() >= timerJacobi.value())
95       cout << "KO : BDC is " <<  timerJacobi.value() / timerBDC.value() << "  times faster than Jacobi" <<endl;
96     else
97       cout << "OK : BDC is " << timerJacobi.value() / timerBDC.value() << "  times faster than Jacobi"  <<endl;
98 
99   }
100   std::cout<< std::endl;
101 }
102 
103 
104 
105 int main(int argc, char* argv[])
106 {
107   std::cout<< std::endl;
108 
109   std::cout<<"On a (Dynamic, Dynamic) (6, 6) Matrix" <<std::endl;
110   bench_svd<Matrix<double,Dynamic,Dynamic> >(Matrix<double,Dynamic,Dynamic>(6, 6));
111 
112   std::cout<<"On a (Dynamic, Dynamic) (32, 32) Matrix" <<std::endl;
113   bench_svd<Matrix<double,Dynamic,Dynamic> >(Matrix<double,Dynamic,Dynamic>(32, 32));
114 
115   //std::cout<<"On a (Dynamic, Dynamic) (128, 128) Matrix" <<std::endl;
116   //bench_svd<Matrix<double,Dynamic,Dynamic> >(Matrix<double,Dynamic,Dynamic>(128, 128));
117 
118   std::cout<<"On a (Dynamic, Dynamic) (160, 160) Matrix" <<std::endl;
119   bench_svd<Matrix<double,Dynamic,Dynamic> >(Matrix<double,Dynamic,Dynamic>(160, 160));
120 
121   std::cout<< "--------------------------------------------------------------------"<< std::endl;
122 
123 }
124