QVision: Qt's Image, Video and Computer Vision Library

00001 /*
00002  *      Copyright (C) 2007, 2008, 2009, 2010, 2011, 2012. PARP Research Group.
00003  *      <http://perception.inf.um.es>
00004  *      University of Murcia, Spain.
00005  *
00006  *      This file is part of the QVision library.
00007  *
00008  *      QVision is free software: you can redistribute it and/or modify
00009  *      it under the terms of the GNU Lesser General Public License as
00010  *      published by the Free Software Foundation, version 3 of the License.
00011  *
00012  *      QVision is distributed in the hope that it will be useful,
00013  *      but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  *      GNU Lesser General Public License for more details.
00016  *
00017  *      You should have received a copy of the GNU Lesser General Public
00018  *      License along with QVision. If not, see <http://www.gnu.org/licenses/>.
00019  */
00020 
00024 
00025 #include <iostream>
00026 
00027 #include <QTime>
00028 
00029 #include <QVApplication>
00030 #include <QVPropertyContainer>
00031 #include <QVMatrix>
00032 #include <QVVector>
00033 
00034 int main(int argc, char *argv[])
00035 {
00036     // QVApplication object:
00037     QVApplication app(argc,argv,"Performance test for QR methods in QVision",false);
00038 
00039     // Container with command line parameters:
00040     QVPropertyContainer arg_container(argv[0]);
00041     arg_container.addProperty<int>("Rows",QVPropertyContainer::inputFlag,200,
00042                                    "Number of rows for the test matrix",1,100000);
00043     arg_container.addProperty<int>("Cols",QVPropertyContainer::inputFlag,100,
00044                                    "Number of columns for the test matrix",1,100000);
00045     arg_container.addProperty<int>("RHS_size",QVPropertyContainer::inputFlag,0,
00046                      ".......................Number of right hand sides to solve for:\n"
00047                      "                                    "
00048                      "If 0, only decompose, if <0 solve and also return decomposition)",-1000,1000);
00049     arg_container.addProperty<int>("n_tests",QVPropertyContainer::inputFlag,1,
00050                                    "Number of tests to average execution time",1,100);
00051     arg_container.addProperty<QString>("method",QVPropertyContainer::inputFlag,"GSL_HOUSEHOLDER_THIN_QR",
00052                      "..................QR method (available methods follow):\n"
00053                      "         GSL_HOUSEHOLDER_THIN_QR | GSL_HOUSEHOLDER_FULL_QR | LAPACK_THIN_DGEQR2 | LAPACK_FULL_DGEQR2");
00054     arg_container.addProperty<bool>("verbose",QVPropertyContainer::inputFlag,false,
00055                                    "If true, print involved matrices and vectors");
00056     arg_container.addProperty<bool>("show_residuals",QVPropertyContainer::inputFlag,false,
00057                                    "If true, print test residuals");
00058 
00059     // Process command line (and check for help or incorrect input parameters):
00060     int ret_value = app.processArguments();
00061     if(ret_value != 1) exit(ret_value);
00062 
00063     // If parameters OK, read possible parameters from command line:
00064     const int Rows = arg_container.getPropertyValue<int>("Rows");
00065     const int Cols = arg_container.getPropertyValue<int>("Cols");
00066     int RHS_size = arg_container.getPropertyValue<int>("RHS_size");
00067     const int n_tests = arg_container.getPropertyValue<int>("n_tests");
00068     const int verbose = arg_container.getPropertyValue<bool>("verbose");
00069     int show_residuals = arg_container.getPropertyValue<bool>("show_residuals");
00070     const QString method = arg_container.getPropertyValue<QString>("method");
00071     TQVQR_Method qr_method;
00072     if(method == "GSL_HOUSEHOLDER_THIN_QR")
00073         qr_method = GSL_HOUSEHOLDER_THIN_QR;
00074     else if(method == "GSL_HOUSEHOLDER_FULL_QR")
00075         qr_method = GSL_HOUSEHOLDER_FULL_QR;
00076     else if(method == "LAPACK_THIN_DGEQR2")
00077         qr_method = LAPACK_THIN_DGEQR2;
00078     else if(method == "LAPACK_FULL_DGEQR2")
00079         qr_method = LAPACK_FULL_DGEQR2;
00080     else {
00081         std::cout << "Incorrect QR method. Use --help to see available methods.\n";
00082         exit(-1);
00083     }
00084 
00085 
00086     std::cout << "Using values: Rows=" << Rows  << " Cols=" << Cols << " RHS_size=" << RHS_size
00087               << " n_tests=" << n_tests << " QR method=" << qPrintable(method) << "\n";
00088 
00089     double total_ms = 0.0;
00090 
00091     for(int i=0;i<n_tests;i++) {
00092 
00093         QVMatrix M = QVMatrix::random(Rows,Cols), Q, R, X, B, X_sol;
00094         QVVector x, b, x_sol;
00095 
00096         // We generate RHS_size random solutions, and the corresponding right hand sides (for testing):
00097         if(qAbs(RHS_size) == 1) {
00098             x_sol = QVVector::random(Cols);
00099             b = M * x_sol;
00100         } else if (qAbs(RHS_size) > 1) {
00101             X_sol = QVMatrix::random(Cols,qAbs(RHS_size));
00102             B = M * X_sol;
00103         }
00104 
00105         QTime t;
00106 
00107         t.start();
00108 
00109         if(RHS_size == 0)
00110             QRDecomposition(M, Q, R, qr_method);
00111         else if(RHS_size == 1)
00112             SolveByQRDecomposition(M, x, b, qr_method);
00113         else if(RHS_size > 1)
00114             SolveByQRDecomposition(M, X, B, qr_method);
00115         else if (RHS_size == -1)
00116             SolveByQRDecomposition(M, x, b, Q, R, qr_method);
00117         else if (RHS_size < -1)
00118             SolveByQRDecomposition(M, X, B, Q, R, qr_method);
00119 
00120         total_ms += t.elapsed();
00121 
00122         if(verbose) {
00123             if(RHS_size <= 0) {
00124                     std::cout << "*****************************************\n";
00125                     std::cout << "M:" << M << "\n";
00126                     std::cout << "Q:" << Q << "\n";
00127                     std::cout << "R:" << R << "\n";
00128                     std::cout << "Q R:" << Q*R << "\n";
00129                     std::cout << "Q^T Q:" << Q.transpose()*Q << "\n";
00130                     std::cout << "*****************************************\n";
00131             }
00132             if(qAbs(RHS_size) == 1) {
00133                     std::cout << "*****************************************\n";
00134                     std::cout << "M:" << M << "\n";
00135                     std::cout << "b:" << b << "\n";
00136                     std::cout << "x:" << x << "\n";
00137                     std::cout << "x_sol:" << x_sol << "\n";
00138                     std::cout << "*****************************************\n";
00139             }
00140             if(qAbs(RHS_size) > 1){
00141                     std::cout << "*****************************************\n";
00142                     std::cout << "M:" << M << "\n";
00143                     std::cout << "B:" << B << "\n";
00144                     std::cout << "X:" << X << "\n";
00145                     std::cout << "X_sol:" << X_sol << "\n";
00146                     std::cout << "*****************************************\n";
00147             }
00148         }
00149 
00150         if(show_residuals) {           
00151             if(RHS_size <= 0) {
00152                 std::cout << "QR residual = " << QRDecompositionResidual(M, Q, R) << "\n";
00153             }
00154             if(qAbs(RHS_size) == 1) {
00155                 std::cout << "Test solve residual = " << SolveResidual(M, x_sol, b) << "\n";
00156                 std::cout << "QR solve residual = " << SolveResidual(M, x, b) << "\n";
00157                 std::cout << "x - x_sol residual = " << (x-x_sol).norm2() << "\n";
00158                 std::cout << "x norm = " << x.norm2() << "\n";
00159                 std::cout << "x_sol norm = " << x_sol.norm2() << "\n";
00160             }
00161             if(qAbs(RHS_size) > 1) {
00162                 std::cout << "Test solve residual = " << SolveResidual(M, X_sol, B) << "\n";
00163                 std::cout << "QR solve residual = " << SolveResidual(M, X, B) << "\n";
00164             }
00165         }
00166     }
00167 
00168     total_ms /= n_tests;
00169 
00170     if(n_tests==1)
00171         std::cout << "Total time: " << total_ms << " ms.\n";
00172     else
00173         std::cout << "Average total time: " << total_ms << " ms.\n";
00174 
00175     std::cout << "Finished.\n";    
00176 }
examples/matrixalgebra-tests/QR-test/QR-test.cpp
