benchEigenSolver.cpp
Go to the documentation of this file.
00001 
00002 // g++ -DNDEBUG -O3 -I.. benchEigenSolver.cpp  -o benchEigenSolver && ./benchEigenSolver
00003 // options:
00004 //  -DBENCH_GMM
00005 //  -DBENCH_GSL -lgsl /usr/lib/libcblas.so.3
00006 //  -DEIGEN_DONT_VECTORIZE
00007 //  -msse2
00008 //  -DREPEAT=100
00009 //  -DTRIES=10
00010 //  -DSCALAR=double
00011 
00012 #include <iostream>
00013 
00014 #include <Eigen/Core>
00015 #include <Eigen/QR>
00016 #include <bench/BenchUtil.h>
00017 using namespace Eigen;
00018 
00019 #ifndef REPEAT
00020 #define REPEAT 1000
00021 #endif
00022 
00023 #ifndef TRIES
00024 #define TRIES 4
00025 #endif
00026 
00027 #ifndef SCALAR
00028 #define SCALAR float
00029 #endif
00030 
00031 typedef SCALAR Scalar;
00032 
00033 template <typename MatrixType>
00034 __attribute__ ((noinline)) void benchEigenSolver(const MatrixType& m)
00035 {
00036   int rows = m.rows();
00037   int cols = m.cols();
00038 
00039   int stdRepeats = std::max(1,int((REPEAT*1000)/(rows*rows*sqrt(rows))));
00040   int saRepeats = stdRepeats * 4;
00041 
00042   typedef typename MatrixType::Scalar Scalar;
00043   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
00044 
00045   MatrixType a = MatrixType::Random(rows,cols);
00046   SquareMatrixType covMat =  a * a.adjoint();
00047 
00048   BenchTimer timerSa, timerStd;
00049 
00050   Scalar acc = 0;
00051   int r = internal::random<int>(0,covMat.rows()-1);
00052   int c = internal::random<int>(0,covMat.cols()-1);
00053   {
00054     SelfAdjointEigenSolver<SquareMatrixType> ei(covMat);
00055     for (int t=0; t<TRIES; ++t)
00056     {
00057       timerSa.start();
00058       for (int k=0; k<saRepeats; ++k)
00059       {
00060         ei.compute(covMat);
00061         acc += ei.eigenvectors().coeff(r,c);
00062       }
00063       timerSa.stop();
00064     }
00065   }
00066 
00067   {
00068     EigenSolver<SquareMatrixType> ei(covMat);
00069     for (int t=0; t<TRIES; ++t)
00070     {
00071       timerStd.start();
00072       for (int k=0; k<stdRepeats; ++k)
00073       {
00074         ei.compute(covMat);
00075         acc += ei.eigenvectors().coeff(r,c);
00076       }
00077       timerStd.stop();
00078     }
00079   }
00080 
00081   if (MatrixType::RowsAtCompileTime==Dynamic)
00082     std::cout << "dyn   ";
00083   else
00084     std::cout << "fixed ";
00085   std::cout << covMat.rows() << " \t"
00086             << timerSa.value() * REPEAT / saRepeats << "s \t"
00087             << timerStd.value() * REPEAT / stdRepeats << "s";
00088 
00089   #ifdef BENCH_GMM
00090   if (MatrixType::RowsAtCompileTime==Dynamic)
00091   {
00092     timerSa.reset();
00093     timerStd.reset();
00094 
00095     gmm::dense_matrix<Scalar> gmmCovMat(covMat.rows(),covMat.cols());
00096     gmm::dense_matrix<Scalar> eigvect(covMat.rows(),covMat.cols());
00097     std::vector<Scalar> eigval(covMat.rows());
00098     eiToGmm(covMat, gmmCovMat);
00099     for (int t=0; t<TRIES; ++t)
00100     {
00101       timerSa.start();
00102       for (int k=0; k<saRepeats; ++k)
00103       {
00104         gmm::symmetric_qr_algorithm(gmmCovMat, eigval, eigvect);
00105         acc += eigvect(r,c);
00106       }
00107       timerSa.stop();
00108     }
00109     // the non-selfadjoint solver does not compute the eigen vectors
00110 //     for (int t=0; t<TRIES; ++t)
00111 //     {
00112 //       timerStd.start();
00113 //       for (int k=0; k<stdRepeats; ++k)
00114 //       {
00115 //         gmm::implicit_qr_algorithm(gmmCovMat, eigval, eigvect);
00116 //         acc += eigvect(r,c);
00117 //       }
00118 //       timerStd.stop();
00119 //     }
00120 
00121     std::cout << " | \t"
00122               << timerSa.value() * REPEAT / saRepeats << "s"
00123               << /*timerStd.value() * REPEAT / stdRepeats << "s"*/ "   na   ";
00124   }
00125   #endif
00126 
00127   #ifdef BENCH_GSL
00128   if (MatrixType::RowsAtCompileTime==Dynamic)
00129   {
00130     timerSa.reset();
00131     timerStd.reset();
00132 
00133     gsl_matrix* gslCovMat = gsl_matrix_alloc(covMat.rows(),covMat.cols());
00134     gsl_matrix* gslCopy = gsl_matrix_alloc(covMat.rows(),covMat.cols());
00135     gsl_matrix* eigvect = gsl_matrix_alloc(covMat.rows(),covMat.cols());
00136     gsl_vector* eigval  = gsl_vector_alloc(covMat.rows());
00137     gsl_eigen_symmv_workspace* eisymm = gsl_eigen_symmv_alloc(covMat.rows());
00138     
00139     gsl_matrix_complex* eigvectz = gsl_matrix_complex_alloc(covMat.rows(),covMat.cols());
00140     gsl_vector_complex* eigvalz  = gsl_vector_complex_alloc(covMat.rows());
00141     gsl_eigen_nonsymmv_workspace* einonsymm = gsl_eigen_nonsymmv_alloc(covMat.rows());
00142     
00143     eiToGsl(covMat, &gslCovMat);
00144     for (int t=0; t<TRIES; ++t)
00145     {
00146       timerSa.start();
00147       for (int k=0; k<saRepeats; ++k)
00148       {
00149         gsl_matrix_memcpy(gslCopy,gslCovMat);
00150         gsl_eigen_symmv(gslCopy, eigval, eigvect, eisymm);
00151         acc += gsl_matrix_get(eigvect,r,c);
00152       }
00153       timerSa.stop();
00154     }
00155     for (int t=0; t<TRIES; ++t)
00156     {
00157       timerStd.start();
00158       for (int k=0; k<stdRepeats; ++k)
00159       {
00160         gsl_matrix_memcpy(gslCopy,gslCovMat);
00161         gsl_eigen_nonsymmv(gslCopy, eigvalz, eigvectz, einonsymm);
00162         acc += GSL_REAL(gsl_matrix_complex_get(eigvectz,r,c));
00163       }
00164       timerStd.stop();
00165     }
00166 
00167     std::cout << " | \t"
00168               << timerSa.value() * REPEAT / saRepeats << "s \t"
00169               << timerStd.value() * REPEAT / stdRepeats << "s";
00170 
00171     gsl_matrix_free(gslCovMat);
00172     gsl_vector_free(gslCopy);
00173     gsl_matrix_free(eigvect);
00174     gsl_vector_free(eigval);
00175     gsl_matrix_complex_free(eigvectz);
00176     gsl_vector_complex_free(eigvalz);
00177     gsl_eigen_symmv_free(eisymm);
00178     gsl_eigen_nonsymmv_free(einonsymm);
00179   }
00180   #endif
00181 
00182   std::cout << "\n";
00183   
00184   // make sure the compiler does not optimize too much
00185   if (acc==123)
00186     std::cout << acc;
00187 }
00188 
00189 int main(int argc, char* argv[])
00190 {
00191   const int dynsizes[] = {4,6,8,12,16,24,32,64,128,256,512,0};
00192   std::cout << "size            selfadjoint       generic";
00193   #ifdef BENCH_GMM
00194   std::cout << "        GMM++          ";
00195   #endif
00196   #ifdef BENCH_GSL
00197   std::cout << "       GSL (double + ATLAS)  ";
00198   #endif
00199   std::cout << "\n";
00200   for (uint i=0; dynsizes[i]>0; ++i)
00201     benchEigenSolver(Matrix<Scalar,Dynamic,Dynamic>(dynsizes[i],dynsizes[i]));
00202 
00203   benchEigenSolver(Matrix<Scalar,2,2>());
00204   benchEigenSolver(Matrix<Scalar,3,3>());
00205   benchEigenSolver(Matrix<Scalar,4,4>());
00206   benchEigenSolver(Matrix<Scalar,6,6>());
00207   benchEigenSolver(Matrix<Scalar,8,8>());
00208   benchEigenSolver(Matrix<Scalar,12,12>());
00209   benchEigenSolver(Matrix<Scalar,16,16>());
00210   return 0;
00211 }
00212 


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