mirror of
https://gitlab.com/libeigen/eigen.git
synced 2026-04-10 11:34:33 +08:00
Apply clang-format to lapack/blas directories
This commit is contained in:
Antonio Sanchez
committed by
Antonio Sánchez
Antonio Sánchez
parent
4eac211e96
commit
186f8205db
@@ -11,59 +11,62 @@
|
||||
#include <Eigen/Cholesky>
|
||||
|
||||
// POTRF computes the Cholesky factorization of a real symmetric positive definite matrix A.
|
||||
EIGEN_LAPACK_FUNC(potrf,(char* uplo, int *n, RealScalar *pa, int *lda, int *info))
|
||||
{
|
||||
EIGEN_LAPACK_FUNC(potrf, (char *uplo, int *n, RealScalar *pa, int *lda, int *info)) {
|
||||
*info = 0;
|
||||
if(UPLO(*uplo)==INVALID) *info = -1;
|
||||
else if(*n<0) *info = -2;
|
||||
else if(*lda<std::max(1,*n)) *info = -4;
|
||||
if(*info!=0)
|
||||
{
|
||||
if (UPLO(*uplo) == INVALID)
|
||||
*info = -1;
|
||||
else if (*n < 0)
|
||||
*info = -2;
|
||||
else if (*lda < std::max(1, *n))
|
||||
*info = -4;
|
||||
if (*info != 0) {
|
||||
int e = -*info;
|
||||
return xerbla_(SCALAR_SUFFIX_UP"POTRF", &e, 6);
|
||||
return xerbla_(SCALAR_SUFFIX_UP "POTRF", &e, 6);
|
||||
}
|
||||
|
||||
Scalar* a = reinterpret_cast<Scalar*>(pa);
|
||||
MatrixType A(a,*n,*n,*lda);
|
||||
Scalar *a = reinterpret_cast<Scalar *>(pa);
|
||||
MatrixType A(a, *n, *n, *lda);
|
||||
int ret;
|
||||
if(UPLO(*uplo)==UP) ret = int(internal::llt_inplace<Scalar, Upper>::blocked(A));
|
||||
else ret = int(internal::llt_inplace<Scalar, Lower>::blocked(A));
|
||||
if (UPLO(*uplo) == UP)
|
||||
ret = int(internal::llt_inplace<Scalar, Upper>::blocked(A));
|
||||
else
|
||||
ret = int(internal::llt_inplace<Scalar, Lower>::blocked(A));
|
||||
|
||||
if (ret >= 0) *info = ret + 1;
|
||||
|
||||
if(ret>=0)
|
||||
*info = ret+1;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
// POTRS solves a system of linear equations A*X = B with a symmetric
|
||||
// positive definite matrix A using the Cholesky factorization
|
||||
// A = U**T*U or A = L*L**T computed by DPOTRF.
|
||||
EIGEN_LAPACK_FUNC(potrs,(char* uplo, int *n, int *nrhs, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, int *info))
|
||||
{
|
||||
EIGEN_LAPACK_FUNC(potrs,
|
||||
(char *uplo, int *n, int *nrhs, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, int *info)) {
|
||||
*info = 0;
|
||||
if(UPLO(*uplo)==INVALID) *info = -1;
|
||||
else if(*n<0) *info = -2;
|
||||
else if(*nrhs<0) *info = -3;
|
||||
else if(*lda<std::max(1,*n)) *info = -5;
|
||||
else if(*ldb<std::max(1,*n)) *info = -7;
|
||||
if(*info!=0)
|
||||
{
|
||||
if (UPLO(*uplo) == INVALID)
|
||||
*info = -1;
|
||||
else if (*n < 0)
|
||||
*info = -2;
|
||||
else if (*nrhs < 0)
|
||||
*info = -3;
|
||||
else if (*lda < std::max(1, *n))
|
||||
*info = -5;
|
||||
else if (*ldb < std::max(1, *n))
|
||||
*info = -7;
|
||||
if (*info != 0) {
|
||||
int e = -*info;
|
||||
return xerbla_(SCALAR_SUFFIX_UP"POTRS", &e, 6);
|
||||
return xerbla_(SCALAR_SUFFIX_UP "POTRS", &e, 6);
|
||||
}
|
||||
|
||||
Scalar* a = reinterpret_cast<Scalar*>(pa);
|
||||
Scalar* b = reinterpret_cast<Scalar*>(pb);
|
||||
MatrixType A(a,*n,*n,*lda);
|
||||
MatrixType B(b,*n,*nrhs,*ldb);
|
||||
Scalar *a = reinterpret_cast<Scalar *>(pa);
|
||||
Scalar *b = reinterpret_cast<Scalar *>(pb);
|
||||
MatrixType A(a, *n, *n, *lda);
|
||||
MatrixType B(b, *n, *nrhs, *ldb);
|
||||
|
||||
if(UPLO(*uplo)==UP)
|
||||
{
|
||||
if (UPLO(*uplo) == UP) {
|
||||
A.triangularView<Upper>().adjoint().solveInPlace(B);
|
||||
A.triangularView<Upper>().solveInPlace(B);
|
||||
}
|
||||
else
|
||||
{
|
||||
} else {
|
||||
A.triangularView<Lower>().solveInPlace(B);
|
||||
A.triangularView<Lower>().adjoint().solveInPlace(B);
|
||||
}
|
||||
|
||||
@@ -11,52 +11,53 @@
|
||||
#include <Eigen/Eigenvalues>
|
||||
|
||||
// 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))
|
||||
{
|
||||
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
|
||||
bool query_size = *lwork==-1;
|
||||
|
||||
bool query_size = *lwork == -1;
|
||||
|
||||
*info = 0;
|
||||
if(*jobz!='N' && *jobz!='V') *info = -1;
|
||||
else if(UPLO(*uplo)==INVALID) *info = -2;
|
||||
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)
|
||||
{
|
||||
if (*jobz != 'N' && *jobz != 'V')
|
||||
*info = -1;
|
||||
else if (UPLO(*uplo) == INVALID)
|
||||
*info = -2;
|
||||
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 e = -*info;
|
||||
return xerbla_(SCALAR_SUFFIX_UP"SYEV ", &e, 6);
|
||||
return xerbla_(SCALAR_SUFFIX_UP "SYEV ", &e, 6);
|
||||
}
|
||||
|
||||
if(query_size)
|
||||
{
|
||||
|
||||
if (query_size) {
|
||||
*lwork = 0;
|
||||
return 0;
|
||||
}
|
||||
|
||||
if(*n==0)
|
||||
return 0;
|
||||
|
||||
PlainMatrixType mat(*n,*n);
|
||||
if(UPLO(*uplo)==UP) mat = matrix(a,*n,*n,*lda).adjoint();
|
||||
else mat = matrix(a,*n,*n,*lda);
|
||||
|
||||
bool computeVectors = *jobz=='V' || *jobz=='v';
|
||||
SelfAdjointEigenSolver<PlainMatrixType> eig(mat,computeVectors?ComputeEigenvectors:EigenvaluesOnly);
|
||||
|
||||
if(eig.info()==NoConvergence)
|
||||
{
|
||||
make_vector(w,*n).setZero();
|
||||
if(computeVectors)
|
||||
matrix(a,*n,*n,*lda).setIdentity();
|
||||
|
||||
if (*n == 0) return 0;
|
||||
|
||||
PlainMatrixType mat(*n, *n);
|
||||
if (UPLO(*uplo) == UP)
|
||||
mat = matrix(a, *n, *n, *lda).adjoint();
|
||||
else
|
||||
mat = matrix(a, *n, *n, *lda);
|
||||
|
||||
bool computeVectors = *jobz == 'V' || *jobz == 'v';
|
||||
SelfAdjointEigenSolver<PlainMatrixType> eig(mat, computeVectors ? ComputeEigenvectors : EigenvaluesOnly);
|
||||
|
||||
if (eig.info() == NoConvergence) {
|
||||
make_vector(w, *n).setZero();
|
||||
if (computeVectors) matrix(a, *n, *n, *lda).setIdentity();
|
||||
//*info = 1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
make_vector(w,*n) = eig.eigenvalues();
|
||||
if(computeVectors)
|
||||
matrix(a,*n,*n,*lda) = eig.eigenvectors();
|
||||
|
||||
|
||||
make_vector(w, *n) = eig.eigenvalues();
|
||||
if (computeVectors) matrix(a, *n, *n, *lda) = eig.eigenvectors();
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
@@ -11,79 +11,74 @@
|
||||
#include <Eigen/LU>
|
||||
|
||||
// computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges
|
||||
EIGEN_LAPACK_FUNC(getrf,(int *m, int *n, RealScalar *pa, int *lda, int *ipiv, int *info))
|
||||
{
|
||||
EIGEN_LAPACK_FUNC(getrf, (int *m, int *n, RealScalar *pa, int *lda, int *ipiv, int *info)) {
|
||||
*info = 0;
|
||||
if(*m<0) *info = -1;
|
||||
else if(*n<0) *info = -2;
|
||||
else if(*lda<std::max(1,*m)) *info = -4;
|
||||
if(*info!=0)
|
||||
{
|
||||
if (*m < 0)
|
||||
*info = -1;
|
||||
else if (*n < 0)
|
||||
*info = -2;
|
||||
else if (*lda < std::max(1, *m))
|
||||
*info = -4;
|
||||
if (*info != 0) {
|
||||
int e = -*info;
|
||||
return xerbla_(SCALAR_SUFFIX_UP"GETRF", &e, 6);
|
||||
return xerbla_(SCALAR_SUFFIX_UP "GETRF", &e, 6);
|
||||
}
|
||||
|
||||
if(*m==0 || *n==0)
|
||||
return 0;
|
||||
if (*m == 0 || *n == 0) return 0;
|
||||
|
||||
Scalar* a = reinterpret_cast<Scalar*>(pa);
|
||||
Scalar *a = reinterpret_cast<Scalar *>(pa);
|
||||
int nb_transpositions;
|
||||
int ret = int(Eigen::internal::partial_lu_impl<Scalar,ColMajor,int>
|
||||
::blocked_lu(*m, *n, a, *lda, ipiv, nb_transpositions));
|
||||
int ret = int(
|
||||
Eigen::internal::partial_lu_impl<Scalar, ColMajor, int>::blocked_lu(*m, *n, a, *lda, ipiv, nb_transpositions));
|
||||
|
||||
for(int i=0; i<std::min(*m,*n); ++i)
|
||||
ipiv[i]++;
|
||||
for (int i = 0; i < std::min(*m, *n); ++i) ipiv[i]++;
|
||||
|
||||
if(ret>=0)
|
||||
*info = ret+1;
|
||||
if (ret >= 0) *info = ret + 1;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
//GETRS solves a system of linear equations
|
||||
// A * X = B or A' * X = B
|
||||
// with a general N-by-N matrix A using the LU factorization computed by GETRF
|
||||
EIGEN_LAPACK_FUNC(getrs,(char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb, int *info))
|
||||
{
|
||||
// GETRS solves a system of linear equations
|
||||
// A * X = B or A' * X = B
|
||||
// with a general N-by-N matrix A using the LU factorization computed by GETRF
|
||||
EIGEN_LAPACK_FUNC(getrs, (char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb,
|
||||
int *info)) {
|
||||
*info = 0;
|
||||
if(OP(*trans)==INVALID) *info = -1;
|
||||
else if(*n<0) *info = -2;
|
||||
else if(*nrhs<0) *info = -3;
|
||||
else if(*lda<std::max(1,*n)) *info = -5;
|
||||
else if(*ldb<std::max(1,*n)) *info = -8;
|
||||
if(*info!=0)
|
||||
{
|
||||
if (OP(*trans) == INVALID)
|
||||
*info = -1;
|
||||
else if (*n < 0)
|
||||
*info = -2;
|
||||
else if (*nrhs < 0)
|
||||
*info = -3;
|
||||
else if (*lda < std::max(1, *n))
|
||||
*info = -5;
|
||||
else if (*ldb < std::max(1, *n))
|
||||
*info = -8;
|
||||
if (*info != 0) {
|
||||
int e = -*info;
|
||||
return xerbla_(SCALAR_SUFFIX_UP"GETRS", &e, 6);
|
||||
return xerbla_(SCALAR_SUFFIX_UP "GETRS", &e, 6);
|
||||
}
|
||||
|
||||
Scalar* a = reinterpret_cast<Scalar*>(pa);
|
||||
Scalar* b = reinterpret_cast<Scalar*>(pb);
|
||||
MatrixType lu(a,*n,*n,*lda);
|
||||
MatrixType B(b,*n,*nrhs,*ldb);
|
||||
Scalar *a = reinterpret_cast<Scalar *>(pa);
|
||||
Scalar *b = reinterpret_cast<Scalar *>(pb);
|
||||
MatrixType lu(a, *n, *n, *lda);
|
||||
MatrixType B(b, *n, *nrhs, *ldb);
|
||||
|
||||
for(int i=0; i<*n; ++i)
|
||||
ipiv[i]--;
|
||||
if(OP(*trans)==NOTR)
|
||||
{
|
||||
B = PivotsType(ipiv,*n) * B;
|
||||
for (int i = 0; i < *n; ++i) ipiv[i]--;
|
||||
if (OP(*trans) == NOTR) {
|
||||
B = PivotsType(ipiv, *n) * B;
|
||||
lu.triangularView<UnitLower>().solveInPlace(B);
|
||||
lu.triangularView<Upper>().solveInPlace(B);
|
||||
}
|
||||
else if(OP(*trans)==TR)
|
||||
{
|
||||
} else if (OP(*trans) == TR) {
|
||||
lu.triangularView<Upper>().transpose().solveInPlace(B);
|
||||
lu.triangularView<UnitLower>().transpose().solveInPlace(B);
|
||||
B = PivotsType(ipiv,*n).transpose() * B;
|
||||
}
|
||||
else if(OP(*trans)==ADJ)
|
||||
{
|
||||
B = PivotsType(ipiv, *n).transpose() * B;
|
||||
} else if (OP(*trans) == ADJ) {
|
||||
lu.triangularView<Upper>().adjoint().solveInPlace(B);
|
||||
lu.triangularView<UnitLower>().adjoint().solveInPlace(B);
|
||||
B = PivotsType(ipiv,*n).transpose() * B;
|
||||
B = PivotsType(ipiv, *n).transpose() * B;
|
||||
}
|
||||
for(int i=0; i<*n; ++i)
|
||||
ipiv[i]++;
|
||||
for (int i = 0; i < *n; ++i) ipiv[i]++;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
203
lapack/svd.inc
203
lapack/svd.inc
@@ -11,128 +11,135 @@
|
||||
#include <Eigen/SVD>
|
||||
|
||||
// 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))
|
||||
{
|
||||
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);
|
||||
|
||||
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)
|
||||
{
|
||||
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);
|
||||
return xerbla_(SCALAR_SUFFIX_UP "GESDD ", &e, 6);
|
||||
}
|
||||
|
||||
if(query_size)
|
||||
{
|
||||
|
||||
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 (*n == 0 || *m == 0) return 0;
|
||||
|
||||
if(*jobz=='A')
|
||||
{
|
||||
matrix(u,*m,*m,*ldu) = svd.matrixU();
|
||||
matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
|
||||
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();
|
||||
}
|
||||
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))
|
||||
{
|
||||
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);
|
||||
|
||||
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)
|
||||
{
|
||||
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);
|
||||
return xerbla_(SCALAR_SUFFIX_UP "GESVD ", &e, 6);
|
||||
}
|
||||
|
||||
if(query_size)
|
||||
{
|
||||
|
||||
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 (*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 (*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();
|
||||
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;
|
||||
}
|
||||
Reference in New Issue
Block a user