// Implementation file for the Matrix class
// Created by:		Michael Edwards
// Date Started:	6/07/04
// Date Finished:
// Compiler Used:	Microsoft Visual ++, Version 6.0 Enterprise Edition

#ifndef MCL_MATRIX_CPP
#define MCL_MATRIX_CPP

#include "Matrix.h"

namespace MCL
{
	template <class NumberType>
	Matrix<NumberType>::Matrix(): column(0), data(0), row(0)		{}

	template <class NumberType>
	Matrix<NumberType>::Matrix(unsigned userRow, unsigned userColumn): column(userColumn), row(userRow)
	// Creates a matrix of row and column with all elements = 0
	// PRECONDITION:  none.
	// POSTCONDITION: a matrix is created with size (row, column) 
	{
		data = new NumberType*[row];
		for (unsigned int index = 0; index < row; ++index)
		{
			data[index] = new NumberType[column];
			memset(data[index], 0, sizeof(NumberType) * column);
		}
	}

	template <class NumberType>
	Matrix<NumberType>::Matrix(unsigned userRow, unsigned userColumn, const NumberType** const userData)
	// Converts a NumberType matrix to a matrix object
	// PRECONDITION:  for all (row, column), userData[row][column] is valid
	// POSTCONDITION: a matrix is created with size (row, column) and filled
	//				  with userData
	{
		Initialize(userRow, userColumn, userData);
	}

	template <class NumberType>
	Matrix<NumberType>::Matrix(unsigned userRow, unsigned userColumn, const NumberType* const userData): column(userColumn), row(userRow)
	// Fills a Matrix with an array of data, userData
	// PRECONDITION:  length of userData >= row * column
	// POSTCONDITION: a matrix is created with size (row, column) and filled
	//				  with userData
	{
		data = CreateMatrix(userRow, userColumn);
		unsigned index = 0;
		for (unsigned int r = 0; r < userRow; ++r, index += column)
			memcpy(data[r], &userData[index], sizeof(NumberType) * column);
	}

	template <class NumberType>
	Matrix<NumberType>::Matrix(const Matrix& rhs): row(rhs.row), column(rhs.column)
	// Copy constructor
	{
		data = CreateMatrix(row, column);
		for (unsigned int r = 0; r < row; ++r)
			memcpy(data[r], &rhs.data[r][0], sizeof(NumberType) * column);
	}

	template <class NumberType>
	Matrix<NumberType>::~Matrix()
	// Default destructor
	{
		DestroyMatrix(row, data);
	}

	template <class NumberType>
	Matrix<NumberType> Matrix<NumberType>::operator + (const Matrix &rhs) const
	// returns the matrix addition of this and rhs
	// PRECONDITION:  matrices must be of same dimension.
	// POSTCONDITION: returns this + rhs
	{
		assert( (row == rhs.row) && (column == rhs.column) );
		Matrix temp(row, column);
		for (unsigned int r = 0; r < row; ++r)
			for (unsigned int c = 0; c < column; ++c)
				temp.data[r][c] = data[r][c] + rhs.data[r][c];
		return temp;
	}

	template <class NumberType>
	Matrix<NumberType> Matrix<NumberType>::operator - (const Matrix &rhs) const
	// returns the matrix substraction of this and rhs
	// PRECONDITION:  matrices must be of same dimension.
	// POSTCONDITION: returns this - rhs
	{
		assert( (row == rhs.row) && (column == rhs.column) );
		Matrix<NumberType> temp(row, column);
		for (unsigned int r = 0; r < row; ++r)
			for (unsigned int c = 0; c < column; ++c)
				temp.data[r][c] = data[r][c] - rhs.data[r][c];
		return temp;
	}

	template <class NumberType>
	Matrix<NumberType> Matrix<NumberType>::operator * (const Matrix &rhs) const
	// returns the cross product of this and rhs
	// NOTE: this throws an exception if this->column != rhs.row
	// PRECONDITION:  this->column == rhs.row
	// POSTCONDITION: returns the cross product of this and rhs
	{
		assert(column == rhs.row);
		Matrix<NumberType> temp(row, column);
		unsigned int r, c, i;
		for (r = 0; r < row; ++r)
			for (c = 0; c < column; ++c)
				for (i = 0; i < column; ++i)
					temp.data[r][c] += data[r][i] * rhs.data[i][c];
		return temp;
	}

