1#ifndef _TCUMATRIX_HPP
2#define _TCUMATRIX_HPP
3/*-------------------------------------------------------------------------
4 * drawElements Quality Program Tester Core
5 * ----------------------------------------
6 *
7 * Copyright 2014 The Android Open Source Project
8 *
9 * Licensed under the Apache License, Version 2.0 (the "License");
10 * you may not use this file except in compliance with the License.
11 * You may obtain a copy of the License at
12 *
13 *      http://www.apache.org/licenses/LICENSE-2.0
14 *
15 * Unless required by applicable law or agreed to in writing, software
16 * distributed under the License is distributed on an "AS IS" BASIS,
17 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18 * See the License for the specific language governing permissions and
19 * limitations under the License.
20 *
21 *//*!
22 * \file
23 * \brief 	Templatized matrix class.
24 *//*--------------------------------------------------------------------*/
25
26#include "tcuDefs.hpp"
27#include "tcuVector.hpp"
28
29namespace tcu
30{
31
32template <typename T, int Size>
33class Array
34{
35public:
36					Array			(void) {}
37					~Array			(void) {}
38
39	inline T		operator[]		(int ndx) const		{ return m_data[ndx]; }
40	inline T&		operator[]		(int ndx)			{ return m_data[ndx]; }
41
42	inline const T*	getPtr			(void) const		{ return m_data; }
43	inline T*		getPtr			(void)				{ return m_data; }
44
45private:
46	T				m_data[Size];
47};
48
49// Templated matrix class.
50template <typename T, int Rows, int Cols>
51class Matrix
52{
53public:
54	typedef	Vector<T, Rows>			Element;
55	typedef	T						Scalar;
56
57	enum
58	{
59		SIZE = Cols,
60		ROWS = Rows,
61		COLS = Cols,
62	};
63
64									Matrix				(void);
65	explicit						Matrix				(const T& src);
66	explicit						Matrix				(const T src[Rows*Cols]);
67									Matrix				(const Vector<T, Rows>& src);
68									Matrix				(const Matrix<T, Rows, Cols>& src);
69									~Matrix				(void);
70
71	Matrix<T, Rows, Cols>&			operator=			(const Matrix<T, Rows, Cols>& src);
72	Matrix<T, Rows, Cols>&			operator*=			(const Matrix<T, Rows, Cols>& src);
73
74	void							setRow				(int rowNdx, const Vector<T, Cols>& vec);
75	void							setColumn			(int colNdx, const Vector<T, Rows>& vec);
76
77	Vector<T, Cols>					getRow				(int ndx) const;
78	Vector<T, Rows>&				getColumn			(int ndx);
79	const Vector<T, Rows>&			getColumn			(int ndx) const;
80
81	Vector<T, Rows>&				operator[]			(int ndx)					{ return getColumn(ndx);	}
82	const Vector<T, Rows>&			operator[]			(int ndx) const				{ return getColumn(ndx);	}
83
84	inline const T&					operator()			(int row, int col) const	{ return m_data[col][row];	}
85	inline T&						operator()			(int row, int col)			{ return m_data[col][row];	}
86
87	Array<T, Rows*Cols>				getRowMajorData		(void) const;
88	Array<T, Rows*Cols>				getColumnMajorData	(void) const;
89
90private:
91	Vector<Vector<T, Rows>, Cols>	m_data;
92};
93
94// Operators.
95
96// Mat * Mat.
97template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
98Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b);
99
100// Mat * Vec (column vector).
101template <typename T, int Rows, int Cols>
102Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec);
103
104// Vec * Mat (row vector).
105template <typename T, int Rows, int Cols>
106Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx);
107
108// Further operations
109
110template <typename T, int Size>
111struct SquareMatrixOps
112{
113	static T						doDeterminant	(const Matrix<T, Size, Size>& mat);
114	static Matrix<T, Size, Size>	doInverse		(const Matrix<T, Size, Size>& mat);
115};
116
117template <typename T>
118struct SquareMatrixOps<T, 2>
119{
120	static T						doDeterminant	(const Matrix<T, 2, 2>& mat);
121	static Matrix<T, 2, 2>			doInverse		(const Matrix<T, 2, 2>& mat);
122};
123
124template <typename T>
125struct SquareMatrixOps<T, 3>
126{
127	static T						doDeterminant	(const Matrix<T, 3, 3>& mat);
128	static Matrix<T, 3, 3>			doInverse		(const Matrix<T, 3, 3>& mat);
129};
130
131template <typename T>
132struct SquareMatrixOps<T, 4>
133{
134	static T						doDeterminant	(const Matrix<T, 4, 4>& mat);
135	static Matrix<T, 4, 4>			doInverse		(const Matrix<T, 4, 4>& mat);
136};
137
138namespace matrix
139{
140
141template <typename T, int Size>
142T determinant (const Matrix<T, Size, Size>& mat)
143{
144	return SquareMatrixOps<T, Size>::doDeterminant(mat);
145}
146
147template <typename T, int Size>
148Matrix<T, Size, Size> inverse (const Matrix<T, Size, Size>& mat)
149{
150	return SquareMatrixOps<T, Size>::doInverse(mat);
151}
152
153} // matrix
154
155// Template implementations.
156
157template <typename T>
158T SquareMatrixOps<T, 2>::doDeterminant (const Matrix<T, 2, 2>& mat)
159{
160	return mat(0,0) * mat(1,1) - mat(1,0) * mat(0,1);
161}
162
163template <typename T>
164T SquareMatrixOps<T, 3>::doDeterminant (const Matrix<T, 3, 3>& mat)
165{
166	return	+ mat(0,0) * mat(1,1) * mat(2,2)
167			+ mat(0,1) * mat(1,2) * mat(2,0)
168			+ mat(0,2) * mat(1,0) * mat(2,1)
169			- mat(0,0) * mat(1,2) * mat(2,1)
170			- mat(0,1) * mat(1,0) * mat(2,2)
171			- mat(0,2) * mat(1,1) * mat(2,0);
172}
173
174template <typename T>
175T SquareMatrixOps<T, 4>::doDeterminant (const Matrix<T, 4, 4>& mat)
176{
177	using matrix::determinant;
178
179	const T minorMatrices[4][3*3] =
180	{
181		{
182			mat(1,1),	mat(2,1),	mat(3,1),
183			mat(1,2),	mat(2,2),	mat(3,2),
184			mat(1,3),	mat(2,3),	mat(3,3),
185		},
186		{
187			mat(1,0),	mat(2,0),	mat(3,0),
188			mat(1,2),	mat(2,2),	mat(3,2),
189			mat(1,3),	mat(2,3),	mat(3,3),
190		},
191		{
192			mat(1,0),	mat(2,0),	mat(3,0),
193			mat(1,1),	mat(2,1),	mat(3,1),
194			mat(1,3),	mat(2,3),	mat(3,3),
195		},
196		{
197			mat(1,0),	mat(2,0),	mat(3,0),
198			mat(1,1),	mat(2,1),	mat(3,1),
199			mat(1,2),	mat(2,2),	mat(3,2),
200		}
201	};
202
203	return	+ mat(0,0) * determinant(Matrix<T, 3, 3>(minorMatrices[0]))
204			- mat(0,1) * determinant(Matrix<T, 3, 3>(minorMatrices[1]))
205			+ mat(0,2) * determinant(Matrix<T, 3, 3>(minorMatrices[2]))
206			- mat(0,3) * determinant(Matrix<T, 3, 3>(minorMatrices[3]));
207}
208
209template <typename T>
210Matrix<T, 2, 2> SquareMatrixOps<T, 2>::doInverse (const Matrix<T, 2, 2>& mat)
211{
212	using matrix::determinant;
213
214	const T			det		= determinant(mat);
215	Matrix<T, 2, 2>	retVal;
216
217	retVal(0, 0) =  mat(1, 1) / det;
218	retVal(0, 1) = -mat(0, 1) / det;
219	retVal(1, 0) = -mat(1, 0) / det;
220	retVal(1, 1) =  mat(0, 0) / det;
221
222	return retVal;
223}
224
225template <typename T>
226Matrix<T, 3, 3> SquareMatrixOps<T, 3>::doInverse (const Matrix<T, 3, 3>& mat)
227{
228	// Blockwise inversion
229	using matrix::inverse;
230
231	const T areaA[2*2] =
232	{
233		mat(0,0),	mat(0,1),
234		mat(1,0),	mat(1,1)
235	};
236	const T areaB[2] =
237	{
238		mat(0,2),
239		mat(1,2),
240	};
241	const T areaC[2] =
242	{
243		mat(2,0),	mat(2,1),
244	};
245	const T areaD[1] =
246	{
247		mat(2,2)
248	};
249	const T nullField[4] = { T(0.0f) };
250
251	const Matrix<T, 2, 2>	invA = inverse(Matrix<T, 2, 2>(areaA));
252	const Matrix<T, 2, 1>	matB =         Matrix<T, 2, 1>(areaB);
253	const Matrix<T, 1, 2>	matC =         Matrix<T, 1, 2>(areaC);
254	const Matrix<T, 1, 1>	matD =         Matrix<T, 1, 1>(areaD);
255
256	const T					schurComplement = T(1.0f) / (matD - matC*invA*matB)(0,0);
257	const Matrix<T, 2, 2>	zeroMat         = Matrix<T, 2, 2>(nullField);
258
259	const Matrix<T, 2, 2>	blockA = invA + invA*matB*schurComplement*matC*invA;
260	const Matrix<T, 2, 1>	blockB = (zeroMat-invA)*matB*schurComplement;
261	const Matrix<T, 1, 2>	blockC = matC*invA*(-schurComplement);
262	const T					blockD = schurComplement;
263
264	const T result[3*3] =
265	{
266		blockA(0,0),	blockA(0,1),	blockB(0,0),
267		blockA(1,0),	blockA(1,1),	blockB(1,0),
268		blockC(0,0),	blockC(0,1),	blockD,
269	};
270
271	return Matrix<T, 3, 3>(result);
272}
273
274template <typename T>
275Matrix<T, 4, 4> SquareMatrixOps<T, 4>::doInverse (const Matrix<T, 4, 4>& mat)
276{
277	// Blockwise inversion
278	using matrix::inverse;
279
280	const T areaA[2*2] =
281	{
282		mat(0,0),	mat(0,1),
283		mat(1,0),	mat(1,1)
284	};
285	const T areaB[2*2] =
286	{
287		mat(0,2),	mat(0,3),
288		mat(1,2),	mat(1,3)
289	};
290	const T areaC[2*2] =
291	{
292		mat(2,0),	mat(2,1),
293		mat(3,0),	mat(3,1)
294	};
295	const T areaD[2*2] =
296	{
297		mat(2,2),	mat(2,3),
298		mat(3,2),	mat(3,3)
299	};
300	const T nullField[4] = { T(0.0f) };
301
302	const Matrix<T, 2, 2> invA = inverse(Matrix<T, 2, 2>(areaA));
303	const Matrix<T, 2, 2> matB =         Matrix<T, 2, 2>(areaB);
304	const Matrix<T, 2, 2> matC =         Matrix<T, 2, 2>(areaC);
305	const Matrix<T, 2, 2> matD =         Matrix<T, 2, 2>(areaD);
306
307	const Matrix<T, 2, 2> schurComplement = inverse(matD - matC*invA*matB);
308	const Matrix<T, 2, 2> zeroMat         = Matrix<T, 2, 2>(nullField);
309
310	const Matrix<T, 2, 2> blockA = invA + invA*matB*schurComplement*matC*invA;
311	const Matrix<T, 2, 2> blockB = (zeroMat-invA)*matB*schurComplement;
312	const Matrix<T, 2, 2> blockC = (zeroMat-schurComplement)*matC*invA;
313	const Matrix<T, 2, 2> blockD = schurComplement;
314
315	const T result[4*4] =
316	{
317		blockA(0,0),	blockA(0,1),	blockB(0,0),	blockB(0,1),
318		blockA(1,0),	blockA(1,1),	blockB(1,0),	blockB(1,1),
319		blockC(0,0),	blockC(0,1),	blockD(0,0),	blockD(0,1),
320		blockC(1,0),	blockC(1,1),	blockD(1,0),	blockD(1,1),
321	};
322
323	return Matrix<T, 4, 4>(result);
324}
325
326// Initialize to identity.
327template <typename T, int Rows, int Cols>
328Matrix<T, Rows, Cols>::Matrix (void)
329{
330	for (int row = 0; row < Rows; row++)
331		for (int col = 0; col < Cols; col++)
332			(*this)(row, col) = (row == col) ? T(1) : T(0);
333}
334
335// Initialize to diagonal matrix.
336template <typename T, int Rows, int Cols>
337Matrix<T, Rows, Cols>::Matrix (const T& src)
338{
339	for (int row = 0; row < Rows; row++)
340		for (int col = 0; col < Cols; col++)
341			(*this)(row, col) = (row == col) ? src : T(0);
342}
343
344// Initialize from data array.
345template <typename T, int Rows, int Cols>
346Matrix<T, Rows, Cols>::Matrix (const T src[Rows*Cols])
347{
348	for (int row = 0; row < Rows; row++)
349		for (int col = 0; col < Cols; col++)
350			(*this)(row, col) = src[row*Cols + col];
351}
352
353// Initialize to diagonal matrix.
354template <typename T, int Rows, int Cols>
355Matrix<T, Rows, Cols>::Matrix (const Vector<T, Rows>& src)
356{
357	DE_STATIC_ASSERT(Rows == Cols);
358	for (int row = 0; row < Rows; row++)
359		for (int col = 0; col < Cols; col++)
360			(*this)(row, col) = (row == col) ? src.m_data[row] : T(0);
361}
362
363// Copy constructor.
364template <typename T, int Rows, int Cols>
365Matrix<T, Rows, Cols>::Matrix (const Matrix<T, Rows, Cols>& src)
366{
367	*this = src;
368}
369
370// Destructor.
371template <typename T, int Rows, int Cols>
372Matrix<T, Rows, Cols>::~Matrix (void)
373{
374}
375
376// Assignment operator.
377template <typename T, int Rows, int Cols>
378Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator= (const Matrix<T, Rows, Cols>& src)
379{
380	for (int row = 0; row < Rows; row++)
381		for (int col = 0; col < Cols; col++)
382			(*this)(row, col) = src(row, col);
383	return *this;
384}
385
386// Multipy and assign op
387template <typename T, int Rows, int Cols>
388Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator*= (const Matrix<T, Rows, Cols>& src)
389{
390	*this = *this * src;
391	return *this;
392}
393
394template <typename T, int Rows, int Cols>
395void Matrix<T, Rows, Cols>::setRow (int rowNdx, const Vector<T, Cols>& vec)
396{
397	for (int col = 0; col < Cols; col++)
398		(*this)(rowNdx, col) = vec.m_data[col];
399}
400
401template <typename T, int Rows, int Cols>
402void Matrix<T, Rows, Cols>::setColumn (int colNdx, const Vector<T, Rows>& vec)
403{
404	m_data[colNdx] = vec;
405}
406
407template <typename T, int Rows, int Cols>
408Vector<T, Cols> Matrix<T, Rows, Cols>::getRow (int rowNdx) const
409{
410	Vector<T, Cols> res;
411	for (int col = 0; col < Cols; col++)
412		res[col] = (*this)(rowNdx, col);
413	return res;
414}
415
416template <typename T, int Rows, int Cols>
417Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx)
418{
419	return m_data[colNdx];
420}
421
422template <typename T, int Rows, int Cols>
423const Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx) const
424{
425	return m_data[colNdx];
426}
427
428template <typename T, int Rows, int Cols>
429Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getColumnMajorData (void) const
430{
431	Array<T, Rows*Cols> a;
432	T* dst = a.getPtr();
433	for (int col = 0; col < Cols; col++)
434		for (int row = 0; row < Rows; row++)
435			*dst++ = (*this)(row, col);
436	return a;
437}
438
439template <typename T, int Rows, int Cols>
440Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getRowMajorData (void) const
441{
442	Array<T, Rows*Cols> a;
443	T* dst = a.getPtr();
444	for (int row = 0; row < Rows; row++)
445		for (int col = 0; col < Cols; col++)
446			*dst++ = (*this)(row, col);
447	return a;
448}
449
450// Multiplication of two matrices.
451template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
452Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b)
453{
454	DE_STATIC_ASSERT(Cols0 == Rows1);
455	Matrix<T, Rows0, Cols1> res;
456	for (int row = 0; row < Rows0; row++)
457	{
458		for (int col = 0; col < Cols1; col++)
459		{
460			T v = T(0);
461			for (int ndx = 0; ndx < Cols0; ndx++)
462				v += a(row,ndx) * b(ndx,col);
463			res(row,col) = v;
464		}
465	}
466	return res;
467}
468
469// Multiply of matrix with column vector.
470template <typename T, int Rows, int Cols>
471Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec)
472{
473	Vector<T, Rows> res;
474	for (int row = 0; row < Rows; row++)
475	{
476		T v = T(0);
477		for (int col = 0; col < Cols; col++)
478			v += mtx(row,col) * vec.m_data[col];
479		res.m_data[row] = v;
480	}
481	return res;
482}
483
484// Multiply of matrix with row vector.
485template <typename T, int Rows, int Cols>
486Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx)
487{
488	Vector<T, Cols> res;
489	for (int col = 0; col < Cols; col++)
490	{
491		T v = T(0);
492		for (int row = 0; row < Rows; row++)
493			v += mtx(row,col) * vec.m_data[row];
494		res.m_data[col] = v;
495	}
496	return res;
497}
498
499// Common typedefs.
500typedef Matrix<float, 2, 2>		Matrix2f;
501typedef Matrix<float, 3, 3>		Matrix3f;
502typedef Matrix<float, 4, 4>		Matrix4f;
503typedef Matrix<double, 2, 2>	Matrix2d;
504typedef Matrix<double, 3, 3>	Matrix3d;
505typedef Matrix<double, 4, 4>	Matrix4d;
506
507// GLSL-style naming \note CxR.
508typedef Matrix2f			Mat2;
509typedef Matrix<float, 3, 2>	Mat2x3;
510typedef Matrix<float, 4, 2>	Mat2x4;
511typedef Matrix<float, 2, 3>	Mat3x2;
512typedef Matrix3f			Mat3;
513typedef Matrix<float, 4, 3>	Mat3x4;
514typedef Matrix<float, 2, 4>	Mat4x2;
515typedef Matrix<float, 3, 4>	Mat4x3;
516typedef Matrix4f			Mat4;
517
518// Matrix-scalar operators.
519
520template <typename T, int Rows, int Cols>
521Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& mtx, T scalar)
522{
523	Matrix<T, Rows, Cols> res;
524	for (int col = 0; col < Cols; col++)
525		for (int row = 0; row < Rows; row++)
526			res(row, col) = mtx(row, col) + scalar;
527	return res;
528}
529
530template <typename T, int Rows, int Cols>
531Matrix<T, Rows, Cols> operator- (const Matrix<T, Rows, Cols>& mtx, T scalar)
532{
533	Matrix<T, Rows, Cols> res;
534	for (int col = 0; col < Cols; col++)
535		for (int row = 0; row < Rows; row++)
536			res(row, col) = mtx(row, col) - scalar;
537	return res;
538}
539
540template <typename T, int Rows, int Cols>
541Matrix<T, Rows, Cols> operator* (const Matrix<T, Rows, Cols>& mtx, T scalar)
542{
543	Matrix<T, Rows, Cols> res;
544	for (int col = 0; col < Cols; col++)
545		for (int row = 0; row < Rows; row++)
546			res(row, col) = mtx(row, col) * scalar;
547	return res;
548}
549
550template <typename T, int Rows, int Cols>
551Matrix<T, Rows, Cols> operator/ (const Matrix<T, Rows, Cols>& mtx, T scalar)
552{
553	Matrix<T, Rows, Cols> res;
554	for (int col = 0; col < Cols; col++)
555		for (int row = 0; row < Rows; row++)
556			res(row, col) = mtx(row, col) / scalar;
557	return res;
558}
559
560// Matrix-matrix component-wise operators.
561
562template <typename T, int Rows, int Cols>
563Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
564{
565	Matrix<T, Rows, Cols> res;
566	for (int col = 0; col < Cols; col++)
567		for (int row = 0; row < Rows; row++)
568			res(row, col) = a(row, col) + b(row, col);
569	return res;
570}
571
572template <typename T, int Rows, int Cols>
573Matrix<T, Rows, Cols> operator- (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
574{
575	Matrix<T, Rows, Cols> res;
576	for (int col = 0; col < Cols; col++)
577		for (int row = 0; row < Rows; row++)
578			res(row, col) = a(row, col) - b(row, col);
579	return res;
580}
581
582template <typename T, int Rows, int Cols>
583Matrix<T, Rows, Cols> operator/ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
584{
585	Matrix<T, Rows, Cols> res;
586	for (int col = 0; col < Cols; col++)
587		for (int row = 0; row < Rows; row++)
588			res(row, col) = a(row, col) / b(row, col);
589	return res;
590}
591
592} // namespace tcu
593
594#endif // _TCUMATRIX_HPP
595