42 #include <visp/vpConfig.h>
44 #ifndef DOXYGEN_SHOULD_SKIP_THIS
46 #include <visp/vpMatrix.h>
47 #include <visp/vpMath.h>
48 #include <visp/vpColVector.h>
49 #include <visp/vpException.h>
50 #include <visp/vpMatrixException.h>
51 #include <visp/vpDebug.h>
65 static double pythag(
double a,
double b)
70 if (absa > absb)
return absa*sqrt(1.0+
vpMath::sqr(absb/absa));
72 else return (std::fabs(absb) <= std::numeric_limits<double>::epsilon() ? 0.0 : absb*sqrt(1.0+
vpMath::sqr(absa/absb)));
78 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
100 #define MAX_ITER_SVD 50
106 double epsilon = 10*std::numeric_limits<double>::epsilon();
108 unsigned int flag,i,its,j,jj,k,l=0,nm=0;
109 double c,f,h,s,x,y,z;
110 double anorm=0.0,g=0.0,scale=0.0;
114 double **a =
new double*[m+1];
115 double **v =
new double*[n+1];
119 double *w =
new double[n+1];
120 for (i=0;i<n;i++) w[i+1] = 0.0;
131 vpERROR_TRACE(
"\n\t\tSVDcmp: You must augment A with extra zero rows") ;
133 "\n\t\tSVDcmp: You must augment A with "
134 "extra zero rows")) ;
136 double* rv1=
new double[n+1];
143 for (k=i;k<=m;k++) scale += fabs(a[k][i]);
145 if ((std::fabs(scale) > epsilon)) {
148 s += a[k][i]*a[k][i];
151 g = -SIGN(sqrt(s),f);
156 for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
158 for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
161 for (k=i;k<=m;k++) a[k][i] *= scale;
166 if (i <= m && i != n) {
167 for (k=l;k<=n;k++) scale += fabs(a[i][k]);
169 if ((std::fabs(scale) > epsilon) ) {
172 s += a[i][k]*a[i][k];
175 g = -SIGN(sqrt(s),f);
178 for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
181 for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
182 for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
185 for (k=l;k<=n;k++) a[i][k] *= scale;
193 if ((std::fabs(g) > epsilon) ) {
195 v[j][i]=(a[i][j]/a[i][l])/g;
197 for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
198 for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
201 for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
211 for (j=l;j<=n;j++) a[i][j]=0.0;
213 if ((std::fabs(g) > epsilon) ) {
217 for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
219 for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
222 for (j=i;j<=m;j++) a[j][i] *= g;
224 for (j=i;j<=m;j++) a[j][i]=0.0;
229 for (its=1;its<=MAX_ITER_SVD;its++) {
234 if ((std::fabs(rv1[l]) <= epsilon) ) {
239 if ((std::fabs(w[nm]) <= epsilon) )
break;
247 if ((std::fabs(f) > epsilon) ) {
267 for (j=1;j<=n;j++) v[j][k]=(-v[j][k]);
271 if (its == MAX_ITER_SVD)
273 for (i=0;i<n;i++) W[i] = w[i+1];
276 std::cout << *
this <<std::endl ;
285 f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
287 f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
289 for (j=l;j<=nm;j++) {
303 for (jj=1;jj<=n;jj++) {
312 if ((std::fabs(z) > epsilon) ) {
319 for (jj=1;jj<=m;jj++) {
331 for (i=0;i<n;i++) W[i] = w[i+1];
365 unsigned int m = this->
rowNum;
366 unsigned int n = this->
colNum;
376 if (std::fabs(w[j]) > std::numeric_limits<double>::epsilon())
378 for (i=0;i<m;i++) s += u[i][j]*b[i];
385 for (jj=0;jj<n;jj++) s += v[j][jj]*tmp[jj];
412 #define TOLERANCE 1.0e-7
415 void svd_internal_use(
double *U,
double *S,
double *V,
416 unsigned int nRow,
unsigned int nCol)
418 unsigned int i, j, k, EstColRank, RotCount, SweepCount, slimit;
419 double eps, e2, tol, vt, p, x0, y0, q, r, c0, s0, d1, d2;
426 e2 = 10.0 * nRow * eps * eps;
430 for (i = 0; i < nCol; i++)
431 for (j = 0; j < nCol; j++) {
432 V[nCol * i + j] = 0.0;
433 V[nCol * i + i] = 1.0;
435 RotCount = EstColRank * (EstColRank - 1) / 2;
436 while (RotCount != 0 && SweepCount <= slimit) {
437 RotCount = EstColRank * (EstColRank - 1) / 2;
439 for (j = 0; j < EstColRank - 1; j++) {
440 for (k = j + 1; k < EstColRank; k++) {
442 for (i = 0; i < nRow; i++) {
443 x0 = U[nCol * i + j];
444 y0 = U[nCol * i + k];
452 if (q <= e2 * S[0] || fabs(p) <= tol * q)
457 vt = sqrt(4 * p * p + r * r);
458 c0 = sqrt(fabs(.5 * (1 + r / vt)));
460 for (i = 0; i < nRow; i++) {
461 d1 = U[nCol * i + j];
462 d2 = U[nCol * i + k];
463 U[nCol * i + j] = d1 * c0 + d2 * s0;
464 U[nCol * i + k] = -d1 * s0 + d2 * c0;
467 for (i = 0; i < nCol; i++) {
468 d1 = V[nCol * i + j];
469 d2 = V[nCol * i + k];
470 V[nCol * i + j] = d1 * c0 + d2 * s0;
471 V[nCol * i + k] = -d1 * s0 + d2 * c0;
478 vt = sqrt(4 * p * p + q * q);
479 s0 = sqrt(fabs(.5 * (1 - q / vt)));
483 for (i = 0; i < nRow; i++) {
484 d1 = U[nCol * i + j];
485 d2 = U[nCol * i + k];
486 U[nCol * i + j] = d1 * c0 + d2 * s0;
487 U[nCol * i + k] = -d1 * s0 + d2 * c0;
490 for (i = 0; i < nCol; i++) {
491 d1 = V[nCol * i + j];
492 d2 = V[nCol * i + k];
493 V[nCol * i + j] = d1 * c0 + d2 * s0;
494 V[nCol * i + k] = -d1 * s0 + d2 * c0;
499 while (EstColRank >= 3 && S[(EstColRank - 1)] <= S[0] * tol + tol * tol)
502 for(i = 0; i < nCol; i++)
504 for(i = 0; i < nCol; i++)
505 for(j = 0; j < nRow; j++)
506 U[nCol * j + i] = U[nCol * j + i] / S[i];
542 #if (VISP_HAVE_OPENCV_VERSION >= 0x020101) // Require opencv >= 2.1.1
543 # include <opencv2/core/core.hpp>
545 int rows = (int)this->
getRows();
546 int cols = (int)this->
getCols();
547 cv::Mat m(rows, cols, CV_64F, this->
data);
548 cv::SVD opencvSVD(m);
549 cv::Mat opencvV = opencvSVD.vt;
550 cv::Mat opencvW = opencvSVD.w;
551 v.
resize((
unsigned int)opencvV.rows, (
unsigned int)opencvV.cols);
552 w.
resize((
unsigned int)(opencvW.rows*opencvW.cols));
554 memcpy(v.
data, opencvV.data, (
size_t)(8*opencvV.rows*opencvV.cols));
556 memcpy(w.
data, opencvW.data, (
size_t)(8*opencvW.rows*opencvW.cols));
557 this->
resize((
unsigned int)opencvSVD.u.rows, (
unsigned int)opencvSVD.u.cols);
558 memcpy(this->
data,opencvSVD.u.data, (
size_t)(8*opencvSVD.u.rows*opencvSVD.u.cols));
563 #ifdef VISP_HAVE_LAPACK
564 extern "C" int dgesdd_(
char *jobz,
int *m,
int *n,
double *a,
int *lda,
double *s,
double *u,
int *ldu,
double *vt,
int *ldvt,
double *work,
int *lwork,
int *iwork,
int *info);
569 int m =
static_cast<int>(this->
getCols()), n = static_cast<int>(this->
getRows()), lda = m, ldu = m, ldvt = std::min(m,n);
575 int* iwork =
new int[8*
static_cast<unsigned int>(std::min(n,m))];
578 double* a =
new double[
static_cast<unsigned int>(lda*n)];
581 double* vt = this->
data;
586 dgesdd_( (
char*)
"S", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork, iwork, &info );
588 work =
new double[
static_cast<unsigned int>(lwork)];
590 dgesdd_( (
char*)
"S", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, iwork, &info );
593 std::cout <<
"The algorithm computing SVD failed to converge." << std::endl;
605 #include <gsl/gsl_linalg.h>
617 gsl_matrix *A = gsl_matrix_alloc(nr, nc) ;
620 for (i=0 ; i < nr ; i++)
623 for (j=0 ; j < nc ; j++)
624 A->data[k+j] = (*
this)[i][j] ;
628 gsl_matrix *V = gsl_matrix_alloc(nc, nc) ;
629 gsl_vector *S = gsl_vector_alloc(nc) ;
630 gsl_vector *work = gsl_vector_alloc(nc) ;
632 gsl_linalg_SV_decomp(A,V,S, work) ;
640 for (i=0 ; i < nr ; i++)
641 for (j=0 ; j < nc ; j++)
642 (*
this)[i][j] = gsl_matrix_get(A,i,j) ;
645 for (i=0 ; i < nc ; i++)
648 for (j=0 ; j < nc ; j++)
649 v[i][j] = V->
data[k+j] ;
652 for (j=0 ; j < nc ; j++)
653 w[j] = gsl_vector_get(S,j) ;
659 gsl_vector_free(work) ;
661 #else //optimisation Anthony 20/03/2008
665 gsl_vector *work = gsl_vector_alloc(nc) ;
696 gsl_linalg_SV_decomp(&A,&V,&S, work) ;
698 gsl_vector_free(work) ;
710 #endif // doxygen should skip this
Definition of the vpMatrix class.
void resize(const unsigned int nrows, const unsigned int ncols, const bool nullify=true)
static Type maximum(const Type &a, const Type &b)
double * data
address of the first element of the data array
double ** rowPtrs
address of the first element of each rows
static double sqr(double x)
vpMatrix transpose() const
unsigned int rowNum
number of rows
Class that provides a data structure for the column vectors as well as a set of operations on these v...
unsigned int getCols() const
Return the number of columns of the matrix.
error that can be emited by the vpMatrix class and its derivates
unsigned int colNum
number of columns
unsigned int getRows() const
Return the number of rows of the matrix.
void resize(const unsigned int i, const bool flagNullify=true)