#ifndef _MATRIX_H_#define _MATRIX_H_#include <stdint.h>#include <vector>template <class T> class Matrix{	public:		Matrix(uint64_t unNumRows, uint64_t unNumCols, long double dbDefaultElement/* = 0.0*/);		~Matrix();		std::vector<std::vector<long double>> GetCore();		void SetCore(std::vector<std::vector<long double>> core);		void GetDimensions(uint64_t & unNumRows, uint64_t & unNumCols);		std::vector<uint8_t> GetDimensions();		uint64_t GetRowSize();		uint64_t GetColSize();		void SetElement(long double dbElement, uint64_t unRow, uint64_t unCol);		long double GetElement(uint64_t unRow, uint64_t unCol);		Matrix operator=(Matrix & rhs);		Matrix operator+(Matrix & rhs);		Matrix operator-(Matrix & rhs);		Matrix operator*(Matrix & rhs);		static Matrix Transpose(Matrix & matrix);	protected:		std::vector<std::vector<long double>> TheMatrix;};#endif

#include "Matrix.h"template <class T>Matrix<T>::Matrix(uint64_t unNumRows, uint64_t unNumCols, long double dbDefaultElement/* = 0.0*/){	for (uint64_t c=0; c<unNumCols; c++)	{		std::vector<long double> NewCol;		for (uint64_t r=0; r<unNumRows; r++)		{			NewCol.push_back(dbDefaultElement);		}		TheMatrix.push_back(NewCol);	}}template <class T>Matrix<T>::~Matrix(){}template <class T>std::vector<std::vector<long double>> Matrix<T>::GetCore(){	return TheMatrix;}template <class T>void Matrix<T>::SetCore(std::vector<std::vector<long double>> core){	TheMatrix = core;}template <class T>Matrix<T> Matrix<T>::operator=(Matrix & rhs){	TheMatrix = rhs.GetCore();	return *this;}template <class T>Matrix<T> Matrix<T>::operator+(Matrix<T> & rhs){	if (rhs.GetDimensions() != this->GetDimensions())	{		Matrix<T> NullMatrix(0, 0);		return NullMatrix;	}	uint64_t unRows, unCols;	GetDimensions(unRows, unCols);	Matrix<T> ResultMatrix(unRows, unCols);	for (uint64_t j=0; j<unCols; j++)	{		for (uint64_t i=0; i<unRows; i++)		{			ResultMatrix.SetElement(this->GetElement(i, j) + rhs.GetElement(i, j), i, j);		}	}	return ResultMatrix;}template <class T>Matrix<T> Matrix<T>::operator-(Matrix & rhs){	if (rhs.GetDimensions() != this->GetDimensions())	{		Matrix NullMatrix(0, 0);		return NullMatrix;	}	uint64_t unRows, unCols;	GetDimensions(unRows, unCols);	Matrix<T> ResultMatrix(unRows, unCols);	for (uint64_t j=0; j<unCols; j++)	{		for (uint64_t i=0; i<unRows; i++)		{			ResultMatrix.SetElement(this->GetElement(i, j) - rhs.GetElement(i, j), i, j);		}	}	return ResultMatrix;}template <class T>Matrix<T> Matrix<T>::operator*(Matrix<T> & rhs){	uint64_t unRowsLeft, unColsLeft, unRowsRight, unColsRight;	rhs.GetDimensions(unRowsRight, unColsRight);	this->GetDimensions(unRowsLeft, unColsLeft);	if (unColsLeft != unRowsRight)	{		Matrix<T> NullMatrix(0, 0);		return NullMatrix;	}	Matrix<T> ResultMatrix(unRowsLeft, unColsRight);	long double dbSum;	for (uint64_t i=0; i<unRowsLeft; i++)	{		for (uint64_t j=0; j<unColsRight; j++)		{			dbSum = 0.0;			for (uint64_t k=0; k<unRowsRight; k++)	// or unColsLeft, they are same			{				dbSum += GetElement(i, k) * rhs.GetElement(k, j);			}			ResultMatrix.SetElement(dbSum, i, j);		}	}	return ResultMatrix;}template <class T>void Matrix<T>::GetDimensions(uint64_t & unNumRows, uint64_t & unNumCols){	unNumCols = TheMatrix.size();	if (unNumCols == 0)	{		unNumRows = 0;	}	else	{		unNumRows = TheMatrix.at(0).size();	}}template <class T>uint64_t Matrix<T>::GetRowSize(){	if (GetColSize() == 0) return 0;	return static_cast<uint64_t>(TheMatrix.at(0).size());}template <class T>uint64_t Matrix<T>::GetColSize(){	return static_cast<uint64_t>(TheMatrix.size());}template <class T>std::vector<uint8_t> Matrix<T>::GetDimensions(){	uint64_t unNumRows, unNumCols;	GetDimensions(unNumRows, unNumCols);	std::vector<uint8_t> Dimensions;	Dimensions.push_back(unNumRows);	Dimensions.push_back(unNumCols);	return Dimensions;}template <class T>void Matrix<T>::SetElement(long double dbElement, uint64_t unRow, uint64_t unCol){	if ((unRow >= GetRowSize()) || (unCol >= GetColSize())) return;	TheMatrix.at(unCol).at(unRow) = dbElement;}template <class T>long double Matrix<T>::GetElement(uint64_t unRow, uint64_t unCol){	if (unRow >= GetRowSize()) return 0.0;	if (unCol >= GetColSize()) return 0.0;	return TheMatrix.at(unCol).at(unRow);}template <class T>Matrix<T> Matrix<T>::Transpose(Matrix & matrix){	uint64_t rows, cols;	matrix.GetDimensions(rows, cols);	Matrix TransposedMatrix(cols, rows);	for (uint64_t c=0; c<rows; c++)	{		for (uint64_t r=0; r<cols; r++)		{			TransposedMatrix.SetElement(matrix.GetElement(r, c), c, r);		}	}	return TransposedMatrix;}

