dense_solvers.cpp (6416B)
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 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 25 void compute(Solver &solver, const MatrixType &A) { 26 solver.compute(A); 27 } 28 29 template<typename Scalar,int 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 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 }