template <typename ftype, unsigned nHeight, unsigned nWidth>
class TMatrix
{
public:
	typedef TMatrix<ftype, nHeight, nWidth> Matrix;
	TMatrix (const ftype* pInit)
	{
		(*this) = pInit;
	}

	TMatrix (float fInit, ...)
	{
		(*this) = &fInit;
	}

	TMatrix ()
	{
	}

	Matrix& operator = (const ftype* pInit)
	{
		memcpy (mx, pInit, sizeof(ftype)*nHeight*nWidth);
		return *this;
	}

	void identity()
	{
		assert (nWidth == nHeight);
		unsigned i;
		memset (mx, 0, sizeof(mx));
		for (i = 0; i < nHeight; ++i)
			mx[i][i] = 1;
	}

	ftype* operator [] (unsigned nRow)
	{
		assert (nRow < nHeight);
		return mx[nRow];
	}

	void invert ()
	{
		// universal code to invert
	}

	const ftype* operator [] (unsigned nRow)const
	{
		assert (nRow < nHeight);
		return mx[nRow];
	}

	template <unsigned l>
	void assignProduct (const TMatrix<ftype, nHeight, l>& mxLeft, const TMatrix<ftype,l, nWidth>& mxRight)
	{
		unsigned i,j,f;
		for (i = 0; i < nHeight; ++i)
			for (j = 0; j < nWidth; ++j)
			{
				ftype fResult = mxLeft[i][0] * mxRight[0][j];
				for (f = 1; f < l; ++f)
					fResult += mxLeft[i][f] * mxRight[f][j];
				mx[i][j] = fResult;
			}
	}



	Matrix& operator += (TMatrix<ftype, nHeight,nWidth>& mxRight)
	{
		unsigned i, j;
		for (i = 0; i < nHeight; ++i)
			for (j = 0; j < nWidth; ++j)
				mx[i][j] += mxRight[i][j];
		return *this;
	}

	Matrix& operator -= (TMatrix<ftype, nHeight,nWidth>& mxRight)
	{
		unsigned i, j;
		for (i = 0; i < nHeight; ++i)
			for (j = 0; j < nWidth; ++j)
				mx[i][j] -= mxRight[i][j];
		return *this;
	}

protected:
	ftype mx[nHeight][nWidth];
};

typedef TMatrix<float,1,3> Matrix13;
typedef TMatrix<float,3,1> Matrix31;
typedef TMatrix<float,3,3> Matrix33;
typedef TMatrix<float,4,4> Matrix44;
typedef TMatrix<float,3,4> Matrix34;
typedef TMatrix<float,4,3> Matrix43;

inline void Matrix33::invert()
{
	// specialized inversion
}

// multiplication of A (k*l) by B (l*m) = result (k*m)
template <typename ftype,unsigned k, unsigned l, unsigned m>
TMatrix<ftype,k,m> operator * (const TMatrix<ftype,k,l>& mxLeft, const TMatrix<ftype,l,m>& mxRight)
{
	TMatrix<ftype, k, m> mxResult;
	mxResult.assignProduct (mxLeft, mxRight);
	return mxResult;
}

Matrix33 operator ^ (const Matrix33& l, const Matrix33& r)
{
	Matrix33 x;
	x[0][0] = l[0][0] * r[0][0] + l[0][1]*r[1][0] + l[0][2]*r[2][0];
	x[0][1] = l[0][0] * r[0][1] + l[0][1]*r[1][1] + l[0][2]*r[2][1];
	x[0][2] = l[0][0] * r[0][2] + l[0][1]*r[1][2] + l[0][2]*r[2][2];
	x[1][0] = l[1][0] * r[0][0] + l[1][1]*r[1][0] + l[1][2]*r[2][0];
	x[1][1] = l[1][0] * r[0][1] + l[1][1]*r[1][1] + l[1][2]*r[2][1];
	x[1][2] = l[1][0] * r[0][2] + l[1][1]*r[1][2] + l[1][2]*r[2][2];
	x[2][0] = l[2][0] * r[0][0] + l[2][1]*r[1][0] + l[2][2]*r[2][0];
	x[2][1] = l[2][0] * r[0][1] + l[2][1]*r[1][1] + l[2][2]*r[2][1];
	x[2][2] = l[2][0] * r[0][2] + l[2][1]*r[1][2] + l[2][2]*r[2][2];
	return x;
}