#include <iostream>#include "Matrix.h"int wmain(int argc, wchar_t * argv[]){	Matrix<long double> mat1(3, 3, 0.0);	mat1.SetElement(11, 0, 0);	system("pause");	return 0;}

// Symbols that begin with an underscore followed by a capital letter// are reserved by the implementation. Although a collision is probably// unlikely, there's really no need for the leading underscore, so you// might as well play it safe and drop it.#ifndef _MATRIX_H_#define _MATRIX_H_#include <stdint.h>#include <vector>template <class T> class Matrix{    public:        Matrix(uint64_t unNumRows, uint64_t unNumCols, long double dbDefaultElement/* = 0.0*/);        ~Matrix();        // I'd recommend either using something like boost::multi_array,        // or just using a single vector (a 1-d array) for storage.        // (You can find plenty of discussions on the forums of why        // 'vector of vector' storage is usually best avoided.)                // Also, note that you're returning a complete copy of the matrix        // data here, which may involve memory allocations, etc. Maybe        // that's what you intend, but I thought I'd point it out just in        // case.        std::vector<std::vector<long double>> GetCore();                // The input container should probably be passed by constant        // reference here to avoid an unnecessary copy.        void SetCore(std::vector<std::vector<long double>> core);        void GetDimensions(uint64_t & unNumRows, uint64_t & unNumCols);                // Returning a resizable container such as 'vector' here is really        // an unnecessary expense. It would make more sense to return a simple        // object such as std::pair, a tuple, or a custom struct/class, or        // a fixed-sized container such as array.        std::vector<uint8_t> GetDimensions();        uint64_t GetRowSize();        uint64_t GetColSize();        // It's fairly typical to overload the () operator for element access        // in a matrix class (which makes for cleaner and more clear code than        // named functions, IMO).        void SetElement(long double dbElement, uint64_t unRow, uint64_t unCol);        long double GetElement(uint64_t unRow, uint64_t unCol);        // Conventional wisdom is that operations such as matrix addition        // should be implemented as binary non-member functions rather than        // as member functions.                // Also, all of these arguments should probably be passed by constant        // reference.        Matrix operator=(Matrix & rhs);        Matrix operator+(Matrix & rhs);        Matrix operator-(Matrix & rhs);        Matrix operator*(Matrix & rhs);        // Again, argument should probably be passed by constant reference.        static Matrix Transpose(Matrix & matrix);    protected:                // Any particular reason this is protected and not private?        std::vector<std::vector<long double>> TheMatrix;};#endif

