• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 #include "tcuArray.hpp"
29 
30 namespace tcu
31 {
32 
33 // Templated matrix class.
34 template <typename T, int Rows, int Cols>
35 class Matrix
36 {
37 public:
38 	typedef	Vector<T, Rows>			Element;
39 	typedef	T						Scalar;
40 
41 	enum
42 	{
43 		SIZE = Cols,
44 		ROWS = Rows,
45 		COLS = Cols,
46 	};
47 
48 									Matrix				(void);
49 	explicit						Matrix				(const T& src);
50 	explicit						Matrix				(const T src[Rows*Cols]);
51 									Matrix				(const Vector<T, Rows>& src);
52 									Matrix				(const Matrix<T, Rows, Cols>& src);
53 									~Matrix				(void);
54 
55 	Matrix<T, Rows, Cols>&			operator=			(const Matrix<T, Rows, Cols>& src);
56 	Matrix<T, Rows, Cols>&			operator*=			(const Matrix<T, Rows, Cols>& src);
57 
58 	void							setRow				(int rowNdx, const Vector<T, Cols>& vec);
59 	void							setColumn			(int colNdx, const Vector<T, Rows>& vec);
60 
61 	Vector<T, Cols>					getRow				(int ndx) const;
62 	Vector<T, Rows>&				getColumn			(int ndx);
63 	const Vector<T, Rows>&			getColumn			(int ndx) const;
64 
operator [](int ndx)65 	Vector<T, Rows>&				operator[]			(int ndx)					{ return getColumn(ndx);	}
operator [](int ndx) const66 	const Vector<T, Rows>&			operator[]			(int ndx) const				{ return getColumn(ndx);	}
67 
operator ()(int row,int col) const68 	inline const T&					operator()			(int row, int col) const	{ return m_data[col][row];	}
operator ()(int row,int col)69 	inline T&						operator()			(int row, int col)			{ return m_data[col][row];	}
70 
71 	Array<T, Rows*Cols>				getRowMajorData		(void) const;
72 	Array<T, Rows*Cols>				getColumnMajorData	(void) const;
73 
74 private:
75 	Vector<Vector<T, Rows>, Cols>	m_data;
76 } DE_WARN_UNUSED_TYPE;
77 
78 // Operators.
79 
80 // Mat * Mat.
81 template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
82 Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b);
83 
84 // Mat * Vec (column vector).
85 template <typename T, int Rows, int Cols>
86 Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec);
87 
88 // Vec * Mat (row vector).
89 template <typename T, int Rows, int Cols>
90 Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx);
91 
92 template <typename T, int Rows, int Cols>
93 bool operator== (const Matrix<T, Rows, Cols>& lhs, const Matrix<T, Rows, Cols>& rhs);
94 
95 template <typename T, int Rows, int Cols>
96 bool operator!= (const Matrix<T, Rows, Cols>& lhs, const Matrix<T, Rows, Cols>& rhs);
97 
98 // Further operations
99 
100 template <typename T, int Size>
101 struct SquareMatrixOps
102 {
103 	static T						doDeterminant	(const Matrix<T, Size, Size>& mat);
104 	static Matrix<T, Size, Size>	doInverse		(const Matrix<T, Size, Size>& mat);
105 };
106 
107 template <typename T>
108 struct SquareMatrixOps<T, 2>
109 {
110 	static T						doDeterminant	(const Matrix<T, 2, 2>& mat);
111 	static Matrix<T, 2, 2>			doInverse		(const Matrix<T, 2, 2>& mat);
112 };
113 
114 template <typename T>
115 struct SquareMatrixOps<T, 3>
116 {
117 	static T						doDeterminant	(const Matrix<T, 3, 3>& mat);
118 	static Matrix<T, 3, 3>			doInverse		(const Matrix<T, 3, 3>& mat);
119 };
120 
121 template <typename T>
122 struct SquareMatrixOps<T, 4>
123 {
124 	static T						doDeterminant	(const Matrix<T, 4, 4>& mat);
125 	static Matrix<T, 4, 4>			doInverse		(const Matrix<T, 4, 4>& mat);
126 };
127 
128 namespace matrix
129 {
130 
131 template <typename T, int Size>
determinant(const Matrix<T,Size,Size> & mat)132 T determinant (const Matrix<T, Size, Size>& mat)
133 {
134 	return SquareMatrixOps<T, Size>::doDeterminant(mat);
135 }
136 
137 template <typename T, int Size>
inverse(const Matrix<T,Size,Size> & mat)138 Matrix<T, Size, Size> inverse (const Matrix<T, Size, Size>& mat)
139 {
140 	return SquareMatrixOps<T, Size>::doInverse(mat);
141 }
142 
143 } // matrix
144 
145 // Template implementations.
146 
147 template <typename T>
doDeterminant(const Matrix<T,2,2> & mat)148 T SquareMatrixOps<T, 2>::doDeterminant (const Matrix<T, 2, 2>& mat)
149 {
150 	return mat(0,0) * mat(1,1) - mat(1,0) * mat(0,1);
151 }
152 
153 template <typename T>
doDeterminant(const Matrix<T,3,3> & mat)154 T SquareMatrixOps<T, 3>::doDeterminant (const Matrix<T, 3, 3>& mat)
155 {
156 	return	+ mat(0,0) * mat(1,1) * mat(2,2)
157 			+ mat(0,1) * mat(1,2) * mat(2,0)
158 			+ mat(0,2) * mat(1,0) * mat(2,1)
159 			- mat(0,0) * mat(1,2) * mat(2,1)
160 			- mat(0,1) * mat(1,0) * mat(2,2)
161 			- mat(0,2) * mat(1,1) * mat(2,0);
162 }
163 
164 template <typename T>
doDeterminant(const Matrix<T,4,4> & mat)165 T SquareMatrixOps<T, 4>::doDeterminant (const Matrix<T, 4, 4>& mat)
166 {
167 	using matrix::determinant;
168 
169 	const T minorMatrices[4][3*3] =
170 	{
171 		{
172 			mat(1,1),	mat(2,1),	mat(3,1),
173 			mat(1,2),	mat(2,2),	mat(3,2),
174 			mat(1,3),	mat(2,3),	mat(3,3),
175 		},
176 		{
177 			mat(1,0),	mat(2,0),	mat(3,0),
178 			mat(1,2),	mat(2,2),	mat(3,2),
179 			mat(1,3),	mat(2,3),	mat(3,3),
180 		},
181 		{
182 			mat(1,0),	mat(2,0),	mat(3,0),
183 			mat(1,1),	mat(2,1),	mat(3,1),
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,1),	mat(2,1),	mat(3,1),
189 			mat(1,2),	mat(2,2),	mat(3,2),
190 		}
191 	};
192 
193 	return	+ mat(0,0) * determinant(Matrix<T, 3, 3>(minorMatrices[0]))
194 			- mat(0,1) * determinant(Matrix<T, 3, 3>(minorMatrices[1]))
195 			+ mat(0,2) * determinant(Matrix<T, 3, 3>(minorMatrices[2]))
196 			- mat(0,3) * determinant(Matrix<T, 3, 3>(minorMatrices[3]));
197 }
198 
199 template <typename T>
doInverse(const Matrix<T,2,2> & mat)200 Matrix<T, 2, 2> SquareMatrixOps<T, 2>::doInverse (const Matrix<T, 2, 2>& mat)
201 {
202 	using matrix::determinant;
203 
204 	const T			det		= determinant(mat);
205 	Matrix<T, 2, 2>	retVal;
206 
207 	retVal(0, 0) =  mat(1, 1) / det;
208 	retVal(0, 1) = -mat(0, 1) / det;
209 	retVal(1, 0) = -mat(1, 0) / det;
210 	retVal(1, 1) =  mat(0, 0) / det;
211 
212 	return retVal;
213 }
214 
215 template <typename T>
doInverse(const Matrix<T,3,3> & mat)216 Matrix<T, 3, 3> SquareMatrixOps<T, 3>::doInverse (const Matrix<T, 3, 3>& mat)
217 {
218 	// Blockwise inversion
219 	using matrix::inverse;
220 
221 	const T areaA[2*2] =
222 	{
223 		mat(0,0),	mat(0,1),
224 		mat(1,0),	mat(1,1)
225 	};
226 	const T areaB[2] =
227 	{
228 		mat(0,2),
229 		mat(1,2),
230 	};
231 	const T areaC[2] =
232 	{
233 		mat(2,0),	mat(2,1),
234 	};
235 	const T areaD[1] =
236 	{
237 		mat(2,2)
238 	};
239 	const T nullField[4] = { T(0.0f) };
240 
241 	const Matrix<T, 2, 2>	invA = inverse(Matrix<T, 2, 2>(areaA));
242 	const Matrix<T, 2, 1>	matB =         Matrix<T, 2, 1>(areaB);
243 	const Matrix<T, 1, 2>	matC =         Matrix<T, 1, 2>(areaC);
244 	const Matrix<T, 1, 1>	matD =         Matrix<T, 1, 1>(areaD);
245 
246 	const T					schurComplement = T(1.0f) / (matD - matC*invA*matB)(0,0);
247 	const Matrix<T, 2, 2>	zeroMat         = Matrix<T, 2, 2>(nullField);
248 
249 	const Matrix<T, 2, 2>	blockA = invA + invA*matB*schurComplement*matC*invA;
250 	const Matrix<T, 2, 1>	blockB = (zeroMat-invA)*matB*schurComplement;
251 	const Matrix<T, 1, 2>	blockC = matC*invA*(-schurComplement);
252 	const T					blockD = schurComplement;
253 
254 	const T result[3*3] =
255 	{
256 		blockA(0,0),	blockA(0,1),	blockB(0,0),
257 		blockA(1,0),	blockA(1,1),	blockB(1,0),
258 		blockC(0,0),	blockC(0,1),	blockD,
259 	};
260 
261 	return Matrix<T, 3, 3>(result);
262 }
263 
264 template <typename T>
doInverse(const Matrix<T,4,4> & mat)265 Matrix<T, 4, 4> SquareMatrixOps<T, 4>::doInverse (const Matrix<T, 4, 4>& mat)
266 {
267 	// Blockwise inversion
268 	using matrix::inverse;
269 
270 	const T areaA[2*2] =
271 	{
272 		mat(0,0),	mat(0,1),
273 		mat(1,0),	mat(1,1)
274 	};
275 	const T areaB[2*2] =
276 	{
277 		mat(0,2),	mat(0,3),
278 		mat(1,2),	mat(1,3)
279 	};
280 	const T areaC[2*2] =
281 	{
282 		mat(2,0),	mat(2,1),
283 		mat(3,0),	mat(3,1)
284 	};
285 	const T areaD[2*2] =
286 	{
287 		mat(2,2),	mat(2,3),
288 		mat(3,2),	mat(3,3)
289 	};
290 	const T nullField[4] = { T(0.0f) };
291 
292 	const Matrix<T, 2, 2> invA = inverse(Matrix<T, 2, 2>(areaA));
293 	const Matrix<T, 2, 2> matB =         Matrix<T, 2, 2>(areaB);
294 	const Matrix<T, 2, 2> matC =         Matrix<T, 2, 2>(areaC);
295 	const Matrix<T, 2, 2> matD =         Matrix<T, 2, 2>(areaD);
296 
297 	const Matrix<T, 2, 2> schurComplement = inverse(matD - matC*invA*matB);
298 	const Matrix<T, 2, 2> zeroMat         = Matrix<T, 2, 2>(nullField);
299 
300 	const Matrix<T, 2, 2> blockA = invA + invA*matB*schurComplement*matC*invA;
301 	const Matrix<T, 2, 2> blockB = (zeroMat-invA)*matB*schurComplement;
302 	const Matrix<T, 2, 2> blockC = (zeroMat-schurComplement)*matC*invA;
303 	const Matrix<T, 2, 2> blockD = schurComplement;
304 
305 	const T result[4*4] =
306 	{
307 		blockA(0,0),	blockA(0,1),	blockB(0,0),	blockB(0,1),
308 		blockA(1,0),	blockA(1,1),	blockB(1,0),	blockB(1,1),
309 		blockC(0,0),	blockC(0,1),	blockD(0,0),	blockD(0,1),
310 		blockC(1,0),	blockC(1,1),	blockD(1,0),	blockD(1,1),
311 	};
312 
313 	return Matrix<T, 4, 4>(result);
314 }
315 
316 // Initialize to identity.
317 template <typename T, int Rows, int Cols>
Matrix(void)318 Matrix<T, Rows, Cols>::Matrix (void)
319 {
320 	for (int row = 0; row < Rows; row++)
321 		for (int col = 0; col < Cols; col++)
322 			(*this)(row, col) = (row == col) ? T(1) : T(0);
323 }
324 
325 // Initialize to diagonal matrix.
326 template <typename T, int Rows, int Cols>
Matrix(const T & src)327 Matrix<T, Rows, Cols>::Matrix (const T& src)
328 {
329 	for (int row = 0; row < Rows; row++)
330 		for (int col = 0; col < Cols; col++)
331 			(*this)(row, col) = (row == col) ? src : T(0);
332 }
333 
334 // Initialize from data array.
335 template <typename T, int Rows, int Cols>
Matrix(const T src[Rows * Cols])336 Matrix<T, Rows, Cols>::Matrix (const T src[Rows*Cols])
337 {
338 	for (int row = 0; row < Rows; row++)
339 		for (int col = 0; col < Cols; col++)
340 			(*this)(row, col) = src[row*Cols + col];
341 }
342 
343 // Initialize to diagonal matrix.
344 template <typename T, int Rows, int Cols>
Matrix(const Vector<T,Rows> & src)345 Matrix<T, Rows, Cols>::Matrix (const Vector<T, Rows>& src)
346 {
347 	DE_STATIC_ASSERT(Rows == Cols);
348 	for (int row = 0; row < Rows; row++)
349 		for (int col = 0; col < Cols; col++)
350 			(*this)(row, col) = (row == col) ? src.m_data[row] : T(0);
351 }
352 
353 // Copy constructor.
354 template <typename T, int Rows, int Cols>
Matrix(const Matrix<T,Rows,Cols> & src)355 Matrix<T, Rows, Cols>::Matrix (const Matrix<T, Rows, Cols>& src)
356 {
357 	*this = src;
358 }
359 
360 // Destructor.
361 template <typename T, int Rows, int Cols>
~Matrix(void)362 Matrix<T, Rows, Cols>::~Matrix (void)
363 {
364 }
365 
366 // Assignment operator.
367 template <typename T, int Rows, int Cols>
operator =(const Matrix<T,Rows,Cols> & src)368 Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator= (const Matrix<T, Rows, Cols>& src)
369 {
370 	for (int row = 0; row < Rows; row++)
371 		for (int col = 0; col < Cols; col++)
372 			(*this)(row, col) = src(row, col);
373 	return *this;
374 }
375 
376 // Multipy and assign op
377 template <typename T, int Rows, int Cols>
operator *=(const Matrix<T,Rows,Cols> & src)378 Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator*= (const Matrix<T, Rows, Cols>& src)
379 {
380 	*this = *this * src;
381 	return *this;
382 }
383 
384 template <typename T, int Rows, int Cols>
setRow(int rowNdx,const Vector<T,Cols> & vec)385 void Matrix<T, Rows, Cols>::setRow (int rowNdx, const Vector<T, Cols>& vec)
386 {
387 	for (int col = 0; col < Cols; col++)
388 		(*this)(rowNdx, col) = vec.m_data[col];
389 }
390 
391 template <typename T, int Rows, int Cols>
setColumn(int colNdx,const Vector<T,Rows> & vec)392 void Matrix<T, Rows, Cols>::setColumn (int colNdx, const Vector<T, Rows>& vec)
393 {
394 	m_data[colNdx] = vec;
395 }
396 
397 template <typename T, int Rows, int Cols>
getRow(int rowNdx) const398 Vector<T, Cols> Matrix<T, Rows, Cols>::getRow (int rowNdx) const
399 {
400 	Vector<T, Cols> res;
401 	for (int col = 0; col < Cols; col++)
402 		res[col] = (*this)(rowNdx, col);
403 	return res;
404 }
405 
406 template <typename T, int Rows, int Cols>
getColumn(int colNdx)407 Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx)
408 {
409 	return m_data[colNdx];
410 }
411 
412 template <typename T, int Rows, int Cols>
getColumn(int colNdx) const413 const Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx) const
414 {
415 	return m_data[colNdx];
416 }
417 
418 template <typename T, int Rows, int Cols>
getColumnMajorData(void) const419 Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getColumnMajorData (void) const
420 {
421 	Array<T, Rows*Cols> a;
422 	T* dst = a.getPtr();
423 	for (int col = 0; col < Cols; col++)
424 		for (int row = 0; row < Rows; row++)
425 			*dst++ = (*this)(row, col);
426 	return a;
427 }
428 
429 template <typename T, int Rows, int Cols>
getRowMajorData(void) const430 Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getRowMajorData (void) const
431 {
432 	Array<T, Rows*Cols> a;
433 	T* dst = a.getPtr();
434 	for (int row = 0; row < Rows; row++)
435 		for (int col = 0; col < Cols; col++)
436 			*dst++ = (*this)(row, col);
437 	return a;
438 }
439 
440 // Multiplication of two matrices.
441 template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
operator *(const Matrix<T,Rows0,Cols0> & a,const Matrix<T,Rows1,Cols1> & b)442 Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b)
443 {
444 	DE_STATIC_ASSERT(Cols0 == Rows1);
445 	Matrix<T, Rows0, Cols1> res;
446 	for (int row = 0; row < Rows0; row++)
447 	{
448 		for (int col = 0; col < Cols1; col++)
449 		{
450 			T v = T(0);
451 			for (int ndx = 0; ndx < Cols0; ndx++)
452 				v += a(row,ndx) * b(ndx,col);
453 			res(row,col) = v;
454 		}
455 	}
456 	return res;
457 }
458 
459 // Multiply of matrix with column vector.
460 template <typename T, int Rows, int Cols>
operator *(const Matrix<T,Rows,Cols> & mtx,const Vector<T,Cols> & vec)461 Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec)
462 {
463 	Vector<T, Rows> res;
464 	for (int row = 0; row < Rows; row++)
465 	{
466 		T v = T(0);
467 		for (int col = 0; col < Cols; col++)
468 			v += mtx(row,col) * vec.m_data[col];
469 		res.m_data[row] = v;
470 	}
471 	return res;
472 }
473 
474 // Multiply of matrix with row vector.
475 template <typename T, int Rows, int Cols>
operator *(const Vector<T,Rows> & vec,const Matrix<T,Rows,Cols> & mtx)476 Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx)
477 {
478 	Vector<T, Cols> res;
479 	for (int col = 0; col < Cols; col++)
480 	{
481 		T v = T(0);
482 		for (int row = 0; row < Rows; row++)
483 			v += mtx(row,col) * vec.m_data[row];
484 		res.m_data[col] = v;
485 	}
486 	return res;
487 }
488 
489 // Common typedefs.
490 typedef Matrix<float, 2, 2>		Matrix2f;
491 typedef Matrix<float, 3, 3>		Matrix3f;
492 typedef Matrix<float, 4, 4>		Matrix4f;
493 typedef Matrix<double, 2, 2>	Matrix2d;
494 typedef Matrix<double, 3, 3>	Matrix3d;
495 typedef Matrix<double, 4, 4>	Matrix4d;
496 
497 // GLSL-style naming \note CxR.
498 typedef Matrix2f			Mat2;
499 typedef Matrix<float, 3, 2>	Mat2x3;
500 typedef Matrix<float, 4, 2>	Mat2x4;
501 typedef Matrix<float, 2, 3>	Mat3x2;
502 typedef Matrix3f			Mat3;
503 typedef Matrix<float, 4, 3>	Mat3x4;
504 typedef Matrix<float, 2, 4>	Mat4x2;
505 typedef Matrix<float, 3, 4>	Mat4x3;
506 typedef Matrix4f			Mat4;
507 
508 // Matrix-scalar operators.
509 
510 template <typename T, int Rows, int Cols>
operator +(const Matrix<T,Rows,Cols> & mtx,T scalar)511 Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& mtx, T scalar)
512 {
513 	Matrix<T, Rows, Cols> res;
514 	for (int col = 0; col < Cols; col++)
515 		for (int row = 0; row < Rows; row++)
516 			res(row, col) = mtx(row, col) + scalar;
517 	return res;
518 }
519 
520 template <typename T, int Rows, int Cols>
operator -(const Matrix<T,Rows,Cols> & mtx,T scalar)521 Matrix<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 
530 template <typename T, int Rows, int Cols>
operator *(const Matrix<T,Rows,Cols> & mtx,T scalar)531 Matrix<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 
540 template <typename T, int Rows, int Cols>
operator /(const Matrix<T,Rows,Cols> & mtx,T scalar)541 Matrix<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 
550 // Matrix-matrix component-wise operators.
551 
552 template <typename T, int Rows, int Cols>
operator +(const Matrix<T,Rows,Cols> & a,const Matrix<T,Rows,Cols> & b)553 Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
554 {
555 	Matrix<T, Rows, Cols> res;
556 	for (int col = 0; col < Cols; col++)
557 		for (int row = 0; row < Rows; row++)
558 			res(row, col) = a(row, col) + b(row, col);
559 	return res;
560 }
561 
562 template <typename T, int Rows, int Cols>
operator -(const Matrix<T,Rows,Cols> & a,const Matrix<T,Rows,Cols> & b)563 Matrix<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 
572 template <typename T, int Rows, int Cols>
operator /(const Matrix<T,Rows,Cols> & a,const Matrix<T,Rows,Cols> & b)573 Matrix<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 
582 template <typename T, int Rows, int Cols>
operator ==(const Matrix<T,Rows,Cols> & lhs,const Matrix<T,Rows,Cols> & rhs)583 bool operator== (const Matrix<T, Rows, Cols>& lhs, const Matrix<T, Rows, Cols>& rhs)
584 {
585 	for (int row = 0; row < Rows; row++)
586 		for (int col = 0; col < Cols; col++)
587 			if (lhs(row, col) != rhs(row, col))
588 				return false;
589 	return true;
590 }
591 
592 template <typename T, int Rows, int Cols>
operator !=(const Matrix<T,Rows,Cols> & lhs,const Matrix<T,Rows,Cols> & rhs)593 bool operator!= (const Matrix<T, Rows, Cols>& lhs, const Matrix<T, Rows, Cols>& rhs)
594 {
595 	return !(lhs == rhs);
596 }
597 
598 } // tcu
599 
600 #endif // _TCUMATRIX_HPP
601