[ublas] [BLAS bindings] help appreciated on using trmv and gbmv

Subject: [ublas] [BLAS bindings] help appreciated on using trmv and gbmv
From: Florent Teichteil (florent.teichteil_at_[hidden])
Date: 2012-10-07 13:03:11

Hi all,

It's me again :)

Has someone ever tried to use the BLAS bindings of trmv and gbmv?
I'm having a hard time to understand why some code returns correct
results, but some other don't. I did not find any example nor
documentation on using these bindings on the web. Either the API is
confusing, or there are some bugs that discourage to use it in
production code. I give here some examples, for which I would greatly
appreciate any help.

Test case I: product of an upper triangular matrix by a vector
I compare the results obtained by gemv and trmv with a (dense)
triangular upper matrix:

matrix<double, column_major> M = zero_matrix<double, column_major>(dim, dim);
for (std::size_t i = 0 ; i < dim ; i++) {
        for (std::size_t j = i ; j < dim ; 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);
noalias(R2) = V;
boost::numeric::bindings::blas::trmv(triangular_adaptor<matrix<double,
column_major>, upper>(M), R2);

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.

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.