#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/ublas/operation.hpp>

using namespace boost::numeric::ublas;
const unsigned int dim = 1000;
const double epsilon = std::numeric_limits<double>::epsilon();
typedef column_major orientation;

template <typename orientation>
void testTriangular() {
	matrix<double, orientation> M = zero_matrix<double, orientation>(dim, dim);
	triangular_matrix<double, upper, orientation> T(dim, dim);
	
	for (std::size_t i = 0 ; i < dim ; i++) {
		for (std::size_t j = i ; j < dim ; j++) {
			double tmp = std::rand();
			M(i, j) = tmp;
			T(i, j) = tmp;
		}
	}
	
	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::ublas::axpy_prod(M, V, R1, true);
	boost::numeric::bindings::blas::gemv(1, M, V, 0, R2);
	std::cout << "gemv: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
	noalias(R2) = V;
	boost::numeric::bindings::blas::trmv(triangular_adaptor<matrix<double, orientation>, upper>(M), R2);
	std::cout << "trmv: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
	//noalias(R2) = V;
	//boost::numeric::bindings::blas::trmv(T, R2);
	std::cout << "trmv with packed triangular matrix: DOES NOT COMPILE" << 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;
	noalias(R2) = V;
	boost::numeric::bindings::blas::trmv(triangular_adaptor<matrix<double, orientation>, lower>(boost::numeric::bindings::trans(M)), R2);
	std::cout << "trmv with trans v1: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
	noalias(R2) = V;
	boost::numeric::bindings::blas::trmv(boost::numeric::bindings::trans(triangular_adaptor<matrix<double, orientation>, upper>(M)), R2);
	std::cout << "trmv with trans v2: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
	//noalias(R2) = V;
	//boost::numeric::bindings::blas::trmv(boost::numeric::bindings::trans(T), R2);
	std::cout << "trmv with trans and packed triangular matrix: DOES NOT COMPILE" << std::endl;
}


template <typename orientation>
void testBanded() {
	matrix<double, orientation> M = zero_matrix<double, orientation>(dim, dim);
	banded_matrix<double, orientation> B(dim, dim, 1, 1);
	
	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++) {
			double tmp = std::rand();
			M(i, j) = tmp;
			B(i, j) = tmp;
		}	
	}
	
	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::ublas::axpy_prod(M, V, R1, true);
	boost::numeric::bindings::blas::gemv(1, M, V, 0, R2);
	std::cout << "gemv: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
	boost::numeric::bindings::blas::gbmv(1, banded_adaptor<matrix<double, orientation> >(M, 1, 1), V, 0, R2);
	std::cout << "gbmv: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
	boost::numeric::bindings::blas::gbmv(1, B, V, 0, R2);
	std::cout << "gbmv with packed banded matrix: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << 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;
	boost::numeric::bindings::blas::gbmv(1, banded_adaptor<matrix<double, orientation> >(boost::numeric::bindings::trans(M), 1, 1), V, 0, R2);
	std::cout << "gbmv with trans v1: " << ((norm_inf(R1 - R2) < epsilon)?("CORRECT"):("WRONG")) << std::endl;
	//boost::numeric::bindings::blas::gbmv(1, boost::numeric::bindings::trans(banded_adaptor<matrix<double, orientation> >(M, 1, 1)), V, 0, R2);
	std::cout << "gbmv with trans v2: DOES NOT COMPILE" << std::endl;
	//boost::numeric::bindings::blas::gbmv(1, boost::numeric::bindings::trans(B), V, 0, R2);
	std::cout << "gbmv with trans and packed banded matrix: DOES NOT COMPILE" << std::endl;
}

int main () {
	std::srand(time(NULL));
	
	std::cout << "FIRST TEST CASE: triangular matrices with column major storage\n" << std::endl;
	testTriangular<column_major>();
	
	std::cout << "\n\nSECOND TEST CASE: triangular matrices with row major storage\n" << std::endl;
	testTriangular<row_major>();
	
	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;
}