|
|
Ublas : |
From: Paul Thompson (p.thompson_at_[hidden])
Date: 2006-10-13 00:39:52
On Friday 13 October 2006 12:29, Nico Galoppo wrote:
> Hi all,
>
> Often, when using the ATLAS bindings such as gesv(A,B) to solve a system A
> X = B, I only want to solve a system with one rhs, ie. A x = b. I find that
> I first have to copy my vector to a matrix and then copy the result back to
> a vector.
>
> Has anyone ever worked on a matrix proxy of a vector, such that this copy
> is not necessary anymore?
>
> Thanks!
>
> --nico
Yes, I achieved this for the LAPACK-ATLAS Cholesky solver. I added a
template specialisation for ublas::vector<T,A> at the point where it becomes
necessary to call traits to get the assumed parameters. So then I used
vector traits to get the necessary terms before calling the C ATLAS function.
It doesn't pretend to represent a vector-as-a-matrix. It does result in some
code duplication. Perhaps a neater solution exists??
Also I only wrote it for boost::numeric::ublas::vector<T,A> not all the
vectors supported by the bindings traits.
This type of mod would be suitable for the LAPACK-ATLAS gesv I suppose.
The modified clapack.hpp is attached. More details below.
Hope it helps!
Paul T.
In more detail:
template <typename SymmA, typename MatrB>
int potrs (SymmA const& a, MatrB& b)
was specialised to:
template <typename SymmA, typename T, typename A>
int potrs (SymmA const& a, boost::numeric::ublas::vector<T,A> & v)
and:
template <typename SymmA, typename MatrB>
int potrs (CBLAS_UPLO const uplo, SymmA const& a, MatrB& b)
was specialised to:
template <typename SymmA, typename T,typename A>
int potrs (CBLAS_UPLO const uplo, SymmA const& a,
boost::numeric::ublas::vector<T,A> & v)
Users can still use the same interface, eg:
cholesky_substitute (SymmA const& a, MatrB& b) { return potrs (a, b); }
still works. The template parameter MatrB ends up as a ublas::vector