benchCholesky.cpp
Go to the documentation of this file.
00001 
00002 // g++ -DNDEBUG -O3 -I.. benchLLT.cpp  -o benchLLT && ./benchLLT
00003 // options:
00004 //  -DBENCH_GSL -lgsl /usr/lib/libcblas.so.3
00005 //  -DEIGEN_DONT_VECTORIZE
00006 //  -msse2
00007 //  -DREPEAT=100
00008 //  -DTRIES=10
00009 //  -DSCALAR=double
00010 
00011 #include <iostream>
00012 
00013 #include <Eigen/Core>
00014 #include <Eigen/Cholesky>
00015 #include <bench/BenchUtil.h>
00016 using namespace Eigen;
00017 
00018 #ifndef REPEAT
00019 #define REPEAT 10000
00020 #endif
00021 
00022 #ifndef TRIES
00023 #define TRIES 10
00024 #endif
00025 
00026 typedef float Scalar;
00027 
00028 template <typename MatrixType>
00029 __attribute__ ((noinline)) void benchLLT(const MatrixType& m)
00030 {
00031   int rows = m.rows();
00032   int cols = m.cols();
00033 
00034   int cost = 0;
00035   for (int j=0; j<rows; ++j)
00036   {
00037     int r = std::max(rows - j -1,0);
00038     cost += 2*(r*j+r+j);
00039   }
00040 
00041   int repeats = (REPEAT*1000)/(rows*rows);
00042 
00043   typedef typename MatrixType::Scalar Scalar;
00044   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
00045 
00046   MatrixType a = MatrixType::Random(rows,cols);
00047   SquareMatrixType covMat =  a * a.adjoint();
00048 
00049   BenchTimer timerNoSqrt, timerSqrt;
00050 
00051   Scalar acc = 0;
00052   int r = internal::random<int>(0,covMat.rows()-1);
00053   int c = internal::random<int>(0,covMat.cols()-1);
00054   for (int t=0; t<TRIES; ++t)
00055   {
00056     timerNoSqrt.start();
00057     for (int k=0; k<repeats; ++k)
00058     {
00059       LDLT<SquareMatrixType> cholnosqrt(covMat);
00060       acc += cholnosqrt.matrixL().coeff(r,c);
00061     }
00062     timerNoSqrt.stop();
00063   }
00064 
00065   for (int t=0; t<TRIES; ++t)
00066   {
00067     timerSqrt.start();
00068     for (int k=0; k<repeats; ++k)
00069     {
00070       LLT<SquareMatrixType> chol(covMat);
00071       acc += chol.matrixL().coeff(r,c);
00072     }
00073     timerSqrt.stop();
00074   }
00075 
00076   if (MatrixType::RowsAtCompileTime==Dynamic)
00077     std::cout << "dyn   ";
00078   else
00079     std::cout << "fixed ";
00080   std::cout << covMat.rows() << " \t"
00081             << (timerNoSqrt.value() * REPEAT) / repeats << "s "
00082             << "(" << 1e-6 * cost*repeats/timerNoSqrt.value() << " MFLOPS)\t"
00083             << (timerSqrt.value() * REPEAT) / repeats << "s "
00084             << "(" << 1e-6 * cost*repeats/timerSqrt.value() << " MFLOPS)\n";
00085 
00086 
00087   #ifdef BENCH_GSL
00088   if (MatrixType::RowsAtCompileTime==Dynamic)
00089   {
00090     timerSqrt.reset();
00091 
00092     gsl_matrix* gslCovMat = gsl_matrix_alloc(covMat.rows(),covMat.cols());
00093     gsl_matrix* gslCopy = gsl_matrix_alloc(covMat.rows(),covMat.cols());
00094 
00095     eiToGsl(covMat, &gslCovMat);
00096     for (int t=0; t<TRIES; ++t)
00097     {
00098       timerSqrt.start();
00099       for (int k=0; k<repeats; ++k)
00100       {
00101         gsl_matrix_memcpy(gslCopy,gslCovMat);
00102         gsl_linalg_cholesky_decomp(gslCopy);
00103         acc += gsl_matrix_get(gslCopy,r,c);
00104       }
00105       timerSqrt.stop();
00106     }
00107 
00108     std::cout << " | \t"
00109               << timerSqrt.value() * REPEAT / repeats << "s";
00110 
00111     gsl_matrix_free(gslCovMat);
00112   }
00113   #endif
00114   std::cout << "\n";
00115   // make sure the compiler does not optimize too much
00116   if (acc==123)
00117     std::cout << acc;
00118 }
00119 
00120 int main(int argc, char* argv[])
00121 {
00122   const int dynsizes[] = {4,6,8,16,24,32,49,64,128,256,512,900,0};
00123   std::cout << "size            no sqrt                           standard";
00124 //   #ifdef BENCH_GSL
00125 //   std::cout << "       GSL (standard + double + ATLAS)  ";
00126 //   #endif
00127   std::cout << "\n";
00128   for (uint i=0; dynsizes[i]>0; ++i)
00129     benchLLT(Matrix<Scalar,Dynamic,Dynamic>(dynsizes[i],dynsizes[i]));
00130 
00131   benchLLT(Matrix<Scalar,2,2>());
00132   benchLLT(Matrix<Scalar,3,3>());
00133   benchLLT(Matrix<Scalar,4,4>());
00134   benchLLT(Matrix<Scalar,5,5>());
00135   benchLLT(Matrix<Scalar,6,6>());
00136   benchLLT(Matrix<Scalar,7,7>());
00137   benchLLT(Matrix<Scalar,8,8>());
00138   benchLLT(Matrix<Scalar,12,12>());
00139   benchLLT(Matrix<Scalar,16,16>());
00140   return 0;
00141 }
00142 


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:30