Home | History | Annotate | Download | only in common
      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 
     65 	Vector<T, Rows>&				operator[]			(int ndx)					{ return getColumn(ndx);	}
     66 	const Vector<T, Rows>&			operator[]			(int ndx) const				{ return getColumn(ndx);	}
     67 
     68 	inline const T&					operator()			(int row, int col) const	{ return m_data[col][row];	}
     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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    362 Matrix<T, Rows, Cols>::~Matrix (void)
    363 {
    364 }
    365 
    366 // Assignment operator.
    367 template <typename T, int Rows, int Cols>
    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>
    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>
    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>
    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>
    398 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>
    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>
    413 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>
    419 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>
    430 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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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>
    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