	template <class NumberType>
	Matrix<NumberType> Matrix<NumberType>::operator * (const NumberType scalar) const
	// does the scalar multiplication of this and a scalar, scalar
	// PRECONDITION:  none.
	// POSTCONDITION: returns the scalar multiplication of this and rhs
	{
		Matrix<NumberType> temp(row, column);
		for (unsigned int r = 0; r < row; ++r)
			for (unsigned int c = 0; c < column; ++c)
				temp.data[r][c] = data[r][c] * scalar;
		return temp;
	}

	template <class NumberType>
	Matrix<NumberType> operator * (const NumberType scalar, const Matrix<NumberType>& rhs)
	{
		Matrix<NumberType> temp(rhs.row, rhs.column);
		for (unsigned int r = 0; r < rhs.row; ++r)
			for (unsigned int c = 0; c < rhs.column; ++c)
				temp.data[r][c] = rhs.data[r][c] * scalar;
		return temp;
	}

	template <class NumberType>
	Matrix<NumberType> Matrix<NumberType>::operator / (const NumberType scalar) const
	// does the scalar division of this and a scalar, scalar
	// PRECONDITION:  scalar != 0.
	// POSTCONDITION: returns the scalar division of this and rhs
	{
		assert(scalar != 0);
		Matrix<NumberType> temp(row, column);
		for (unsigned int r = 0; r < row; ++r)
			for (unsigned int c = 0; c < column; ++c)
				temp.data[r][c] = data[r][c] / scalar;
		return temp;
	}

	template <class NumberType>
	Matrix<NumberType>& Matrix<NumberType>::operator = (const Matrix &rhs)
	// assignment operator for the matrix class
	// PRECONDITION:  none.
	// POSTCONDITION: this = rhs
	{
		if ((row != rhs.row) || (column != rhs.column))
		{
			DestroyMatrix(row, data);
			data = CreateMatrix(rhs.row, rhs.column);
			row = rhs.row;
			column = rhs.column;
		}

		for (unsigned r = 0; r < row; ++r)
			memcpy(data[r], rhs.data[r], sizeof(NumberType) * column);
				
		return *this;
	}

	template <class NumberType>
	Matrix<NumberType>& Matrix<NumberType>::operator -= (const Matrix &rhs)
	// returns the matrix substraction of this and rhs
	// NOTE: this and rhs must be the same size
	// PRECONDITION:  matrices must be of same dimension.
	// POSTCONDITION: returns this - rhs
	{
		assert( (row == rhs.row) && (column == rhs.column) );
		for (unsigned int r = 0; r < row; ++r)
			for (unsigned int c = 0; c < column; ++c)
				data[r][c] -= rhs.data[r][c];
		return *this;
	}

	template <class NumberType>
	Matrix<NumberType>& Matrix<NumberType>::operator += (const Matrix &rhs)
	// returns the matrix addition of this and rhs
	// NOTE: this and rhs must be the same size
	// PRECONDITION:  matrices must be of same dimension.
	// POSTCONDITION: returns this + rhs
	{
		assert( (row == rhs.row) && (column == rhs.column) );
		for (unsigned int r = 0; r < row; ++r)
			for (unsigned int c = 0; c < column; ++c)
				data[r][c] += rhs.data[r][c];
		return *this;
	}

	template <class NumberType>
	Matrix<NumberType>& Matrix<NumberType>::operator *= (const NumberType scalar)
	// does the scalar multiplication of this and a scalar, scalar
	// PRECONDITION:  none.
	// POSTCONDITION: returns the scalar multiplication of this and rhs
	{
		for (unsigned int r = 0; r < row; ++r)
			for (unsigned int c = 0; c < column; ++c)
				data[r][c] *= scalar;
		return *this;
	}

	template <class NumberType>
	Matrix<NumberType>& Matrix<NumberType>::operator /= (const NumberType scalar)
	// does the scalar division of this and a scalar, scalar
	// PRECONDITION:  scalar != 0.
	// POSTCONDITION: returns the scalar division of this and rhs
	{
		assert(scalar != 0);
		for (unsigned int r = 0; r < row; ++r)
			for (unsigned int c = 0; c < column; ++c)
				data[r][c] /= scalar;
		return *this;
	}

	template <class NumberType>
	const NumberType* const Matrix<NumberType>::operator[](unsigned int k) const
	// grabs member pk from the coordinate set
	// PRECONDITION:  k: [0, row)
	// POSTCONDITION: row k is returned.
	{
		assert(k < row);
		return data[k];
	}

	/*
	NumberType* Matrix::operator[](unsigned k)
	// grabs member pk from the coordinate set
	// PRECONDITION:  k: [0, row)
	// POSTCONDITION: row k is returned.
	{
		assert(k < row);
		return data[k];
	}
	*/

