Visual Servoing Platform version 3.6.0
Loading...
Searching...
No Matches
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>
62typedef MKL_INT integer;
63#else
64#if defined(VISP_HAVE_LAPACK_BUILT_IN)
65typedef long int integer;
66#else
67typedef int integer;
68#endif
69extern "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
78SVD related functions
79
80---------------------------------------------------------------------*/
81
82#if defined(VISP_HAVE_OPENCV) // Require opencv >= 2.1.1
83
151{
152 int rows = (int)this->getRows();
153 int cols = (int)this->getCols();
154 cv::Mat m(rows, cols, CV_64F, this->data);
155 cv::SVD opencvSVD(m);
156 cv::Mat opencvV = opencvSVD.vt;
157 cv::Mat opencvW = opencvSVD.w;
158 V.resize((unsigned int)opencvV.rows, (unsigned int)opencvV.cols);
159 w.resize((unsigned int)(opencvW.rows * opencvW.cols));
160
161 memcpy(V.data, opencvV.data, (size_t)(8 * opencvV.rows * opencvV.cols));
162 V = V.transpose();
163 memcpy(w.data, opencvW.data, (size_t)(8 * opencvW.rows * opencvW.cols));
164 this->resize((unsigned int)opencvSVD.u.rows, (unsigned int)opencvSVD.u.cols);
165 memcpy(this->data, opencvSVD.u.data, (size_t)(8 * opencvSVD.u.rows * opencvSVD.u.cols));
166}
167
168#endif
169
170#if defined(VISP_HAVE_LAPACK)
238{
239#ifdef VISP_HAVE_GSL
240 {
241 // GSL cannot consider M < N. In that case we transpose input matrix
242 vpMatrix U;
243 unsigned int nc = getCols();
244 unsigned int nr = getRows();
245
246 if (rowNum < colNum) {
247 U = this->transpose();
248 nc = getRows();
249 nr = getCols();
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 } else {
267 A.data = this->data;
268 }
269 A.owner = 0;
270 A.block = 0;
271
272 gsl_matrix V_;
273 V_.size1 = nc;
274 V_.size2 = nc;
275 V_.tda = V_.size2;
276 V_.data = V.data;
277 V_.owner = 0;
278 V_.block = 0;
279
280 gsl_vector S;
281 S.size = nc;
282 S.stride = 1;
283 S.data = w.data;
284 S.owner = 0;
285 S.block = 0;
286
287 gsl_linalg_SV_decomp(&A, &V_, &S, work);
288
289 if (rowNum < colNum) {
290 (*this) = V.transpose();
291 V = U;
292 }
293
294 gsl_vector_free(work);
295 }
296#else
297 {
298 w.resize(this->getCols());
299 V.resize(this->getCols(), this->getCols());
300
301 integer m = (integer)(this->getCols());
302 integer n = (integer)(this->getRows());
303 integer lda = m;
304 integer ldu = m;
305 integer ldvt = (std::min)(m, n);
306 integer info, lwork;
307
308 double wkopt;
309 double *work;
310
311 integer *iwork = new integer[8 * static_cast<integer>((std::min)(n, m))];
312
313 double *s = w.data;
314 double *a = new double[static_cast<unsigned int>(lda * n)];
315 memcpy(a, this->data, this->getRows() * this->getCols() * sizeof(double));
316 double *u = V.data;
317 double *vt = this->data;
318
319 lwork = -1;
320 dgesdd_((char *)"S", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork, iwork, &info);
321 lwork = (int)wkopt;
322 work = new double[static_cast<unsigned int>(lwork)];
323
324 dgesdd_((char *)"S", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, iwork, &info);
325
326 if (info > 0) {
327 throw(vpMatrixException(vpMatrixException::fatalError, "The algorithm computing SVD failed to converge."));
328 }
329
330 V = V.transpose();
331 delete[] work;
332 delete[] iwork;
333 delete[] a;
334 }
335#endif
336}
337#endif
338
339#ifdef VISP_HAVE_EIGEN3
407{
408 w.resize(this->getCols());
409 V.resize(this->getCols(), this->getCols());
410
411 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > M(this->data, this->getRows(),
412 this->getCols());
413
414 Eigen::JacobiSVD<Eigen::MatrixXd> svd(M, Eigen::ComputeThinU | Eigen::ComputeThinV);
415
416 Eigen::Map<Eigen::VectorXd> w_(w.data, w.size());
417 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > V_(V.data, V.getRows(),
418 V.getCols());
419 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > U_(this->data, this->getRows(),
420 this->getCols());
421 w_ = svd.singularValues();
422 V_ = svd.matrixV();
423 U_ = svd.matrixU();
424}
425#endif
unsigned int getCols() const
Definition vpArray2D.h:280
double * data
Address of the first element of the data array.
Definition vpArray2D.h:144
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition vpArray2D.h:305
unsigned int rowNum
Number of rows in the array.
Definition vpArray2D.h:134
unsigned int size() const
Return the number of elements of the 2D array.
Definition vpArray2D.h:292
unsigned int getRows() const
Definition vpArray2D.h:290
unsigned int colNum
Number of columns in the array.
Definition vpArray2D.h:136
Implementation of column vector and the associated operations.
void resize(unsigned int i, bool flagNullify=true)
@ 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:152
void svdLapack(vpColVector &w, vpMatrix &V)
void svd(vpColVector &w, vpMatrix &V)
void svdOpenCV(vpColVector &w, vpMatrix &V)
void svdEigen3(vpColVector &w, vpMatrix &V)
vpMatrix transpose() const
Definition vpMatrix.cpp:468