|
|
Ublas : |
Subject: Re: [ublas] Eigenvalues & -vectors
From: Thomas Klimpel (Thomas.Klimpel_at_[hidden])
Date: 2010-07-04 18:26:54
Kraus Philipp wrote:
> I have looked at the example for the dsyev call under[snip]
> Can I call boost syev in the same way?
I probably misunderstand this question. You can use exactly the same argument list as in this MKL example when you replace "dsyev" by "LAPACK_DSYEV", but what would be the point of using the bindings then?
> I want to calculate for an arbitrary square matrix the
> eigenvectors and -values (like the Matlab command [SNIP]).
> Do I understand the steps in the correct order?
> 1. I create a lower bounded matrix of my sqaured matrix
> 2. I call the boost command syev with the bounded matrix
One simple way to do it would be the following code:
#include "boost/numeric/bindings/std.hpp"
#include "boost/numeric/bindings/ublas.hpp"
#include "boost/numeric/bindings/lapack/driver/syev.hpp"
//...
namespace ublas = boost::numeric::ublas;
namespace bindings = boost::numeric::bindings;
typedef ublas::matrix<double, ublas::column_major> matrix_type;
//...
matrix_type A(n,n);
// fill the lower triangle of A with the correct values
//...
std::vector<double> w(n);
bindings::lapack::syev('V', bindings::lower(A), w); // writing "Vectors" instead of 'V' won't work here, sorry for that :(
> Can I use LAPACK commands for creating the lower bounded matrix? How?
I guess you mean the "bindings::lower(A)" part in the above call. An alternative may be
typedef ublas::symmetric_adaptor<matrix_type, ublas::lower> symmetric_type ;
symmetric_type AL (A);
bindings::lapack::syev('V', AL, w);
or if you prefer to use the adaptors from the bindings and want to be really explicit
bindings::lapack::syev('V', bindings::symm(bindings::lower(A)), w);
but I'm not sure whether this would pass const correctness checks.
> What do I use the workspace in the syev function call?
You can use "bindings::minimal_workspace", "bindings::optimal_workspace" (or just omit the workspace argument) or "bindings::workspace(your_userdefined_workspace_array)". It's possible that the "bindings::workspace(your_userdefined_workspace_array)" call currently doesn't pass const correctness checks, but that should be relatively easy to fix if that's important.
Regards,
Thomas