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>
16using namespace Eigen;
17
18#ifndef REPEAT
19#define REPEAT 10000
20#endif
21
22#ifndef TRIES
23#define TRIES 10
24#endif
25
26typedef float Scalar;
27
28template <typename MatrixType>
29__attribute__ ((noinline)) void benchLLT(const MatrixType& m)
30{
31  int rows = m.rows();
32  int cols = m.cols();
33
34  int 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.value() * REPEAT) / repeats << "s "
82            << "(" << 1e-6 * cost*repeats/timerNoSqrt.value() << " MFLOPS)\t"
83            << (timerSqrt.value() * REPEAT) / repeats << "s "
84            << "(" << 1e-6 * cost*repeats/timerSqrt.value() << " MFLOPS)\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
120int main(int argc, char* argv[])
121{
122  const int dynsizes[] = {4,6,8,16,24,32,49,64,128,256,512,900,0};
123  std::cout << "size            no sqrt                           standard";
124//   #ifdef BENCH_GSL
125//   std::cout << "       GSL (standard + double + ATLAS)  ";
126//   #endif
127  std::cout << "\n";
128  for (uint 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