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 
     29 namespace tcu
     30 {
     31 
     32 template <typename T, int Size>
     33 class Array
     34 {
     35 public:
     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 
     45 private:
     46 	T				m_data[Size];
     47 };
     48 
     49 // Templated matrix class.
     50 template <typename T, int Rows, int Cols>
     51 class Matrix
     52 {
     53 public:
     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 
     90 private:
     91 	Vector<Vector<T, Rows>, Cols>	m_data;
     92 };
     93 
     94 // Operators.
     95 
     96 // Mat * Mat.
     97 template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
     98 Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b);
     99 
    100 // Mat * Vec (column vector).
    101 template <typename T, int Rows, int Cols>
    102 Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec);
    103 
    104 // Vec * Mat (row vector).
    105 template <typename T, int Rows, int Cols>
    106 Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx);
    107 
    108 // Further operations
    109 
    110 template <typename T, int Size>
    111 struct 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 
    117 template <typename T>
    118 struct 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 
    124 template <typename T>
    125 struct 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 
    131 template <typename T>
    132 struct 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 
    138 namespace matrix
    139 {
    140 
    141 template <typename T, int Size>
    142 T determinant (const Matrix<T, Size, Size>& mat)
    143 {
    144 	return SquareMatrixOps<T, Size>::doDeterminant(mat);
    145 }
    146 
    147 template <typename T, int Size>
    148 Matrix<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 
    157 template <typename T>
    158 T 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 
    163 template <typename T>
    164 T 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 
    174 template <typename T>
    175 T 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 
    209 template <typename T>
    210 Matrix<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 
    225 template <typename T>
    226 Matrix<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 
    274 template <typename T>
    275 Matrix<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.
    327 template <typename T, int Rows, int Cols>
    328 Matrix<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.
    336 template <typename T, int Rows, int Cols>
    337 Matrix<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.
    345 template <typename T, int Rows, int Cols>
    346 Matrix<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.
    354 template <typename T, int Rows, int Cols>
    355 Matrix<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.
    364 template <typename T, int Rows, int Cols>
    365 Matrix<T, Rows, Cols>::Matrix (const Matrix<T, Rows, Cols>& src)
    366 {
    367 	*this = src;
    368 }
    369 
    370 // Destructor.
    371 template <typename T, int Rows, int Cols>
    372 Matrix<T, Rows, Cols>::~Matrix (void)
    373 {
    374 }
    375 
    376 // Assignment operator.
    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 	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
    387 template <typename T, int Rows, int Cols>
    388 Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator*= (const Matrix<T, Rows, Cols>& src)
    389 {
    390 	*this = *this * src;
    391 	return *this;
    392 }
    393 
    394 template <typename T, int Rows, int Cols>
    395 void 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 
    401 template <typename T, int Rows, int Cols>
    402 void Matrix<T, Rows, Cols>::setColumn (int colNdx, const Vector<T, Rows>& vec)
    403 {
    404 	m_data[colNdx] = vec;
    405 }
    406 
    407 template <typename T, int Rows, int Cols>
    408 Vector<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 
    416 template <typename T, int Rows, int Cols>
    417 Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx)
    418 {
    419 	return m_data[colNdx];
    420 }
    421 
    422 template <typename T, int Rows, int Cols>
    423 const Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx) const
    424 {
    425 	return m_data[colNdx];
    426 }
    427 
    428 template <typename T, int Rows, int Cols>
    429 Array<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 
    439 template <typename T, int Rows, int Cols>
    440 Array<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.
    451 template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
    452 Matrix<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.
    470 template <typename T, int Rows, int Cols>
    471 Vector<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.
    485 template <typename T, int Rows, int Cols>
    486 Vector<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.
    500 typedef Matrix<float, 2, 2>		Matrix2f;
    501 typedef Matrix<float, 3, 3>		Matrix3f;
    502 typedef Matrix<float, 4, 4>		Matrix4f;
    503 typedef Matrix<double, 2, 2>	Matrix2d;
    504 typedef Matrix<double, 3, 3>	Matrix3d;
    505 typedef Matrix<double, 4, 4>	Matrix4d;
    506 
    507 // GLSL-style naming \note CxR.
    508 typedef Matrix2f			Mat2;
    509 typedef Matrix<float, 3, 2>	Mat2x3;
    510 typedef Matrix<float, 4, 2>	Mat2x4;
    511 typedef Matrix<float, 2, 3>	Mat3x2;
    512 typedef Matrix3f			Mat3;
    513 typedef Matrix<float, 4, 3>	Mat3x4;
    514 typedef Matrix<float, 2, 4>	Mat4x2;
    515 typedef Matrix<float, 3, 4>	Mat4x3;
    516 typedef Matrix4f			Mat4;
    517 
    518 // Matrix-scalar operators.
    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 template <typename T, int Rows, int Cols>
    551 Matrix<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 
    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 Matrix<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