|
|
Ublas : |
Subject: Re: [ublas] [BLAS bindings] help appreciated on using trmv and gbmv
From: Oswin Krause (Oswin.Krause_at_[hidden])
Date: 2012-10-07 14:57:54
Hi,
is this the real code used? because in the matrix initialisation loop
you have
> M(i, j) = = std::rand();
which looks wrong to me :).
>
> OK, the previous code works: I get R1==R2. But not if I change the
> storage to row_major. I can understand it, since BLAS requires
> column_major storage, so the bindings API should prevent from using
> row_storage by producing a compilation error, which it doesn't.
You are using the new bindings, or? In this case, this is actually
checked, by lines like that:
"BOOST_STATIC_ASSERT( (is_same<Order, tag::column_major>::value) );"
so maybe there is an error in the traits, and you should not even get
the computation of R1 running right now.
Aside from that: i don't think that the assert is needed, because a
row_major matrix can always be interpreted as a transposed column major
matrix. at least for my own implementations of the bindings (there was a
bit of overlap with the implementation of this one) this trick worked
like a charm.
I don't see any obvious errors here, maybe someone else can shd some
light on it?
>
> Now, I try to compute the product of the transpose of an upper
> triangular matrix by a vector:
>
> boost::numeric::bindings::blas::gemv(1,
> boost::numeric::bindings::trans(M), V, 0, R1);
> noalias(R2) = V;
>
> boost::numeric::bindings::blas::trmv(boost::numeric::bindings::trans(triangular_adaptor<matrix<double,
> column_major>, upper>(M)), R2);
>
> This time, I get R1 != R2
> I also tried the following option, that does not produce the correct
> answer either:
>
> boost::numeric::bindings::blas::trmv(triangular_adaptor<matrix<double,
> column_major>, lower>(boost::numeric::bindings::trans(M)), R2);
>
> Could someone explain me why the 1st code works, but not the second
> one? How can we compute the product of the transpose of an upper
> triangular matrix by a vector?
>
> Test case II: product of a banded matrix by a vector
> I compare the results obtained by gemv and gbmv with a (dense) banded
> matrix:
>
> matrix<double, column_major> dM = zero_matrix<double,
> column_major>(dim, dim);
> for (std::size_t i = 0 ; i < dim ; i++) {
> for (std::size_t j = std::max((std::size_t) 0,i-1) ; j <=
> std::min(dim-1,i+1) ; j++) {
> M(i, j) = std::rand();
> }
> }
> boost::numeric::ublas::vector<double> V(dim);
> for (std::size_t i = 0 ; i < dim ; i++) {
> V(i) = std::rand();
> }
> boost::numeric::ublas::vector<double> R1(dim);
> boost::numeric::ublas::vector<double> R2(dim);
> boost::numeric::bindings::blas::gemv(1, M, V, 0, R1);
> boost::numeric::bindings::blas::gbmv(1, banded_adaptor<matrix<double,
> column_major> >(M, 1, 1), V, 0, R2);
>
> Unfortunately, I get R1 != R2
> I checked that R1 gets the correct answer, computing
> voost::numeric::ublas::axpy_prod(M, V, R2, true), for which I get
> R1==R2 this time.
>
> So my conclusion is: either I'm misusing the API or there are bugs in
> these binding functions.
> I would greatly appreciate any help or examples showing me how to use
> these bindings with triangular or banded matrices. I'm really
> confused.
>
> Thanks!
> Florent
>
> ps: subsidiary question on other tests: why can we use gbmv with
> boost::numeric::ublas::banded_matrix (even if it produces also
> incorrect results), whereas using trmv with
> boost::numeric::ublas::triangular_matrix does not even compile? What
> is the rationale behind it? According to netlib's BLAS documentation,
> the storage of ublas' banded and triangular matrices is the same as
> netlib's ones, so both options should compile and give correct
> results.
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: Oswin.Krause_at_[hidden]