#include <complex>
#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);

	logical_t* dummy_selctg(0);
	fortran_int_t sdim;
	out_vector_type dummy_alphar(n);
	out_vector_type dummy_alphai(n);
	out_vector_type dummy_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
		dummy_alphar,	// alphar
		dummy_alphai,	// alphai
		dummy_beta,		// beta
		Q,				// VSR
		Z				// VSR
	);

	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;

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