#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;
}