eigen/lapack/svd.cpp

// 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>


// 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);

  if(*jobz=='A')
  {
    BDCSVD<PlainMatrixType, ComputeFullU|ComputeFullV> svd(mat);
    make_vector(s,diag_size) = svd.singularValues().head(diag_size);
    matrix(u,*m,*m,*ldu)     = svd.matrixU();
    matrix(vt,*n,*n,*ldvt)   = svd.matrixV().adjoint();
  }
  else if(*jobz=='S')
  {
    BDCSVD<PlainMatrixType, ComputeThinU|ComputeThinV> svd(mat);
    make_vector(s,diag_size)      = svd.singularValues().head(diag_size);
    matrix(u,*m,diag_size,*ldu)   = svd.matrixU();
    matrix(vt,diag_size,*n,*ldvt) = svd.matrixV().adjoint();
  }
  else if(*jobz=='O' && *m>=*n)
  {
    BDCSVD<PlainMatrixType, ComputeThinU|ComputeThinV> svd(mat);
    make_vector(s,diag_size) = svd.singularValues().head(diag_size);
    matrix(a,*m,*n,*lda)     = svd.matrixU();
    matrix(vt,*n,*n,*ldvt)   = svd.matrixV().adjoint();
  }
  else if(*jobz=='O')
  {
    BDCSVD<PlainMatrixType, ComputeThinU|ComputeThinV> svd(mat);
    make_vector(s,diag_size)    = svd.singularValues().head(diag_size);
    matrix(u,*m,*m,*ldu)        = svd.matrixU();
    matrix(a,diag_size,*n,*lda) = svd.matrixV().adjoint();
  }
  else 
  {
    BDCSVD<PlainMatrixType> svd(mat);
    make_vector(s,diag_size)    = svd.singularValues().head(diag_size);
  }
    
  return 0;
}


template<typename MatrixType, int Options>
void gesvdAssignmentHelper(MatrixType& mat, char* jobu, char* jobv, int* m, int* n, int diag_size, Scalar* a, int* lda, RealScalar* s, Scalar* u, int* ldu, Scalar* vt, int* ldvt)
{
  JacobiSVD<MatrixType, Options> svd(mat);
  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();
  }
}

template<typename MatrixType, int Options, typename ...Args>
void gesvdSetVOptions(MatrixType& mat, char* jobu, char* jobv, Args... args)
{
  if (*jobv=='A')
  {
    gesvdAssignmentHelper<MatrixType, Options | ComputeFullV>(mat, jobu, jobv, args...);
  }
  else if (*jobv=='S' || *jobv=='O')
  {
    gesvdAssignmentHelper<MatrixType, Options | ComputeThinV>(mat, jobu, jobv, args...);
  }
  else 
  {
    gesvdAssignmentHelper<MatrixType, Options>(mat, jobu, jobv, args...);
  }
}


template<typename MatrixType, typename ...Args>
void gesvdSetUOptions(MatrixType& mat, char* jobu, char* jobv, Args... args)
{
  if (*jobu=='A')
  {
    gesvdSetVOptions<MatrixType, ComputeFullU>(mat, jobu, jobv, args...);
  }
  else if (*jobu=='S' || *jobu=='O')
  {
    gesvdSetVOptions<MatrixType, ComputeThinU>(mat, jobu, jobv, args...);
  }
  else 
  {
    gesvdSetVOptions<MatrixType, 0>(mat, jobu, jobv, args...);
  }
}

// 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);

  gesvdSetUOptions<PlainMatrixType>(mat, jobu, jobv, m, n, diag_size, a, lda, s, u, ldu, vt, ldvt); 
  
  return 0;
}