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>
benchEigenSolver(const MatrixType & m)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
main(int argc,char * argv[])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