#include <complex>
#include <boost/numeric/bindings/lapack/computational/tgsen.hpp>
#include <boost/numeric/bindings/lapack/driver/gges.hpp>
#include <boost/numeric/bindings/ublas.hpp>
#include <boost/numeric/ublas/expression_types.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <iostream>
#include <vector>

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

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

	const std::size_t n = 2;

	in_matrix_type A(n,n);
	A(0,0) = 1; A(0,1) = 2;
	A(1,0) = 3; A(1,1) = 4;

	in_matrix_type B(n,n);
	B(0,0) = 5; B(0,1) = 6;
	B(1,0) = 7; B(1,1) = 8;
	
	out_matrix_type S(A);
	out_matrix_type T(B);
	out_matrix_type Q(n,n);
	out_matrix_type Z(n,n);

	::external_fp dummy_selctg(0);
	::fortran_int_t sdim;
	out_vector_type alphar(n);
	out_vector_type alphai(n);
	out_vector_type beta(n);

	// Find S,T,Q,Z s.t.: A = Q*S*Z' and B = Q*T*Z'
	bindings::lapack::gges(
		'V',			// jobvsl
		'V',			// jobvsr
		'N',			// sort
		dummy_selctg,	// selctg
		S,				// on entry A; on exit S
		T,				// on entry B; on exit T
		sdim,			// sdim
		alphar,			// alphar
		alphai,			// alphai
		beta,			// beta
		Q,				// VSR
		Z				// VSR
	);

	std::cout << "==> GGES" << std::endl;
	std::cout << "A = " << A << std::endl;
	std::cout << "B = " << B << std::endl;
	std::cout << "S = " << S << std::endl;
	std::cout << "T = " << T << std::endl;
	std::cout << "Q = " << Q << std::endl;
	std::cout << "Z = " << Z << std::endl;
	std::cout << "alphar = " << alphar << std::endl;
	std::cout << "alphai = " << alphai << std::endl;
	std::cout << "beta = " << beta << std::endl;

	ublas::vector< ::fortran_bool_t > select(n, 1); // Select ALL eigenvalues ... just a test
	out_vector_type sorted_alphar(n, 0);
	out_vector_type sorted_alphai(n, 0);
	out_vector_type sorted_beta(n, 0);
	::fortran_int_t m(0);
	out_value_type pl(0);
	out_value_type pr(0);
	out_vector_type dif(2, 0);

	bindings::lapack::tgsen(
		0,				// jobi
		1,				// wantq
		1,				// wantz
		select,			// select
		S,				// on entry S; on exit the reordered S
		T,				// on entry T; on exit the reordered T
		sorted_alphar,	// alphar
		sorted_alphai,	// alphai
		sorted_beta,	// beta
		Q,				// on entry Q; on exit the reordered Q
		Z,				// on entry Z; on exit the reordered Z
		m,				// m
		pl,				// pl
		pr,				// pr
		dif				// dif
	);

	std::cout << "==> TGSEN" << std::endl;
	std::cout << "reordered S = " << S << std::endl;
	std::cout << "reordered T = " << T << std::endl;
	std::cout << "reordered Q = " << Q << std::endl;
	std::cout << "reordered Z = " << Z << std::endl;
	std::cout << "reorderd alphar = " << sorted_alphar << std::endl;
	std::cout << "reorderd alphai = " << sorted_alphai << std::endl;
	std::cout << "reorderd beta = " << sorted_beta << std::endl;

	std::cout << "==> CHECK FOR CORRECTNESS" << std::endl;

	in_matrix_type X;
	
	X = ublas::prod(Q,S);
	X = ublas::prod(X, ublas::trans(Z));
	std::cout << "A == Q*S*Z' = " << X << std::endl;

	X = ublas::prod(Q,T);
	X = ublas::prod(X, ublas::trans(Z));
	std::cout << "B == Q*T*Z' = " << X << std::endl;
}