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:
Array(void)36 Array (void) {}
~Array(void)37 ~Array (void) {}
38
operator [](int ndx) const39 inline T operator[] (int ndx) const { return m_data[ndx]; }
operator [](int ndx)40 inline T& operator[] (int ndx) { return m_data[ndx]; }
41
getPtr(void) const42 inline const T* getPtr (void) const { return m_data; }
getPtr(void)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
operator [](int ndx)81 Vector<T, Rows>& operator[] (int ndx) { return getColumn(ndx); }
operator [](int ndx) const82 const Vector<T, Rows>& operator[] (int ndx) const { return getColumn(ndx); }
83
operator ()(int row,int col) const84 inline const T& operator() (int row, int col) const { return m_data[col][row]; }
operator ()(int row,int col)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>
determinant(const Matrix<T,Size,Size> & mat)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>
inverse(const Matrix<T,Size,Size> & mat)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>
doDeterminant(const Matrix<T,2,2> & mat)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>
doDeterminant(const Matrix<T,3,3> & mat)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>
doDeterminant(const Matrix<T,4,4> & mat)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>
doInverse(const Matrix<T,2,2> & mat)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>
doInverse(const Matrix<T,3,3> & mat)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>
doInverse(const Matrix<T,4,4> & mat)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>
Matrix(void)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>
Matrix(const T & src)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>
Matrix(const T src[Rows * 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>
Matrix(const Vector<T,Rows> & src)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>
Matrix(const Matrix<T,Rows,Cols> & src)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>
~Matrix(void)372 Matrix<T, Rows, Cols>::~Matrix (void)
373 {
374 }
375
376 // Assignment operator.
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 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>
operator *=(const Matrix<T,Rows,Cols> & src)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>
setRow(int rowNdx,const Vector<T,Cols> & vec)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>
setColumn(int colNdx,const Vector<T,Rows> & vec)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>
getRow(int rowNdx) const408 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>
getColumn(int colNdx)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>
getColumn(int colNdx) const423 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>
getColumnMajorData(void) const429 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>
getRowMajorData(void) const440 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>
operator *(const Matrix<T,Rows0,Cols0> & a,const Matrix<T,Rows1,Cols1> & b)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>
operator *(const Matrix<T,Rows,Cols> & mtx,const Vector<T,Cols> & vec)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>
operator *(const Vector<T,Rows> & vec,const Matrix<T,Rows,Cols> & mtx)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>
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 template <typename T, int Rows, int Cols>
operator /(const Matrix<T,Rows,Cols> & mtx,T scalar)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>
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> & a,const Matrix<T,Rows,Cols> & b)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