Re: [ublas] lapack-bindings: geev

From: Andreas Klöckner (kloeckner_at_[hidden])
Date: 2005-04-19 02:54:55

Hi Franz,

Franz Müller schrieb:
> I'm using VC++.NET 7.1 and the newest lapack-bindings from CVS.
> 
> I tried to get the eigenvalues of a matrix, but unfortunately didn't
> succeed. I'm calling geev using 
> lapack::geev(a, Eigenvals, &leftEV, &Eigenvecs, 
> 		lapack::optimal_workspace());
> 
> It compiles fine, but then on execution there's a runtime-error saying
> "** On entry to DGEEV , parameter number 13 had an illegal value"
> and I'll receive only zeros as eigenvalues.

First of all, which case do you end up with? You may have seen in
geev.hpp that the dispatch routine distinguishes real_case, mixed_case
and complex_case. It would be helpful to know which of those it is that
is faulty.

Parameter 13, if I'm counting right, is LWORK, the length of the
workspace. Could you please find out whether the error is caused by the
workspace query (the first call to geev_backend, with lwork==-1) or the
second one? It would also be helpful to know the value of lwork after
the workspace query, just to see if something is horribly off.

> Am I doing something wrong? Unfortunately I couldn't figure out what's wrong
> on my own.

Without seeing your code, well, ... ;)

> If anybody has a working example, could you please send it to me? (By the
> way, I'd appreciate any example or additional doc using the lapack-bindings,
> not only for general eigenvalues)

This code works for me, but it might not help you that much. It's from a
BPL binding to UBlas that I'll release in a while.

8< ---------------------------------------------------------------------
template <typename ValueType>
static python::object eigenvectorsWrapper(unsigned get_left_vectors,
                                          unsigned get_right_vectors,
                                          const ublas::matrix<ValueType> &a)
{
  typedef ublas::matrix<ValueType> mat;
  typedef ublas::matrix<ValueType, ublas::column_major> fortran_mat;
  typedef ublas::vector<typename helpers::complexify<ValueType>::type >
eval_vec;
  typedef ublas::matrix<typename helpers::complexify<ValueType>::type,
ublas::column_major> evec_mat;

  int const n = a.size1();
  fortran_mat a_copy(a);
  std::auto_ptr<eval_vec> w(new eval_vec(a.size1()));

  std::auto_ptr<evec_mat> vl, vr;
  if (get_left_vectors)
    vl = std::auto_ptr<evec_mat>(new evec_mat(n, n));
  if (get_right_vectors)
    vr = std::auto_ptr<evec_mat>(new evec_mat(n, n));

  int ierr = boost::numeric::bindings::lapack::geev(a_copy, *w,
vl.get(), vr.get(),
                                                    boost::numeric::bindings::lapack::optimal_workspace());

  if (ierr < 0)
    throw std::runtime_error("invalid argument to geev");
  else if (ierr > 0)
    throw std::runtime_error("no convergence for given matrix");

  typedef ublas::matrix<typename helpers::complexify<ValueType>::type> cmat;
  typedef ublas::vector<ValueType> vec;

  if (get_left_vectors)
  {
    if (get_right_vectors)
      return python::make_tuple(python::object(w.release()),
                                python::object(new cmat(*vl)),
                                python::object(new cmat(*vr)));
    else
      return python::make_tuple(python::object(w.release()),
                                python::object(new cmat(*vl)),
                                python::object());
  }
  else
  {
    if (get_right_vectors)
      return python::make_tuple(python::object(w.release()),
                                python::object(),
                                python::object(new cmat(*vr)));
    else
      return python::make_tuple(python::object(w.release()),
                                python::object(),
                                python::object());
  }
}
8< ---------------------------------------------------------------------

> One last question: is anybody actively working on the bindings? Any plan to
> add them to boost (non-sandbox)?

I'll leave that for someone else to answer, noting that at least I'm an
active *user* who has contributed a little (notably, geev.hpp is my fault).

Andreas

-- 
--------------------------------------+-----------------------------
Dipl.-Math. techn. Andreas Klöckner   | FON [phone]
Institut für Angewandte Mathematik II | FAX [phone]
Universität Karlsruhe                 | EMAIL kloeckner_at_[hidden]
Englerstraße 2                        |
76128 Karlsruhe                       | RAUM 022
Germany                               | Kollegiengebäude Mathematik
--------------------------------------+-----------------------------