#include <boost/numeric/bindings/lapack/driver/geev.hpp>
#include <boost/numeric/bindings/ublas.hpp>
#include <boost/numeric/ublas/fwd.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <complex>
#include <iostream>

namespace ublas = boost::numeric::ublas;
namespace bindings = boost::numeric::bindings;

namespace detail {

template <typename MatrixExprT, typename OutRealVectorT, typename OutImagVectorT, typename OutLeftMatrixT, typename OutRightMatrixT>
//[NOT COMPILE]
//void eigen_impl(ublas::matrix_expression<MatrixExprT> const& A, OutRealVectorT& rw, OutImagVectorT& iw, OutLeftMatrixT& LV, OutRightMatrixT& RV, ublas::column_major_tag)
//void eigen_impl(ublas::matrix_expression<MatrixExprT>& A, OutRealVectorT& rw, OutImagVectorT& iw, OutLeftMatrixT& LV, OutRightMatrixT& RV, ublas::column_major_tag)
//void eigen_impl(MatrixExprT const& A, OutRealVectorT& rw, OutImagVectorT& iw, OutLeftMatrixT& LV, OutRightMatrixT& RV, ublas::column_major_tag)
//[/NOT COMPILE]
void eigen_impl(MatrixExprT& A, OutRealVectorT& rw, OutImagVectorT& iw, OutLeftMatrixT& LV, OutRightMatrixT& RV, ublas::column_major_tag)
{
	char jobvl;
	char jobvr;

	jobvr = jobvl = 'V';

	bindings::lapack::geev(jobvl, jobvr, A, rw, iw, LV, RV);
}

template <typename MatrixExprT, typename OutRealVectorT, typename OutImagVectorT, typename OutLeftMatrixT, typename OutRightMatrixT>
//[NOT COMPILE]
//void eigen_impl(ublas::matrix_expression<MatrixExprT> const& A, OutRealVectorT& rw, OutImagVectorT& iw, OutLeftMatrixT& LV, OutRightMatrixT& RV, ublas::row_major_tag)
//void eigen_impl(ublas::matrix_expression<MatrixExprT>& A, OutRealVectorT& rw, OutImagVectorT& iw, OutLeftMatrixT& LV, OutRightMatrixT& RV, ublas::row_major_tag)
//void eigen_impl(MatrixExprT const& A, OutRealVectorT& rw, OutImagVectorT& iw, OutLeftMatrixT& LV, OutRightMatrixT& RV, ublas::row_major_tag)
//[/NOT COMPILE]
void eigen_impl(MatrixExprT& A, OutRealVectorT& rw, OutImagVectorT& iw, OutLeftMatrixT& LV, OutRightMatrixT& RV, ublas::row_major_tag)
{
	// use a temp column-major matrix
}

} // Namespace detail


template <typename MatrixExprT, typename OutRealVectorT, typename OutImagVectorT, typename OutLeftMatrixT, typename OutRightMatrixT>
//[NOT COMPILE]
//void eigen(ublas::matrix_expression<MatrixExprT> const& A, OutRealVectorT& rv, OutImagVectorT& iv, OutLeftMatrixT& LV, OutRightMatrixT& RV)
//void eigen(ublas::matrix_expression<MatrixExprT>& A, OutRealVectorT& rv, OutImagVectorT& iv, OutLeftMatrixT& LV, OutRightMatrixT& RV)
//void eigen(MatrixExprT const& A, OutRealVectorT& rv, OutImagVectorT& iv, OutLeftMatrixT& LV, OutRightMatrixT& RV)
//[/NOT COMPILE]
void eigen(MatrixExprT& A, OutRealVectorT& rv, OutImagVectorT& iv, OutLeftMatrixT& LV, OutRightMatrixT& RV)
{
	typedef typename ublas::matrix_traits<MatrixExprT>::orientation_category orientation_category;

    detail::eigen_impl(A, rv, iv, LV, RV, orientation_category());
}

int main()
{
	typedef double value_type;
	typedef ublas::matrix<value_type, ublas::column_major> in_matrix_type;
	typedef ublas::matrix<value_type, ublas::column_major> out_matrix_type;
	typedef ublas::vector<value_type> out_vector_type;

	std::size_t n = 5;

	in_matrix_type A(n,n);

	A(0,0) =  1.96; A(0,1) = -6.49; A(0,2) = -0.47; A(0,3) = -7.20; A(0,4) = -0.65;
	A(1,0) = -6.49; A(1,1) =  3.80; A(1,2) = -6.39; A(1,3) =  1.50; A(1,4) = -6.34;
	A(2,0) = -0.47; A(2,1) = -6.39; A(2,2) =  4.17; A(2,3) = -1.51; A(2,4) =  2.67;
	A(3,0) = -7.20; A(3,1) =  1.50; A(3,2) = -1.51; A(3,3) =  5.70; A(3,4) =  1.80;
	A(4,0) = -0.65; A(4,1) = -6.34; A(4,2) =  2.67; A(4,3) =  1.80; A(4,4) = -7.10;

	out_matrix_type LV(n,n);
	out_matrix_type RV(n,n);
	out_vector_type wr(n);
	out_vector_type wi(n);

	eigen(A, wr, wi, LV, RV);

	std::cout << "Input Matrix: " << A << std::endl;
	std::cout << "Eigenvalues: ";
	for (std::size_t i = 0; i < wr.size(); ++i)
	{
		std::cout << " " << std::complex<value_type>(wr(i), wi(i));
	}
	std::cout << std::endl;
	std::cout << "Left Eigenvectors (stored columwise): " << LV << std::endl;
	std::cout << "Right Eigenvectors (stored columwise): " << RV << std::endl;
}