#ifndef _TCUMATRIX_HPP
#define _TCUMATRIX_HPP
/*-------------------------------------------------------------------------
 * drawElements Quality Program Tester Core
 * ----------------------------------------
 *
 * Copyright 2014 The Android Open Source Project
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 *//*!
 * \file
 * \brief Templatized matrix class.
 *//*--------------------------------------------------------------------*/

#include "tcuDefs.hpp"
#include "tcuVector.hpp"
#include "tcuArray.hpp"

namespace tcu
{

// Templated matrix class.
template <typename T, int Rows, int Cols>
class Matrix
{
public:
	typedef	Vector<T, Rows>			Element;
	typedef	T						Scalar;

	enum
	{
		SIZE = Cols,
		ROWS = Rows,
		COLS = Cols,
	};

									Matrix				(void);
	explicit						Matrix				(const T& src);
	explicit						Matrix				(const T src[Rows*Cols]);
									Matrix				(const Vector<T, Rows>& src);
									Matrix				(const Matrix<T, Rows, Cols>& src);
									~Matrix				(void);

	Matrix<T, Rows, Cols>&			operator=			(const Matrix<T, Rows, Cols>& src);
	Matrix<T, Rows, Cols>&			operator*=			(const Matrix<T, Rows, Cols>& src);

	void							setRow				(int rowNdx, const Vector<T, Cols>& vec);
	void							setColumn			(int colNdx, const Vector<T, Rows>& vec);

	Vector<T, Cols>					getRow				(int ndx) const;
	Vector<T, Rows>&				getColumn			(int ndx);
	const Vector<T, Rows>&			getColumn			(int ndx) const;

	Vector<T, Rows>&				operator[]			(int ndx)					{ return getColumn(ndx);	}
	const Vector<T, Rows>&			operator[]			(int ndx) const				{ return getColumn(ndx);	}

	inline const T&					operator()			(int row, int col) const	{ return m_data[col][row];	}
	inline T&						operator()			(int row, int col)			{ return m_data[col][row];	}

	Array<T, Rows*Cols>				getRowMajorData		(void) const;
	Array<T, Rows*Cols>				getColumnMajorData	(void) const;

private:
	Vector<Vector<T, Rows>, Cols>	m_data;
} DE_WARN_UNUSED_TYPE;

// Operators.

// Mat * Mat.
template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b);

// Mat * Vec (column vector).
template <typename T, int Rows, int Cols>
Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec);

// Vec * Mat (row vector).
template <typename T, int Rows, int Cols>
Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx);

// Further operations

template <typename T, int Size>
struct SquareMatrixOps
{
	static T						doDeterminant	(const Matrix<T, Size, Size>& mat);
	static Matrix<T, Size, Size>	doInverse		(const Matrix<T, Size, Size>& mat);
};

template <typename T>
struct SquareMatrixOps<T, 2>
{
	static T						doDeterminant	(const Matrix<T, 2, 2>& mat);
	static Matrix<T, 2, 2>			doInverse		(const Matrix<T, 2, 2>& mat);
};

template <typename T>
struct SquareMatrixOps<T, 3>
{
	static T						doDeterminant	(const Matrix<T, 3, 3>& mat);
	static Matrix<T, 3, 3>			doInverse		(const Matrix<T, 3, 3>& mat);
};

template <typename T>
struct SquareMatrixOps<T, 4>
{
	static T						doDeterminant	(const Matrix<T, 4, 4>& mat);
	static Matrix<T, 4, 4>			doInverse		(const Matrix<T, 4, 4>& mat);
};

namespace matrix
{

template <typename T, int Size>
T determinant (const Matrix<T, Size, Size>& mat)
{
	return SquareMatrixOps<T, Size>::doDeterminant(mat);
}

template <typename T, int Size>
Matrix<T, Size, Size> inverse (const Matrix<T, Size, Size>& mat)
{
	return SquareMatrixOps<T, Size>::doInverse(mat);
}

} // matrix

// Template implementations.

template <typename T>
T SquareMatrixOps<T, 2>::doDeterminant (const Matrix<T, 2, 2>& mat)
{
	return mat(0,0) * mat(1,1) - mat(1,0) * mat(0,1);
}

template <typename T>
T SquareMatrixOps<T, 3>::doDeterminant (const Matrix<T, 3, 3>& mat)
{
	return	+ mat(0,0) * mat(1,1) * mat(2,2)
			+ mat(0,1) * mat(1,2) * mat(2,0)
			+ mat(0,2) * mat(1,0) * mat(2,1)
			- mat(0,0) * mat(1,2) * mat(2,1)
			- mat(0,1) * mat(1,0) * mat(2,2)
			- mat(0,2) * mat(1,1) * mat(2,0);
}

