1 
2 // g++ -DNDEBUG -O3 -I.. benchEigenSolver.cpp  -o benchEigenSolver && ./benchEigenSolver
3 // options:
4 //  -DBENCH_GMM
5 //  -DBENCH_GSL -lgsl /usr/lib/libcblas.so.3
6 //  -DEIGEN_DONT_VECTORIZE
7 //  -msse2
8 //  -DREPEAT=100
9 //  -DTRIES=10
10 //  -DSCALAR=double
11 
12 #include <iostream>
13 
14 #include <Eigen/Core>
15 #include <Eigen/QR>
16 #include <bench/BenchUtil.h>
17 using namespace Eigen;
18 
19 #ifndef REPEAT
20 #define REPEAT 1000
21 #endif
22 
23 #ifndef TRIES
24 #define TRIES 4
25 #endif
26 
27 #ifndef SCALAR
28 #define SCALAR float
29 #endif
30 
31 typedef SCALAR Scalar;
32 
33 template <typename MatrixType>
34 __attribute__ ((noinline)) void benchEigenSolver(const MatrixType& m)
35 {
36   int rows = m.rows();
37   int cols = m.cols();
38 
39   int stdRepeats = std::max(1,int((REPEAT*1000)/(rows*rows*sqrt(rows))));
40   int saRepeats = stdRepeats * 4;
41 
42   typedef typename MatrixType::Scalar Scalar;
43   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
44 
45   MatrixType a = MatrixType::Random(rows,cols);
46   SquareMatrixType covMat =  a * a.adjoint();
47 
48   BenchTimer timerSa, timerStd;
49 
50   Scalar acc = 0;
51   int r = internal::random<int>(0,covMat.rows()-1);
52   int c = internal::random<int>(0,covMat.cols()-1);
53   {
54     SelfAdjointEigenSolver<SquareMatrixType> ei(covMat);
55     for (int t=0; t<TRIES; ++t)
56     {
57       timerSa.start();
58       for (int k=0; k<saRepeats; ++k)
59       {
60         ei.compute(covMat);
61         acc += ei.eigenvectors().coeff(r,c);
62       }
63       timerSa.stop();
64     }
65   }
66 
67   {
68     EigenSolver<SquareMatrixType> ei(covMat);
69     for (int t=0; t<TRIES; ++t)
70     {
71       timerStd.start();
72       for (int k=0; k<stdRepeats; ++k)
73       {
74         ei.compute(covMat);
75         acc += ei.eigenvectors().coeff(r,c);
76       }
77       timerStd.stop();
78     }
79   }
80 
81   if (MatrixType::RowsAtCompileTime==Dynamic)
82     std::cout << "dyn   ";
83   else
84     std::cout << "fixed ";
85   std::cout << covMat.rows() << " \t"
86             << timerSa.value() * REPEAT / saRepeats << "s \t"
87             << timerStd.value() * REPEAT / stdRepeats << "s";
88 
89   #ifdef BENCH_GMM
90   if (MatrixType::RowsAtCompileTime==Dynamic)
91   {
92     timerSa.reset();
93     timerStd.reset();
94 
95     gmm::dense_matrix<Scalar> gmmCovMat(covMat.rows(),covMat.cols());
96     gmm::dense_matrix<Scalar> eigvect(covMat.rows(),covMat.cols());
97     std::vector<Scalar> eigval(covMat.rows());
98     eiToGmm(covMat, gmmCovMat);
99     for (int t=0; t<TRIES; ++t)
100     {
101       timerSa.start();
102       for (int k=0; k<saRepeats; ++k)
103       {
104         gmm::symmetric_qr_algorithm(gmmCovMat, eigval, eigvect);
105         acc += eigvect(r,c);
106       }
107       timerSa.stop();
108     }
109     // the non-selfadjoint solver does not compute the eigen vectors
110 //     for (int t=0; t<TRIES; ++t)
111 //     {
112 //       timerStd.start();
113 //       for (int k=0; k<stdRepeats; ++k)
114 //       {
115 //         gmm::implicit_qr_algorithm(gmmCovMat, eigval, eigvect);
116 //         acc += eigvect(r,c);
117 //       }
118 //       timerStd.stop();
119 //     }
120 
121     std::cout << " | \t"
122               << timerSa.value() * REPEAT / saRepeats << "s"
123               << /*timerStd.value() * REPEAT / stdRepeats << "s"*/ "   na   ";
124   }
125   #endif
126 
127   #ifdef BENCH_GSL
128   if (MatrixType::RowsAtCompileTime==Dynamic)
129   {
130     timerSa.reset();
131     timerStd.reset();
132 
133     gsl_matrix* gslCovMat = gsl_matrix_alloc(covMat.rows(),covMat.cols());
134     gsl_matrix* gslCopy = gsl_matrix_alloc(covMat.rows(),covMat.cols());
135     gsl_matrix* eigvect = gsl_matrix_alloc(covMat.rows(),covMat.cols());
136     gsl_vector* eigval  = gsl_vector_alloc(covMat.rows());
137     gsl_eigen_symmv_workspace* eisymm = gsl_eigen_symmv_alloc(covMat.rows());
138 
139     gsl_matrix_complex* eigvectz = gsl_matrix_complex_alloc(covMat.rows(),covMat.cols());
140     gsl_vector_complex* eigvalz  = gsl_vector_complex_alloc(covMat.rows());
141     gsl_eigen_nonsymmv_workspace* einonsymm = gsl_eigen_nonsymmv_alloc(covMat.rows());
142 
143     eiToGsl(covMat, &gslCovMat);
144     for (int t=0; t<TRIES; ++t)
145     {
146       timerSa.start();
147       for (int k=0; k<saRepeats; ++k)
148       {
149         gsl_matrix_memcpy(gslCopy,gslCovMat);
150         gsl_eigen_symmv(gslCopy, eigval, eigvect, eisymm);
151         acc += gsl_matrix_get(eigvect,r,c);
152       }
153       timerSa.stop();
154     }
155     for (int t=0; t<TRIES; ++t)
156     {
157       timerStd.start();
158       for (int k=0; k<stdRepeats; ++k)
159       {
160         gsl_matrix_memcpy(gslCopy,gslCovMat);
161         gsl_eigen_nonsymmv(gslCopy, eigvalz, eigvectz, einonsymm);
162         acc += GSL_REAL(gsl_matrix_complex_get(eigvectz,r,c));
163       }
164       timerStd.stop();
165     }
166 
167     std::cout << " | \t"
168               << timerSa.value() * REPEAT / saRepeats << "s \t"
169               << timerStd.value() * REPEAT / stdRepeats << "s";
170 
171     gsl_matrix_free(gslCovMat);
172     gsl_vector_free(gslCopy);
173     gsl_matrix_free(eigvect);
174     gsl_vector_free(eigval);
175     gsl_matrix_complex_free(eigvectz);
176     gsl_vector_complex_free(eigvalz);
177     gsl_eigen_symmv_free(eisymm);
178     gsl_eigen_nonsymmv_free(einonsymm);
179   }
180   #endif
181 
182   std::cout << "\n";
183 
184   // make sure the compiler does not optimize too much
185   if (acc==123)
186     std::cout << acc;
187 }
188 
189 int main(int argc, char* argv[])
190 {
191   const int dynsizes[] = {4,6,8,12,16,24,32,64,128,256,512,0};
192   std::cout << "size            selfadjoint       generic";
193   #ifdef BENCH_GMM
194   std::cout << "        GMM++          ";
195   #endif
196   #ifdef BENCH_GSL
197   std::cout << "       GSL (double + ATLAS)  ";
198   #endif
199   std::cout << "\n";
200   for (uint i=0; dynsizes[i]>0; ++i)
201     benchEigenSolver(Matrix<Scalar,Dynamic,Dynamic>(dynsizes[i],dynsizes[i]));
202 
203   benchEigenSolver(Matrix<Scalar,2,2>());
204   benchEigenSolver(Matrix<Scalar,3,3>());
205   benchEigenSolver(Matrix<Scalar,4,4>());
206   benchEigenSolver(Matrix<Scalar,6,6>());
207   benchEigenSolver(Matrix<Scalar,8,8>());
208   benchEigenSolver(Matrix<Scalar,12,12>());
209   benchEigenSolver(Matrix<Scalar,16,16>());
210   return 0;
211 }
212 
213