Add lapack interface to JacobiSVD and BDCSVD
diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h
index 7b39285..b49a18c 100644
--- a/Eigen/src/Core/util/Memory.h
+++ b/Eigen/src/Core/util/Memory.h
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2008-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2009 Kenneth Riddile <kfriddile@yahoo.com>
// Copyright (C) 2010 Hauke Heibel <hauke.heibel@gmail.com>
diff --git a/Eigen/src/SparseCore/CompressedStorage.h b/Eigen/src/SparseCore/CompressedStorage.h
index 2845e44..e09ad97 100644
--- a/Eigen/src/SparseCore/CompressedStorage.h
+++ b/Eigen/src/SparseCore/CompressedStorage.h
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
diff --git a/lapack/complex_double.cpp b/lapack/complex_double.cpp
index 424d2b8..c9c5752 100644
--- a/lapack/complex_double.cpp
+++ b/lapack/complex_double.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -15,3 +15,4 @@
#include "cholesky.cpp"
#include "lu.cpp"
+#include "svd.cpp"
diff --git a/lapack/complex_single.cpp b/lapack/complex_single.cpp
index c0b2d29..6d11b26 100644
--- a/lapack/complex_single.cpp
+++ b/lapack/complex_single.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -15,3 +15,4 @@
#include "cholesky.cpp"
#include "lu.cpp"
+#include "svd.cpp"
diff --git a/lapack/double.cpp b/lapack/double.cpp
index d86549e..ea78bb6 100644
--- a/lapack/double.cpp
+++ b/lapack/double.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -15,3 +15,4 @@
#include "cholesky.cpp"
#include "lu.cpp"
#include "eigenvalues.cpp"
+#include "svd.cpp"
diff --git a/lapack/eigenvalues.cpp b/lapack/eigenvalues.cpp
index 6141032..921c515 100644
--- a/lapack/eigenvalues.cpp
+++ b/lapack/eigenvalues.cpp
@@ -7,10 +7,10 @@
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-#include "common.h"
+#include "lapack_common.h"
#include <Eigen/Eigenvalues>
-// computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges
+// computes eigen values and vectors of a general N-by-N matrix A
EIGEN_LAPACK_FUNC(syev,(char *jobz, char *uplo, int* n, Scalar* a, int *lda, Scalar* w, Scalar* /*work*/, int* lwork, int *info))
{
// TODO exploit the work buffer
@@ -22,24 +22,7 @@
else if(*n<0) *info = -3;
else if(*lda<std::max(1,*n)) *info = -5;
else if((!query_size) && *lwork<std::max(1,3**n-1)) *info = -8;
-
-// if(*info==0)
-// {
-// int nb = ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 )
-// LWKOPT = MAX( 1, ( NB+2 )*N )
-// WORK( 1 ) = LWKOPT
-// *
-// IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY )
-// $ INFO = -8
-// END IF
-// *
-// IF( INFO.NE.0 ) THEN
-// CALL XERBLA( 'SSYEV ', -INFO )
-// RETURN
-// ELSE IF( LQUERY ) THEN
-// RETURN
-// END IF
-
+
if(*info!=0)
{
int e = -*info;
diff --git a/lapack/lapack_common.h b/lapack/lapack_common.h
index e558c14..a935987 100644
--- a/lapack/lapack_common.h
+++ b/lapack/lapack_common.h
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2010-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2010-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -18,6 +18,11 @@
typedef Eigen::Map<Eigen::Transpositions<Eigen::Dynamic,Eigen::Dynamic,int> > PivotsType;
+#if ISCOMPLEX
+#define EIGEN_LAPACK_ARG_IF_COMPLEX(X) X,
+#else
+#define EIGEN_LAPACK_ARG_IF_COMPLEX(X)
+#endif
#endif // EIGEN_LAPACK_COMMON_H
diff --git a/lapack/single.cpp b/lapack/single.cpp
index a64ed44..c7da3ef 100644
--- a/lapack/single.cpp
+++ b/lapack/single.cpp
@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
-// Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -15,3 +15,4 @@
#include "cholesky.cpp"
#include "lu.cpp"
#include "eigenvalues.cpp"
+#include "svd.cpp"
diff --git a/lapack/svd.cpp b/lapack/svd.cpp
new file mode 100644
index 0000000..ecac3ba
--- /dev/null
+++ b/lapack/svd.cpp
@@ -0,0 +1,138 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "lapack_common.h"
+#include <Eigen/SVD>
+#include <unsupported/Eigen/BDCSVD>
+
+// computes the singular values/vectors a general M-by-N matrix A using divide-and-conquer
+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,
+ EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar */*rwork*/) int * /*iwork*/, int *info))
+{
+ // TODO exploit the work buffer
+ bool query_size = *lwork==-1;
+ int diag_size = (std::min)(*m,*n);
+
+ *info = 0;
+ if(*jobz!='A' && *jobz!='S' && *jobz!='O' && *jobz!='N') *info = -1;
+ else if(*m<0) *info = -2;
+ else if(*n<0) *info = -3;
+ else if(*lda<std::max(1,*m)) *info = -5;
+ else if(*lda<std::max(1,*m)) *info = -8;
+ else if(*ldu <1 || (*jobz=='A' && *ldu <*m)
+ || (*jobz=='O' && *m<*n && *ldu<*m)) *info = -8;
+ else if(*ldvt<1 || (*jobz=='A' && *ldvt<*n)
+ || (*jobz=='S' && *ldvt<diag_size)
+ || (*jobz=='O' && *m>=*n && *ldvt<*n)) *info = -10;
+
+ if(*info!=0)
+ {
+ int e = -*info;
+ return xerbla_(SCALAR_SUFFIX_UP"GESDD ", &e, 6);
+ }
+
+ if(query_size)
+ {
+ *lwork = 0;
+ return 0;
+ }
+
+ if(*n==0 || *m==0)
+ return 0;
+
+ PlainMatrixType mat(*m,*n);
+ mat = matrix(a,*m,*n,*lda);
+
+ int option = *jobz=='A' ? ComputeFullU|ComputeFullV
+ : *jobz=='S' ? ComputeThinU|ComputeThinV
+ : *jobz=='O' ? ComputeThinU|ComputeThinV
+ : 0;
+
+ BDCSVD<PlainMatrixType> svd(mat,option);
+
+ make_vector(s,diag_size) = svd.singularValues().head(diag_size);
+
+ if(*jobz=='A')
+ {
+ matrix(u,*m,*m,*ldu) = svd.matrixU();
+ matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
+ }
+ else if(*jobz=='S')
+ {
+ matrix(u,*m,diag_size,*ldu) = svd.matrixU();
+ matrix(vt,diag_size,*n,*ldvt) = svd.matrixV().adjoint();
+ }
+ else if(*jobz=='O' && *m>=*n)
+ {
+ matrix(a,*m,*n,*lda) = svd.matrixU();
+ matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
+ }
+ else if(*jobz=='O')
+ {
+ matrix(u,*m,*m,*ldu) = svd.matrixU();
+ matrix(a,diag_size,*n,*lda) = svd.matrixV().adjoint();
+ }
+
+ return 0;
+}
+
+// computes the singular values/vectors a general M-by-N matrix A using two sided jacobi algorithm
+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,
+ EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar */*rwork*/) int *info))
+{
+ // TODO exploit the work buffer
+ bool query_size = *lwork==-1;
+ int diag_size = (std::min)(*m,*n);
+
+ *info = 0;
+ if( *jobu!='A' && *jobu!='S' && *jobu!='O' && *jobu!='N') *info = -1;
+ else if((*jobv!='A' && *jobv!='S' && *jobv!='O' && *jobv!='N')
+ || (*jobu=='O' && *jobv=='O')) *info = -2;
+ else if(*m<0) *info = -3;
+ else if(*n<0) *info = -4;
+ else if(*lda<std::max(1,*m)) *info = -6;
+ else if(*ldu <1 || ((*jobu=='A' || *jobu=='S') && *ldu<*m)) *info = -9;
+ else if(*ldvt<1 || (*jobv=='A' && *ldvt<*n)
+ || (*jobv=='S' && *ldvt<diag_size)) *info = -11;
+
+ if(*info!=0)
+ {
+ int e = -*info;
+ return xerbla_(SCALAR_SUFFIX_UP"GESVD ", &e, 6);
+ }
+
+ if(query_size)
+ {
+ *lwork = 0;
+ return 0;
+ }
+
+ if(*n==0 || *m==0)
+ return 0;
+
+ PlainMatrixType mat(*m,*n);
+ mat = matrix(a,*m,*n,*lda);
+
+ int option = (*jobu=='A' ? ComputeFullU : *jobu=='S' || *jobu=='O' ? ComputeThinU : 0)
+ | (*jobv=='A' ? ComputeFullV : *jobv=='S' || *jobv=='O' ? ComputeThinV : 0);
+
+ JacobiSVD<PlainMatrixType> svd(mat,option);
+
+ make_vector(s,diag_size) = svd.singularValues().head(diag_size);
+
+ if(*jobu=='A') matrix(u,*m,*m,*ldu) = svd.matrixU();
+ else if(*jobu=='S') matrix(u,*m,diag_size,*ldu) = svd.matrixU();
+ else if(*jobu=='O') matrix(a,*m,diag_size,*lda) = svd.matrixU();
+
+ if(*jobv=='A') matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
+ else if(*jobv=='S') matrix(vt,diag_size,*n,*ldvt) = svd.matrixV().adjoint();
+ else if(*jobv=='O') matrix(a,diag_size,*n,*lda) = svd.matrixV().adjoint();
+
+ return 0;
+}
