source: pacpussensors/trunk/Vislab/lib3dv-1.2.0/lib3dv/eigen/bench/benchEigenSolver.cpp

Last change on this file was 136, checked in by ldecherf, 8 years ago

Doc

File size: 5.7 KB
Line 
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
Note: See TracBrowser for help on using the repository browser.