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>
17using 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
31typedef SCALAR Scalar;
32
33template <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
189int 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