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 (at) 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, 86 EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar */*rwork*/) int *info)) 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 } 139
