#include <itpp/itbase.h>
using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

//#define NDEBUG

mat matmul ( mat &A, mat &B ) {
	mat C ( A.rows(), B.cols() );
	int i, j, k;
	double sum;

	for ( i = 0; i < A.rows(); i++ ) {
		for ( j = 0; j < A.cols(); j++ ) {
			sum = 0.0;
			for ( k = 0; k < A.cols(); k++ ) {
				sum += A._elem ( i, k ) * B._elem ( k, j );
			}
			C ( i, j ) = sum;
		}
	}
	return C;
}

void matmul2 ( int n,  double *A, double *B, double *C ) {
	int i, j, k;
	double sum;

	for ( i = 0; i < n; i++ ) {
		for ( j = 0; j < n; j++ ) {
			sum = 0.0;
			for ( k = 0; k < n; k++ ) {
				sum += A [ i*n+k ] * B [ k*n+j ];
			}
			C[ i*n+j] = sum;
		}
	}
//	return C;
}

int main() {
	Real_Timer tt;
	vec exec_times ( 4 );
	vec exec_times_b ( 4 );
	vec exec_times_c ( 4 );

	mat A;
	mat B;
	mat C;

	vec vn = "5 50 200 500";
	int n;

	for ( int i = 0; i < vn.length(); i++ ) {
		n = vn ( i );
		A = randu ( n, n );
		B = randu ( n, n );

		tt.tic();
		for ( int ii = 0; ii < 10; ii++ ) {
			C = matmul ( A, B );
		}
		exec_times ( i ) = tt.toc();

		tt.tic();
		for ( int ii = 0; ii < 10; ii++ ) {
			C = A * B;
		}
		exec_times_b ( i ) = tt.toc();

		C = zeros ( n, n );
		tt.tic();
		for ( int ii = 0; ii < 10; ii++ ) {
			matmul2 ( n, A._data(), B._data(), C._data() );
		}
		exec_times_c ( i ) = tt.toc();
	}
	cout << exec_times << endl;
	cout << exec_times_b << endl;
	cout << exec_times_c << endl;

	it_file itf ( "blas_test.it" );
	itf << Name ( "exec_times" ) << exec_times;
	itf << Name ( "exec_times_b" ) << exec_times_b;
	itf << Name ( "exec_times_c" ) << exec_times_c;

	return 0;
}