svd.cpp
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #include "lapack_common.h"
11 #include <Eigen/SVD>
12 
13 // computes the singular values/vectors a general M-by-N matrix A using divide-and-conquer
14 EIGEN_LAPACK_FUNC(gesdd,(char *jobz, int *m, int* n, Scalar* a, int *lda, RealScalar *s, Scalar *u, int *ldu, Scalar *vt, int *ldvt, Scalar* /*work*/, int* lwork,
15  EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar */*rwork*/) int * /*iwork*/, int *info))
16 {
17  // TODO exploit the work buffer
18  bool query_size = *lwork==-1;
19  int diag_size = (std::min)(*m,*n);
20 
21  *info = 0;
22  if(*jobz!='A' && *jobz!='S' && *jobz!='O' && *jobz!='N') *info = -1;
23  else if(*m<0) *info = -2;
24  else if(*n<0) *info = -3;
25  else if(*lda<std::max(1,*m)) *info = -5;
26  else if(*lda<std::max(1,*m)) *info = -8;
27  else if(*ldu <1 || (*jobz=='A' && *ldu <*m)
28  || (*jobz=='O' && *m<*n && *ldu<*m)) *info = -8;
29  else if(*ldvt<1 || (*jobz=='A' && *ldvt<*n)
30  || (*jobz=='S' && *ldvt<diag_size)
31  || (*jobz=='O' && *m>=*n && *ldvt<*n)) *info = -10;
32 
33  if(*info!=0)
34  {
35  int e = -*info;
36  return xerbla_(SCALAR_SUFFIX_UP"GESDD ", &e, 6);
37  }
38 
39  if(query_size)
40  {
41  *lwork = 0;
42  return 0;
43  }
44 
45  if(*n==0 || *m==0)
46  return 0;
47 
48  PlainMatrixType mat(*m,*n);
49  mat = matrix(a,*m,*n,*lda);
50 
51  int option = *jobz=='A' ? ComputeFullU|ComputeFullV
52  : *jobz=='S' ? ComputeThinU|ComputeThinV
53  : *jobz=='O' ? ComputeThinU|ComputeThinV
54  : 0;
55 
56  BDCSVD<PlainMatrixType> svd(mat,option);
57 
58  make_vector(s,diag_size) = svd.singularValues().head(diag_size);
59 
60  if(*jobz=='A')
61  {
62  matrix(u,*m,*m,*ldu) = svd.matrixU();
63  matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
64  }
65  else if(*jobz=='S')
66  {
67  matrix(u,*m,diag_size,*ldu) = svd.matrixU();
68  matrix(vt,diag_size,*n,*ldvt) = svd.matrixV().adjoint();
69  }
70  else if(*jobz=='O' && *m>=*n)
71  {
72  matrix(a,*m,*n,*lda) = svd.matrixU();
73  matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
74  }
75  else if(*jobz=='O')
76  {
77  matrix(u,*m,*m,*ldu) = svd.matrixU();
78  matrix(a,diag_size,*n,*lda) = svd.matrixV().adjoint();
79  }
80 
81  return 0;
82 }
83 
84 // computes the singular values/vectors a general M-by-N matrix A using two sided jacobi algorithm
85 EIGEN_LAPACK_FUNC(gesvd,(char *jobu, char *jobv, int *m, int* n, Scalar* a, int *lda, RealScalar *s, Scalar *u, int *ldu, Scalar *vt, int *ldvt, Scalar* /*work*/, int* lwork,
87 {
88  // TODO exploit the work buffer
89  bool query_size = *lwork==-1;
90  int diag_size = (std::min)(*m,*n);
91 
92  *info = 0;
93  if( *jobu!='A' && *jobu!='S' && *jobu!='O' && *jobu!='N') *info = -1;
94  else if((*jobv!='A' && *jobv!='S' && *jobv!='O' && *jobv!='N')
95  || (*jobu=='O' && *jobv=='O')) *info = -2;
96  else if(*m<0) *info = -3;
97  else if(*n<0) *info = -4;
98  else if(*lda<std::max(1,*m)) *info = -6;
99  else if(*ldu <1 || ((*jobu=='A' || *jobu=='S') && *ldu<*m)) *info = -9;
100  else if(*ldvt<1 || (*jobv=='A' && *ldvt<*n)
101  || (*jobv=='S' && *ldvt<diag_size)) *info = -11;
102 
103  if(*info!=0)
104  {
105  int e = -*info;
106  return xerbla_(SCALAR_SUFFIX_UP"GESVD ", &e, 6);
107  }
108 
109  if(query_size)
110  {
111  *lwork = 0;
112  return 0;
113  }
114 
115  if(*n==0 || *m==0)
116  return 0;
117 
118  PlainMatrixType mat(*m,*n);
119  mat = matrix(a,*m,*n,*lda);
120 
121  int option = (*jobu=='A' ? ComputeFullU : *jobu=='S' || *jobu=='O' ? ComputeThinU : 0)
122  | (*jobv=='A' ? ComputeFullV : *jobv=='S' || *jobv=='O' ? ComputeThinV : 0);
123 
124  JacobiSVD<PlainMatrixType> svd(mat,option);
125 
126  make_vector(s,diag_size) = svd.singularValues().head(diag_size);
127  {
128  if(*jobu=='A') matrix(u,*m,*m,*ldu) = svd.matrixU();
129  else if(*jobu=='S') matrix(u,*m,diag_size,*ldu) = svd.matrixU();
130  else if(*jobu=='O') matrix(a,*m,diag_size,*lda) = svd.matrixU();
131  }
132  {
133  if(*jobv=='A') matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
134  else if(*jobv=='S') matrix(vt,diag_size,*n,*ldvt) = svd.matrixV().adjoint();
135  else if(*jobv=='O') matrix(a,diag_size,*n,*lda) = svd.matrixV().adjoint();
136  }
137  return 0;
138 }
#define SCALAR_SUFFIX_UP
Matrix3f m
SCALAR Scalar
Definition: bench_gemm.cpp:46
cout<< "Here is the matrix m:"<< endl<< m<< endl;JacobiSVD< MatrixXf > svd(m, ComputeThinU|ComputeThinV)
#define max(a, b)
Definition: datatypes.h:20
#define EIGEN_LAPACK_ARG_IF_COMPLEX(X)
Definition: lapack_common.h:25
#define min(a, b)
Definition: datatypes.h:19
const MatrixUType & matrixU() const
Definition: SVDBase.h:101
int n
else if n * info
EIGEN_LAPACK_FUNC(gesdd,(char *jobz, int *m, int *n, Scalar *a, int *lda, RealScalar *s, Scalar *u, int *ldu, Scalar *vt, int *ldvt, Scalar *, int *lwork, EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar *) int *, int *info))
Definition: svd.cpp:14
EIGEN_WEAK_LINKING int xerbla_(const char *msg, int *info, int)
Definition: xerbla.cpp:15
Array< double, 1, 3 > e(1./3., 0.5, 2.)
RealScalar s
* lda
Definition: eigenvalues.cpp:59
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
class Bidiagonal Divide and Conquer SVD
Map< Matrix< T, Dynamic, 1 >, 0, InnerStride< Dynamic > > make_vector(T *data, int size, int incr)
Two-sided Jacobi SVD decomposition of a rectangular matrix.
const SingularValuesType & singularValues() const
Definition: SVDBase.h:129
const MatrixVType & matrixV() const
Definition: SVDBase.h:117
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
The matrix class, also used for vectors and row-vectors.


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:36:29