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  *      published by the Free Software Foundation, version 3 of the License.
00010  *
00011  *      QVision is distributed in the hope that it will be useful,
00012  *      but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  *      GNU Lesser General Public License for more details.
00015  *
00016  *      You should have received a copy of the GNU Lesser General Public
00017  *      License along with QVision. If not, see <http://www.gnu.org/licenses/>.
00018  */
00019 
00023 
00024 #ifdef QVMATRIXALGEBRA_AVAILABLE
00025 
00026 #define CHECK_MKL_DSS_ERROR \
00027     if (error != MKL_DSS_SUCCESS) { \
00028         std::cout << "Solver returned error code " << error << std::endl; \
00029         exit(1); \
00030     }
00031 
00032 #ifndef QMATRIXALGEBRA_H
00033 #define QMATRIXALGEBRA_H
00034 
00035 #include <QV3DPointF>
00036 #include <qvmath/qvmatrix.h>
00037 #include <QVSparseBlockMatrix>
00038 
00039 // Substituted default values.
00040 //      #define DEFAULT_TQVSVD_METHOD                   GSL_THIN_DECOMP_MOD
00041 //      #define DEFAULT_TQVSV_METHOD                    LAPACK_ONLY_DGESVD
00042 //      #define DEFAULT_TQVCHOLESKY_METHOD              GSL_CHOLESKY
00043 //      #define DEFAULT_TQVQR_METHOD                    GSL_HOUSEHOLDER_THIN_QR
00044 //      #define DEFAULT_TQVLU_METHOD                    GSL_LU
00045 //      #define DEFAULT_TQVEIGENDECOMPOSITION_METHOD    GSL_EIGENSYMM
00046 //      #define DEFAULT_TQVEIGENVALUES_METHOD           GSL_EIGENSYMM_ONLY
00047 
00048 // Set the default values for the matrix algebra functions to use GSL or LAPACK, whichever is available.
00049 #ifdef LAPACK_AVAILABLE
00050         #define DEFAULT_TQVSVD_METHOD                   LAPACK_FULL_DGESVD              
00051         #define DEFAULT_TQVSV_METHOD                    LAPACK_ONLY_DGESVD
00052         #define DEFAULT_TQVCHOLESKY_METHOD              LAPACK_CHOLESKY_DPOTRF
00053         #define DEFAULT_TQVQR_METHOD                    LAPACK_FULL_DGEQR2
00054         #define DEFAULT_TQVLU_METHOD                    LAPACK_DGETRF
00055         #define DEFAULT_TQVEIGENDECOMPOSITION_METHOD    LAPACK_DSYEV
00056         #define DEFAULT_TQVEIGENVALUES_METHOD           LAPACK_DSYEV_ONLY
00057 #elif GSL_AVAILABLE
00058         #define DEFAULT_TQVSVD_METHOD                   GSL_THIN_DECOMP_MOD
00059         #define DEFAULT_TQVSV_METHOD                    GSL_ONLY_DECOMP
00060         #define DEFAULT_TQVCHOLESKY_METHOD              GSL_CHOLESKY
00061         #define DEFAULT_TQVQR_METHOD                    GSL_HOUSEHOLDER_THIN_QR
00062         #define DEFAULT_TQVLU_METHOD                    GSL_LU
00063         #define DEFAULT_TQVEIGENDECOMPOSITION_METHOD    GSL_EIGENSYMM
00064         #define DEFAULT_TQVEIGENVALUES_METHOD           GSL_EIGENSYMM_ONLY
00065 #endif
00066 
00067 #ifdef MKL_AVAILABLE
00068         #define DEFAULT_TQVSPARSESOLVE_METHOD           QVMKL_DSS
00069 #else
00070         #define DEFAULT_TQVSPARSESOLVE_METHOD           QV_DENSE
00071 #endif
00072 
00075 typedef enum {
00076     GSL_THIN_DECOMP_MOD,    
00077     GSL_THIN_DECOMP,        
00078     GSL_THIN_DECOMP_JACOBI, 
00079     LAPACK_FULL_DGESVD,     
00080     LAPACK_FULL_DGESDD,     
00081     LAPACK_THIN_DGESVD,     
00082     LAPACK_THIN_DGESDD      
00083 } TQVSVD_Method;
00084 
00087 typedef enum {
00088     GSL_ONLY_DECOMP_MOD,    
00089     GSL_ONLY_DECOMP,        
00090     GSL_ONLY_DECOMP_JACOBI, 
00091     LAPACK_ONLY_DGESVD,     
00092     LAPACK_ONLY_DGESDD      
00093 } TQVSV_Method;
00094 
00097 typedef enum {
00098     GSL_CHOLESKY,           
00099     LAPACK_CHOLESKY_DPOTRF  
00100 } TQVCholesky_Method;
00101 
00104 typedef enum {
00105     GSL_HOUSEHOLDER_THIN_QR,     
00106     GSL_HOUSEHOLDER_FULL_QR,     
00107     LAPACK_THIN_DGEQR2,           
00108     LAPACK_FULL_DGEQR2           
00109 } TQVQR_Method;
00110 
00113 typedef enum {
00114     GSL_LU,                 
00115     LAPACK_DGETRF           
00116 } TQVLU_Method;
00117 
00120 typedef enum {
00121     GSL_EIGENSYMM,         
00122     LAPACK_DSYEV           
00123 } TQVEigenDecomposition_Method;
00124 
00127 typedef enum {
00128     GSL_EIGENSYMM_ONLY,    
00129     LAPACK_DSYEV_ONLY     
00130 } TQVEigenValues_Method;
00131 
00134 typedef enum {
00136     QVMKL_DSS = 0,
00138     QVMKL_ISS = 1,
00140     QVCHOLMOD_DSS = 2,
00142     QV_SCG = 3,
00144     QV_BJPCG = 4,
00146     QV_DENSE = 5
00147 } TQVSparseSolve_Method;
00148 
00149 
00174 void singularValueDecomposition(const QVMatrix &M, QVMatrix &U, QVVector &s, QVMatrix &V, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00175 
00176 #ifndef DOXYGEN_IGNORE_THIS
00177 // This function is deprecated, and should be substituted by the former one.
00178 void SingularValueDecomposition(const QVMatrix &M, QVMatrix &U, QVMatrix &S, QVMatrix &V, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00179 #endif
00180 
00198 void solveFromSingularValueDecomposition(const QVMatrix &U, const QVVector &s, const QVMatrix &V, QVMatrix &X, const QVMatrix &B);
00199 
00217 void solveFromSingularValueDecomposition(const QVMatrix &U, const QVVector &s, const QVMatrix &V, QVVector &x, const QVVector &b);
00218 
00243 void solveBySingularValueDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, const TQVSVD_Method method  = DEFAULT_TQVSVD_METHOD);
00244 
00267 void solveBySingularValueDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, const TQVSVD_Method method  = DEFAULT_TQVSVD_METHOD);
00268 
00293 void solveBySingularValueDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B,
00294                                        QVMatrix &U, QVVector &s, QVMatrix &V, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00295 
00296 
00321 void solveBySingularValueDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b,
00322                                        QVMatrix &U, QVVector &s, QVMatrix &V, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
00323 
00338 double singularValueDecompositionResidual(const QVMatrix &M, const QVMatrix &U, const QVVector &s, const QVMatrix &V);
00339 
00353 void singularValues(const QVMatrix &M, QVVector &s, const TQVSV_Method method = DEFAULT_TQVSV_METHOD);
00354 
00367 double singularValuesResidual(const QVMatrix &M, const QVVector &s);
00368 
00369 
00390 void CholeskyDecomposition(const QVMatrix &M, QVMatrix &L, const TQVCholesky_Method method = DEFAULT_TQVCHOLESKY_METHOD);
00391 
00405 void solveFromCholeskyDecomposition(const QVMatrix &L, QVMatrix &X, const QVMatrix &B);
00406 
00419 void solveFromCholeskyDecomposition(const QVMatrix &L, QVVector &x, const QVVector &b);
00420 
00436 void solveByCholeskyDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, const TQVCholesky_Method method = DEFAULT_TQVCHOLESKY_METHOD);
00437 
00451 void solveByCholeskyDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, const TQVCholesky_Method method = DEFAULT_TQVCHOLESKY_METHOD);
00452 
00469 void solveByCholeskyDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, QVMatrix &L, const TQVCholesky_Method method = DEFAULT_TQVCHOLESKY_METHOD);
00470 
00487 void solveByCholeskyDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, QVMatrix &L, const TQVCholesky_Method method = DEFAULT_TQVCHOLESKY_METHOD);
00488 
00501 double CholeskyDecompositionResidual(const QVMatrix &M, const QVMatrix &L);
00502 
00503 
00527 void LUDecomposition(const QVMatrix &M, QVMatrix &P, QVMatrix &L, QVMatrix &U, const TQVLU_Method method = DEFAULT_TQVLU_METHOD);
00528 
00546 void solveFromLUDecomposition(const QVMatrix &P, const QVMatrix &L, const QVMatrix &U, QVMatrix &X, const QVMatrix &B);
00547 
00562 void solveFromLUDecomposition(const QVMatrix &P, const QVMatrix &L, const QVMatrix &U, QVVector &x, const QVVector &b);
00563 
00579 void solveByLUDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, const TQVLU_Method method = DEFAULT_TQVLU_METHOD);
00580 
00596 void solveByLUDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, const TQVLU_Method method = DEFAULT_TQVLU_METHOD);
00597 
00618 void solveByLUDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, QVMatrix &P, QVMatrix &L, QVMatrix &U, const TQVLU_Method method = DEFAULT_TQVLU_METHOD);
00619 
00639 void solveByLUDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b,QVMatrix &P, QVMatrix &L, QVMatrix &U, const TQVLU_Method method = DEFAULT_TQVLU_METHOD);
00640 
00655 double LUDecompositionResidual(const QVMatrix &M, const QVMatrix &P, const QVMatrix &L, const QVMatrix &U);
00656 
00657 
00680 void QRDecomposition(const QVMatrix &M, QVMatrix &Q, QVMatrix &R, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
00681 
00695 double QRDecompositionResidual(const QVMatrix &M, const QVMatrix &Q, const QVMatrix &R);
00696 
00721 void QLDecomposition(const QVMatrix &M, QVMatrix &Q, QVMatrix &L, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
00722 
00736 double QLDecompositionResidual(const QVMatrix &M, const QVMatrix &Q, const QVMatrix &L);
00737 
00762 void RQDecomposition(const QVMatrix &M, QVMatrix &R, QVMatrix &Q, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
00763 
00777 double RQDecompositionResidual(const QVMatrix &M, const QVMatrix &R, const QVMatrix &Q);
00778 
00779 
00804 void LQDecomposition(const QVMatrix &M, QVMatrix &L, QVMatrix &Q, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
00805 
00819 double LQDecompositionResidual(const QVMatrix &M, const QVMatrix &L, const QVMatrix &Q);
00820 
00837 void solveFromQRDecomposition(const QVMatrix &Q, const QVMatrix &R, QVMatrix &X, const QVMatrix &B);
00838 
00854 void solveFromQRDecomposition(const QVMatrix &Q, const QVMatrix &R, QVVector &x, const QVVector &b);
00855 
00874 void solveByQRDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
00875 
00892 void solveByQRDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
00893 
00914 void solveByQRDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B,
00915                             QVMatrix &Q, QVMatrix &R, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
00916 
00917 
00938 void solveByQRDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b,
00939                             QVMatrix &Q, QVMatrix &R, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
00940 
00964 void eigenDecomposition(const QVMatrix &M, QVVector &lambda, QVMatrix &Q, const TQVEigenDecomposition_Method method = DEFAULT_TQVEIGENDECOMPOSITION_METHOD);
00965 
00979 double eigenDecompositionResidual(const QVMatrix &M, const QVVector &lambda, const QVMatrix &Q);
00980 
00994 void eigenValues(const QVMatrix &M, QVVector &lambda, const TQVEigenValues_Method method = DEFAULT_TQVEIGENVALUES_METHOD);
00995 
01008 double eigenValuesResidual(const QVMatrix &M, const QVVector &lambda);
01009 
01010 
01022 QVMatrix pseudoInverse(const QVMatrix &M, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD, const double epsilon = 1.0E-10);
01023 
01038 double determinant(const QVMatrix &M, const TQVLU_Method method = DEFAULT_TQVLU_METHOD);
01039 
01056 void solveHomogeneous(const QVMatrix &A, QVector<double> &x, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
01057 
01073 double solveResidual(const QVMatrix &M, const QVMatrix &X, const QVMatrix &B);
01074 
01090 double solveResidual(const QVMatrix &M, const QVVector &x, const QVVector &b);
01091 
01092 #ifndef DOXYGEN_IGNORE_THIS
01093 extern int dummy; // For optional output parameter of next SparseSolve function:
01094 #endif
01095 
01123 double sparseSolve(const QVSparseBlockMatrix &M, QVVector &x, const QVVector &b,
01124                  const bool isSymmetric = false, const bool isPosDefinite = false,
01125                  const TQVSparseSolve_Method method = DEFAULT_TQVSPARSESOLVE_METHOD, const bool start_from_x = false,
01126                  const bool iters_or_resid = true, const int iters = 0, const double resid = 1.0E-10,
01127                  int &final_iter_count = dummy);
01128 
01150 bool solveHomogeneous(const QVSparseBlockMatrix &A, QVVector &x, const int maxIterations = 10, const double minRelativeError = 0.0, const TQVSparseSolve_Method method=QVMKL_DSS);
01151 
01152 #ifndef DOXYGEN_IGNORE_THIS
01153 void cold_start_mkl_initialization(const int size = 2);
01154 
01155 double solveConjugateGradient(const QVSparseBlockMatrix &A, QVVector &x, const QVVector &b, const int maxIters, const int minIters = 0, const double minAbsoluteError = 0.0);
01156 
01157 double sparseSolve(const MKLPardisoSparseFormat &pardisomatrix, QVVector &x, const QVVector &b,
01158                  const bool isSymmetric = false, const bool isPosDefinite = false,
01159                  const TQVSparseSolve_Method method = DEFAULT_TQVSPARSESOLVE_METHOD, const bool start_from_x = false,
01160                  const bool iters_or_resid = true, const int iters = 0, const double resid = 1.0E-10,
01161                  int &final_iter_count = dummy);
01162 
01163 // Internal use only. Fast function for homogeneous solving.
01164 #ifdef GSL_AVAILABLE
01165 void solveHomogeneousEig(const QVMatrix &M, QVVector &result);
01166 #endif
01167 
01169 void solveLinear(const QVMatrix &A, QVVector &x, const QVVector &b);
01170 
01172 void solveLinear(const QVMatrix &A, QVMatrix &X, const QVMatrix &B);
01173 
01175 void solveOverDetermined(const QVMatrix &A, QVMatrix &X, const QVMatrix &B);
01176 
01178 void solveHomogeneousLinear(const QVMatrix &A, QVector<double> &x);
01179 
01181 double homogLineFromMoments(double x,double y,double xx,double xy,double yy,double &a,double &b,double &c);
01182 
01184 QVVector regressionLine(const QVMatrix &points);
01185 
01186 // Deprecated.
01187 void SingularValueDecomposition(const QVMatrix &M, QVMatrix &U, QVVector &s, QVMatrix &V, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
01188 void SolveFromSingularValueDecomposition(const QVMatrix &U, const QVVector &s, const QVMatrix &V, QVMatrix &X, const QVMatrix &B);
01189 void SolveFromSingularValueDecomposition(const QVMatrix &U, const QVVector &s, const QVMatrix &V, QVVector &x, const QVVector &b);
01190 void SolveBySingularValueDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, const TQVSVD_Method method  = DEFAULT_TQVSVD_METHOD);
01191 void SolveBySingularValueDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, const TQVSVD_Method method  = DEFAULT_TQVSVD_METHOD);
01192 void SolveBySingularValueDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, QVMatrix &U, QVVector &s, QVMatrix &V, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
01193 void SolveBySingularValueDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, QVMatrix &U, QVVector &s, QVMatrix &V, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
01194 double SingularValueDecompositionResidual(const QVMatrix &M, const QVMatrix &U, const QVVector &s, const QVMatrix &V);
01195 void SingularValues(const QVMatrix &M, QVVector &s, const TQVSV_Method method = DEFAULT_TQVSV_METHOD);
01196 double SingularValuesResidual(const QVMatrix &M, const QVVector &s);
01197 void SolveFromCholeskyDecomposition(const QVMatrix &L, QVMatrix &X, const QVMatrix &B);
01198 void SolveFromCholeskyDecomposition(const QVMatrix &L, QVVector &x, const QVVector &b);
01199 void SolveByCholeskyDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, const TQVCholesky_Method method = DEFAULT_TQVCHOLESKY_METHOD);
01200 void SolveByCholeskyDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, const TQVCholesky_Method method = DEFAULT_TQVCHOLESKY_METHOD);
01201 void SolveByCholeskyDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, QVMatrix &L, const TQVCholesky_Method method = DEFAULT_TQVCHOLESKY_METHOD);
01202 void SolveByCholeskyDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, QVMatrix &L, const TQVCholesky_Method method = DEFAULT_TQVCHOLESKY_METHOD);
01203 void SolveFromLUDecomposition(const QVMatrix &P, const QVMatrix &L, const QVMatrix &U, QVMatrix &X, const QVMatrix &B);
01204 void SolveFromLUDecomposition(const QVMatrix &P, const QVMatrix &L, const QVMatrix &U, QVVector &x, const QVVector &b);
01205 void SolveByLUDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, const TQVLU_Method method = DEFAULT_TQVLU_METHOD);
01206 void SolveByLUDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, const TQVLU_Method method = DEFAULT_TQVLU_METHOD);
01207 void SolveByLUDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, QVMatrix &P, QVMatrix &L, QVMatrix &U, const TQVLU_Method method = DEFAULT_TQVLU_METHOD);
01208 void SolveByLUDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b,QVMatrix &P, QVMatrix &L, QVMatrix &U, const TQVLU_Method method = DEFAULT_TQVLU_METHOD);
01209 void SolveFromQRDecomposition(const QVMatrix &Q, const QVMatrix &R, QVMatrix &X, const QVMatrix &B);
01210 void SolveFromQRDecomposition(const QVMatrix &Q, const QVMatrix &R, QVVector &x, const QVVector &b);
01211 void SolveByQRDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
01212 void SolveByQRDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
01213 void SolveByQRDecomposition(const QVMatrix &M, QVMatrix &X, const QVMatrix &B, QVMatrix &Q, QVMatrix &R, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
01214 void SolveByQRDecomposition(const QVMatrix &M, QVVector &x, const QVVector &b, QVMatrix &Q, QVMatrix &R, const TQVQR_Method method = DEFAULT_TQVQR_METHOD);
01215 void EigenDecomposition(const QVMatrix &M, QVVector &lambda, QVMatrix &Q, const TQVEigenDecomposition_Method method = DEFAULT_TQVEIGENDECOMPOSITION_METHOD);
01216 double EigenDecompositionResidual(const QVMatrix &M, const QVVector &lambda, const QVMatrix &Q);
01217 void EigenValues(const QVMatrix &M, QVVector &lambda, const TQVEigenValues_Method method = DEFAULT_TQVEIGENVALUES_METHOD);
01218 double EigenValuesResidual(const QVMatrix &M, const QVVector &lambda);
01219 void SolveHomogeneous(const QVMatrix &A, QVector<double> &x, const TQVSVD_Method method = DEFAULT_TQVSVD_METHOD);
01220 double SolveResidual(const QVMatrix &M, const QVMatrix &X, const QVMatrix &B);
01221 double SolveResidual(const QVMatrix &M, const QVVector &x, const QVVector &b);
01222 double SparseSolve(const QVSparseBlockMatrix M, QVVector &x, const QVVector b,
01223                  const bool isSymmetric = false, const bool isPosDefinite = false,
01224                  const TQVSparseSolve_Method method = DEFAULT_TQVSPARSESOLVE_METHOD, const bool start_from_x = false,
01225                  const bool iters_or_resid = true, const int iters = 0, const double resid = 1.0E-10,
01226                  int &final_iter_count = dummy);
01227 void SolveHomogeneous(const QVSparseBlockMatrix &A, QVVector &x, const int maxIterations = 10, const double minRelativeError = 0.0, const TQVSparseSolve_Method method=QVMKL_DSS);
01228 #endif
01229 
01230 #endif // QMATRIXALGEBRA_H
01231 
01232 #endif // QVMATRIXALGEBRA_AVAILABLE
src/qvmath/qvmatrixalgebra.h
