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 // Further operations
     93 
     94 template <typename T, int Size>
     95 struct SquareMatrixOps
     96 {
     97 	static T						doDeterminant	(const Matrix<T, Size, Size>& mat);
     98 	static Matrix<T, Size, Size>	doInverse		(const Matrix<T, Size, Size>& mat);
     99 };
    100 
    101 template <typename T>
    102 struct SquareMatrixOps<T, 2>
    103 {
    104 	static T						doDeterminant	(const Matrix<T, 2, 2>& mat);
    105 	static Matrix<T, 2, 2>			doInverse		(const Matrix<T, 2, 2>& mat);
    106 };
    107 
    108 template <typename T>
    109 struct SquareMatrixOps<T, 3>
    110 {
    111 	static T						doDeterminant	(const Matrix<T, 3, 3>& mat);
    112 	static Matrix<T, 3, 3>			doInverse		(const Matrix<T, 3, 3>& mat);
    113 };
    114 
    115 template <typename T>
    116 struct SquareMatrixOps<T, 4>
    117 {
    118 	static T						doDeterminant	(const Matrix<T, 4, 4>& mat);
    119 	static Matrix<T, 4, 4>			doInverse		(const Matrix<T, 4, 4>& mat);
    120 };
    121 
    122 namespace matrix
    123 {
    124 
    125 template <typename T, int Size>
    126 T determinant (const Matrix<T, Size, Size>& mat)
    127 {
    128 	return SquareMatrixOps<T, Size>::doDeterminant(mat);
    129 }
    130 
    131 template <typename T, int Size>
    132 Matrix<T, Size, Size> inverse (const Matrix<T, Size, Size>& mat)
    133 {
    134 	return SquareMatrixOps<T, Size>::doInverse(mat);
    135 }
    136 
    137 } // matrix
    138 
    139 // Template implementations.
    140 
    141 template <typename T>
    142 T SquareMatrixOps<T, 2>::doDeterminant (const Matrix<T, 2, 2>& mat)
    143 {
    144 	return mat(0,0) * mat(1,1) - mat(1,0) * mat(0,1);
    145 }
    146 
    147 template <typename T>
    148 T SquareMatrixOps<T, 3>::doDeterminant (const Matrix<T, 3, 3>& mat)
    149 {
    150 	return	+ mat(0,0) * mat(1,1) * mat(2,2)
    151 			+ mat(0,1) * mat(1,2) * mat(2,0)
    152 			+ mat(0,2) * mat(1,0) * mat(2,1)
    153 			- mat(0,0) * mat(1,2) * mat(2,1)
    154 			- mat(0,1) * mat(1,0) * mat(2,2)
    155 			- mat(0,2) * mat(1,1) * mat(2,0);
    156 }
    157 
    158 template <typename T>
    159 T SquareMatrixOps<T, 4>::doDeterminant (const Matrix<T, 4, 4>& mat)
    160 {
    161 	using matrix::determinant;
    162 
    163 	const T minorMatrices[4][3*3] =
    164 	{
    165 		{
    166 			mat(1,1),	mat(2,1),	mat(3,1),
    167 			mat(1,2),	mat(2,2),	mat(3,2),
    168 			mat(1,3),	mat(2,3),	mat(3,3),
    169 		},
    170 		{
    171 			mat(1,0),	mat(2,0),	mat(3,0),
    172 			mat(1,2),	mat(2,2),	mat(3,2),
    173 			mat(1,3),	mat(2,3),	mat(3,3),
    174 		},
    175 		{
    176 			mat(1,0),	mat(2,0),	mat(3,0),
    177 			mat(1,1),	mat(2,1),	mat(3,1),
    178 			mat(1,3),	mat(2,3),	mat(3,3),
    179 		},
    180 		{
    181 			mat(1,0),	mat(2,0),	mat(3,0),
    182 			mat(1,1),	mat(2,1),	mat(3,1),
    183 			mat(1,2),	mat(2,2),	mat(3,2),
    184 		}
    185 	};
    186 
    187 	return	+ mat(0,0) * determinant(Matrix<T, 3, 3>(minorMatrices[0]))
    188 			- mat(0,1) * determinant(Matrix<T, 3, 3>(minorMatrices[1]))
    189 			+ mat(0,2) * determinant(Matrix<T, 3, 3>(minorMatrices[2]))
    190 			- mat(0,3) * determinant(Matrix<T, 3, 3>(minorMatrices[3]));
    191 }
    192 
    193 template <typename T>
    194 Matrix<T, 2, 2> SquareMatrixOps<T, 2>::doInverse (const Matrix<T, 2, 2>& mat)
    195 {
    196 	using matrix::determinant;
    197 
    198 	const T			det		= determinant(mat);
    199 	Matrix<T, 2, 2>	retVal;
    200 
    201 	retVal(0, 0) =  mat(1, 1) / det;
    202 	retVal(0, 1) = -mat(0, 1) / det;
    203 	retVal(1, 0) = -mat(1, 0) / det;
    204 	retVal(1, 1) =  mat(0, 0) / det;
    205 
    206 	return retVal;
    207 }
    208 
    209 template <typename T>
    210 Matrix<T, 3, 3> SquareMatrixOps<T, 3>::doInverse (const Matrix<T, 3, 3>& mat)
    211 {
    212 	// Blockwise inversion
    213 	using matrix::inverse;
    214 
    215 	const T areaA[2*2] =
    216 	{
    217 		mat(0,0),	mat(0,1),
    218 		mat(1,0),	mat(1,1)
    219 	};
    220 	const T areaB[2] =
    221 	{
    222 		mat(0,2),
    223 		mat(1,2),
    224 	};
    225 	const T areaC[2] =
    226 	{
    227 		mat(2,0),	mat(2,1),
    228 	};
    229 	const T areaD[1] =
    230 	{
    231 		mat(2,2)
    232 	};
    233 	const T nullField[4] = { T(0.0f) };
    234 
    235 	const Matrix<T, 2, 2>	invA = inverse(Matrix<T, 2, 2>(areaA));
    236 	const Matrix<T, 2, 1>	matB =         Matrix<T, 2, 1>(areaB);
    237 	const Matrix<T, 1, 2>	matC =         Matrix<T, 1, 2>(areaC);
    238 	const Matrix<T, 1, 1>	matD =         Matrix<T, 1, 1>(areaD);
    239 
    240 	const T					schurComplement = T(1.0f) / (matD - matC*invA*matB)(0,0);
    241 	const Matrix<T, 2, 2>	zeroMat         = Matrix<T, 2, 2>(nullField);
    242 
    243 	const Matrix<T, 2, 2>	blockA = invA + invA*matB*schurComplement*matC*invA;
    244 	const Matrix<T, 2, 1>	blockB = (zeroMat-invA)*matB*schurComplement;
    245 	const Matrix<T, 1, 2>	blockC = matC*invA*(-schurComplement);
    246 	const T					blockD = schurComplement;
    247 
    248 	const T result[3*3] =
    249 	{
    250 		blockA(0,0),	blockA(0,1),	blockB(0,0),
    251 		blockA(1,0),	blockA(1,1),	blockB(1,0),
    252 		blockC(0,0),	blockC(0,1),	blockD,
    253 	};
    254 
    255 	return Matrix<T, 3, 3>(result);
    256 }
    257 
    258 template <typename T>
    259 Matrix<T, 4, 4> SquareMatrixOps<T, 4>::doInverse (const Matrix<T, 4, 4>& mat)
    260 {
    261 	// Blockwise inversion
    262 	using matrix::inverse;
    263 
    264 	const T areaA[2*2] =
    265 	{
    266 		mat(0,0),	mat(0,1),
    267 		mat(1,0),	mat(1,1)
    268 	};
    269 	const T areaB[2*2] =
    270 	{
    271 		mat(0,2),	mat(0,3),
    272 		mat(1,2),	mat(1,3)
    273 	};
    274 	const T areaC[2*2] =
    275 	{
    276 		mat(2,0),	mat(2,1),
    277 		mat(3,0),	mat(3,1)
    278 	};
    279 	const T areaD[2*2] =
    280 	{
    281 		mat(2,2),	mat(2,3),
    282 		mat(3,2),	mat(3,3)
    283 	};
    284 	const T nullField[4] = { T(0.0f) };
    285 
    286 	const Matrix<T, 2, 2> invA = inverse(Matrix<T, 2, 2>(areaA));
    287 	const Matrix<T, 2, 2> matB =         Matrix<T, 2, 2>(areaB);
    288 	const Matrix<T, 2, 2> matC =         Matrix<T, 2, 2>(areaC);
    289 	const Matrix<T, 2, 2> matD =         Matrix<T, 2, 2>(areaD);
    290 
    291 	const Matrix<T, 2, 2> schurComplement = inverse(matD - matC*invA*matB);
    292 	const Matrix<T, 2, 2> zeroMat         = Matrix<T, 2, 2>(nullField);
    293 
    294 	const Matrix<T, 2, 2> blockA = invA + invA*matB*schurComplement*matC*invA;
    295 	const Matrix<T, 2, 2> blockB = (zeroMat-invA)*matB*schurComplement;
    296 	const Matrix<T, 2, 2> blockC = (zeroMat-schurComplement)*matC*invA;
    297 	const Matrix<T, 2, 2> blockD = schurComplement;
    298 
    299 	const T result[4*4] =
    300 	{
    301 		blockA(0,0),	blockA(0,1),	blockB(0,0),	blockB(0,1),
    302 		blockA(1,0),	blockA(1,1),	blockB(1,0),	blockB(1,1),
    303 		blockC(0,0),	blockC(0,1),	blockD(0,0),	blockD(0,1),
    304 		blockC(1,0),	blockC(1,1),	blockD(1,0),	blockD(1,1),
    305 	};
    306 
    307 	return Matrix<T, 4, 4>(result);
    308 }
    309 
    310 // Initialize to identity.
    311 template <typename T, int Rows, int Cols>
    312 Matrix<T, Rows, Cols>::Matrix (void)
    313 {
    314 	for (int row = 0; row < Rows; row++)
    315 		for (int col = 0; col < Cols; col++)
    316 			(*this)(row, col) = (row == col) ? T(1) : T(0);
    317 }
    318 
    319 // Initialize to diagonal matrix.
    320 template <typename T, int Rows, int Cols>
    321 Matrix<T, Rows, Cols>::Matrix (const T& src)
    322 {
    323 	for (int row = 0; row < Rows; row++)
    324 		for (int col = 0; col < Cols; col++)
    325 			(*this)(row, col) = (row == col) ? src : T(0);
    326 }
    327 
    328 // Initialize from data array.
    329 template <typename T, int Rows, int Cols>
    330 Matrix<T, Rows, Cols>::Matrix (const T src[Rows*Cols])
    331 {
    332 	for (int row = 0; row < Rows; row++)
    333 		for (int col = 0; col < Cols; col++)
    334 			(*this)(row, col) = src[row*Cols + col];
    335 }
    336 
    337 // Initialize to diagonal matrix.
    338 template <typename T, int Rows, int Cols>
    339 Matrix<T, Rows, Cols>::Matrix (const Vector<T, Rows>& src)
    340 {
    341 	DE_STATIC_ASSERT(Rows == Cols);
    342 	for (int row = 0; row < Rows; row++)
    343 		for (int col = 0; col < Cols; col++)
    344 			(*this)(row, col) = (row == col) ? src.m_data[row] : T(0);
    345 }
    346 
    347 // Copy constructor.
    348 template <typename T, int Rows, int Cols>
    349 Matrix<T, Rows, Cols>::Matrix (const Matrix<T, Rows, Cols>& src)
    350 {
    351 	*this = src;
    352 }
    353 
    354 // Destructor.
    355 template <typename T, int Rows, int Cols>
    356 Matrix<T, Rows, Cols>::~Matrix (void)
    357 {
    358 }
    359 
    360 // Assignment operator.
    361 template <typename T, int Rows, int Cols>
    362 Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator= (const Matrix<T, Rows, Cols>& src)
    363 {
    364 	for (int row = 0; row < Rows; row++)
    365 		for (int col = 0; col < Cols; col++)
    366 			(*this)(row, col) = src(row, col);
    367 	return *this;
    368 }
    369 
    370 // Multipy and assign op
    371 template <typename T, int Rows, int Cols>
    372 Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator*= (const Matrix<T, Rows, Cols>& src)
    373 {
    374 	*this = *this * src;
    375 	return *this;
    376 }
    377 
    378 template <typename T, int Rows, int Cols>
    379 void Matrix<T, Rows, Cols>::setRow (int rowNdx, const Vector<T, Cols>& vec)
    380 {
    381 	for (int col = 0; col < Cols; col++)
    382 		(*this)(rowNdx, col) = vec.m_data[col];
    383 }
    384 
    385 template <typename T, int Rows, int Cols>
    386 void Matrix<T, Rows, Cols>::setColumn (int colNdx, const Vector<T, Rows>& vec)
    387 {
    388 	m_data[colNdx] = vec;
    389 }
    390 
    391 template <typename T, int Rows, int Cols>
    392 Vector<T, Cols> Matrix<T, Rows, Cols>::getRow (int rowNdx) const
    393 {
    394 	Vector<T, Cols> res;
    395 	for (int col = 0; col < Cols; col++)
    396 		res[col] = (*this)(rowNdx, col);
    397 	return res;
    398 }
    399 
    400 template <typename T, int Rows, int Cols>
    401 Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx)
    402 {
    403 	return m_data[colNdx];
    404 }
    405 
    406 template <typename T, int Rows, int Cols>
    407 const Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx) const
    408 {
    409 	return m_data[colNdx];
    410 }
    411 
    412 template <typename T, int Rows, int Cols>
    413 Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getColumnMajorData (void) const
    414 {
    415 	Array<T, Rows*Cols> a;
    416 	T* dst = a.getPtr();
    417 	for (int col = 0; col < Cols; col++)
    418 		for (int row = 0; row < Rows; row++)
    419 			*dst++ = (*this)(row, col);
    420 	return a;
    421 }
    422 
    423 template <typename T, int Rows, int Cols>
    424 Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getRowMajorData (void) const
    425 {
    426 	Array<T, Rows*Cols> a;
    427 	T* dst = a.getPtr();
    428 	for (int row = 0; row < Rows; row++)
    429 		for (int col = 0; col < Cols; col++)
    430 			*dst++ = (*this)(row, col);
    431 	return a;
    432 }
    433 
    434 // Multiplication of two matrices.
    435 template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
    436 Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b)
    437 {
    438 	DE_STATIC_ASSERT(Cols0 == Rows1);
    439 	Matrix<T, Rows0, Cols1> res;
    440 	for (int row = 0; row < Rows0; row++)
    441 	{
    442 		for (int col = 0; col < Cols1; col++)
    443 		{
    444 			T v = T(0);
    445 			for (int ndx = 0; ndx < Cols0; ndx++)
    446 				v += a(row,ndx) * b(ndx,col);
    447 			res(row,col) = v;
    448 		}
    449 	}
    450 	return res;
    451 }
    452 
    453 // Multiply of matrix with column vector.
    454 template <typename T, int Rows, int Cols>
    455 Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec)
    456 {
    457 	Vector<T, Rows> res;
    458 	for (int row = 0; row < Rows; row++)
    459 	{
    460 		T v = T(0);
    461 		for (int col = 0; col < Cols; col++)
    462 			v += mtx(row,col) * vec.m_data[col];
    463 		res.m_data[row] = v;
    464 	}
    465 	return res;
    466 }
    467 
    468 // Multiply of matrix with row vector.
    469 template <typename T, int Rows, int Cols>
    470 Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx)
    471 {
    472 	Vector<T, Cols> res;
    473 	for (int col = 0; col < Cols; col++)
    474 	{
    475 		T v = T(0);
    476 		for (int row = 0; row < Rows; row++)
    477 			v += mtx(row,col) * vec.m_data[row];
    478 		res.m_data[col] = v;
    479 	}
    480 	return res;
    481 }
    482 
    483 // Common typedefs.
    484 typedef Matrix<float, 2, 2>		Matrix2f;
    485 typedef Matrix<float, 3, 3>		Matrix3f;
    486 typedef Matrix<float, 4, 4>		Matrix4f;
    487 typedef Matrix<double, 2, 2>	Matrix2d;
    488 typedef Matrix<double, 3, 3>	Matrix3d;
    489 typedef Matrix<double, 4, 4>	Matrix4d;
    490 
    491 // GLSL-style naming \note CxR.
    492 typedef Matrix2f			Mat2;
    493 typedef Matrix<float, 3, 2>	Mat2x3;
    494 typedef Matrix<float, 4, 2>	Mat2x4;
    495 typedef Matrix<float, 2, 3>	Mat3x2;
    496 typedef Matrix3f			Mat3;
    497 typedef Matrix<float, 4, 3>	Mat3x4;
    498 typedef Matrix<float, 2, 4>	Mat4x2;
    499 typedef Matrix<float, 3, 4>	Mat4x3;
    500 typedef Matrix4f			Mat4;
    501 
    502 // Matrix-scalar operators.
    503 
    504 template <typename T, int Rows, int Cols>
    505 Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& mtx, T scalar)
    506 {
    507 	Matrix<T, Rows, Cols> res;
    508 	for (int col = 0; col < Cols; col++)
    509 		for (int row = 0; row < Rows; row++)
    510 			res(row, col) = mtx(row, col) + scalar;
    511 	return res;
    512 }
    513 
    514 template <typename T, int Rows, int Cols>
    515 Matrix<T, Rows, Cols> operator- (const Matrix<T, Rows, Cols>& mtx, T scalar)
    516 {
    517 	Matrix<T, Rows, Cols> res;
    518 	for (int col = 0; col < Cols; col++)
    519 		for (int row = 0; row < Rows; row++)
    520 			res(row, col) = mtx(row, col) - scalar;
    521 	return res;
    522 }
    523 
    524 template <typename T, int Rows, int Cols>
    525 Matrix<T, Rows, Cols> operator* (const Matrix<T, Rows, Cols>& mtx, T scalar)
    526 {
    527 	Matrix<T, Rows, Cols> res;
    528 	for (int col = 0; col < Cols; col++)
    529 		for (int row = 0; row < Rows; row++)
    530 			res(row, col) = mtx(row, col) * scalar;
    531 	return res;
    532 }
    533 
    534 template <typename T, int Rows, int Cols>
    535 Matrix<T, Rows, Cols> operator/ (const Matrix<T, Rows, Cols>& mtx, T scalar)
    536 {
    537 	Matrix<T, Rows, Cols> res;
    538 	for (int col = 0; col < Cols; col++)
    539 		for (int row = 0; row < Rows; row++)
    540 			res(row, col) = mtx(row, col) / scalar;
    541 	return res;
    542 }
    543 
    544 // Matrix-matrix component-wise operators.
    545 
    546 template <typename T, int Rows, int Cols>
    547 Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
    548 {
    549 	Matrix<T, Rows, Cols> res;
    550 	for (int col = 0; col < Cols; col++)
    551 		for (int row = 0; row < Rows; row++)
    552 			res(row, col) = a(row, col) + b(row, col);
    553 	return res;
    554 }
    555 
    556 template <typename T, int Rows, int Cols>
    557 Matrix<T, Rows, Cols> operator- (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
    558 {
    559 	Matrix<T, Rows, Cols> res;
    560 	for (int col = 0; col < Cols; col++)
    561 		for (int row = 0; row < Rows; row++)
    562 			res(row, col) = a(row, col) - b(row, col);
    563 	return res;
    564 }
    565 
    566 template <typename T, int Rows, int Cols>
    567 Matrix<T, Rows, Cols> operator/ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
    568 {
    569 	Matrix<T, Rows, Cols> res;
    570 	for (int col = 0; col < Cols; col++)
    571 		for (int row = 0; row < Rows; row++)
    572 			res(row, col) = a(row, col) / b(row, col);
    573 	return res;
    574 }
    575 
    576 } // tcu
    577 
    578 #endif // _TCUMATRIX_HPP
    579