#ifndef _Matrix
#define _Matrix

/*

Templated two-dimensional numerical array which
looks like a conventional C multiarray.
Storage corresponds to C (row-major) ordering.
Elements are accessed via A(i,j) notation.

Array assignment is by reference (i.e. shallow assignment).
That is, B = A implies that the A and B point to the
same array, so modifications to the elements of A
will be reflected in B. If an independent copy
is required, then B = A.Copy() can be used.  Note
that this facilitates returning arrays from functions
without relying on compiler optimizations to eliminate
extensive data copying.

This class employs its own garbage collection via
the use of reference counts.  That is, whenever
an internal array storage no longer has any references
to it, it is destoryed.

Taken from Template Numerical Toolkit (TNT)

*/

template <class T>
class Matrix
{
	private:

    T*  p;
	int m;				// Number of rows in the Matrix
    int n;				// Number of columns in the Matrix
    int *RefCount;		// Check for uneeded arrays

	void Initialize();
    void Destroy();

	public:

   	string Name;		// Name of the Matrix

	       Matrix();
	       Matrix(int m, int n);
	       Matrix(int m, int n, const T &a);
	       Matrix(int m, int n, T *a);
	inline Matrix(const Matrix &A);
	inline Matrix & operator=(const Matrix &A);
	inline Matrix & Ref(const Matrix &A);
	       Matrix Copy() const;
	inline T& operator[](int i);
	inline const T& operator[](int i) const;
	inline T& operator()(int i, int j);
	inline const T& operator()(int i, int j) const;
	inline int Rows() const;
	inline int Cols() const;
	inline Vector<int> Size() const;
	inline T* Begin();
	inline const T* Begin() const;
	inline T* End();
	inline int GetRefCount() const;
           ~Matrix();

};

/* Public Funcitons */

/*

Copy constructor

Array data is NOT copied, but shared.
Thus, in Matrix B(A), subsequent changes to A will
be reflected in B.  For an indepent copy of A, use
Matrix B(A.Copy()), or B = A.Copy(), instead.

*/

template <class T>
Matrix<T>::Matrix(const Matrix<T> &A)
{
	p = A.p;
	m = A.m;
	n = A.n;
	RefCount = A.RefCount;
	Name = A.Name;
	(*RefCount)++;
}

/*

Constructor 1:

Create a new (m x n) array, WITHOUT initializing array elements.
This version avoids the O(m*n) initialization overhead and
is used just before manual assignment.

@param m:	The first (row) dimension of the new Matrix.
@param n:	The second (column) dimension of the new Matrix.

*/

template <class T>
Matrix<T>::Matrix(int m, int n)
{
	p = 0;
	this->m = m;
	this->n = n;
	RefCount = 0;
	Name = "NoName";

	if(m < 0 || n < 0)
	{
		Print("Error in Matrix::Matrix(int m, int n)");
		Print("m < 0 or n < 0");
		Print("m",m);
		Print("n",n);
		Print("HaltExecution");
		exit(0);
	}

	if(m > 0 && n > 0) Initialize();
	RefCount = new int;
	*RefCount = 1;
}

/*

Constructor 2:

Create a new (m x n) array, initializing array elements to
constant specified by argument.

@param m:	The first (row) dimension of the new Matrix.
@param n:	The second (column) dimension of the new Matrix.
@param a:	The constant value to set all elements of the new array to.

*/

template <class T>
Matrix<T>::Matrix(int m, int n, const T &a)
{
	p = 0;
	this->m = m;
	this->n = n;
	RefCount = 0;
	Name = "NoName";

	if(m < 0 || n < 0)
	{
		cout << "Error in Matrix::Matrix(int m, int n, T a)\n";
		cout << "m < 0 or n < 0\n";
		cout << "m: " << m << "\n";
		cout << "n: " << n << "\n";
		cout << "Halting Execution\n";
		exit(0);
	}

	if(m > 0 && n > 0) Initialize();
	for(int i=0; i<m*n; i++) p[i] = a;
	RefCount = new int;
	*RefCount = 1;
}

/*

Constructor 3:

Create a new (m x n) array,  as a view of an existing one-dimensional
array stored in <b>C order</b>, i.e. right-most dimension varying fastest.
(Often referred to as "row-major" ordering.)
Note that the storage for this pre-existing array will
never be garbage collected by the Array2D class.

@param m the first (row) dimension of the new Matrix.
@param n the second (column) dimension of the new Matrix.
@param a the one dimensional C array to use as data storage for the array.

*/

template <class T>
Matrix<T>::Matrix(int m, int n, T *a)
{
	p = a;
	this->m = m;
	this->n = n;
	RefCount = 0;
	Name = "NoName";
	RefCount = new int;
	*RefCount = 2; // This avoid destroying original data.

}

/* ... */
template <class T>
inline T& Matrix<T>::operator[](int i)
{
	return p[i];
}

/* ... */
template <class T>
inline const T& Matrix<T>::operator[](int i) const
{
	return p[i];
}

