Re: [BLAS bindings] help appreciated on using trmv and gbmv

"Oswin Krause" <[email protected]>
Newsgroups gmane.comp.lib.boost.ublas
Message-ID <[email protected]>
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
> [email protected]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: [email protected]

_______________________________________________
ublas mailing list
[email protected]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [email protected]
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.