template <typename T>
T SquareMatrixOps<T, 4>::doDeterminant (const Matrix<T, 4, 4>& mat)
{
	using matrix::determinant;

	const T minorMatrices[4][3*3] =
	{
		{
			mat(1,1),	mat(2,1),	mat(3,1),
			mat(1,2),	mat(2,2),	mat(3,2),
			mat(1,3),	mat(2,3),	mat(3,3),
		},
		{
			mat(1,0),	mat(2,0),	mat(3,0),
			mat(1,2),	mat(2,2),	mat(3,2),
			mat(1,3),	mat(2,3),	mat(3,3),
		},
		{
			mat(1,0),	mat(2,0),	mat(3,0),
			mat(1,1),	mat(2,1),	mat(3,1),
			mat(1,3),	mat(2,3),	mat(3,3),
		},
		{
			mat(1,0),	mat(2,0),	mat(3,0),
			mat(1,1),	mat(2,1),	mat(3,1),
			mat(1,2),	mat(2,2),	mat(3,2),
		}
	};

	return	+ mat(0,0) * determinant(Matrix<T, 3, 3>(minorMatrices[0]))
			- mat(0,1) * determinant(Matrix<T, 3, 3>(minorMatrices[1]))
			+ mat(0,2) * determinant(Matrix<T, 3, 3>(minorMatrices[2]))
			- mat(0,3) * determinant(Matrix<T, 3, 3>(minorMatrices[3]));
}

template <typename T>
Matrix<T, 2, 2> SquareMatrixOps<T, 2>::doInverse (const Matrix<T, 2, 2>& mat)
{
	using matrix::determinant;

	const T			det		= determinant(mat);
	Matrix<T, 2, 2>	retVal;

	retVal(0, 0) =  mat(1, 1) / det;
	retVal(0, 1) = -mat(0, 1) / det;
	retVal(1, 0) = -mat(1, 0) / det;
	retVal(1, 1) =  mat(0, 0) / det;

	return retVal;
}

template <typename T>
Matrix<T, 3, 3> SquareMatrixOps<T, 3>::doInverse (const Matrix<T, 3, 3>& mat)
{
	// Blockwise inversion
	using matrix::inverse;

	const T areaA[2*2] =
	{
		mat(0,0),	mat(0,1),
		mat(1,0),	mat(1,1)
	};
	const T areaB[2] =
	{
		mat(0,2),
		mat(1,2),
	};
	const T areaC[2] =
	{
		mat(2,0),	mat(2,1),
	};
	const T areaD[1] =
	{
		mat(2,2)
	};
	const T nullField[4] = { T(0.0f) };

	const Matrix<T, 2, 2>	invA = inverse(Matrix<T, 2, 2>(areaA));
	const Matrix<T, 2, 1>	matB =         Matrix<T, 2, 1>(areaB);
	const Matrix<T, 1, 2>	matC =         Matrix<T, 1, 2>(areaC);
	const Matrix<T, 1, 1>	matD =         Matrix<T, 1, 1>(areaD);

	const T					schurComplement = T(1.0f) / (matD - matC*invA*matB)(0,0);
	const Matrix<T, 2, 2>	zeroMat         = Matrix<T, 2, 2>(nullField);

	const Matrix<T, 2, 2>	blockA = invA + invA*matB*schurComplement*matC*invA;
	const Matrix<T, 2, 1>	blockB = (zeroMat-invA)*matB*schurComplement;
	const Matrix<T, 1, 2>	blockC = matC*invA*(-schurComplement);
	const T					blockD = schurComplement;

	const T result[3*3] =
	{
		blockA(0,0),	blockA(0,1),	blockB(0,0),
		blockA(1,0),	blockA(1,1),	blockB(1,0),
		blockC(0,0),	blockC(0,1),	blockD,
	};

	return Matrix<T, 3, 3>(result);
}

