Re: [eigen] complex SVD?

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Eigen's JacobiSVD does two-sided Jacobi rotations, which is a slow but
very accurate svd algorithm. That's why it's a lot slower than LAPACK.

Note that Intel recently contributed a MKL/BLAS/LAPACK backend, so why
don't you use that. Gael could confirm whether the Intel LAPACK
backend is fully generic LAPACK compatible.

changeset:   4410:a1f843f013c9
user:        karturov@xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
date:        Mon Dec 05 14:52:21 2011 +0700
summary:     Intel(R) MKL support added.

Documentation is here:

http://eigen.tuxfamily.org/dox-devel/TopicUsingIntelMKL.html

Benoit

2011/12/20 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?
>
>
>