#include "Matrix.h"template <class T>Matrix<T>::Matrix(uint64_t unNumRows, uint64_t unNumCols, long double dbDefaultElement/* = 0.0*/){    // Again, a different form of storage would probably be preferable,    // but note that this:        /*    for (uint64_t c=0; c<unNumCols; c++)    {        std::vector<long double> NewCol;        for (uint64_t r=0; r<unNumRows; r++)        {            NewCol.push_back(dbDefaultElement);        }        TheMatrix.push_back(NewCol);    }    */        // Can be more succinctly (and efficiently) written as (untested):    TheMatrix.resize(unNumCols, std::vector(unNumRows, dbDefaultElement));}// Don't implement empty non-virtual destructors (the compiler-generated// destructor will suffice in this case).template <class T>Matrix<T>::~Matrix(){}template <class T>std::vector<std::vector<long double>> Matrix<T>::GetCore(){    return TheMatrix;}template <class T>void Matrix<T>::SetCore(std::vector<std::vector<long double>> core){    TheMatrix = core;}// Shouldn't need to implement this either (the compiler should generate// an equivalent assignment operator, IIRC).template <class T>Matrix<T> Matrix<T>::operator=(Matrix & rhs){    TheMatrix = rhs.GetCore();    return *this;}template <class T>Matrix<T> Matrix<T>::operator+(Matrix<T> & rhs){    if (rhs.GetDimensions() != this->GetDimensions())    {        Matrix<T> NullMatrix(0, 0);        return NullMatrix;    }    uint64_t unRows, unCols;    GetDimensions(unRows, unCols);    Matrix<T> ResultMatrix(unRows, unCols);    for (uint64_t j=0; j<unCols; j++)    {        for (uint64_t i=0; i<unRows; i++)        {            ResultMatrix.SetElement(this->GetElement(i, j) + rhs.GetElement(i, j), i, j);        }    }    return ResultMatrix;}template <class T>Matrix<T> Matrix<T>::operator-(Matrix & rhs){    if (rhs.GetDimensions() != this->GetDimensions())    {        Matrix NullMatrix(0, 0);        return NullMatrix;    }    uint64_t unRows, unCols;    GetDimensions(unRows, unCols);    Matrix<T> ResultMatrix(unRows, unCols);    for (uint64_t j=0; j<unCols; j++)    {        for (uint64_t i=0; i<unRows; i++)        {            ResultMatrix.SetElement(this->GetElement(i, j) - rhs.GetElement(i, j), i, j);        }    }    return ResultMatrix;}template <class T>Matrix<T> Matrix<T>::operator*(Matrix<T> & rhs){    uint64_t unRowsLeft, unColsLeft, unRowsRight, unColsRight;    rhs.GetDimensions(unRowsRight, unColsRight);    this->GetDimensions(unRowsLeft, unColsLeft);    if (unColsLeft != unRowsRight)    {        Matrix<T> NullMatrix(0, 0);        return NullMatrix;    }    Matrix<T> ResultMatrix(unRowsLeft, unColsRight);    long double dbSum;    for (uint64_t i=0; i<unRowsLeft; i++)    {        for (uint64_t j=0; j<unColsRight; j++)        {            dbSum = 0.0;            for (uint64_t k=0; k<unRowsRight; k++)    // or unColsLeft, they are same            {                dbSum += GetElement(i, k) * rhs.GetElement(k, j);            }            ResultMatrix.SetElement(dbSum, i, j);        }    }    return ResultMatrix;}template <class T>void Matrix<T>::GetDimensions(uint64_t & unNumRows, uint64_t & unNumCols){    unNumCols = TheMatrix.size();    if (unNumCols == 0)    {        unNumRows = 0;    }    else    {        unNumRows = TheMatrix.at(0).size();    }}template <class T>uint64_t Matrix<T>::GetRowSize(){    if (GetColSize() == 0) return 0;    return static_cast<uint64_t>(TheMatrix.at(0).size());}template <class T>uint64_t Matrix<T>::GetColSize(){    return static_cast<uint64_t>(TheMatrix.size());}template <class T>std::vector<uint8_t> Matrix<T>::GetDimensions(){    uint64_t unNumRows, unNumCols;    GetDimensions(unNumRows, unNumCols);    std::vector<uint8_t> Dimensions;    Dimensions.push_back(unNumRows);    Dimensions.push_back(unNumCols);    return Dimensions;}template <class T>void Matrix<T>::SetElement(long double dbElement, uint64_t unRow, uint64_t unCol){    // at() will always throw on an invalid index (IIRC), which means    // that either there's no point in returning on an invalid index,    // or there's no point in using at().        // Also, silently failing (both here and in GetElement()) probably    // isn't the best choice here. Better, IMO, to let an exception be    // thrown (if you're going to use at(), that is) so that the user will    // know what happened.        // Finally, keep in mind that using at() rather than operator[]() can    // incur a performance penalty, which could be measurable if you're    // accessing elements frequently.    if ((unRow >= GetRowSize()) || (unCol >= GetColSize())) return;    TheMatrix.at(unCol).at(unRow) = dbElement;}template <class T>long double Matrix<T>::GetElement(uint64_t unRow, uint64_t unCol){    if (unRow >= GetRowSize()) return 0.0;    if (unCol >= GetColSize()) return 0.0;    return TheMatrix.at(unCol).at(unRow);}template <class T>Matrix<T> Matrix<T>::Transpose(Matrix & matrix){    uint64_t rows, cols;    matrix.GetDimensions(rows, cols);    Matrix TransposedMatrix(cols, rows);    for (uint64_t c=0; c<rows; c++)    {        for (uint64_t r=0; r<cols; r++)        {            TransposedMatrix.SetElement(matrix.GetElement(r, c), c, r);        }    }    return TransposedMatrix;}