template <typename T>
Matrix<T, 4, 4> SquareMatrixOps<T, 4>::doInverse (const Matrix<T, 4, 4>& mat)
{
	// Blockwise inversion
	using matrix::inverse;

	const T areaA[2*2] =
	{
		mat(0,0),	mat(0,1),
		mat(1,0),	mat(1,1)
	};
	const T areaB[2*2] =
	{
		mat(0,2),	mat(0,3),
		mat(1,2),	mat(1,3)
	};
	const T areaC[2*2] =
	{
		mat(2,0),	mat(2,1),
		mat(3,0),	mat(3,1)
	};
	const T areaD[2*2] =
	{
		mat(2,2),	mat(2,3),
		mat(3,2),	mat(3,3)
	};
	const T nullField[4] = { T(0.0f) };

	const Matrix<T, 2, 2> invA = inverse(Matrix<T, 2, 2>(areaA));
	const Matrix<T, 2, 2> matB =         Matrix<T, 2, 2>(areaB);
	const Matrix<T, 2, 2> matC =         Matrix<T, 2, 2>(areaC);
	const Matrix<T, 2, 2> matD =         Matrix<T, 2, 2>(areaD);

	const Matrix<T, 2, 2> schurComplement = inverse(matD - matC*invA*matB);
	const Matrix<T, 2, 2> zeroMat         = Matrix<T, 2, 2>(nullField);

	const Matrix<T, 2, 2> blockA = invA + invA*matB*schurComplement*matC*invA;
	const Matrix<T, 2, 2> blockB = (zeroMat-invA)*matB*schurComplement;
	const Matrix<T, 2, 2> blockC = (zeroMat-schurComplement)*matC*invA;
	const Matrix<T, 2, 2> blockD = schurComplement;

	const T result[4*4] =
	{
		blockA(0,0),	blockA(0,1),	blockB(0,0),	blockB(0,1),
		blockA(1,0),	blockA(1,1),	blockB(1,0),	blockB(1,1),
		blockC(0,0),	blockC(0,1),	blockD(0,0),	blockD(0,1),
		blockC(1,0),	blockC(1,1),	blockD(1,0),	blockD(1,1),
	};

	return Matrix<T, 4, 4>(result);
}

// Initialize to identity.
template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols>::Matrix (void)
{
	for (int row = 0; row < Rows; row++)
		for (int col = 0; col < Cols; col++)
			(*this)(row, col) = (row == col) ? T(1) : T(0);
}

// Initialize to diagonal matrix.
template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols>::Matrix (const T& src)
{
	for (int row = 0; row < Rows; row++)
		for (int col = 0; col < Cols; col++)
			(*this)(row, col) = (row == col) ? src : T(0);
}

// Initialize from data array.
template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols>::Matrix (const T src[Rows*Cols])
{
	for (int row = 0; row < Rows; row++)
		for (int col = 0; col < Cols; col++)
			(*this)(row, col) = src[row*Cols + col];
}

// Initialize to diagonal matrix.
template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols>::Matrix (const Vector<T, Rows>& src)
{
	DE_STATIC_ASSERT(Rows == Cols);
	for (int row = 0; row < Rows; row++)
		for (int col = 0; col < Cols; col++)
			(*this)(row, col) = (row == col) ? src.m_data[row] : T(0);
}

// Copy constructor.
template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols>::Matrix (const Matrix<T, Rows, Cols>& src)
{
	*this = src;
}

// Destructor.
template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols>::~Matrix (void)
{
}

// Assignment operator.
template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator= (const Matrix<T, Rows, Cols>& src)
{
	for (int row = 0; row < Rows; row++)
		for (int col = 0; col < Cols; col++)
			(*this)(row, col) = src(row, col);
	return *this;
}

// Multipy and assign op
template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator*= (const Matrix<T, Rows, Cols>& src)
{
	*this = *this * src;
	return *this;
}

template <typename T, int Rows, int Cols>
void Matrix<T, Rows, Cols>::setRow (int rowNdx, const Vector<T, Cols>& vec)
{
	for (int col = 0; col < Cols; col++)
		(*this)(rowNdx, col) = vec.m_data[col];
}

template <typename T, int Rows, int Cols>
void Matrix<T, Rows, Cols>::setColumn (int colNdx, const Vector<T, Rows>& vec)
{
	m_data[colNdx] = vec;
}

template <typename T, int Rows, int Cols>
Vector<T, Cols> Matrix<T, Rows, Cols>::getRow (int rowNdx) const
{
	Vector<T, Cols> res;
	for (int col = 0; col < Cols; col++)
		res[col] = (*this)(rowNdx, col);
	return res;
}

template <typename T, int Rows, int Cols>
Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx)
{
	return m_data[colNdx];
}

template <typename T, int Rows, int Cols>
const Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx) const
{
	return m_data[colNdx];
}

template <typename T, int Rows, int Cols>
Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getColumnMajorData (void) const
{
	Array<T, Rows*Cols> a;
	T* dst = a.getPtr();
	for (int col = 0; col < Cols; col++)
		for (int row = 0; row < Rows; row++)
			*dst++ = (*this)(row, col);
	return a;
}

