$search
00001 00009 /***************************************************************************** 00010 ** Includes 00011 *****************************************************************************/ 00012 00013 #include <iostream> 00014 #include <ecl/linear_algebra.hpp> 00015 #include <ecl/threads/priority.hpp> 00016 #include <ecl/time/stopwatch.hpp> 00017 #include <ecl/time/timestamp.hpp> 00018 #include <ecl/exceptions/standard_exception.hpp> 00019 #include <ecl/formatters.hpp> 00020 00021 /***************************************************************************** 00022 ** Using 00023 *****************************************************************************/ 00024 00025 using ecl::StandardException; 00026 using ecl::StopWatch; 00027 using ecl::TimeStamp; 00028 using Eigen::Matrix3f; 00029 using Eigen::MatrixXf; 00030 00031 /***************************************************************************** 00032 ** Main 00033 *****************************************************************************/ 00034 00035 int main(int argc, char **argv) { 00036 00037 const unsigned int repeats = 5; 00038 try { 00039 ecl::set_priority(ecl::RealTimePriority4); 00040 } catch ( StandardException &e ) { 00041 // dont worry about it. 00042 } 00043 StopWatch stopwatch; 00044 TimeStamp times[13]; 00045 00046 MatrixXf A = MatrixXf::Random(100,100); 00047 MatrixXf b = MatrixXf::Identity(100,100); 00048 MatrixXf x(100,1); 00049 00050 // get rid of caching effects. 00051 for ( unsigned int i = 0; i < repeats; ++i ) { 00052 x = A.colPivHouseholderQr().solve(b); 00053 } 00054 00055 times[0].stamp(0,0); 00056 for ( unsigned int i = 0; i < repeats; ++i ) { 00057 A = MatrixXf::Random(100,100); 00058 stopwatch.restart(); 00059 x = A.colPivHouseholderQr().solve(b); 00060 times[0] += stopwatch.split(); 00061 } 00062 times[1].stamp(0,0); 00063 for ( unsigned int i = 0; i < repeats; ++i ) { 00064 A = MatrixXf::Random(100,100); 00065 stopwatch.restart(); 00066 x = A.fullPivHouseholderQr().solve(b); 00067 times[1] += stopwatch.split(); 00068 } 00069 times[2].stamp(0,0); 00070 for ( unsigned int i = 0; i < repeats; ++i ) { 00071 A = MatrixXf::Random(100,100); 00072 stopwatch.restart(); 00073 x = A.partialPivLu().solve(b); 00074 times[2] += stopwatch.split(); 00075 } 00076 times[3].stamp(0,0); 00077 for ( unsigned int i = 0; i < repeats; ++i ) { 00078 A = MatrixXf::Random(100,100); 00079 stopwatch.restart(); 00080 x = A.fullPivLu().solve(b); 00081 times[3] += stopwatch.split(); 00082 } 00083 00084 ecl::Format<double> format(15,9,ecl::RightAlign); 00085 std::cout << std::endl; 00086 std::cout << "************** Eigen3 Decompositions ***************" << std::endl; 00087 std::cout << std::endl; 00088 std::cout << "Ax=b where A = Rand(100,100) and b = Identity(100x100)" << std::endl; 00089 std::cout << "i.e. solving the inverse (results averaged over " << repeats << " runs)" << std::endl; 00090 std::cout << std::endl; 00091 std::cout << "ColPivHouseHolder : " << format(static_cast<double>(times[0])/static_cast<double>(repeats)) << std::endl; 00092 std::cout << "FullPivHouseHolder : " << format(static_cast<double>(times[1])/static_cast<double>(repeats)) << std::endl; 00093 std::cout << "PartialPivLu : " << format(static_cast<double>(times[2])/static_cast<double>(repeats)) << std::endl; 00094 std::cout << "FullPivLu : " << format(static_cast<double>(times[3])/static_cast<double>(repeats)) << std::endl; 00095 std::cout << std::endl; 00096 00097 b = MatrixXf::Random(100,1); 00098 times[4].stamp(0,0); 00099 for ( unsigned int i = 0; i < repeats; ++i ) { 00100 A = MatrixXf::Random(100,100); 00101 stopwatch.restart(); 00102 x = A.colPivHouseholderQr().solve(b); 00103 times[4] += stopwatch.split(); 00104 } 00105 times[5].stamp(0,0); 00106 for ( unsigned int i = 0; i < repeats; ++i ) { 00107 A = MatrixXf::Random(100,100); 00108 stopwatch.restart(); 00109 x = A.fullPivHouseholderQr().solve(b); 00110 times[5] += stopwatch.split(); 00111 } 00112 times[6].stamp(0,0); 00113 for ( unsigned int i = 0; i < repeats; ++i ) { 00114 A = MatrixXf::Random(100,100); 00115 stopwatch.restart(); 00116 x = A.partialPivLu().solve(b); 00117 times[6] += stopwatch.split(); 00118 } 00119 times[7].stamp(0,0); 00120 for ( unsigned int i = 0; i < repeats; ++i ) { 00121 A = MatrixXf::Random(100,100); 00122 stopwatch.restart(); 00123 x = A.fullPivLu().solve(b); 00124 times[7] += stopwatch.split(); 00125 } 00126 00127 std::cout << "Ax=b where A = Rand(100,100) and b = Random(100x1)" << std::endl; 00128 std::cout << "i.e. solving linear equations (results averaged over " << repeats << " runs)" << std::endl; 00129 std::cout << std::endl; 00130 std::cout << "ColPivHouseHolder : " << format(static_cast<double>(times[4])/static_cast<double>(repeats)) << std::endl; 00131 std::cout << "FullPivHouseHolder : " << format(static_cast<double>(times[5])/static_cast<double>(repeats)) << std::endl; 00132 std::cout << "PartialPivLu : " << format(static_cast<double>(times[6])/static_cast<double>(repeats)) << std::endl; 00133 std::cout << "FullPivLu : " << format(static_cast<double>(times[7])/static_cast<double>(repeats)) << std::endl; 00134 std::cout << std::endl; 00135 00136 b = MatrixXf::Identity(100,100); 00137 times[8].stamp(0,0); 00138 for ( unsigned int i = 0; i < repeats; ++i ) { 00139 MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous 00140 A = D.transpose()*D; // because this will only be positive definite if D is invertible 00141 stopwatch.restart(); 00142 x = A.colPivHouseholderQr().solve(b); 00143 times[8] += stopwatch.split(); 00144 } 00145 times[9].stamp(0,0); 00146 for ( unsigned int i = 0; i < repeats; ++i ) { 00147 MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous 00148 A = D.transpose()*D; // because this will only be positive definite if D is invertible 00149 stopwatch.restart(); 00150 x = A.fullPivHouseholderQr().solve(b); 00151 times[9] += stopwatch.split(); 00152 } 00153 times[10].stamp(0,0); 00154 for ( unsigned int i = 0; i < repeats; ++i ) { 00155 MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous 00156 A = D.transpose()*D; // because this will only be positive definite if D is invertible 00157 stopwatch.restart(); 00158 x = A.partialPivLu().solve(b); 00159 times[10] += stopwatch.split(); 00160 } 00161 times[11].stamp(0,0); 00162 for ( unsigned int i = 0; i < repeats; ++i ) { 00163 MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous 00164 A = D.transpose()*D; // because this will only be positive definite if D is invertible 00165 stopwatch.restart(); 00166 x = A.fullPivLu().solve(b); 00167 times[11] += stopwatch.split(); 00168 } 00169 times[12].stamp(0,0); 00170 for ( unsigned int i = 0; i < repeats; ++i ) { 00171 MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous 00172 A = D.transpose()*D; // because this will only be positive definite if D is invertible 00173 stopwatch.restart(); 00174 x = A.llt().solve(b); 00175 // A.llt().solveInPlace(b); // only marginally faster 00176 times[12] += stopwatch.split(); 00177 } 00178 times[13].stamp(0,0); 00179 for ( unsigned int i = 0; i < repeats; ++i ) { 00180 MatrixXf D = MatrixXf::Random(100,100); // This is actually a bit dangerous 00181 A = D.transpose()*D; // because this will only be positive definite if D is invertible 00182 stopwatch.restart(); 00183 x = A.ldlt().solve(b); 00184 // A.ldlt().solveInPlace(b); // only marginally faster 00185 times[13] += stopwatch.split(); 00186 } 00187 00188 std::cout << "Ax=b where A = RandPosDef(100,100) and b = Identity(100x100)" << std::endl; 00189 std::cout << "(results averaged over " << repeats << " runs)" << std::endl; 00190 std::cout << std::endl; 00191 std::cout << "ColPivHouseHolder : " << format(static_cast<double>(times[8])/static_cast<double>(repeats)) << std::endl; 00192 std::cout << "FullPivHouseHolder : " << format(static_cast<double>(times[9])/static_cast<double>(repeats)) << std::endl; 00193 std::cout << "PartialPivLu : " << format(static_cast<double>(times[10])/static_cast<double>(repeats)) << std::endl; 00194 std::cout << "FullPivLu : " << format(static_cast<double>(times[11])/static_cast<double>(repeats)) << std::endl; 00195 std::cout << "llt : " << format(static_cast<double>(times[12])/static_cast<double>(repeats)) << std::endl; 00196 std::cout << "ldlt : " << format(static_cast<double>(times[13])/static_cast<double>(repeats)) << std::endl; 00197 std::cout << std::endl; 00198 00199 return 0; 00200 }