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>
16 using namespace Eigen;
17
18 #ifndef REPEAT
19 #define REPEAT 10000
20 #endif
21
22 #ifndef TRIES
23 #define TRIES 10
24 #endif
25
26 typedef float Scalar;
27
28 template <typename MatrixType>
benchLLT(const MatrixType & m)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
main(int argc,char * argv[])120 int 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