Re: [ublas] [bindings] How do I call the new gees? Why is it different from the old one?

Subject: Re: [ublas] [bindings] How do I call the new gees? Why is it different from the old one?
From: Paul Leopardi (paul.leopardi_at_[hidden])
Date: 2010-04-19 15:16:15

On Tuesday 20 April 2010 04:10:45 Rutger ter Borg wrote:
> thanks for all the feedback so far! I am wondering: isn't geev more suited
> for your job (or syev, or .. )? You don't seem to be using the Schur stuff.
> 
> I missed the point that you're calling an overload which is suited for
> complex types only. I'll see if I can include an static assert / add
> comments about that to the overloads.
> 
> The overload available for real types in MatrixA is:
> 
> gees( const char jobvs, const char sort, logical_t* select, MatrixA& a,
>         int_t& sdim, VectorWR& wr, VectorWI& wi, MatrixVS& vs );
> 
> I do remember having seen something in the older bindings (interlace?) that
> does stuff with complex/real vectors, but I am not sure if it was applied
> in this case.
> 
> Admitted, the current binding interface do not yet contain that many
> defaulting overloads. They do accept much more data types. And if it
> compiles, it will most likely work, too. :-) The idea is that some high-
> level interface can do all kinds of cool stuff with the bindings. Think of
> libraries such as uBLAS, GLAS, MTL, etc.. Adding high-level functions to
> the bindings is a lot of work; getting the bindings as robust as possible
> is a first priority for me (and getting it to a RFC-status).

Hi Rutger,
Thanks. My code now works, after I split up lambda into real and imaginary 
parts, and reassembled them.

I originally did want to use the real Schur form, which is why I called gees, 
but now I'm really only interested in the complex eigenvalues of the real 
matrix A. The calling routine uses std::arg to get the arguments of the 
eigenvalues, so it's somewhat painful to have to scan two arrays to reassemble 
the complex values for this purpose. 

I may also have a look at using geev, tomorrow, but I'm not hopeful about it, 
since I think I need to allocate left and right matrices for eigenvectors, 
which will not be used anyway. 

It's now past 5am here in Canberra, and I'm going to catch some sleep. 

The old bindings worked a treat. The new ones are a pain (*almost* as 
complicated to call as the original LAPACK routines themselves :-)

Thanks for all your help.
Best, Paul

---
#if defined(_GLUCAT_USE_BINDINGS_V1)
#include <boost/numeric/bindings/lapack/workspace.hpp>
#include <boost/numeric/bindings/lapack/gees.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#endif
#if defined(_GLUCAT_USE_BINDINGS)
#include <boost/numeric/bindings/lapack/driver/gees.hpp>
#include <boost/numeric/bindings/ublas.hpp>
#endif
[snip...]
  /// Eigenvalues of a matrix
  template< typename Matrix_T >
  ublas::vector< std::complex<double> >
  eigenvalues(const Matrix_T& val)
  {
    typedef std::complex<double> complex_t;
    typedef typename ublas::vector< complex_t > complex_vector_t;
    typedef typename Matrix_T::size_type matrix_index_t;
    const matrix_index_t dim = val.size1();
    complex_vector_t lambda = complex_vector_t(dim);
#if defined(_GLUCAT_USE_BINDINGS_V1) || defined(_GLUCAT_USE_BINDINGS)
    namespace lapack = boost::numeric::bindings::lapack;
    typedef typename ublas::matrix<double, ublas::column_major>
                     lapack_matrix_t;
    lapack_matrix_t T = to_lapack(val);
    lapack_matrix_t V = T;
#endif
#if defined(_GLUCAT_USE_BINDINGS_V1)
    lapack::gees (T, lambda, V, lapack::optimal_workspace() );
#elif defined(_GLUCAT_USE_BINDINGS)
    typedef typename ublas::vector< double > vector_t;
    vector_t real_lambda = vector_t(dim);
    vector_t imag_lambda = vector_t(dim);
    fortran_int_t sdim = 0;
    lapack::gees ('N','N',(logical_t*)0,T,sdim,real_lambda,imag_lambda,V);
    typedef typename vector_t::size_type vector_index_t;
    lambda.clear();
    for (vector_index_t k=0; k!= dim; ++k)
      lambda[k] = complex_t(real_lambda[k], imag_lambda[k]);
#else
    lambda.clear();
#endif
    return lambda;
  }