template <typename T, int Rows, int Cols>
Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getRowMajorData (void) const
{
	Array<T, Rows*Cols> a;
	T* dst = a.getPtr();
	for (int row = 0; row < Rows; row++)
		for (int col = 0; col < Cols; col++)
			*dst++ = (*this)(row, col);
	return a;
}

// Multiplication of two matrices.
template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b)
{
	DE_STATIC_ASSERT(Cols0 == Rows1);
	Matrix<T, Rows0, Cols1> res;
	for (int row = 0; row < Rows0; row++)
	{
		for (int col = 0; col < Cols1; col++)
		{
			T v = T(0);
			for (int ndx = 0; ndx < Cols0; ndx++)
				v += a(row,ndx) * b(ndx,col);
			res(row,col) = v;
		}
	}
	return res;
}

// Multiply of matrix with column vector.
template <typename T, int Rows, int Cols>
Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec)
{
	Vector<T, Rows> res;
	for (int row = 0; row < Rows; row++)
	{
		T v = T(0);
		for (int col = 0; col < Cols; col++)
			v += mtx(row,col) * vec.m_data[col];
		res.m_data[row] = v;
	}
	return res;
}

// Multiply of matrix with row vector.
template <typename T, int Rows, int Cols>
Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx)
{
	Vector<T, Cols> res;
	for (int col = 0; col < Cols; col++)
	{
		T v = T(0);
		for (int row = 0; row < Rows; row++)
			v += mtx(row,col) * vec.m_data[row];
		res.m_data[col] = v;
	}
	return res;
}

// Common typedefs.
typedef Matrix<float, 2, 2>		Matrix2f;
typedef Matrix<float, 3, 3>		Matrix3f;
typedef Matrix<float, 4, 4>		Matrix4f;
typedef Matrix<double, 2, 2>	Matrix2d;
typedef Matrix<double, 3, 3>	Matrix3d;
typedef Matrix<double, 4, 4>	Matrix4d;

// GLSL-style naming \note CxR.
typedef Matrix2f			Mat2;
typedef Matrix<float, 3, 2>	Mat2x3;
typedef Matrix<float, 4, 2>	Mat2x4;
typedef Matrix<float, 2, 3>	Mat3x2;
typedef Matrix3f			Mat3;
typedef Matrix<float, 4, 3>	Mat3x4;
typedef Matrix<float, 2, 4>	Mat4x2;
typedef Matrix<float, 3, 4>	Mat4x3;
typedef Matrix4f			Mat4;

// Matrix-scalar operators.

template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& mtx, T scalar)
{
	Matrix<T, Rows, Cols> res;
	for (int col = 0; col < Cols; col++)
		for (int row = 0; row < Rows; row++)
			res(row, col) = mtx(row, col) + scalar;
	return res;
}

template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols> operator- (const Matrix<T, Rows, Cols>& mtx, T scalar)
{
	Matrix<T, Rows, Cols> res;
	for (int col = 0; col < Cols; col++)
		for (int row = 0; row < Rows; row++)
			res(row, col) = mtx(row, col) - scalar;
	return res;
}

template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols> operator* (const Matrix<T, Rows, Cols>& mtx, T scalar)
{
	Matrix<T, Rows, Cols> res;
	for (int col = 0; col < Cols; col++)
		for (int row = 0; row < Rows; row++)
			res(row, col) = mtx(row, col) * scalar;
	return res;
}

template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols> operator/ (const Matrix<T, Rows, Cols>& mtx, T scalar)
{
	Matrix<T, Rows, Cols> res;
	for (int col = 0; col < Cols; col++)
		for (int row = 0; row < Rows; row++)
			res(row, col) = mtx(row, col) / scalar;
	return res;
}

// Matrix-matrix component-wise operators.

template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
{
	Matrix<T, Rows, Cols> res;
	for (int col = 0; col < Cols; col++)
		for (int row = 0; row < Rows; row++)
			res(row, col) = a(row, col) + b(row, col);
	return res;
}

template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols> operator- (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
{
	Matrix<T, Rows, Cols> res;
	for (int col = 0; col < Cols; col++)
		for (int row = 0; row < Rows; row++)
			res(row, col) = a(row, col) - b(row, col);
	return res;
}

template <typename T, int Rows, int Cols>
Matrix<T, Rows, Cols> operator/ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
{
	Matrix<T, Rows, Cols> res;
	for (int col = 0; col < Cols; col++)
		for (int row = 0; row < Rows; row++)
			res(row, col) = a(row, col) / b(row, col);
	return res;
}

} // tcu

#endif // _TCUMATRIX_HPP