#define _CheckBounds

#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <ctime>
using namespace std;

#include "ToString.cpp"
#include "Timer.cpp"
#include "Vector.cpp"
#include "Matrix.cpp"
#include "PrintOp.cpp"


/* Loop has wrong order */
void Prod1(Vector<double> &y, const Matrix<double> &A, const Vector<double> &x);
void Prod1(Vector<double> &y, const Matrix<double> &A, const Vector<double> &x)
{
	int m = A.Rows();
	int n = A.Cols();

	for(int i=0; i<m; i++) y(i) = 0.0;
	for(int j=0; j<n; j++) for(int i=0; i<m; i++) y(i) += A(i,j) * x(j);
}

/* Ordering fixed */
void Prod2(Vector<double> &y, const Matrix<double> &A, const Vector<double> &x);
void Prod2(Vector<double> &y, const Matrix<double> &A, const Vector<double> &x)
{
	int m = A.Rows();
	int n = A.Cols();

	for(int i=0; i<m; i++) y(i) = 0.0;
	for(int i=0; i<m; i++) for(int j=0; j<n; j++) y(i) += A(i,j) * x(j);
}

/* No overhead */
void Prod3(Vector<double> &y, const Matrix<double> &A, const Vector<double> &x);
void Prod3(Vector<double> &y, const Matrix<double> &A, const Vector<double> &x)
{
	int m = A.Rows();
	int n = A.Cols();

	for(int i=0; i<m; i++) y(i) = 0.0;
	int Curr = 0;
	for(int i=0; i<m; i++) for(int j=0; j<n; j++) y[i] += A[Curr++] * x[j];
}

/* Poiters */
void Prod4(Vector<double> &y, const Matrix<double> &A, const Vector<double> &x);
void Prod4(Vector<double> &y, const Matrix<double> &A, const Vector<double> &x)
{
	int m = A.Rows();
	int n = A.Cols();

	double *y1 = y.Begin();
	const double *A1 = A.Begin();

	for(int i=0; i<m; i++)
	{
		*y1 = 0.0;
		double *x1 = x.Begin();
		for(int j=0; j<n; j++) *y1 += *A1++ * *x1++;
		*y1++;
	}
}

/* Loop unrolling and pointers */
void Prod5(Vector<double> &y, const Matrix<double> &A, const Vector<double> &x);
void Prod5(Vector<double> &y, const Matrix<double> &A, const Vector<double> &x)
{
	int m = A.Rows();
	int n = A.Cols();

	double *y1 = y.Begin();
	const double *A1 = A.Begin();

	for(int i=0; i<m; i++)
	{
		*y1 = 0.0;
		double *x1 = x.Begin();
		for(int j=0; j<n/5; j++)
		{
			*y1 +=   A1[0] * x1[0]
			       + A1[1] * x1[1]
			       + A1[2] * x1[2]
			       + A1[3] * x1[3]
			       + A1[4] * x1[4];
			A1+=5;
			x1+=5;
		}
		*y1++;
	}
}


int main()
{
	cout << "Begin\n";

	int m = 5000;
	int n = 5000;
	Vector<double> x(n,3.0);
	Matrix<double> A(m,n,0.001);


	Vector<double> y1(m);
	Vector<double> y2(m);
	Vector<double> y3(m);
	Vector<double> y4(m);
	Vector<double> y5(m);
	Vector<double> y6(m);

	Timer tm1;
	tm1.PrintLevel = 0;
	tm1.Start();
	Prod1(y1,A,x);
	cout << "Cpu Time Used, Prod1: " << tm1.Stop() << "\n";
	cout << "y1(0): " << y1(0) << "\n";

	Timer tm2;
	tm2.PrintLevel = 0;
	tm2.Start();
	Prod2(y2,A,x);
	cout << "Cpu Time Used, Prod2: " << tm2.Stop() << "\n";
	cout << "y2(0): " << y2(0) << "\n";

	Timer tm3;
	tm3.PrintLevel = 0;
	tm3.Start();
	Prod3(y3,A,x);
	cout << "Cpu Time Used, Prod3: " << tm3.Stop() << "\n";
	cout << "y3(0): " << y3(0) << "\n";

	Timer tm4;
	tm4.PrintLevel = 0;
	tm4.Start();
	Prod4(y4,A,x);
	cout << "Cpu Time Used, Prod4: " << tm4.Stop() << "\n";
	cout << "y4(0): " << y4(0) << "\n";

	Timer tm5;
	tm5.PrintLevel = 0;
	tm5.Start();
	Prod5(y5,A,x);
	cout << "Cpu Time Used, Prod5: " << tm5.Stop() << "\n";
	cout << "y5(0): " << y5(0) << "\n";

	cout << "End";
	return 0;
}