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 Singular Value Decomposition 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_THIN_DECOMP_MOD",
00052                      ".............SVD/SV method (available methods follow):\n"
00053                      "                GSL_THIN_DECOMP_MOD | GSL_THIN_DECOMP | GSL_THIN_DECOMP_JACOBI |\n"
00054                      "                LAPACK_FULL_DGESVD | LAPACK_FULL_DGESDD | LAPACK_THIN_DGESVD | \n"
00055                      "                LAPACK_FULL_DGESDD ....................................(for SVD)\n"
00056                      "                GSL_ONLY_DECOMP_MOD | GSL_ONLY_DECOMP | GSL_ONLY_DECOMP_JACOBI |\n"
00057                      "                LAPACK_ONLY_DGESVD | LAPACK_ONLY_DGESDD ...........(for SV only)");
00058     arg_container.addProperty<bool>("verbose",QVPropertyContainer::inputFlag,false,
00059                                    "If true, print involved matrices and vectors");
00060     arg_container.addProperty<bool>("show_residuals",QVPropertyContainer::inputFlag,false,
00061                                    "If true, print test residuals");
00062 
00063     // Process command line (and check for help or incorrect input parameters):
00064     int ret_value = app.processArguments();
00065     if(ret_value != 1) exit(ret_value);
00066 
00067     // If parameters OK, read possible parameters from command line:
00068     const int Rows = arg_container.getPropertyValue<int>("Rows");
00069     const int Cols = arg_container.getPropertyValue<int>("Cols");
00070     int RHS_size = arg_container.getPropertyValue<int>("RHS_size");
00071     const int n_tests = arg_container.getPropertyValue<int>("n_tests");
00072     const int verbose = arg_container.getPropertyValue<bool>("verbose");
00073     int show_residuals = arg_container.getPropertyValue<bool>("show_residuals");
00074     const QString method = arg_container.getPropertyValue<QString>("method");
00075     TQVSVD_Method svd_method;
00076     TQVSV_Method sv_method;
00077     bool only_sv = false;
00078     if(method == "GSL_THIN_DECOMP_MOD")
00079         svd_method = GSL_THIN_DECOMP_MOD;
00080     else if(method == "GSL_THIN_DECOMP")
00081         svd_method = GSL_THIN_DECOMP;
00082     else if(method == "GSL_THIN_DECOMP_JACOBI")
00083         svd_method = GSL_THIN_DECOMP_JACOBI;
00084     else if(method == "LAPACK_FULL_DGESVD")
00085         svd_method = LAPACK_FULL_DGESVD;
00086     else if(method == "LAPACK_FULL_DGESDD")
00087         svd_method = LAPACK_FULL_DGESDD;
00088     else if(method == "LAPACK_THIN_DGESVD")
00089         svd_method = LAPACK_THIN_DGESVD;
00090     else if(method == "LAPACK_THIN_DGESDD")
00091         svd_method = LAPACK_THIN_DGESDD;
00092     else if(method == "LAPACK_ONLY_DGESVD") {
00093         sv_method = LAPACK_ONLY_DGESVD;
00094         only_sv = true;
00095     } else if(method == "GSL_ONLY_DECOMP_MOD") {
00096         sv_method = GSL_ONLY_DECOMP_MOD;
00097         only_sv = true;
00098     } else if(method == "GSL_ONLY_DECOMP") {
00099         sv_method = GSL_ONLY_DECOMP;
00100         only_sv = true;
00101     } else if(method == "GSL_ONLY_DECOMP_JACOBI") {
00102         sv_method = GSL_ONLY_DECOMP_JACOBI;
00103         only_sv = true;
00104     } else if(method == "LAPACK_ONLY_DGESVD") {
00105         sv_method = LAPACK_ONLY_DGESVD;
00106         only_sv = true;
00107     } else if(method == "LAPACK_ONLY_DGESDD") {
00108         sv_method = LAPACK_ONLY_DGESDD;
00109         only_sv = true;
00110     } else {
00111         std::cout << "Incorrect SVD/SV method. Use --help to see available methods.\n";
00112         exit(-1);
00113     }
00114 
00115     if(only_sv) {
00116         // Only singular values asked for; no solving or residuals computed:
00117         if(RHS_size != 0) {
00118             std::cout << "Only singular values asked for; RHS_size forced 0.\n";
00119             RHS_size = 0;
00120         }
00121     }
00122 
00123     std::cout << "Using values: Rows=" << Rows  << " Cols=" << Cols << " RHS_size=" << RHS_size
00124               << " n_tests=" << n_tests << " SVD method=" << qPrintable(method) << "\n";
00125 
00126     double total_ms = 0.0;
00127 
00128     for(int i=0;i<n_tests;i++) {
00129 
00130         QVMatrix M = QVMatrix::random(Rows,Cols), U, V, X, B, X_sol;
00131         QVVector s, x, b, x_sol;
00132 
00133         if(not only_sv) {
00134             // We generate RHS_size random solutions, and the corresponding right hand sides (for testing):
00135             if(qAbs(RHS_size) == 1) {
00136                 x_sol = QVVector::random(Cols);
00137                 b = M * x_sol;
00138             } else if (qAbs(RHS_size) > 1) {
00139                 X_sol = QVMatrix::random(Cols,qAbs(RHS_size));
00140                 B = M * X_sol;
00141             }
00142         }
00143 
00144         QTime t;
00145 
00146         t.start();
00147 
00148         if(only_sv)
00149             SingularValues(M, s, sv_method);
00150         else if(RHS_size == 0)
00151             SingularValueDecomposition(M, U, s, V, svd_method);
00152         else if(RHS_size == 1)
00153             SolveBySingularValueDecomposition(M, x, b, svd_method);
00154         else if(RHS_size > 1)
00155             SolveBySingularValueDecomposition(M, X, B, svd_method);
00156         else if (RHS_size == -1)
00157             SolveBySingularValueDecomposition(M, x, b, U, s, V, svd_method);
00158         else if (RHS_size < -1)
00159             SolveBySingularValueDecomposition(M, X, B, U, s, V, svd_method);
00160 
00161          total_ms += t.elapsed();
00162 
00163         if(verbose) {
00164             if(only_sv) {
00165                 std::cout << "*****************************************\n";
00166                 std::cout << "M:" << M << "\n";
00167                 std::cout << "s:" << s << "\n";
00168                 std::cout << "*****************************************\n";
00169             } else {
00170                 if(RHS_size <= 0) {
00171                     std::cout << "*****************************************\n";
00172                     std::cout << "M:" << M << "\n";
00173                     std::cout << "U:" << U << "\n";
00174                     std::cout << "s:" << s << "\n";
00175                     std::cout << "V:" << V << "\n";
00176                     QVMatrix S = QVMatrix::diagonal(s);
00177                     const int m = U.getCols(), n = V.getCols();
00178                     if (m>n)
00179                         S = S & QVMatrix(m-n,n,0.0);
00180                     else if (m<n)
00181                         S = S | QVMatrix(m,n-m,0.0);
00182                     std::cout << "U diag(s) V^T:" << U*S*V.transpose() << "\n";
00183                     std::cout << "U^T U:" << U.transpose()*U << "\n";
00184                     std::cout << "V^T V:" << V.transpose()*V << "\n";
00185                     std::cout << "*****************************************\n";
00186                 }
00187                 if(qAbs(RHS_size) == 1) {
00188                     std::cout << "*****************************************\n";
00189                     std::cout << "M:" << M << "\n";
00190                     std::cout << "b:" << b << "\n";
00191                     std::cout << "x:" << x << "\n";
00192                     std::cout << "x_sol:" << x_sol << "\n";
00193                     std::cout << "*****************************************\n";
00194                 }
00195                 if(qAbs(RHS_size) > 1){
00196                     std::cout << "*****************************************\n";
00197                     std::cout << "M:" << M << "\n";
00198                     std::cout << "B:" << B << "\n";
00199                     std::cout << "X:" << X << "\n";
00200                     std::cout << "X_sol:" << X_sol << "\n";
00201                     std::cout << "*****************************************\n";
00202                 }
00203             }
00204         }        
00205 
00206         if(show_residuals) {
00207             if(only_sv) {
00208                 std::cout << "SV residual = " << SingularValuesResidual(M, s) << "\n";
00209             } else {
00210                 if(RHS_size <= 0) {
00211                     std::cout << "SVD residual = " << SingularValueDecompositionResidual(M, U, s, V) << "\n";
00212                 }
00213                 if(qAbs(RHS_size) == 1) {
00214                     std::cout << "Test solve residual = " << SolveResidual(M, x_sol, b) << "\n";
00215                     std::cout << "SVD solve residual = " << SolveResidual(M, x, b) << "\n";
00216                     std::cout << "x - x_sol residual = " << (x-x_sol).norm2() << "\n";
00217                     std::cout << "x norm = " << x.norm2() << "\n";
00218                     std::cout << "x_sol norm = " << x_sol.norm2() << "\n";
00219                 }
00220                 if(qAbs(RHS_size) > 1) {
00221                     std::cout << "Test solve residual = " << SolveResidual(M, X_sol, B) << "\n";
00222                     std::cout << "SVD solve residual = " << SolveResidual(M, X, B) << "\n";
00223                 }
00224             }
00225         }
00226     }
00227 
00228     total_ms /= n_tests;
00229 
00230     if(n_tests==1)
00231         std::cout << "Total time: " << total_ms << " ms.\n";
00232     else
00233         std::cout << "Average total time: " << total_ms << " ms.\n";
00234 
00235     std::cout << "Finished.\n";    
00236 }
examples/matrixalgebra-tests/SVD-test/SVD-test.cpp
