Matrix<T> Matrix<T>::Cat(int dimension, Vector<T>& m1, Matrix<T>& m2)
{
	return Matrix<T>::Cat(dimension, (Matrix<T>)m1, m2);
}

template< class T >
Matrix<T> Matrix<T>::Cat(int dimension, Vector<T>& m1, Vector<T>& m2)
{
	return Matrix<T>::Cat(dimension, (Matrix<T>)m1, (Matrix<T>)m2);
}


template< class T >
Matrix<T> Matrix<T>::Ones(int side)
{
	Matrix<T> temp(side,side,1);
	return temp;
}


template< class T >
Matrix<T> Matrix<T>::Ones(int rows, int cols)
{
	Matrix<T> temp(rows,cols,1);
	return temp;
}

template< class T >
Matrix<T> Matrix<T>::Zeros(int side)
{
	Matrix<T> temp(side,side,0);
	return temp;
}


template< class T >
Matrix<T> Matrix<T>::Zeros(int rows, int cols)
{
	Matrix<T> temp(rows,cols,0);
	return temp;
}







template< class T >
bool Matrix<T>::IsSquare(Matrix<T>& m)
{
	if(m.xDim == m.yDim)
	{
		return true;
	}
	else
	{
		return false;
	}
}





/// \brief are these matrices (and in this order) are compatible for matrix multiplication
template< class T >
bool Matrix<T>::IsM2MCompatible(Matrix<T>& m1, Matrix<T>& m2)
{
	if(m1.xDim == m2.yDim)
	{
		return true;
	}
	else
	{
		return false;
	}
}


/// \brief Matrix multiplication. Also used to overload & operator.
template< class T >
Matrix<T> Matrix<T>::MMultiply(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsM2MCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible for matrix multiplication!");
	}

	Matrix<T> temp(m1.yDim, m2.xDim);
	temp.Init(0);
	for(int i=0; i<m1.yDim; i++)
	{
		for(int j=0; j<m2.xDim; j++)
		{
			for(int k = 0; k<m1.xDim; k++)
			{
				temp[j][i] += m1[k][i]*m2[j][k]; 
			}
		}
	}
	
	return temp;
}


/// \brief Matrix to Vector multiplication. Vector is taken as a column vector.
template< class T >
Matrix<T> Matrix<T>::MMultiply(Matrix<T>& m1, Vector<T>& m2)
{
	if(m1.xDim != m2.length)
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible for matrix multiplication!");
	}

	Matrix<T> temp(m1.yDim, 1);
	temp.Init(0);
	for(int i=0; i<m1.yDim; i++)
	{
		for(int j=0; j<m2.length; j++)
		{
			temp.data[i] += m1[j][i]*m2[j]; 
		}
	}
	
	return temp;
}

/// \brief Vector to Matrix multiplication. Vector is taken as a row vector.
template< class T >
Matrix<T> Matrix<T>::MMultiply(Vector<T>& m1, Matrix<T>& m2)
{
	if(m2.Rows() != m1.length)
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible for matrix multiplication!");
	}

	Matrix<T> temp(1, m2.Columns());
	temp.Init(0);
	for(int i=0; i<m2.Columns(); i++)
	{
		for(int j=0; j<m1.length; j++)
		{
				temp.data[i] += m2[i][j]*m1[j]; 
		}
	}
	
	return temp;
}


/// \brief Transpose of a matrix. Allocates a new Matrix.
template< class T >
Matrix<T> Matrix<T>::Transpose(Matrix<T>& m)
{
	Matrix<T> temp(m.Columns(), m.Rows());
	for(int i=0; i<m.Rows(); i++)
	{
		for(int j=0; j<m.Columns(); j++)
		{
			temp[i][j] = m[j][i]; 
		}
	}
	return temp;
}


/// \brief Transpose of a matrix. Does not allocate a new Matrix. Creates transpose inline.
template< class T >
Matrix<T>& Matrix<T>::Transpose()
{
	if(!Matrix<T>::IsSquare(*this))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix is not square!");
	}
	
	for(int i=0; i<yDim; i++)
	{
		for(int j=0; j<i; j++)
		{
			T temp = columns[j][i];
			columns[j][i] = columns[i][j];
			columns[i][j] = temp;
		}
	}
	return *this;

}


/// \brief Check if two matrices have the same sizes.
template< class T >
bool Matrix<T>::IsCompatible(Matrix<T>& m1, Matrix<T>& m2)
{
	if(m1.xDim != m2.xDim || m1.yDim != m2.yDim)
	{
		return false;
	}
	else
	{
		return true;
	}

}


// //////////////////
// Boolean Operations...
// //////////////////

/// \brief Elementwise AND operator. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::And(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] != 0 && m2.data[i] != 0) ? 1 : 0;
	}
	
	return temp;
}


/// \brief Elementwise OR operator. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Or(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] == 0 && m2.data[i] == 0) ? 0 : 1;
	}
	
	return temp;
}

/// \brief Elementwise < (less than) operator. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Lt(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] < m2.data[i]) ? 1 : 0;
	}
	
	return temp;
}

/// \brief Elementwise > (greater than) operator. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Gt(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] > m2.data[i]) ? 1 : 0;
	}
	
	return temp;
}

/// \brief Elementwise <= (less than or equal) operator. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Le(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] <= m2.data[i]) ? 1 : 0;
	}
	
	return temp;
}

