#include <boost/numeric/ublas/matrix.hpp>
#include <boost/version.hpp>
#include <iostream>

#define BOOST_TEST_MODULE "ublas matrix test suite"
#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>

#define MY_EXPAND(x) x
#define MY_TRACE(x) std::cerr << MY_EXPAND(x) << std::endl;

static const double TOL(1.0e-5);


BOOST_AUTO_TEST_CASE( boost_version )
{
	MY_TRACE( "Using Boost: " << (BOOST_VERSION / 100000) << "." << (BOOST_VERSION / 100 % 1000) << "." << (BOOST_VERSION % 100) );
}


BOOST_AUTO_TEST_CASE( matrix_scalar_prod )
{
    MY_TRACE( "Test Case: <<matrix-scalar prod>>" );

    typedef double real_type;
    typedef boost::numeric::ublas::matrix<real_type> matrix_type;

    matrix_type A(5,4);

    A(0,0) = 0.555950; A(0,1) = 0.274690; A(0,2) = 0.540605; A(0,3) = 0.798938;
    A(1,0) = 0.108929; A(1,1) = 0.830123; A(1,2) = 0.891726; A(1,3) = 0.895283;
    A(2,0) = 0.948014; A(2,1) = 0.973234; A(2,2) = 0.216504; A(2,3) = 0.883152;
    A(3,0) = 0.023787; A(3,1) = 0.675382; A(3,2) = 0.231751; A(3,3) = 0.450332;
    A(4,0) = 1.023787; A(4,1) = 1.675382; A(4,2) = 1.231751; A(4,3) = 1.450332;

    real_type b(0.5);

    matrix_type C(5,4);

    C = b*A;

    for (matrix_type::const_iterator1 row_it = C.begin1(); row_it != C.end1(); ++row_it)
    {
        for (matrix_type::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it)
        {
            matrix_type::size_type row(col_it.index1());
            matrix_type::size_type col(col_it.index2());

            MY_TRACE( "C(" << row << "," << col << ") = " << *col_it << " ==> " << (b*A(row,col)) );

            BOOST_CHECK( std::fabs(*col_it - (b*A(row,col))) <= TOL );
        }
    }
}


BOOST_AUTO_TEST_CASE( matrix_matrix_sum )
{
	typedef double real_type;
	typedef boost::numeric::ublas::matrix<real_type> matrix_type;

	matrix_type A(5,4);

	A(0,0) = 0.555950; A(0,1) = 0.274690; A(0,2) = 0.540605; A(0,3) = 0.798938;
	A(1,0) = 0.108929; A(1,1) = 0.830123; A(1,2) = 0.891726; A(1,3) = 0.895283;
	A(2,0) = 0.948014; A(2,1) = 0.973234; A(2,2) = 0.216504; A(2,3) = 0.883152;
	A(3,0) = 0.023787; A(3,1) = 0.675382; A(3,2) = 0.231751; A(3,3) = 0.450332;
	A(4,0) = 1.023787; A(4,1) = 1.675382; A(4,2) = 1.231751; A(4,3) = 1.450332;

	matrix_type B(5,4);

	B(0,0) = 0.555950; B(0,1) = 0.274690; B(0,2) = 0.540605; B(0,3) = 0.798938;
	B(1,0) = 0.108929; B(1,1) = 0.830123; B(1,2) = 0.891726; B(1,3) = 0.895283;
	B(2,0) = 0.948014; B(2,1) = 0.973234; B(2,2) = 0.216504; B(2,3) = 0.883152;
	B(3,0) = 0.023787; B(3,1) = 0.675382; B(3,2) = 0.231751; B(3,3) = 0.450332;
	B(4,0) = 1.023787; B(4,1) = 1.675382; B(4,2) = 1.231751; B(4,3) = 1.450332;

	matrix_type C(5,4);

	C = A + B;

	for (matrix_type::const_iterator1 row_it = C.begin1(); row_it != C.end1(); ++row_it)
	{
		for (matrix_type::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it)
		{
			matrix_type::size_type row(col_it.index1());
			matrix_type::size_type col(col_it.index2());

			MY_TRACE( "C(" << row << "," << col << ") = " << A(row,col) << " + " << B(row,col) << " ==> " << (A(row,col) + B(row,col)) << " = " << *col_it );
			BOOST_CHECK( std::fabs(*col_it - (A(row,col)+B(row,col))) <= TOL );
		}
	}
}


BOOST_AUTO_TEST_CASE( matrix_matrix_prod )
{
    MY_TRACE( "Test Case: <<matrix-matrix prod>>" );

    typedef double real_type;
    typedef boost::numeric::ublas::matrix<real_type> matrix_type;

    matrix_type A(5,4);

    A(0,0) = 0.555950; A(0,1) = 0.274690; A(0,2) = 0.540605; A(0,3) = 0.798938;
    A(1,0) = 0.108929; A(1,1) = 0.830123; A(1,2) = 0.891726; A(1,3) = 0.895283;
    A(2,0) = 0.948014; A(2,1) = 0.973234; A(2,2) = 0.216504; A(2,3) = 0.883152;
    A(3,0) = 0.023787; A(3,1) = 0.675382; A(3,2) = 0.231751; A(3,3) = 0.450332;
    A(4,0) = 1.023787; A(4,1) = 1.675382; A(4,2) = 1.231751; A(4,3) = 1.450332;

    matrix_type B(4,5);

    B(0,0) = 0.555950; B(0,1) = 0.108929; B(0,2) = 0.948014; B(0,3) = 0.023787; B(0,4) = 1.023787;
    B(1,0) = 0.274690; B(1,1) = 0.830123; B(1,2) = 0.973234; B(1,3) = 0.675382; B(1,4) = 1.675382;
    B(2,0) = 0.540605; B(2,1) = 0.891726; B(2,2) = 0.216504; B(2,3) = 0.231751; B(2,4) = 1.231751;
    B(3,0) = 0.798938; B(3,1) = 0.895283; B(3,2) = 0.883152; B(3,3) = 0.450332; B(3,4) = 1.450332;

    matrix_type T(A.size1(),B.size2());

    for (matrix_type::iterator1 row_it = T.begin1(); row_it != T.end1(); ++row_it)
    {
        for (matrix_type::iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it)
        {
            for (matrix_type::size_type i = 0; i < A.size2(); ++i)
            {
                    *col_it += A(col_it.index1(), i) * B(i, col_it.index2());
            }
        }
    }

    matrix_type C(A.size1(),B.size2());

    C = boost::numeric::ublas::prod(A, B);

    for (matrix_type::const_iterator1 row_it = C.begin1(); row_it != C.end1(); ++row_it)
    {
        for (matrix_type::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it)
        {
            matrix_type::size_type row(col_it.index1());
            matrix_type::size_type col(col_it.index2());

            MY_TRACE( "C(" << row << "," << col << ") = " << *col_it << " ==> " << T(row,col) );

            BOOST_CHECK( std::fabs(*col_it - T(row,col)) <= TOL );
        }
    }
}