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>
benchLLT(const MatrixType & m)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
main(int argc,char * argv[])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
142