Re: [ublas] [bindings] [lapack] geev and the type of the input matrix

Subject: Re: [ublas] [bindings] [lapack] geev and the type of the input matrix
From: Paul Leopardi (paul.leopardi_at_[hidden])
Date: 2010-07-05 22:01:50

Hi,
On Tuesday 06 July 2010 06:45:31 Thomas Klimpel wrote:
> Marco Guazzone wrote:
> > However, there are situations where the input matrix has only real
> > numbers and the output eigen{vectors,values} are complex.
> 
> Yes, and they will occur in complex conjugate pairs then.
> 
> > So, since (in addition to a complex matrix) a non-complex matrix may
> > also give complex eigen{vectors,values},  why is the check on the
> > value type not done on the output types (i.e., VectorW, MatrixVL,
> > MatrixVL)?
> 
> Because this is an automatically generated straightforward binding against
> the corresponding lapack routine, without any further modifications for
> greater convenience.
> 
> > IMHO, this would avoid the user (e.g., me ;) ) to copy its working
> > (possibly non-complex) matrix variable into a temp complex matrix
> > before using lapack::geev.
> 
> You shouldn't do this. Just use the non-complex version of geev with your
> non-complex matrix, and decode the complex conjugate pairs in the output
> correctly. In case you just need the eigenvalues, all you need to do is
> create one vector of complex numbers from two vectors of real numbers
> without any special case handling. The decoding of the eigenvectors is
> less nice, and will also double the memory consumption.

In GluCat ( http://glucat.sf.net ) there is a very similar situation, but with 
gees instead of geev. The file glucat/matrix_imp.h has (simplified):

  /// 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);
    lambda.clear();

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

    lambda.clear();
    for (vector_t::size_type  k=0; k!= dim; ++k)
      lambda[k] = complex_t(real_lambda[k], imag_lambda[k]);
    return lambda;
  }