[bindings] gbmv error
Rutger ter Borg <[email protected]>
| Newsgroups | gmane.comp.lib.boost.ublas |
|---|---|
| Message-ID | <[email protected]> |
Hello all, after some further testing of the gbmv blas routine, I have encountered something that is a bug, but I am a bit puzzled by it. The tests are doing basic multiplying, see below. Matrix is 10x10 with contents shown, vector V is a unit vector. I'm puzzled, because BLAS doesn't return an error, the numbers are really close to what they should be, and the errors only occur if column_major banded matrices are used, and all arguments to the backend seem to make sense. See the ublas/bindings comparison below 'WRONG'. I've attached the test case, and I've committed the necessary upgrades to bindings::trans. Any ideas? Cheers, Rutger THIRD TEST CASE: banded matrices with column major storage row 0 [10](0,0,0,0,0,0,0,0,0,0) row 1 [10](1,1,1,0,0,0,0,0,0,0) row 2 [10](0,1,1,1,0,0,0,0,0,0) row 3 [10](0,0,1,1,1,0,0,0,0,0) row 4 [10](0,0,0,1,1,1,0,0,0,0) row 5 [10](0,0,0,0,1,1,1,0,0,0) row 6 [10](0,0,0,0,0,1,1,1,0,0) row 7 [10](0,0,0,0,0,0,1,1,1,0) row 8 [10](0,0,0,0,0,0,0,1,1,1) row 9 [10](0,0,0,0,0,0,0,0,1,1) gemv: CORRECT ublas: [10](0,3,3,3,3,3,3,3,3,2) bindings: [10](0,3,3,3,3,3,3,3,3,2) gbmv: WRONG ublas: [10](0,3,3,3,3,3,3,3,3,2) bindings: [10](1,3,3,2,3,3,3,3,3,2) T R A N S P O S E gemv with trans: CORRECT ublas: [10](1,2,3,3,3,3,3,3,3,2) bindings: [10](1,2,3,3,3,3,3,3,3,2) gbmv with trans: WRONG ublas: [10](1,2,3,3,3,3,3,3,3,2) bindings: [10](1,3,3,2,3,3,3,3,3,2) FOURTH TEST CASE: banded matrices with row major storage row 0 [10](0,0,0,0,0,0,0,0,0,0) row 1 [10](1,1,1,0,0,0,0,0,0,0) row 2 [10](0,1,1,1,0,0,0,0,0,0) row 3 [10](0,0,1,1,1,0,0,0,0,0) row 4 [10](0,0,0,1,1,1,0,0,0,0) row 5 [10](0,0,0,0,1,1,1,0,0,0) row 6 [10](0,0,0,0,0,1,1,1,0,0) row 7 [10](0,0,0,0,0,0,1,1,1,0) row 8 [10](0,0,0,0,0,0,0,1,1,1) row 9 [10](0,0,0,0,0,0,0,0,1,1) gemv: CORRECT ublas: [10](0,3,3,3,3,3,3,3,3,2) bindings: [10](0,3,3,3,3,3,3,3,3,2) gbmv: CORRECT ublas: [10](0,3,3,3,3,3,3,3,3,2) bindings: [10](0,3,3,3,3,3,3,3,3,2) T R A N S P O S E gemv with trans: CORRECT ublas: [10](1,2,3,3,3,3,3,3,3,2) bindings: [10](1,2,3,3,3,3,3,3,3,2) gbmv with trans: CORRECT ublas: [10](1,2,3,3,3,3,3,3,3,2) bindings: [10](1,2,3,3,3,3,3,3,3,2) _______________________________________________ ublas mailing list [email protected] http://lists.boost.org/mailman/listinfo.cgi/ublas Sent to: [email protected]
test_gbmv.cpp
(text/x-c++src, 3.1 KB)
#include <cstdlib>
#include <iostream>
#include <limits>
#include <boost/numeric/bindings/ublas.hpp>
#include <boost/numeric/bindings/blas.hpp>
#include <boost/numeric/bindings/trans.hpp>
#include <boost/numeric/bindings/upper.hpp>
#include <boost/numeric/bindings/unit_upper.hpp>
#include <boost/numeric/bindings/lower.hpp>
#include <boost/numeric/bindings/unit_lower.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/bindings/io.hpp>
#include <boost/numeric/bindings/noop.hpp>
using namespace boost::numeric::ublas;
const unsigned int dim = 10;
const double epsilon = 10 * std::numeric_limits<double>::epsilon();
template <typename orientation>
void testBanded() {
matrix<double, orientation> M = zero_matrix<double, orientation>(dim, dim);
banded_matrix<double, orientation> B(dim, dim, 1, 1);
B = zero_matrix<double, orientation>(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<unsigned int>(dim-1,i+1) ; j++) {
double tmp = 1.0;
M(i, j) = tmp;
B(i, j) = tmp;
}
}
for ( int i =0; i<dim; ++i ) {
std::cout << "row " << i << " " << row( B, i ) << std::endl;
}
boost::numeric::ublas::vector<double> V(dim);
for (std::size_t i = 0 ; i < dim ; i++) {
V(i) = 1.0;
}
boost::numeric::ublas::vector<double> R1(dim);
boost::numeric::ublas::vector<double> R2(dim);
boost::numeric::ublas::axpy_prod(M, V, R1, true);
boost::numeric::bindings::blas::gemv(1, M, V, 0, R2);
std::cout << "\ngemv: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
std::cout << "ublas: " << R1 << std::endl;
std::cout << "bindings: " << R2 << std::endl;
boost::numeric::ublas::axpy_prod(B, V, R1, true);
boost::numeric::bindings::blas::gbmv(1, B, V, 0, R2);
std::cout << "\ngbmv: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
std::cout << "ublas: " << R1 << std::endl;
std::cout << "bindings: " << R2 << std::endl;
std::cout << std::endl << "T R A N S P O S E " << std::endl;
boost::numeric::ublas::axpy_prod(boost::numeric::ublas::trans(M), V, R1, true);
boost::numeric::bindings::blas::gemv(1, boost::numeric::bindings::trans(M), V, 0, R2);
std::cout << "\ngemv with trans: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
std::cout << "ublas: " << R1 << std::endl;
std::cout << "bindings: " << R2 << std::endl;
boost::numeric::ublas::axpy_prod(boost::numeric::ublas::trans(B), V, R1, true);
boost::numeric::bindings::blas::gbmv(1, boost::numeric::bindings::trans(B), V, 0, R2);
std::cout << "\ngbmv with trans: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
std::cout << "ublas: " << R1 << std::endl;
std::cout << "bindings: " << R2 << std::endl;
}
int main () {
std::srand(time(NULL));
std::cout << "\n\nTHIRD TEST CASE: banded matrices with column major storage\n" << std::endl;
testBanded<column_major>();
std::cout << "\n\nFOURTH TEST CASE: banded matrices with row major storage\n" << std::endl;
testBanded<row_major>();
return 0;
}