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
