#include <boost/numeric/bindings/ublas.hpp>
#include <boost/numeric/bindings/blas.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>



using namespace boost::numeric::ublas;




template< typename Orientation > 
void run_test() {
    
int dim = 7;    

matrix<double, Orientation> M = zero_matrix<double, Orientation>(dim, dim);
for (std::size_t i = 0 ; i < dim ; i++) {
        for (std::size_t j = i ; j < dim ; j++) {
                M(i, j) = double(std::rand()) / 1e9;
        }
}

boost::numeric::ublas::vector<double> V(dim);
for (std::size_t i = 0 ; i < dim ; i++) {
        V(i) = double(std::rand()) / 1e9;
}
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,
Orientation>, upper>(M), R2);
    
std::cout << prod( M, V ) << std::endl;
std::cout << R1 << std::endl;
std::cout << R2 << std::endl;
        
}







int main() {
    run_test< column_major >();
    run_test< row_major>();
    
}