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

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

	// Question 1: Is there a neater syntax to set all elements of a matrix to zero?
	matrix<std::complex<double>, column_major> A4(zero_matrix<std::complex<double>> (4, 4));
	A4(0,0)=2;
	A4(0,1)=4; 
	A4(1,0)=1; 
	A4(1,1)=3; 
	A4(2,0)=0; 
	A4(2,3)=std::complex<double>(7.1,-3.88); 
	A4(3,0)=0; 
	A4(3,1)=3.2; 

	matrix<std::complex<double>, column_major> B4(zero_matrix<std::complex<double>> (4, 4));
	B4(0,0)=11;
	B4(0,1)=3; 
	B4(1,0)=1.2; 
	B4(1,1)=-1; 
	B4(2,0)=1; 
	B4(2,3)=std::complex<double>(9,1); 
	B4(3,0)=-3.338; 
	B4(3,1)=3.2; 

	// Question 2: Is there a better way to pretty-print a matrix for debugging?
	std::cout << std::endl << std::endl << "A4 = " << A4 << std::endl;
	std::cout << std::endl << std::endl << "B4 = " << B4 << std::endl;
	
	matrix<std::complex<double>, column_major> vecs4 (4,4);
	vector<std::complex<double>> w4(4); 

	// Question 3: Is this syntax the "correct" way to call geev?
	geev(A4,w4,(matrix<std::complex<double>, column_major>*)NULL,&vecs4,boost::numeric::bindings::lapack::optimal_workspace());

	std::cout <<
		"eigenvals (w4) = " << w4 << std::endl <<
		"eigenvecs (vecs4) = " << vecs4 << std::endl;

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

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

	std::cout << "vecs4 * A4 = " << matrix<std::complex<double>, column_major>(prod(vecs4,A4)) << std::endl;

	// Question 4: Why does this next line fail?
	//std::cout << "vecs4 * A4 = " << prod<matrix<std::complex<double>, column_major>>(vecs4,A4) << std::endl;

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

	return 0;
}