/* Access elements of the Matrix. Use _BoundsCheck for bounds-checking */
template <class T>
inline T& Matrix<T>::operator()(int i, int j)
{
	#ifdef _CheckBounds
		if(i < 0 || i >= m || j < 0 || j >= n)
		{
			cout << "Error in Matrix\n";
			cout << "Index Out of Bounds\n";
			cout << "Name: " << Name << "\n";
			cout << "i: " << i << "\n";
			cout << "j: " << j << "\n";
			cout << "m: " << m << "\n";
			cout << "n: " << n << "\n";
			cout << "Halting Execution\n";
			exit(0);
		}
	#endif

	return p[n*i+j];
}

/* Access elements of the Matrix. Use _BoundsCheck for bounds-checking */
template <class T>
inline const T& Matrix<T>::operator()(int i, int j) const
{
	#ifdef _CheckBounds
		if(i < 0 || i >= m || j < 0 || j >= n)
		{
			cout << "Error in Matrix" << "\n";
			cout << "Index Out of Bounds" << "\n";
			cout << "Name: " << Name << "\n";
			cout << "i: " << i << "\n";
			cout << "j: " << j << "\n";
			cout << "m: " << m << "\n";
			cout << "n: " << n << "\n";
			cout << "Halting Execution\n";
			exit(0);
		}
	#endif

	return p[n*i+j];
}

/*

Create a new copy of existing Matrix.
Used in B = A.copy() or in the construction
of B, e.g. Matrix B(A.copy()), to create a
new array that does not share data.

*/

template <class T>
Matrix<T> Matrix<T>::Copy() const
{
	Matrix A(m,n);
	for(int i=0; i<m*n; i++) A[i] = p[i];
	return A;
}

/*

Create a reference (shallow assignment) to another existing array.
In B.ref(A), B and A shared the same data and subsequent changes
to the array elements of one will be reflected in the other.

<p>
This is what operator= calls, and B = A and B.ref(A) are equivalent
operations.

@return The new referenced array: in B.ref(A), it returns B.

*/

template <class T>
Matrix<T> & Matrix<T>::Ref(const Matrix<T> &A)
{
	if(this != &A)
	{
		(*RefCount)--;
		if(*RefCount < 1) Destroy();
		m = A.m;
		n = A.n;
		p = A.p;
		Name = A.Name;
		RefCount = A.RefCount;
		(*RefCount)++ ;
	}

	return *this;
}

/* B = A is shorthand notation for B.Ref(A) */
template <class T>
Matrix<T> & Matrix<T>::operator=(const Matrix<T> &A)
{
	return Ref(A);
}

/* @return the number of rows in the Matrix */
template <class T>
inline int Matrix<T>::Rows() const
{
	return m;
}

/* @return the number of columns in the Matrix */
template <class T>
inline int Matrix<T>::Cols() const
{
	return n;
}

/* Return the Size of a Matrix in a 2-element Vector */
template <class T>
inline Vector<int> Matrix<T>::Size() const
{
	Vector<int> sz(2);
	sz[0] = m;
	sz[1] = n;
	return sz;
}

/* Returns location of first element, i.e. A[0][0] (mutable) */
template <class T>
T* Matrix<T>::Begin()
{
	return &p[0];
}

/* Returns location of first element, i.e. A[0][0] (mutable) */
template <class T>
const T* Matrix<T>::Begin() const
{
	return &p[0];
}

/* Returns location of first element, i.e. A[0][0] (mutable) */
template <class T>
T* Matrix<T>::End()
{
	return &p[0] + m * n;
}

/*

@return the number of arrays that share the same storage area
as this one.  (Must be at least one.)

*/

template <class T>
inline int Matrix<T>::GetRefCount() const
{
	return *RefCount;
}

/* Destructor */
template <class T>
Matrix<T>::~Matrix()
{
	(*RefCount)--;
	if(*RefCount < 1) Destroy();
}

/* Private internal functions */

/* Allocate memory for the Matrix */
template <class T>
void Matrix<T>::Initialize()
{
	try
	{
		p = new T[m*n];
	}
	catch(bad_alloc ba)
	{
		cout << "Error in Matrix\n";
		cout << "Memory Allocation Failed\n";
		cout << "m: " << m << "\n";
		cout << "n: " << n << "\n";
		cout << "Memory per Entry: " << sizeof(T) << "\n";
		cout << "Memory Needed in Megabytes: " << m*n*sizeof(T)/1e6 << "\n";
		cout << "Halting Execution\n";
		exit(0);
	}
}

/* Delete the Matrix from memory */
template <class T>
void Matrix<T>::Destroy()
{
	if(p != 0)
	{
		delete[] (p);
	}
	if(RefCount != 0) delete RefCount;
}

/* Create a null (0x0) array */
template <class T>
Matrix<T>::Matrix()
{
	p = 0;
	m = 0;
	n = 0;
	Name = "NoName";
	RefCount = new int;
	*RefCount = 1;
}

#endif