|
|
Ublas : |
Subject: Re: [ublas] [bindings] [lapack] geev and the type of the input matrix
From: Marco Guazzone (marco.guazzone_at_[hidden])
Date: 2010-07-06 02:32:43
On Tue, Jul 6, 2010 at 4:01 AM, Paul Leopardi
<paul.leopardi_at_[hidden]> wrote:
> 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;
> }
Thomas and Paul,
Thank you very much for the hints!
Cheers,
-- Marco