Boost logo

Ublas :

Subject: Re: [ublas] Schur Decomposition - is Bindings/LAPACK the simplest solution?
From: Thomas Klimpel (Thomas.Klimpel_at_[hidden])
Date: 2009-06-17 16:49:37


Paul C. Leopardi wrote:
> I'd like to use the Schur decomposition in GluCat. Since GluCat already uses
> uBLAS, a "simple" way to do this would be to use Bindings with LAPACK.

So I guess the LAPACK routine you need is "gees" (sgees, dgees, cgees or zgees). Special versions are only available for matrices in Hessenberg form ("hseqr" and friends...), but the code takes the structure of the matrix into account by explicitly checking for zeros, and the result would have no special structure anyway. It's possible however that this has created a dilemma for some of the libraries that offer optimized implementations of LAPACK, because the explicit checking for zero's somehow prevents efficient exploitation of level 3 BLAS routines.

> 1. What is the status of Bindings? Is it now an official part of Boost?
> If not, where would I obtain the code?

There exists now a new "automatically" generated version at

https://svn.boost.org/svn/boost/sandbox/numeric_bindings

that will hopefully become a part of Boost someday. What's currently lacking in this version is test coverage.

The old hand written version is now available at

https://svn.boost.org/svn/boost/sandbox/numeric_bindings-v1

> 2. Is the Bindings/LAPACK code able to operate on compressed
> matrices as well as dense ones?

LAPACK works for dense (ge...) and banded (gb...) matrices. Symmetric and hermitian matrices can be dense (sy..., he...), packed (sp..., hp...) or banded+packed (sb..., hb...). In principle, these types are also supported by the Bindings.

For the solution of linear systems, compressed matrices can be handled by umfpack or mumps. (There also seem to exist bindings to SuperLU, but these are not part of the Bindings library). For eigenvalue problems, the only option to exploit detailed matrix structure seems to be to reduce it to Hessenberg form or symmetric real triangular form yourself (somehow exploiting the matrix structure...), and use the LAPACK routines for Hessenberg (hs...) or symmetric real triangular (st...) matrices.

Regards,
Thomas