/// \brief Elementwise >= (greater than or equal) operator. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Ge(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] >= m2.data[i]) ? 1 : 0;
	}
	
	return temp;
}

/// \brief Elementwise == operator. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Eq(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] == m2.data[i]) ? 1 : 0;
	}
	
	return temp;
}

/// \brief Elementwise != (not equal) operator. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Ne(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] != m2.data[i]) ? 1 : 0;
	}
	
	return temp;
}




// ////////////////
// Boolean Operations with value types...
// ////////////////


/// \brief AND operator between Matrix elements and a value. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::And(Matrix<T>& m1, T v)
{
	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] != 0 && v != 0) ? 1 : 0;
	}
	
	return temp;
}


/// \brief OR operator between Matrix elements and a value. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Or(Matrix<T>& m1, T v)
{
	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] == 0 && v == 0) ? 0 : 1;
	}
	
	return temp;
}


/// \brief < (Less than) operator between Matrix elements and a value. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Lt(Matrix<T>& m1, T v)
{
	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] < v) ? 1 : 0;
	}
	
	return temp;
}


/// \brief > (Greater than) operator between Matrix elements and a value. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Gt(Matrix<T>& m1, T v)
{
	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] > v) ? 1 : 0;
	}
	
	return temp;
}


/// \brief <= (Less than or equal) operator between Matrix elements and a value. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Le(Matrix<T>& m1, T v)
{
	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] <= v) ? 1 : 0;
	}
	
	return temp;
}


/// \brief <= (Greater than or equal) operator between Matrix elements and a value. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Ge(Matrix<T>& m1, T v)
{
	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] >= v) ? 1 : 0;
	}
	
	return temp;
}


/// \brief == operator between Matrix elements and a value. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Eq(Matrix<T>& m1, T v)
{
	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] == v) ? 1 : 0;
	}
	
	return temp;
}


/// \brief != (Not equal) operator between Matrix elements and a value. Returns 1 or 0 for each element.
template< class T >
Matrix<int> Matrix<T>::Ne(Matrix<T>& m1, T v)
{
	Matrix<int> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.Length();i++)
	{
		temp.Data()[i] = (m1.data[i] != v) ? 1 : 0;
	}
	
	return temp;
}











// /////////////////
// Elementwise matrix arithmetic
// ////////////////


template< class T >
Matrix<T> Matrix<T>::Add(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<T> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.length;i++)
	{
		temp.data[i] = m1.data[i] + m2.data[i];
	}
	
	return temp;
}


template< class T >
Matrix<T> Matrix<T>::Subtract(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<T> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.length;i++)
	{
		temp.data[i] = m1.data[i] - m2.data[i];
	}
	
	return temp;
}

template< class T >
Matrix<T> Matrix<T>::Multiply(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<T> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.length;i++)
	{
		temp.data[i] = m1.data[i] * m2.data[i];
	}
	
	return temp;
}

template< class T >
Matrix<T> Matrix<T>::Divide(Matrix<T>& m1, Matrix<T>& m2)
{
	if(!Matrix<T>::IsCompatible(m1, m2))
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Matrix sizes are not compatible!");
	}

	Matrix<T> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.length;i++)
	{
		//FW: catch division exception instead...
		if(m2.data[i] != 0)
		{
			temp.data[i] = m1.data[i] / m2.data[i];
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Divide by zero in matrix division!");
		}
	}
	
	return temp;
}

// ////////
// Arithmetic operations between a matrix and a value.
// ////////


template< class T >
Matrix<T> Matrix<T>::Add(Matrix<T>& m1, T v2)
{
	Matrix<T> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.length;i++)
	{
		temp.data[i] = m1.data[i] + v2;
	}
	return temp;
}

template< class T >
Matrix<T> Matrix<T>::Subtract(Matrix<T>& m1, T v2)
{
	Matrix<T> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.length;i++)
	{
		temp.data[i] = m1.data[i] - v2;
	}
	return temp;
}

template< class T >
Matrix<T> Matrix<T>::Subtract(T v2, Matrix<T>& m1)
{
	Matrix<T> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.length;i++)
	{
		temp.data[i] = v2 - m1.data[i];
	}
	return temp;
}


template< class T >
Matrix<T> Matrix<T>::Multiply(Matrix<T>& m1, T v2)
{
	Matrix<T> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.length;i++)
	{
		temp.data[i] = m1.data[i] * v2;
	}
	return temp;
}

template< class T >
Matrix<T> Matrix<T>::Divide(Matrix<T>& m1, T v2)
{
	if(v2 == 0)
	{
		cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
		Utility::RunTimeError("Divide by zero in matrix by value division!");
	}

	Matrix<T> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.length;i++)
	{
		temp.data[i] = m1.data[i] / v2;
	}
	return temp;
}

template< class T >
Matrix<T> Matrix<T>::Divide(T v2, Matrix<T>& m1)
{
	Matrix<T> temp(m1.yDim, m1.xDim);
	for(int i=0;i<temp.length;i++)
	{
		if(m1.data[i] != 0)
		{
			temp.data[i] = v2 / m1.data[i];
		}
		else
		{
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
			Utility::RunTimeError("Divide by zero in value by matrix division!");
		}
	}
	return temp;
}


// //////////
// Elementwise matrix arithmetic with "this" matrix
// /////////