// test_prod.cpp: Exploring syntax for matrix-matrix multiplication in uBLAS

#define BIND_FORTRAN_LOWERCASE_UNDERSCORE

#include<boost/numeric/bindings/traits/ublas_matrix.hpp> 
#include<boost/numeric/bindings/traits/ublas_vector.hpp> 
#include <boost/numeric/ublas/matrix.hpp> 
#include <boost/numeric/ublas/vector.hpp> 
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/bindings/lapack/geev.hpp> 

using namespace boost::numeric::bindings::lapack; 
using namespace boost::numeric::ublas; 

typedef matrix<std::complex<double>, column_major> cmat;
typedef vector<std::complex<double>> cvec;
typedef std::complex<double> cdouble;

void prettyprint(const std::string strName, const cmat &A)
{
	std::cout << strName << " has size [" << A.size1() << "," << A.size2() << "]:" << std::endl;

	for (int rowdx=0; rowdx < 5 && rowdx < A.size1(); rowdx++)
	{
		for (int coldx=0; coldx < 4 && coldx < A.size2(); coldx++)
		{
			if (A(rowdx,coldx).imag()==0)
				printf("%0.8g\t",real(A(rowdx,coldx)));
			else
				printf("%0.8g%+0.8gi\t",real(A(rowdx,coldx)),imag(A(rowdx,coldx)));
		}

		if (A.size1()>=5)
		{
			printf("...\n");
		}
		else
		{
			printf("\n");
		}
	}

	if (A.size2()>=5)
	{
		printf("...\n");
	}
	else
	{
		printf("\n");
	}
}

void prettyprint(const std::string strName, const cvec &A)
{
	std::cout << strName << " has size [" << A.size() << "]:" << std::endl;

	for (int coldx=0; coldx < 4 && coldx < A.size(); coldx++)
	{
		if (A(coldx).imag()==0)
			printf("%0.8g\t",real(A(coldx)));
		else
			printf("%0.8g%+0.8gi\t",real(A(coldx)),imag(A(coldx)));
	}

	if (A.size()>=5)
	{
		printf("...\n\n");
	}
	else
	{
		printf("\n\n");
	}
}

int main(int argc, char* argv[])
{
	// Now try to compute some matrix products

	cmat A4(4,4);
	A4.clear();
	A4(0,0)=2;
	A4(0,1)=4; 
	A4(1,0)=1; 
	A4(1,1)=3; 
	A4(2,0)=0; 
	A4(2,3)=cdouble(7.1,-3.88); 
	A4(3,0)=0; 
	A4(3,1)=3.2; 

	cmat B4(4, 4);
	B4.clear();
	B4(0,0)=11;
	B4(0,1)=3; 
	B4(1,0)=1.2; 
	B4(1,1)=-1; 
	B4(2,0)=1; 
	B4(2,3)=cdouble(9,1); 
	B4(3,0)=-3.338; 
	B4(3,1)=3.2; 

	// Question 1: Is this what you meant by a "manipulator"?
	prettyprint("A4", A4);
	prettyprint("B4", B4);

	cmat vecs4 (4,4);
	cvec w4(4); 

	geev(A4,w4,(cmat*) NULL,&vecs4,optimal_workspace());

	prettyprint("eigenvals (w4)",w4);
	prettyprint("eigenvecs (vecs4)",vecs4);

	// Now test multiplying matrices together:
	std::cout << "Testing matrix products: " << std::endl;

	std::cout << "vecs4 * A4 = " << prod(vecs4,A4) << std::endl << std::endl;

	std::cout << "vecs4 * A4 = " << cmat(prod(vecs4,A4)) << std::endl << std::endl;

	// Question 2: Why does this next line fail?
	//std::cout << "vecs4 * A4 = " << prod<cmat>(vecs4,A4) << std::endl;
	std::cout << "vecs4 * A4 = " << prod<cmat>(herm(vecs4),A4) << std::endl;

	// Question 3: Why can I not do prod(prod(...)) like in th next line?
	//std::cout << "vecs4 * A4 * vecs4 = " << prod(cmat(prod(vecs4,A4)),vecs4) << std::endl;

	// but this one works?
	std::cout << "vecs4 * A4 * vecs4 = " << prod(prod<cmat>(vecs4,A4),vecs4) << std::endl;
	
	return 0;
}