1 #include <iostream>
2 #include "BenchTimer.h"
3 #include <Eigen/Dense>
4 #include <map>
5 #include <vector>
6 #include <string>
7 #include <sstream>
8 using namespace Eigen;
9
10 std::map<std::string,Array<float,1,8,DontAlign|RowMajor> > results;
11 std::vector<std::string> labels;
12 std::vector<Array2i> sizes;
13
14 template<typename Solver,typename MatrixType>
15 EIGEN_DONT_INLINE
compute_norm_equation(Solver & solver,const MatrixType & A)16 void compute_norm_equation(Solver &solver, const MatrixType &A) {
17 if(A.rows()!=A.cols())
18 solver.compute(A.transpose()*A);
19 else
20 solver.compute(A);
21 }
22
23 template<typename Solver,typename MatrixType>
24 EIGEN_DONT_INLINE
compute(Solver & solver,const MatrixType & A)25 void compute(Solver &solver, const MatrixType &A) {
26 solver.compute(A);
27 }
28
29 template<typename Scalar,int Size>
bench(int id,int rows,int size=Size)30 void bench(int id, int rows, int size = Size)
31 {
32 typedef Matrix<Scalar,Dynamic,Size> Mat;
33 typedef Matrix<Scalar,Dynamic,Dynamic> MatDyn;
34 typedef Matrix<Scalar,Size,Size> MatSquare;
35 Mat A(rows,size);
36 A.setRandom();
37 if(rows==size)
38 A = A*A.adjoint();
39 BenchTimer t_llt, t_ldlt, t_lu, t_fplu, t_qr, t_cpqr, t_cod, t_fpqr, t_jsvd, t_bdcsvd;
40
41 int svd_opt = ComputeThinU|ComputeThinV;
42
43 int tries = 5;
44 int rep = 1000/size;
45 if(rep==0) rep = 1;
46 // rep = rep*rep;
47
48 LLT<MatSquare> llt(size);
49 LDLT<MatSquare> ldlt(size);
50 PartialPivLU<MatSquare> lu(size);
51 FullPivLU<MatSquare> fplu(size,size);
52 HouseholderQR<Mat> qr(A.rows(),A.cols());
53 ColPivHouseholderQR<Mat> cpqr(A.rows(),A.cols());
54 CompleteOrthogonalDecomposition<Mat> cod(A.rows(),A.cols());
55 FullPivHouseholderQR<Mat> fpqr(A.rows(),A.cols());
56 JacobiSVD<MatDyn> jsvd(A.rows(),A.cols());
57 BDCSVD<MatDyn> bdcsvd(A.rows(),A.cols());
58
59 BENCH(t_llt, tries, rep, compute_norm_equation(llt,A));
60 BENCH(t_ldlt, tries, rep, compute_norm_equation(ldlt,A));
61 BENCH(t_lu, tries, rep, compute_norm_equation(lu,A));
62 if(size<=1000)
63 BENCH(t_fplu, tries, rep, compute_norm_equation(fplu,A));
64 BENCH(t_qr, tries, rep, compute(qr,A));
65 BENCH(t_cpqr, tries, rep, compute(cpqr,A));
66 BENCH(t_cod, tries, rep, compute(cod,A));
67 if(size*rows<=10000000)
68 BENCH(t_fpqr, tries, rep, compute(fpqr,A));
69 if(size<500) // JacobiSVD is really too slow for too large matrices
70 BENCH(t_jsvd, tries, rep, jsvd.compute(A,svd_opt));
71 // if(size*rows<=20000000)
72 BENCH(t_bdcsvd, tries, rep, bdcsvd.compute(A,svd_opt));
73
74 results["LLT"][id] = t_llt.best();
75 results["LDLT"][id] = t_ldlt.best();
76 results["PartialPivLU"][id] = t_lu.best();
77 results["FullPivLU"][id] = t_fplu.best();
78 results["HouseholderQR"][id] = t_qr.best();
79 results["ColPivHouseholderQR"][id] = t_cpqr.best();
80 results["CompleteOrthogonalDecomposition"][id] = t_cod.best();
81 results["FullPivHouseholderQR"][id] = t_fpqr.best();
82 results["JacobiSVD"][id] = t_jsvd.best();
83 results["BDCSVD"][id] = t_bdcsvd.best();
84 }
85
86
main()87 int main()
88 {
89 labels.push_back("LLT");
90 labels.push_back("LDLT");
91 labels.push_back("PartialPivLU");
92 labels.push_back("FullPivLU");
93 labels.push_back("HouseholderQR");
94 labels.push_back("ColPivHouseholderQR");
95 labels.push_back("CompleteOrthogonalDecomposition");
96 labels.push_back("FullPivHouseholderQR");
97 labels.push_back("JacobiSVD");
98 labels.push_back("BDCSVD");
99
100 for(int i=0; i<labels.size(); ++i)
101 results[labels[i]].fill(-1);
102
103 const int small = 8;
104 sizes.push_back(Array2i(small,small));
105 sizes.push_back(Array2i(100,100));
106 sizes.push_back(Array2i(1000,1000));
107 sizes.push_back(Array2i(4000,4000));
108 sizes.push_back(Array2i(10000,small));
109 sizes.push_back(Array2i(10000,100));
110 sizes.push_back(Array2i(10000,1000));
111 sizes.push_back(Array2i(10000,4000));
112
113 using namespace std;
114
115 for(int k=0; k<sizes.size(); ++k)
116 {
117 cout << sizes[k](0) << "x" << sizes[k](1) << "...\n";
118 bench<float,Dynamic>(k,sizes[k](0),sizes[k](1));
119 }
120
121 cout.width(32);
122 cout << "solver/size";
123 cout << " ";
124 for(int k=0; k<sizes.size(); ++k)
125 {
126 std::stringstream ss;
127 ss << sizes[k](0) << "x" << sizes[k](1);
128 cout.width(10); cout << ss.str(); cout << " ";
129 }
130 cout << endl;
131
132
133 for(int i=0; i<labels.size(); ++i)
134 {
135 cout.width(32); cout << labels[i]; cout << " ";
136 ArrayXf r = (results[labels[i]]*100000.f).floor()/100.f;
137 for(int k=0; k<sizes.size(); ++k)
138 {
139 cout.width(10);
140 if(r(k)>=1e6) cout << "-";
141 else cout << r(k);
142 cout << " ";
143 }
144 cout << endl;
145 }
146
147 // HTML output
148 cout << "<table class=\"manual\">" << endl;
149 cout << "<tr><th>solver/size</th>" << endl;
150 for(int k=0; k<sizes.size(); ++k)
151 cout << " <th>" << sizes[k](0) << "x" << sizes[k](1) << "</th>";
152 cout << "</tr>" << endl;
153 for(int i=0; i<labels.size(); ++i)
154 {
155 cout << "<tr";
156 if(i%2==1) cout << " class=\"alt\"";
157 cout << "><td>" << labels[i] << "</td>";
158 ArrayXf r = (results[labels[i]]*100000.f).floor()/100.f;
159 for(int k=0; k<sizes.size(); ++k)
160 {
161 if(r(k)>=1e6) cout << "<td>-</td>";
162 else
163 {
164 cout << "<td>" << r(k);
165 if(i>0)
166 cout << " (x" << numext::round(10.f*results[labels[i]](k)/results["LLT"](k))/10.f << ")";
167 if(i<4 && sizes[k](0)!=sizes[k](1))
168 cout << " <sup><a href=\"#note_ls\">*</a></sup>";
169 cout << "</td>";
170 }
171 }
172 cout << "</tr>" << endl;
173 }
174 cout << "</table>" << endl;
175
176 // cout << "LLT (ms) " << (results["LLT"]*1000.).format(fmt) << "\n";
177 // cout << "LDLT (%) " << (results["LDLT"]/results["LLT"]).format(fmt) << "\n";
178 // cout << "PartialPivLU (%) " << (results["PartialPivLU"]/results["LLT"]).format(fmt) << "\n";
179 // cout << "FullPivLU (%) " << (results["FullPivLU"]/results["LLT"]).format(fmt) << "\n";
180 // cout << "HouseholderQR (%) " << (results["HouseholderQR"]/results["LLT"]).format(fmt) << "\n";
181 // cout << "ColPivHouseholderQR (%) " << (results["ColPivHouseholderQR"]/results["LLT"]).format(fmt) << "\n";
182 // cout << "CompleteOrthogonalDecomposition (%) " << (results["CompleteOrthogonalDecomposition"]/results["LLT"]).format(fmt) << "\n";
183 // cout << "FullPivHouseholderQR (%) " << (results["FullPivHouseholderQR"]/results["LLT"]).format(fmt) << "\n";
184 // cout << "JacobiSVD (%) " << (results["JacobiSVD"]/results["LLT"]).format(fmt) << "\n";
185 // cout << "BDCSVD (%) " << (results["BDCSVD"]/results["LLT"]).format(fmt) << "\n";
186 }
187