Home | History | Annotate | Download | only in bench
      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