template <class T>Matrix<T>& Matrix<T>::operator -= (const Matrix & rhs){    // add code to subtract rhs from *this    // ...    return *this;}template <class T>const Matrix<T> Matrix<T>::operator - (const Matrix & rhs){    // use copy constructor to create a matrix identical to *this    Matrix<T> result( *this );    // and use opertor -= on the copy to perform the subtraction and return the    // result    return result -= rhs;}

#ifndef MATRIX_H#define MATRIX_H#include <stdint.h>#include <vector>template <class T> class Matrix{	public:		Matrix(uint64_t unNumRows = 0, uint64_t unNumCols = 0);		void GetDimensions(uint64_t & unNumRows, uint64_t & unNumCols);		std::pair<uint64_t, uint64_t> GetDimensions();		void SetDimensions(uint64_t unNumRows, uint64_t unNumCols);		uint64_t GetRowSize();		uint64_t GetColSize();		void SetElement(T dbElement, uint64_t unRow, uint64_t unCol);		T GetElement(uint64_t unRow, uint64_t unCol);		//Matrix operator=(const Matrix & rhs);	// Compiler generate this automatically		Matrix operator+(const Matrix & rhs);		Matrix operator-(const Matrix & rhs);		Matrix operator*(const Matrix & rhs);		Matrix & operator+=(const Matrix & rhs);		Matrix & operator-=(const Matrix & rhs);		Matrix & operator*=(const Matrix & rhs);		static Matrix Transpose	(const Matrix & matrix);		static Matrix Multiply	(const Matrix & LeftMatrix, const Matrix & RightMatrix);		static Matrix Add	(const Matrix & LeftMatrix, const Matrix & RightMatrix);		static Matrix Subtract	(const Matrix & LeftMatrix, const Matrix & RightMatrix);		static Matrix Negate	(const Matrix & matrix);		// TO DO:		static bool IsSquare(const Matrix & matrix);		static bool IsFullRowRank(const Matrix & matrix);		static bool IsFullColRank(const Matrix & matrix);		// TO DO:		static uint64_t GetRowRank(const Matrix & matrix);		static uint64_t GetColRank(const Matrix & matrix);		// These will be transferred into SquareMatrix class:		//static bool IsIdentity(const SquareMatrix & matrix);		//static bool IsDiagonal(const SquareMatrix & matrix);		//static bool IsUpperTriangular(SquareMatrix Matrix & matrix);		//static bool IsLowerTriangular(SquareMatrix Matrix & matrix);		//static bool IsFullRank(const SquareMatrix & matrix);		//static T GetDeterminant(const SquareMatrix & matrix);		//static <std::vector<T>> GetEigenValues(const SquareMatrix & matrix);		//static <std::vector<std::vector<T>>> GetEigenVectors(const SquareMatrix & matrix);		//static SquareMatrix GetInverse(const SquareMatrix & matrix);		//static SquareMatrix GetDiagonalizedForm(const SquareMatrix & matrix);	protected:		std::vector<T> TheMatrix;		uint64_t m_unRowSize;		uint64_t m_unColSize;};#endif//////////////////////////////////////////////////////////////////////////////template <class T>Matrix<T>::Matrix(uint64_t unNumRows/* = 0*/, uint64_t unNumCols/* = 0*/){	SetDimensions(unNumRows, unNumCols);}template <class T>Matrix<T> Matrix<T>::operator+(const Matrix<T> & rhs){	return Add(*this, rhs);}template <class T>Matrix<T> Matrix<T>::operator-(const Matrix<T> & rhs){	return Subtract(*this, rhs);}template <class T>Matrix<T> Matrix<T>::operator*(const Matrix<T> & rhs){	return Multiply(*this, rhs);}template <class T>Matrix<T> & Matrix<T>::operator+=(const Matrix & rhs){	*this = Add(*this, rhs);	return *this;}template <class T>Matrix<T> & Matrix<T>::operator-=(const Matrix & rhs){	*this = Subtract(*this, rhs);	return *this;}template <class T>Matrix<T> & Matrix<T>::operator*=(const Matrix & rhs){	*this = Multiply(*this, rhs);	return *this;}template <class T>void Matrix<T>::GetDimensions(uint64_t & unNumRows, uint64_t & unNumCols){	unNumRows = m_unRowSize;	unNumCols = m_unColSize;}template <class T>void Matrix<T>::SetDimensions(uint64_t unNumRows, uint64_t unNumCols){	m_unRowSize = unNumRows;	m_unColSize = unNumCols;	TheMatrix.resize(m_unRowSize * m_unColSize);}template <class T>uint64_t Matrix<T>::GetRowSize(){	return m_unRowSize;}template <class T>uint64_t Matrix<T>::GetColSize(){	return m_unColSize;}template <class T>std::pair<uint64_t, uint64_t> Matrix<T>::GetDimensions(){	std::pair<uint64_t, uint64_t> Dimensions;	Dimensions.first = GetRowSize();	Dimensions.second = GetColSize();	return Dimensions;}template <class T>void Matrix<T>::SetElement(T Element, uint64_t unRow, uint64_t unCol){	try	{		TheMatrix.at(m_unColSize * unRow + unCol) = Element;	}	catch(std::out_of_range & e)	{		// Do nothing.	}}template <class T>T Matrix<T>::GetElement(uint64_t unRow, uint64_t unCol){	try	{		return TheMatrix.at(m_unColSize * unRow + unCol);	}	catch(std::out_of_range & e)	{		return static_cast<T>(0);	}}template <class T>Matrix<T> Matrix<T>::Transpose(const Matrix & matrix){	uint64_t rows, cols;	const_cast<Matrix<T>&>(matrix).GetDimensions(rows, cols);	Matrix TransposedMatrix(cols, rows);	for (uint64_t c=0; c<rows; c++)	{		for (uint64_t r=0; r<cols; r++)		{			TransposedMatrix.SetElement(const_cast<Matrix<T>&>(matrix).GetElement(r, c), c, r);		}	}	return TransposedMatrix;}template <class T>Matrix<T> Matrix<T>::Multiply(const Matrix & LeftMatrix, const Matrix & RightMatrix){	uint64_t unRowsLeft, unColsLeft, unRowsRight, unColsRight;	const_cast<Matrix<T>&>(RightMatrix).GetDimensions(unRowsRight, unColsRight);	const_cast<Matrix<T>&>(LeftMatrix).GetDimensions(unRowsLeft, unColsLeft);	if (unColsLeft != unRowsRight)	{		Matrix<T> NullMatrix(0, 0);		return NullMatrix;	}	Matrix<T> ResultMatrix(unRowsLeft, unColsRight);	T dbSum;	for (uint64_t i=0; i<unRowsLeft; i++)	{		for (uint64_t j=0; j<unColsRight; j++)		{			dbSum = 0.0;			for (uint64_t k=0; k<unRowsRight; k++)	// or unColsLeft, they are same			{				dbSum += const_cast<Matrix<T>&>(LeftMatrix).GetElement(i, k)				       * const_cast<Matrix<T>&>(RightMatrix).GetElement(k, j);			}			ResultMatrix.SetElement(dbSum, i, j);		}	}	return ResultMatrix;}template <class T>Matrix<T> Matrix<T>::Add(const Matrix & LeftMatrix, const Matrix & RightMatrix){	if (const_cast<Matrix<T>&>(RightMatrix).GetDimensions()		!= const_cast<Matrix<T>&>(LeftMatrix).GetDimensions())	{		Matrix<T> NullMatrix(0, 0);		return NullMatrix;	}	uint64_t unRows, unCols;	const_cast<Matrix<T>&>(LeftMatrix).GetDimensions(unRows, unCols);	Matrix<T> ResultMatrix(unRows, unCols);	for (uint64_t i=0; i<unRows; i++)	{		for (uint64_t j=0; j<unCols; j++)		{			ResultMatrix.SetElement(const_cast<Matrix<T>&>(LeftMatrix).GetElement(i, j)				+ const_cast<Matrix<T>&>(RightMatrix).GetElement(i, j), i, j);		}	}	return ResultMatrix;}template <class T>Matrix<T> Matrix<T>::Subtract(const Matrix & LeftMatrix, const Matrix & RightMatrix){	return Add(RightMatrix, Negate(LeftMatrix));}template <class T>Matrix<T> Matrix<T>::Negate(const Matrix & matrix){	uint64_t unNumRows, unNumCols;	const_cast<Matrix<T>&>(matrix).GetDimensions(unNumRows, unNumCols);	Matrix<T> NegatedMatrix(unNumRows, unNumCols);	for (uint64_t i=0; i<unNumRows; i++)	{		for (uint64_t j=0; j<unNumCols; j++)		{			NegatedMatrix.SetElement(-const_cast<Matrix<T>&>(matrix).GetElement(i, j), i, j);		}	}	return NegatedMatrix;}