eigen3_decompositions.cpp
Go to the documentation of this file.
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 }


ecl_core_apps
Author(s): Daniel Stonier (d.stonier@gmail.com)
autogenerated on Thu Jan 2 2014 11:13:30