benchEigenSolver.cpp
Go to the documentation of this file.
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;
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  {
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  {
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 
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Eigen::SelfAdjointEigenSolver::compute
EIGEN_DEVICE_FUNC SelfAdjointEigenSolver & compute(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
Scalar
SCALAR Scalar
Definition: benchEigenSolver.cpp:31
Eigen::EigenSolver::compute
EigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
Eigen::BenchTimer::value
double value(int TIMER=CPU_TIMER) const
Definition: BenchTimer.h:104
BenchUtil.h
Eigen::SelfAdjointEigenSolver::eigenvectors
const EIGEN_DEVICE_FUNC EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:277
main
int main(int argc, char *argv[])
Definition: benchEigenSolver.cpp:189
Eigen::SelfAdjointEigenSolver
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Definition: SelfAdjointEigenSolver.h:76
Eigen::BenchTimer::start
void start()
Definition: BenchTimer.h:81
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:22
Eigen::BenchTimer
Definition: BenchTimer.h:59
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
Eigen::EigenSolver::eigenvectors
EigenvectorsType eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: EigenSolver.h:345
Eigen::BenchTimer::reset
void reset()
Definition: BenchTimer.h:75
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
TRIES
#define TRIES
Definition: benchEigenSolver.cpp:24
REPEAT
#define REPEAT
Definition: benchEigenSolver.cpp:20
Eigen::PlainObjectBase::coeff
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:152
Eigen::EigenSolver
Computes eigenvalues and eigenvectors of general matrices.
Definition: EigenSolver.h:64
Eigen::BenchTimer::stop
void stop()
Definition: BenchTimer.h:86
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
eiToGmm
void eiToGmm(const EigenSparseMatrix &src, GmmSparse &dst)
Definition: BenchSparseUtil.h:74
__attribute__
__attribute__((noinline)) void benchEigenSolver(const MatrixType &m)
Definition: benchEigenSolver.cpp:34
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
align_3::t
Point2 t(10, 10)
max
#define max(a, b)
Definition: datatypes.h:20
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
SCALAR
#define SCALAR
Definition: benchEigenSolver.cpp:28


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:00