cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

benchCholesky.cpp (3548B)

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