framework/common/tcuMatrix.hpp

   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