Boost logo

Ublas :

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

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

    namespace lapack = boost::numeric::bindings::lapack;
    typedef typename ublas::matrix<double, ublas::column_major>

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

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