Home | History | Annotate | Download | only in bench 
      1 
      2 // g++ -DNDEBUG -O3 -I.. benchLLT.cpp  -o benchLLT && ./benchLLT
      3 // options:
      4 //  -DBENCH_GSL -lgsl /usr/lib/libcblas.so.3
      5 //  -DEIGEN_DONT_VECTORIZE
      6 //  -msse2
      7 //  -DREPEAT=100
      8 //  -DTRIES=10
      9 //  -DSCALAR=double
     10 
     11 #include <iostream>
     12 
     13 #include <Eigen/Core>
     14 #include <Eigen/Cholesky>
     15 #include <bench/BenchUtil.h>
     16 using namespace Eigen;
     17 
     18 #ifndef REPEAT
     19 #define REPEAT 10000
     20 #endif
     21 
     22 #ifndef TRIES
     23 #define TRIES 10
     24 #endif
     25 
     26 typedef float Scalar;
     27 
     28 template <typename MatrixType>
     29 __attribute__ ((noinline)) void benchLLT(const MatrixType& m)
     30 {
     31   int rows = m.rows();
     32   int cols = m.cols();
     33 
     34   double cost = 0;
     35   for (int j=0; j<rows; ++j)
     36   {
     37     int r = std::max(rows - j -1,0);
     38     cost += 2*(r*j+r+j);
     39   }
     40 
     41   int repeats = (REPEAT*1000)/(rows*rows);
     42 
     43   typedef typename MatrixType::Scalar Scalar;
     44   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
     45 
     46   MatrixType a = MatrixType::Random(rows,cols);
     47   SquareMatrixType covMat =  a * a.adjoint();
     48 
     49   BenchTimer timerNoSqrt, timerSqrt;
     50 
     51   Scalar acc = 0;
     52   int r = internal::random<int>(0,covMat.rows()-1);
     53   int c = internal::random<int>(0,covMat.cols()-1);
     54   for (int t=0; t<TRIES; ++t)
     55   {
     56     timerNoSqrt.start();
     57     for (int k=0; k<repeats; ++k)
     58     {
     59       LDLT<SquareMatrixType> cholnosqrt(covMat);
     60       acc += cholnosqrt.matrixL().coeff(r,c);
     61     }
     62     timerNoSqrt.stop();
     63   }
     64 
     65   for (int t=0; t<TRIES; ++t)
     66   {
     67     timerSqrt.start();
     68     for (int k=0; k<repeats; ++k)
     69     {
     70       LLT<SquareMatrixType> chol(covMat);
     71       acc += chol.matrixL().coeff(r,c);
     72     }
     73     timerSqrt.stop();
     74   }
     75 
     76   if (MatrixType::RowsAtCompileTime==Dynamic)
     77     std::cout << "dyn   ";
     78   else
     79     std::cout << "fixed ";
     80   std::cout << covMat.rows() << " \t"
     81             << (timerNoSqrt.best()) / repeats << "s "
     82             << "(" << 1e-9 * cost*repeats/timerNoSqrt.best() << " GFLOPS)\t"
     83             << (timerSqrt.best()) / repeats << "s "
     84             << "(" << 1e-9 * cost*repeats/timerSqrt.best() << " GFLOPS)\n";
     85 
     86 
     87   #ifdef BENCH_GSL
     88   if (MatrixType::RowsAtCompileTime==Dynamic)
     89   {
     90     timerSqrt.reset();
     91 
     92     gsl_matrix* gslCovMat = gsl_matrix_alloc(covMat.rows(),covMat.cols());
     93     gsl_matrix* gslCopy = gsl_matrix_alloc(covMat.rows(),covMat.cols());
     94 
     95     eiToGsl(covMat, &gslCovMat);
     96     for (int t=0; t<TRIES; ++t)
     97     {
     98       timerSqrt.start();
     99       for (int k=0; k<repeats; ++k)
    100       {
    101         gsl_matrix_memcpy(gslCopy,gslCovMat);
    102         gsl_linalg_cholesky_decomp(gslCopy);
    103         acc += gsl_matrix_get(gslCopy,r,c);
    104       }
    105       timerSqrt.stop();
    106     }
    107 
    108     std::cout << " | \t"
    109               << timerSqrt.value() * REPEAT / repeats << "s";
    110 
    111     gsl_matrix_free(gslCovMat);
    112   }
    113   #endif
    114   std::cout << "\n";
    115   // make sure the compiler does not optimize too much
    116   if (acc==123)
    117     std::cout << acc;
    118 }
    119 
    120 int main(int argc, char* argv[])
    121 {
    122   const int dynsizes[] = {4,6,8,16,24,32,49,64,128,256,512,900,1500,0};
    123   std::cout << "size            LDLT                            LLT";
    124 //   #ifdef BENCH_GSL
    125 //   std::cout << "       GSL (standard + double + ATLAS)  ";
    126 //   #endif
    127   std::cout << "\n";
    128   for (int i=0; dynsizes[i]>0; ++i)
    129     benchLLT(Matrix<Scalar,Dynamic,Dynamic>(dynsizes[i],dynsizes[i]));
    130 
    131   benchLLT(Matrix<Scalar,2,2>());
    132   benchLLT(Matrix<Scalar,3,3>());
    133   benchLLT(Matrix<Scalar,4,4>());
    134   benchLLT(Matrix<Scalar,5,5>());
    135   benchLLT(Matrix<Scalar,6,6>());
    136   benchLLT(Matrix<Scalar,7,7>());
    137   benchLLT(Matrix<Scalar,8,8>());
    138   benchLLT(Matrix<Scalar,12,12>());
    139   benchLLT(Matrix<Scalar,16,16>());
    140   return 0;
    141 }
    142 
    143