Visual Servoing Platform  version 3.6.1 under development (2024-04-18)
vpMatrix_svd.cpp
1 /****************************************************************************
2  *
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2023 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 https://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 *****************************************************************************/
35 
36 #include <visp3/core/vpColVector.h>
37 #include <visp3/core/vpConfig.h>
38 #include <visp3/core/vpDebug.h>
39 #include <visp3/core/vpException.h>
40 #include <visp3/core/vpMath.h>
41 #include <visp3/core/vpMatrix.h>
42 #include <visp3/core/vpMatrixException.h>
43 
44 #include <cmath> // std::fabs
45 #include <iostream>
46 #include <limits> // numeric_limits
47 
48 #ifdef VISP_HAVE_EIGEN3
49 #include <Eigen/SVD>
50 #endif
51 
52 #if defined(VISP_HAVE_OPENCV) // Require opencv >= 2.1.1
53 #include <opencv2/core/core.hpp>
54 #endif
55 
56 #ifdef VISP_HAVE_LAPACK
57 #ifdef VISP_HAVE_GSL
58 #include <gsl/gsl_linalg.h>
59 #endif
60 #ifdef VISP_HAVE_MKL
61 #include <mkl.h>
62 typedef MKL_INT integer;
63 #else
64 #if defined(VISP_HAVE_LAPACK_BUILT_IN)
65 typedef long int integer;
66 #else
67 typedef int integer;
68 #endif
69 extern "C" int dgesdd_(char *jobz, integer *m, integer *n, double *a, integer *lda, double *s, double *u, integer *ldu,
70  double *vt, integer *ldvt, double *work, integer *lwork, integer *iwork, integer *info);
71 
72 #include <stdio.h>
73 #include <string.h>
74 #endif
75 #endif
76 /*---------------------------------------------------------------------
77 
78 SVD related functions
79 
80 ---------------------------------------------------------------------*/
81 
82 #if defined(VISP_HAVE_OPENCV) // Require opencv >= 2.1.1
83 
150 {
151  int rows = (int)this->getRows();
152  int cols = (int)this->getCols();
153  cv::Mat m(rows, cols, CV_64F, this->data);
154  cv::SVD opencvSVD(m);
155  cv::Mat opencvV = opencvSVD.vt;
156  cv::Mat opencvW = opencvSVD.w;
157  V.resize((unsigned int)opencvV.rows, (unsigned int)opencvV.cols);
158  w.resize((unsigned int)(opencvW.rows * opencvW.cols));
159 
160  memcpy(V.data, opencvV.data, (size_t)(8 * opencvV.rows * opencvV.cols));
161  V = V.transpose();
162  memcpy(w.data, opencvW.data, (size_t)(8 * opencvW.rows * opencvW.cols));
163  this->resize((unsigned int)opencvSVD.u.rows, (unsigned int)opencvSVD.u.cols);
164  memcpy(this->data, opencvSVD.u.data, (size_t)(8 * opencvSVD.u.rows * opencvSVD.u.cols));
165 }
166 
167 #endif
168 
169 #if defined(VISP_HAVE_LAPACK)
237 {
238 #ifdef VISP_HAVE_GSL
239  {
240  // GSL cannot consider M < N. In that case we transpose input matrix
241  vpMatrix U;
242  unsigned int nc = getCols();
243  unsigned int nr = getRows();
244 
245  if (rowNum < colNum) {
246  U = this->transpose();
247  nc = getRows();
248  nr = getCols();
249  }
250  else {
251  nc = getCols();
252  nr = getRows();
253  }
254 
255  w.resize(nc);
256  V.resize(nc, nc);
257 
258  gsl_vector *work = gsl_vector_alloc(nc);
259 
260  gsl_matrix A;
261  A.size1 = nr;
262  A.size2 = nc;
263  A.tda = A.size2;
264  if (rowNum < colNum) {
265  A.data = U.data;
266  }
267  else {
268  A.data = this->data;
269  }
270  A.owner = 0;
271  A.block = 0;
272 
273  gsl_matrix V_;
274  V_.size1 = nc;
275  V_.size2 = nc;
276  V_.tda = V_.size2;
277  V_.data = V.data;
278  V_.owner = 0;
279  V_.block = 0;
280 
281  gsl_vector S;
282  S.size = nc;
283  S.stride = 1;
284  S.data = w.data;
285  S.owner = 0;
286  S.block = 0;
287 
288  gsl_linalg_SV_decomp(&A, &V_, &S, work);
289 
290  if (rowNum < colNum) {
291  (*this) = V;
292  V = U;
293  }
294 
295  gsl_vector_free(work);
296  }
297 #else
298  {
299  vpMatrix U;
300  unsigned int nc = getCols();
301  unsigned int nr = getRows();
302 
303  if (rowNum < colNum) {
304  U = this->transpose();
305  nc = getRows();
306  nr = getCols();
307  }
308  else {
309  nc = getCols();
310  nr = getRows();
311  }
312 
313  w.resize(nc);
314  V.resize(nc, nc);
315 
316  double *a = new double[static_cast<unsigned int>(nr * nc)];
317  if (rowNum < colNum) {
318  memcpy(a, U.data, this->getRows() * this->getCols() * sizeof(double));
319  }
320  else {
321  memcpy(a, this->data, this->getRows() * this->getCols() * sizeof(double));
322  }
323 
324  integer m = (integer)(nc);
325  integer n = (integer)(nr);
326  integer lda = nc;
327  integer ldu = nc;
328  integer ldvt = std::min<integer>(nr, nc);
329  integer info, lwork;
330 
331  double wkopt;
332  double *work;
333 
334  integer *iwork = new integer[8 * static_cast<integer>(std::min<integer>(nr, nc))];
335 
336  double *s = w.data;
337  double *u = V.data;
338  double *vt;
339  if (rowNum < colNum) {
340  vt = U.data;
341  }
342  else {
343  vt = this->data;
344  }
345 
346  lwork = -1;
347  dgesdd_((char *)"S", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork, iwork, &info);
348  lwork = (int)wkopt;
349  work = new double[static_cast<unsigned int>(lwork)];
350 
351  dgesdd_((char *)"S", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, iwork, &info);
352 
353  if (info > 0) {
354  throw(vpMatrixException(vpMatrixException::fatalError, "The algorithm computing SVD failed to converge."));
355  }
356 
357  if (rowNum < colNum) {
358  (*this) = V.transpose();
359  V = U;
360  }
361  else {
362  V = V.transpose();
363  }
364  delete[] work;
365  delete[] iwork;
366  delete[] a;
367  }
368 #endif
369 }
370 #endif
371 
372 #ifdef VISP_HAVE_EIGEN3
440 {
441  if (rowNum < colNum) {
442  w.resize(rowNum);
443  V.resize(colNum, rowNum);
444  }
445  else {
446  w.resize(colNum);
447  V.resize(colNum, colNum);
448  }
449 
450  Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > M(this->data, this->getRows(),
451  this->getCols());
452  Eigen::JacobiSVD<Eigen::MatrixXd> svd(M, Eigen::ComputeThinU | Eigen::ComputeThinV);
453 
454  Eigen::Map<Eigen::VectorXd> w_(w.data, w.size());
455 
456  if (rowNum < colNum) {
457  this->resize(rowNum, rowNum);
458  }
459  else {
460  this->resize(rowNum, colNum);
461  }
462  Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > V_(V.data, V.getRows(), V.getCols());
463  Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > U_(this->data, rowNum, colNum);
464  w_ = svd.singularValues();
465  V_ = svd.matrixV();
466  U_ = svd.matrixU();
467 }
468 #endif
unsigned int getCols() const
Definition: vpArray2D.h:314
double * data
Address of the first element of the data array.
Definition: vpArray2D.h:139
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition: vpArray2D.h:339
unsigned int rowNum
Number of rows in the array.
Definition: vpArray2D.h:129
unsigned int size() const
Return the number of elements of the 2D array.
Definition: vpArray2D.h:326
unsigned int getRows() const
Definition: vpArray2D.h:324
unsigned int colNum
Number of columns in the array.
Definition: vpArray2D.h:131
Implementation of column vector and the associated operations.
Definition: vpColVector.h:163
void resize(unsigned int i, bool flagNullify=true)
Definition: vpColVector.h:1056
@ fatalError
Fatal error.
Definition: vpException.h:84
error that can be emitted by the vpMatrix class and its derivatives
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:146
void svdLapack(vpColVector &w, vpMatrix &V)
void svd(vpColVector &w, vpMatrix &V)
Definition: vpMatrix.cpp:2132
void svdOpenCV(vpColVector &w, vpMatrix &V)
void svdEigen3(vpColVector &w, vpMatrix &V)
vpMatrix transpose() const
Definition: vpMatrix.cpp:472