Visual Servoing Platform  version 3.2.1 under development (2019-05-26)
vpMatrix_svd.cpp
1 /****************************************************************************
2  *
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2019 by Inria. All rights reserved.
5  *
6  * This software is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  * See the file LICENSE.txt at the root directory of this source
11  * distribution for additional information about the GNU GPL.
12  *
13  * For using ViSP with software that can not be combined with the GNU
14  * GPL, please contact Inria about acquiring a ViSP Professional
15  * Edition License.
16  *
17  * See http://visp.inria.fr for more information.
18  *
19  * This software was developed at:
20  * Inria Rennes - Bretagne Atlantique
21  * Campus Universitaire de Beaulieu
22  * 35042 Rennes Cedex
23  * France
24  *
25  * If you have questions regarding the use of this file, please contact
26  * Inria at visp@inria.fr
27  *
28  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30  *
31  * Description:
32  * Matrix SVD decomposition.
33  *
34  * Authors:
35  * Eric Marchand
36  *
37  *****************************************************************************/
38 
39 #include <visp3/core/vpColVector.h>
40 #include <visp3/core/vpConfig.h>
41 #include <visp3/core/vpDebug.h>
42 #include <visp3/core/vpException.h>
43 #include <visp3/core/vpMath.h>
44 #include <visp3/core/vpMatrix.h>
45 #include <visp3/core/vpMatrixException.h>
46 
47 #include <cmath> // std::fabs
48 #include <iostream>
49 #include <limits> // numeric_limits
50 
51 #ifdef VISP_HAVE_EIGEN3
52 #include <Eigen/SVD>
53 #endif
54 
55 #ifdef VISP_HAVE_GSL
56 #include <gsl/gsl_linalg.h>
57 #endif
58 
59 #if (VISP_HAVE_OPENCV_VERSION >= 0x020101) // Require opencv >= 2.1.1
60 #include <opencv2/core/core.hpp>
61 #endif
62 
63 #ifdef VISP_HAVE_LAPACK
64 # ifdef VISP_HAVE_MKL
65 #include <mkl.h>
66 typedef MKL_INT integer;
67 # else
68 # if defined(VISP_HAVE_LAPACK_BUILT_IN)
69 typedef long int integer;
70 # else
71 typedef int integer;
72 # endif
73 extern "C" int dgesdd_(char *jobz, integer *m, integer *n, double *a, integer *lda, double *s, double *u, integer *ldu,
74  double *vt, integer *ldvt, double *work, integer *lwork, integer *iwork, integer *info);
75 
76 #include <stdio.h>
77 #include <string.h>
78 # endif
79 #endif
80 /*---------------------------------------------------------------------
81 
82 SVD related functions
83 
84 ---------------------------------------------------------------------*/
85 
86 #if (VISP_HAVE_OPENCV_VERSION >= 0x020101) // Require opencv >= 2.1.1
87 
155 {
156  int rows = (int)this->getRows();
157  int cols = (int)this->getCols();
158  cv::Mat m(rows, cols, CV_64F, this->data);
159  cv::SVD opencvSVD(m);
160  cv::Mat opencvV = opencvSVD.vt;
161  cv::Mat opencvW = opencvSVD.w;
162  V.resize((unsigned int)opencvV.rows, (unsigned int)opencvV.cols);
163  w.resize((unsigned int)(opencvW.rows * opencvW.cols));
164 
165  memcpy(V.data, opencvV.data, (size_t)(8 * opencvV.rows * opencvV.cols));
166  V = V.transpose();
167  memcpy(w.data, opencvW.data, (size_t)(8 * opencvW.rows * opencvW.cols));
168  this->resize((unsigned int)opencvSVD.u.rows, (unsigned int)opencvSVD.u.cols);
169  memcpy(this->data, opencvSVD.u.data, (size_t)(8 * opencvSVD.u.rows * opencvSVD.u.cols));
170 }
171 
172 #endif
173 
174 #if defined(VISP_HAVE_LAPACK)
175 
242 {
243  w.resize(this->getCols());
244  V.resize(this->getCols(), this->getCols());
245 
246  integer m = (integer)(this->getCols());
247  integer n = (integer)(this->getRows());
248  integer lda = m;
249  integer ldu = m;
250  integer ldvt = (std::min)(m, n);
251  integer info, lwork;
252 
253  double wkopt;
254  double *work;
255 
256  integer *iwork = new integer[8 * static_cast<integer>((std::min)(n, m))];
257 
258  double *s = w.data;
259  double *a = new double[static_cast<unsigned int>(lda * n)];
260  memcpy(a, this->data, this->getRows() * this->getCols() * sizeof(double));
261  double *u = V.data;
262  double *vt = this->data;
263 
264  lwork = -1;
265  dgesdd_((char *)"S", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork, iwork, &info);
266  lwork = (int)wkopt;
267  work = new double[static_cast<unsigned int>(lwork)];
268 
269  dgesdd_((char *)"S", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, iwork, &info);
270 
271  if (info > 0) {
272  throw(vpMatrixException(vpMatrixException::fatalError, "The algorithm computing SVD failed to converge."));
273  }
274 
275  V = V.transpose();
276  delete[] work;
277  delete[] iwork;
278  delete[] a;
279 }
280 #endif
281 
282 #ifdef VISP_HAVE_GSL
283 
351 {
352  w.resize(this->getCols());
353  V.resize(this->getCols(), this->getCols());
354 
355  unsigned int nc = getCols();
356  unsigned int nr = getRows();
357  gsl_vector *work = gsl_vector_alloc(nc);
358 
359  gsl_matrix A;
360  A.size1 = nr;
361  A.size2 = nc;
362  A.tda = A.size2;
363  A.data = this->data;
364  A.owner = 0;
365  A.block = 0;
366 
367  gsl_matrix V_;
368  V_.size1 = nc;
369  V_.size2 = nc;
370  V_.tda = V_.size2;
371  V_.data = V.data;
372  V_.owner = 0;
373  V_.block = 0;
374 
375  gsl_vector S;
376  S.size = nc;
377  S.stride = 1;
378  S.data = w.data;
379  S.owner = 0;
380  S.block = 0;
381 
382  gsl_linalg_SV_decomp(&A, &V_, &S, work);
383 
384  gsl_vector_free(work);
385 }
386 #endif // # #GSL
387 
388 #ifdef VISP_HAVE_EIGEN3
389 
456 {
457  w.resize(this->getCols());
458  V.resize(this->getCols(), this->getCols());
459 
460  Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > M(this->data, this->getRows(),
461  this->getCols());
462 
463  Eigen::JacobiSVD<Eigen::MatrixXd> svd(M, Eigen::ComputeThinU | Eigen::ComputeThinV);
464 
465  Eigen::Map<Eigen::VectorXd> w_(w.data, w.size());
466  Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > V_(V.data, V.getRows(),
467  V.getCols());
468  Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > U_(this->data, this->getRows(),
469  this->getCols());
470  w_ = svd.singularValues();
471  V_ = svd.matrixV();
472  U_ = svd.matrixU();
473 }
474 #endif
void svd(vpColVector &w, vpMatrix &V)
Definition: vpMatrix.cpp:1881
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:164
void resize(const unsigned int nrows, const unsigned int ncols, const bool flagNullify=true, const bool recopy_=true)
Definition: vpArray2D.h:305
double * data
Address of the first element of the data array.
Definition: vpArray2D.h:145
unsigned int size() const
Return the number of elements of the 2D array.
Definition: vpArray2D.h:291
unsigned int getCols() const
Definition: vpArray2D.h:279
void svdLapack(vpColVector &w, vpMatrix &V)
vpMatrix transpose() const
Definition: vpMatrix.cpp:394
void svdOpenCV(vpColVector &w, vpMatrix &V)
unsigned int getRows() const
Definition: vpArray2D.h:289
void svdGsl(vpColVector &w, vpMatrix &V)
Implementation of column vector and the associated operations.
Definition: vpColVector.h:130
error that can be emited by the vpMatrix class and its derivates
void svdEigen3(vpColVector &w, vpMatrix &V)
void resize(const unsigned int i, const bool flagNullify=true)
Definition: vpColVector.h:310