	template <class NumberType>
	NumberType** Matrix<NumberType>::CreateMatrix(unsigned int row, unsigned int column) const
	// creates a matrix that is row by column
	// PRECONDITION:  none.
	// POSTCONDITION: the matrix is returned
	{
		NumberType** temp;
		temp = new NumberType*[row];
		for (unsigned int index = 0; index < row; ++index)
			temp[index] = new NumberType[column];
		return temp;
	}

	template <class NumberType>
	void Matrix<NumberType>::DestroyMatrix(unsigned int row, NumberType** myMatrix)
	// destroys a matrix pointed to by myMatrix
	// PRECONDITION:
	// POSTCONDITION:
	{
		for (unsigned int r = 0; r < row; ++r)
			delete [] myMatrix[r];
		delete [] myMatrix;
	}

	/*
	NumberType Matrix::Determinant() const
	// Finds the determinant of the matrix
	// PRECONDITION:  row == column
	// POSTCONDITION:
	{
		assert(row == column);

		return 0;
	}
	*/

	template <class NumberType>
	Matrix<NumberType> Matrix<NumberType>::GenerateIdentity(unsigned int userRow, 
															unsigned int userColumn) const
	// creates an identity matrix of size [userRow, userColumn]
	// PRECONDITION:
	// POSTCONDITION:
	{
		Matrix temp(userRow, userColumn);
		for (unsigned int i = 0; i < ((userRow < userColumn) ? userRow: userColumn); ++i)
			temp.data[i][i] = 1;
		return temp;
	}

	template <class NumberType>
	void Matrix<NumberType>::Initialize(unsigned int userRow, unsigned int userColumn, 
										const NumberType** userData)
	// Fills a Matrix with an array of data, userData
	// PRECONDITION:  length of userData >= row * column
	// POSTCONDITION: a matrix is created with size (row, column) and filled
	//				  with userData
	{
		row = userRow;
		column = userColumn;
		data = CreateMatrix(userRow, userColumn);
		for (unsigned int r = 0; r < userRow; ++r)
			memcpy(data[r], &userData[r][0], sizeof(NumberType) * column);
	}

	template <class NumberType>
	Matrix<NumberType> Matrix<NumberType>::Inverse() const
	// returns an inverted form of the matrix
	// PRECONDITION:  the current matrix is square
	// POSTCONDITION: returns the inverted form of the current matrix
	{
		assert(row == column);
		Matrix<NumberType> temp(*this);
		Matrix<NumberType> inverse = GenerateIdentity(3, 3);
		Matrix<NumberType> lhs;

		for (unsigned int c = 0; c < column; ++c)
		{
			lhs = GenerateIdentity(row, column);
			if (temp.data[c][c] != 0)
				lhs.data[c][c] = 1 / temp.data[c][c];
			else 
				return Matrix(row, column);
			temp	= lhs * temp;
			inverse = lhs * inverse;

			for (unsigned int r = 0; r < row; ++r)
				if (r != c)
				{
					lhs = GenerateIdentity(row, column);
					lhs.data[r][c] = -1 * temp.data[r][c];
					temp	= lhs * temp;
					inverse = lhs * inverse;
				}
		}
		return inverse;
	}

	template <class NumberType>
	Matrix<NumberType> Matrix<NumberType>::Reduce() const
	// returns the Matrix in Row-Echelon Form 
	// PRECONDITION:  none.
	// POSTCONDITION: returns the reduced matrix of this
	{
		Matrix<NumberType> temp(*this);

		Matrix<NumberType> lhs;
		for (unsigned int c = 0; (c < column) && (c < row); ++c)
		{
			lhs = GenerateIdentity(row, column);
			if (temp.data[c][c] != 0)
				lhs.data[c][c] = 1 / temp.data[c][c];
			else 
				return temp;
			temp = lhs * temp;

			for (unsigned int r = 0; r < row; ++r)
				if (r != c)
				{
					lhs = GenerateIdentity(row, column);
					lhs.data[r][c] = -1 * temp.data[r][c];
					temp = lhs * temp;
				}
		}
		return temp;
	}

	template <class NumberType>
	Matrix<NumberType> Matrix<NumberType>::Transpose() const
	// returns the transposed matrix of this
	// PRECONDITION:  none.
	// POSTCONDITION: returns the transposition of the current matrix
	{
		Matrix<NumberType> temp(column, row);
		for (unsigned int r = 0; r < row; ++r)
			for (unsigned int c = 0; c < column; ++c)
				temp.data[c][r] = data[r][c];
		return temp;
	}
}

#endif