eigen3_decompositions.cpp
Go to the documentation of this file.
1 
9 /*****************************************************************************
10 ** Includes
11 *****************************************************************************/
12 
13 #include <iostream>
14 #include <ecl/linear_algebra.hpp>
15 #include <ecl/threads/priority.hpp>
16 #include <ecl/time/stopwatch.hpp>
17 #include <ecl/time/timestamp.hpp>
18 #include <ecl/exceptions/standard_exception.hpp>
19 #include <ecl/formatters.hpp>
20 
21 /*****************************************************************************
22 ** Using
23 *****************************************************************************/
24 
25 using ecl::StandardException;
26 using ecl::StopWatch;
27 using ecl::TimeStamp;
28 using Eigen::Matrix3f;
29 using Eigen::MatrixXf;
30 
31 /*****************************************************************************
32 ** Main
33 *****************************************************************************/
34 
35 int main(int argc, char **argv) {
36 
37  const unsigned int repeats = 5;
38  try {
39  ecl::set_priority(ecl::RealTimePriority4);
40  } catch ( StandardException &e ) {
41  // dont worry about it.
42  }
43  StopWatch stopwatch;
44  TimeStamp times[13];
45 
46  MatrixXf A = MatrixXf::Random(100,100);
47  MatrixXf b = MatrixXf::Identity(100,100);
48  MatrixXf x(100,1);
49 
50  // get rid of caching effects.
51  for ( unsigned int i = 0; i < repeats; ++i ) {
52  x = A.colPivHouseholderQr().solve(b);
53  }
54 
55  times[0].stamp(0,0);
56  for ( unsigned int i = 0; i < repeats; ++i ) {
57  A = MatrixXf::Random(100,100);
58  stopwatch.restart();
59  x = A.colPivHouseholderQr().solve(b);
60  times[0] += stopwatch.split();
61  }
62  times[1].stamp(0,0);
63  for ( unsigned int i = 0; i < repeats; ++i ) {
64  A = MatrixXf::Random(100,100);
65  stopwatch.restart();
66  x = A.fullPivHouseholderQr().solve(b);
67  times[1] += stopwatch.split();
68  }
69  times[2].stamp(0,0);
70  for ( unsigned int i = 0; i < repeats; ++i ) {
71  A = MatrixXf::Random(100,100);
72  stopwatch.restart();
73  x = A.partialPivLu().solve(b);
74  times[2] += stopwatch.split();
75  }
76  times[3].stamp(0,0);
77  for ( unsigned int i = 0; i < repeats; ++i ) {
78  A = MatrixXf::Random(100,100);
79  stopwatch.restart();
80  x = A.fullPivLu().solve(b);
81  times[3] += stopwatch.split();
82  }
83 
84  ecl::Format<double> format(15,9,ecl::RightAlign);
85  std::cout << std::endl;
86  std::cout << "************** Eigen3 Decompositions ***************" << std::endl;
87  std::cout << std::endl;
88  std::cout << "Ax=b where A = Rand(100,100) and b = Identity(100x100)" << std::endl;
89  std::cout << "i.e. solving the inverse (results averaged over " << repeats << " runs)" << std::endl;
90  std::cout << std::endl;
91  std::cout << "ColPivHouseHolder : " << format(static_cast<double>(times[0])/static_cast<double>(repeats)) << std::endl;
92  std::cout << "FullPivHouseHolder : " << format(static_cast<double>(times[1])/static_cast<double>(repeats)) << std::endl;
93  std::cout << "PartialPivLu : " << format(static_cast<double>(times[2])/static_cast<double>(repeats)) << std::endl;
94  std::cout << "FullPivLu : " << format(static_cast<double>(times[3])/static_cast<double>(repeats)) << std::endl;
95  std::cout << std::endl;
96 
97  b = MatrixXf::Random(100,1);
98  times[4].stamp(0,0);
99  for ( unsigned int i = 0; i < repeats; ++i ) {
100  A = MatrixXf::Random(100,100);
101  stopwatch.restart();
102  x = A.colPivHouseholderQr().solve(b);
103  times[4] += stopwatch.split();
104  }
105  times[5].stamp(0,0);
106  for ( unsigned int i = 0; i < repeats; ++i ) {
107  A = MatrixXf::Random(100,100);
108  stopwatch.restart();
109  x = A.fullPivHouseholderQr().solve(b);
110  times[5] += stopwatch.split();
111  }
112  times[6].stamp(0,0);
113  for ( unsigned int i = 0; i < repeats; ++i ) {
114  A = MatrixXf::Random(100,100);
115  stopwatch.restart();
116  x = A.partialPivLu().solve(b);
117  times[6] += stopwatch.split();
118  }
119  times[7].stamp(0,0);
120  for ( unsigned int i = 0; i < repeats; ++i ) {
121  A = MatrixXf::Random(100,100);
122  stopwatch.restart();
123  x = A.fullPivLu().solve(b);
124  times[7] += stopwatch.split();
125  }
126 
127  std::cout << "Ax=b where A = Rand(100,100) and b = Random(100x1)" << std::endl;
128  std::cout << "i.e. solving linear equations (results averaged over " << repeats << " runs)" << std::endl;
129  std::cout << std::endl;
130  std::cout << "ColPivHouseHolder : " << format(static_cast<double>(times[4])/static_cast<double>(repeats)) << std::endl;
131  std::cout << "FullPivHouseHolder : " << format(static_cast<double>(times[5])/static_cast<double>(repeats)) << std::endl;
132  std::cout << "PartialPivLu : " << format(static_cast<double>(times[6])/static_cast<double>(repeats)) << std::endl;
133  std::cout << "FullPivLu : " << format(static_cast<double>(times[7])/static_cast<double>(repeats)) << std::endl;
134  std::cout << std::endl;
135 
136  b = MatrixXf::Identity(100,100);
137  times[8].stamp(0,0);
138  for ( unsigned int i = 0; i < repeats; ++i ) {
139  MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous
140  A = D.transpose()*D; // because this will only be positive definite if D is invertible
141  stopwatch.restart();
142  x = A.colPivHouseholderQr().solve(b);
143  times[8] += stopwatch.split();
144  }
145  times[9].stamp(0,0);
146  for ( unsigned int i = 0; i < repeats; ++i ) {
147  MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous
148  A = D.transpose()*D; // because this will only be positive definite if D is invertible
149  stopwatch.restart();
150  x = A.fullPivHouseholderQr().solve(b);
151  times[9] += stopwatch.split();
152  }
153  times[10].stamp(0,0);
154  for ( unsigned int i = 0; i < repeats; ++i ) {
155  MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous
156  A = D.transpose()*D; // because this will only be positive definite if D is invertible
157  stopwatch.restart();
158  x = A.partialPivLu().solve(b);
159  times[10] += stopwatch.split();
160  }
161  times[11].stamp(0,0);
162  for ( unsigned int i = 0; i < repeats; ++i ) {
163  MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous
164  A = D.transpose()*D; // because this will only be positive definite if D is invertible
165  stopwatch.restart();
166  x = A.fullPivLu().solve(b);
167  times[11] += stopwatch.split();
168  }
169  times[12].stamp(0,0);
170  for ( unsigned int i = 0; i < repeats; ++i ) {
171  MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous
172  A = D.transpose()*D; // because this will only be positive definite if D is invertible
173  stopwatch.restart();
174  x = A.llt().solve(b);
175 // A.llt().solveInPlace(b); // only marginally faster
176  times[12] += stopwatch.split();
177  }
178  times[13].stamp(0,0);
179  for ( unsigned int i = 0; i < repeats; ++i ) {
180  MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous
181  A = D.transpose()*D; // because this will only be positive definite if D is invertible
182  stopwatch.restart();
183  x = A.ldlt().solve(b);
184 // A.ldlt().solveInPlace(b); // only marginally faster
185  times[13] += stopwatch.split();
186  }
187 
188  std::cout << "Ax=b where A = RandPosDef(100,100) and b = Identity(100x100)" << std::endl;
189  std::cout << "(results averaged over " << repeats << " runs)" << std::endl;
190  std::cout << std::endl;
191  std::cout << "ColPivHouseHolder : " << format(static_cast<double>(times[8])/static_cast<double>(repeats)) << std::endl;
192  std::cout << "FullPivHouseHolder : " << format(static_cast<double>(times[9])/static_cast<double>(repeats)) << std::endl;
193  std::cout << "PartialPivLu : " << format(static_cast<double>(times[10])/static_cast<double>(repeats)) << std::endl;
194  std::cout << "FullPivLu : " << format(static_cast<double>(times[11])/static_cast<double>(repeats)) << std::endl;
195  std::cout << "llt : " << format(static_cast<double>(times[12])/static_cast<double>(repeats)) << std::endl;
196  std::cout << "ldlt : " << format(static_cast<double>(times[13])/static_cast<double>(repeats)) << std::endl;
197  std::cout << std::endl;
198 
199  return 0;
200 }
ecl::Format< double >
formatters.hpp
main
int main(int argc, char **argv)
Definition: eigen3_decompositions.cpp:35
ecl::RightAlign
RightAlign
linear_algebra.hpp
ecl::StandardException
stopwatch.hpp
standard_exception.hpp
priority.hpp
ecl::RealTimePriority4
RealTimePriority4
timestamp.hpp

ecl_core_apps
Author(s): Daniel Stonier
autogenerated on Wed Mar 2 2022 00:16:52