Re: [eigen] complex SVD?

[Mark Borgerding <mark@xxxxxxxxxxxxxx>]

Last night when I thought Eigen couldn't handle complex SVD (based on 
apparently outdated comments in SVD and JacobiSvd.h),
I held my nose and cra**ed out some C++ code to use LAPACK's CGESDD 
function for complex<float> SVD.
The lapack routines are all in Fortran with no C/C++ declarations and 
thus have some ugly calling conventions

extern "C" {
typedef void* fcomplex_ptr;
    void cgesdd_( const char* jobz, const int* m,const int* n, 
fcomplex_ptr a,
            const int* lda, float* s, fcomplex_ptr u,const int* ldu, 
fcomplex_ptr vt,const int* ldvt,
            fcomplex_ptr work, const int* lwork, float* rwork, int* 
iwork, int* info );
#define cgesdd cgesdd_
}
.... lots of initialization ...
 cgesdd( "S", &m, &n, &_a[0], &lda, &_s[0], &_u[0], &ldu, &_vt[0], 
&ldvt, &_work[0], &lwork, &_rwork[0], &_iwork[0], &info );
....


The Good News:

This morning I saw Gael's note and was overjoyed that I could ditch the 
ugly code in favor of much more elegant Eigen code:

        _svd = JacobiSVD<typename 
MatType::PlainObject>(_rows,_cols,ComputeThinU|ComputeThinV);
        ...
        memcpy(_A.data(),inblock.data(),inblock.bytes() );
        _svd.compute(_A);
        _singVals = _svd.singularValues();


The Bad News:

The elegant code consistently runs much slower than LAPACK. :(
E.g. ~5x longer to do a 128x1000 complex<float> SVD  (see system details 
below)

Both are computing full set of singular values and the thinned subset of 
the left&right singular vectors  -- (This is the "S" argument for 
cgesdd; ComputeThinU|ComputeThinV for JacobiSvd) .

Both are initialized outside the loop and then called several times. 
(calling cgesdd with lwork=-1 to find size, reusing JacobiSvd object )

I double checked that Eigen is using SIMD ( i.e. lots of code in the 
assembly like 'mulps  %xmm7,%xmm6')

EIGEN_NO_DEBUG is defined

Is there anything else I could be doing wrong that would cause such a 
big speed gap?

-- Mark

P.S.
System details:
ubuntu 11.10 x86_64
Core i7 CPU       Q 820  @ 1.73GHz
g++ 4.6.1
eigen 3.0.3
liblapack-dev    3.3.1-1



On 12/20/2011 02:45 AM, Gael Guennebaud wrote:
> JacobiSVD does work for complex matrices.
>
> gael
>
> On Tue, Dec 20, 2011 at 3:08 AM, Mark Borgerding<mark@xxxxxxxxxxxxxx>  wrote:
> > Is there an SVD decomposition for complex matrices?