1 //=============================================================================
2 //
3 // Copyright (c) Kitware, Inc.
4 // All rights reserved.
5 // See LICENSE.txt for details.
6 //
7 // This software is distributed WITHOUT ANY WARRANTY; without even
8 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
9 // PURPOSE. See the above copyright notice for more information.
10 //
11 // Copyright 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
12 // Copyright 2015 UT-Battelle, LLC.
13 // Copyright 2015 Los Alamos National Security.
14 //
15 // Under the terms of Contract DE-NA0003525 with NTESS,
16 // the U.S. Government retains certain rights in this software.
17 // Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
18 // Laboratory (LANL), the U.S. Government retains certain rights in
19 // this software.
20 //
21 //=============================================================================
22 #ifndef vtk_m_Matrix_h
23 #define vtk_m_Matrix_h
24
25 #include <vtkm/Assert.h>
26 #include <vtkm/Math.h>
27 #include <vtkm/TypeTraits.h>
28 #include <vtkm/Types.h>
29 #include <vtkm/VecTraits.h>
30
31 namespace vtkm
32 {
33
34 /// \brief Basic Matrix type.
35 ///
36 /// The Matrix class holds a small two dimensional array for simple linear
37 /// algebra and vector operations. VTK-m provides several Matrix-based
38 /// operations to assist in visualization computations.
39 ///
40 /// A Matrix is not intended to hold very large arrays. Rather, they are a
41 /// per-thread data structure to hold information like geometric transforms and
42 /// tensors.
43 ///
44 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
45 class Matrix
46 {
47 public:
48 using ComponentType = T;
49 static constexpr vtkm::IdComponent NUM_ROWS = NumRow;
50 static constexpr vtkm::IdComponent NUM_COLUMNS = NumCol;
51
52 VTKM_EXEC_CONT
Matrix()53 Matrix() {}
54
55 VTKM_EXEC_CONT
Matrix(const ComponentType & value)56 explicit Matrix(const ComponentType& value)
57 : Components(vtkm::Vec<ComponentType, NUM_COLUMNS>(value))
58 {
59 }
60
61 /// Brackets are used to reference a matrix like a 2D array (i.e.
62 /// matrix[row][column]).
63 ///
64 VTKM_EXEC_CONT
65 const vtkm::Vec<ComponentType, NUM_COLUMNS>& operator[](vtkm::IdComponent rowIndex) const
66 {
67 VTKM_ASSERT(rowIndex >= 0);
68 VTKM_ASSERT(rowIndex < NUM_ROWS);
69 return this->Components[rowIndex];
70 }
71
72 /// Brackets are used to referens a matrix like a 2D array i.e.
73 /// matrix[row][column].
74 ///
75 VTKM_EXEC_CONT
76 vtkm::Vec<ComponentType, NUM_COLUMNS>& operator[](vtkm::IdComponent rowIndex)
77 {
78 VTKM_ASSERT(rowIndex >= 0);
79 VTKM_ASSERT(rowIndex < NUM_ROWS);
80 return this->Components[rowIndex];
81 }
82
83 /// Parentheses are used to reference a matrix using mathematical tuple
84 /// notation i.e. matrix(row,column).
85 ///
86 VTKM_EXEC_CONT
operator()87 const ComponentType& operator()(vtkm::IdComponent rowIndex, vtkm::IdComponent colIndex) const
88 {
89 VTKM_ASSERT(rowIndex >= 0);
90 VTKM_ASSERT(rowIndex < NUM_ROWS);
91 VTKM_ASSERT(colIndex >= 0);
92 VTKM_ASSERT(colIndex < NUM_COLUMNS);
93 return (*this)[rowIndex][colIndex];
94 }
95
96 /// Parentheses are used to reference a matrix using mathematical tuple
97 /// notation i.e. matrix(row,column).
98 ///
99 VTKM_EXEC_CONT
operator()100 ComponentType& operator()(vtkm::IdComponent rowIndex, vtkm::IdComponent colIndex)
101 {
102 VTKM_ASSERT(rowIndex >= 0);
103 VTKM_ASSERT(rowIndex < NUM_ROWS);
104 VTKM_ASSERT(colIndex >= 0);
105 VTKM_ASSERT(colIndex < NUM_COLUMNS);
106 return (*this)[rowIndex][colIndex];
107 }
108
109 private:
110 vtkm::Vec<vtkm::Vec<ComponentType, NUM_COLUMNS>, NUM_ROWS> Components;
111 };
112
113 /// Returns a tuple containing the given row (indexed from 0) of the given
114 /// matrix.
115 ///
116 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
MatrixGetRow(const vtkm::Matrix<T,NumRow,NumCol> & matrix,vtkm::IdComponent rowIndex)117 VTKM_EXEC_CONT const vtkm::Vec<T, NumCol>& MatrixGetRow(
118 const vtkm::Matrix<T, NumRow, NumCol>& matrix,
119 vtkm::IdComponent rowIndex)
120 {
121 return matrix[rowIndex];
122 }
123
124 /// Returns a tuple containing the given column (indexed from 0) of the given
125 /// matrix. Might not be as efficient as the \c MatrixGetRow function.
126 ///
127 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
MatrixGetColumn(const vtkm::Matrix<T,NumRow,NumCol> & matrix,vtkm::IdComponent columnIndex)128 VTKM_EXEC_CONT vtkm::Vec<T, NumRow> MatrixGetColumn(const vtkm::Matrix<T, NumRow, NumCol>& matrix,
129 vtkm::IdComponent columnIndex)
130 {
131 vtkm::Vec<T, NumRow> columnValues;
132 for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
133 {
134 columnValues[rowIndex] = matrix(rowIndex, columnIndex);
135 }
136 return columnValues;
137 }
138
139 /// Convenience function for setting a row of a matrix.
140 ///
141 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
MatrixSetRow(vtkm::Matrix<T,NumRow,NumCol> & matrix,vtkm::IdComponent rowIndex,const vtkm::Vec<T,NumCol> & rowValues)142 VTKM_EXEC_CONT void MatrixSetRow(vtkm::Matrix<T, NumRow, NumCol>& matrix,
143 vtkm::IdComponent rowIndex,
144 const vtkm::Vec<T, NumCol>& rowValues)
145 {
146 matrix[rowIndex] = rowValues;
147 }
148
149 /// Convenience function for setting a column of a matrix.
150 ///
151 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
MatrixSetColumn(vtkm::Matrix<T,NumRow,NumCol> & matrix,vtkm::IdComponent columnIndex,const vtkm::Vec<T,NumRow> & columnValues)152 VTKM_EXEC_CONT void MatrixSetColumn(vtkm::Matrix<T, NumRow, NumCol>& matrix,
153 vtkm::IdComponent columnIndex,
154 const vtkm::Vec<T, NumRow>& columnValues)
155 {
156 for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
157 {
158 matrix(rowIndex, columnIndex) = columnValues[rowIndex];
159 }
160 }
161
162 /// Standard matrix multiplication.
163 ///
164 template <typename T,
165 vtkm::IdComponent NumRow,
166 vtkm::IdComponent NumCol,
167 vtkm::IdComponent NumInternal>
MatrixMultiply(const vtkm::Matrix<T,NumRow,NumInternal> & leftFactor,const vtkm::Matrix<T,NumInternal,NumCol> & rightFactor)168 VTKM_EXEC_CONT vtkm::Matrix<T, NumRow, NumCol> MatrixMultiply(
169 const vtkm::Matrix<T, NumRow, NumInternal>& leftFactor,
170 const vtkm::Matrix<T, NumInternal, NumCol>& rightFactor)
171 {
172 vtkm::Matrix<T, NumRow, NumCol> result;
173 for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
174 {
175 for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
176 {
177 T sum = T(leftFactor(rowIndex, 0) * rightFactor(0, colIndex));
178 for (vtkm::IdComponent internalIndex = 1; internalIndex < NumInternal; internalIndex++)
179 {
180 sum = T(sum + (leftFactor(rowIndex, internalIndex) * rightFactor(internalIndex, colIndex)));
181 }
182 result(rowIndex, colIndex) = sum;
183 }
184 }
185 return result;
186 }
187
188 /// Standard matrix-vector multiplication.
189 ///
190 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
MatrixMultiply(const vtkm::Matrix<T,NumRow,NumCol> & leftFactor,const vtkm::Vec<T,NumCol> & rightFactor)191 VTKM_EXEC_CONT vtkm::Vec<T, NumRow> MatrixMultiply(
192 const vtkm::Matrix<T, NumRow, NumCol>& leftFactor,
193 const vtkm::Vec<T, NumCol>& rightFactor)
194 {
195 vtkm::Vec<T, NumRow> product;
196 for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
197 {
198 product[rowIndex] = vtkm::Dot(vtkm::MatrixGetRow(leftFactor, rowIndex), rightFactor);
199 }
200 return product;
201 }
202
203 /// Standard vector-matrix multiplication
204 ///
205 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
MatrixMultiply(const vtkm::Vec<T,NumRow> & leftFactor,const vtkm::Matrix<T,NumRow,NumCol> & rightFactor)206 VTKM_EXEC_CONT vtkm::Vec<T, NumCol> MatrixMultiply(
207 const vtkm::Vec<T, NumRow>& leftFactor,
208 const vtkm::Matrix<T, NumRow, NumCol>& rightFactor)
209 {
210 vtkm::Vec<T, NumCol> product;
211 for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
212 {
213 product[colIndex] = vtkm::Dot(leftFactor, vtkm::MatrixGetColumn(rightFactor, colIndex));
214 }
215 return product;
216 }
217
218 /// Returns the identity matrix.
219 ///
220 template <typename T, vtkm::IdComponent Size>
MatrixIdentity()221 VTKM_EXEC_CONT vtkm::Matrix<T, Size, Size> MatrixIdentity()
222 {
223 vtkm::Matrix<T, Size, Size> result(T(0));
224 for (vtkm::IdComponent index = 0; index < Size; index++)
225 {
226 result(index, index) = T(1.0);
227 }
228 return result;
229 }
230
231 /// Fills the given matrix with the identity matrix.
232 ///
233 template <typename T, vtkm::IdComponent Size>
MatrixIdentity(vtkm::Matrix<T,Size,Size> & matrix)234 VTKM_EXEC_CONT void MatrixIdentity(vtkm::Matrix<T, Size, Size>& matrix)
235 {
236 matrix = vtkm::MatrixIdentity<T, Size>();
237 }
238
239 /// Returns the transpose of the given matrix.
240 ///
241 template <typename T, vtkm::IdComponent NumRows, vtkm::IdComponent NumCols>
MatrixTranspose(const vtkm::Matrix<T,NumRows,NumCols> & matrix)242 VTKM_EXEC_CONT vtkm::Matrix<T, NumCols, NumRows> MatrixTranspose(
243 const vtkm::Matrix<T, NumRows, NumCols>& matrix)
244 {
245 vtkm::Matrix<T, NumCols, NumRows> result;
246 for (vtkm::IdComponent index = 0; index < NumRows; index++)
247 {
248 vtkm::MatrixSetColumn(result, index, vtkm::MatrixGetRow(matrix, index));
249 #ifdef VTKM_ICC
250 // For reasons I do not really understand, the Intel compiler with with
251 // optimization on is sometimes removing this for loop. It appears that the
252 // optimizer sometimes does not recognize that the MatrixSetColumn function
253 // has side effects. I cannot fathom any reason for this other than a bug in
254 // the compiler, but unfortunately I do not know a reliable way to
255 // demonstrate the problem.
256 __asm__("");
257 #endif
258 }
259 return result;
260 }
261
262 namespace detail
263 {
264
265 // Used with MatrixLUPFactor.
266 template <typename T, vtkm::IdComponent Size>
MatrixLUPFactorFindPivot(vtkm::Matrix<T,Size,Size> & A,vtkm::Vec<vtkm::IdComponent,Size> & permutation,vtkm::IdComponent topCornerIndex,T & inversionParity,bool & valid)267 VTKM_EXEC_CONT void MatrixLUPFactorFindPivot(vtkm::Matrix<T, Size, Size>& A,
268 vtkm::Vec<vtkm::IdComponent, Size>& permutation,
269 vtkm::IdComponent topCornerIndex,
270 T& inversionParity,
271 bool& valid)
272 {
273 vtkm::IdComponent maxRowIndex = topCornerIndex;
274 T maxValue = vtkm::Abs(A(maxRowIndex, topCornerIndex));
275 for (vtkm::IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++)
276 {
277 T compareValue = vtkm::Abs(A(rowIndex, topCornerIndex));
278 if (maxValue < compareValue)
279 {
280 maxValue = compareValue;
281 maxRowIndex = rowIndex;
282 }
283 }
284
285 if (maxValue < vtkm::Epsilon<T>())
286 {
287 valid = false;
288 }
289
290 if (maxRowIndex != topCornerIndex)
291 {
292 // Swap rows in matrix.
293 vtkm::Vec<T, Size> maxRow = vtkm::MatrixGetRow(A, maxRowIndex);
294 vtkm::MatrixSetRow(A, maxRowIndex, vtkm::MatrixGetRow(A, topCornerIndex));
295 vtkm::MatrixSetRow(A, topCornerIndex, maxRow);
296
297 // Record change in permutation matrix.
298 vtkm::IdComponent maxOriginalRowIndex = permutation[maxRowIndex];
299 permutation[maxRowIndex] = permutation[topCornerIndex];
300 permutation[topCornerIndex] = maxOriginalRowIndex;
301
302 // Keep track of inversion parity.
303 inversionParity = -inversionParity;
304 }
305 }
306
307 // Used with MatrixLUPFactor
308 template <typename T, vtkm::IdComponent Size>
MatrixLUPFactorFindUpperTriangleElements(vtkm::Matrix<T,Size,Size> & A,vtkm::IdComponent topCornerIndex)309 VTKM_EXEC_CONT void MatrixLUPFactorFindUpperTriangleElements(vtkm::Matrix<T, Size, Size>& A,
310 vtkm::IdComponent topCornerIndex)
311 {
312 // Compute values for upper triangle on row topCornerIndex
313 for (vtkm::IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
314 {
315 A(topCornerIndex, colIndex) /= A(topCornerIndex, topCornerIndex);
316 }
317
318 // Update the rest of the matrix for calculations on subsequent rows
319 for (vtkm::IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++)
320 {
321 for (vtkm::IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
322 {
323 A(rowIndex, colIndex) -= A(rowIndex, topCornerIndex) * A(topCornerIndex, colIndex);
324 }
325 }
326 }
327
328 /// Performs an LUP-factorization on the given matrix using Crout's method. The
329 /// LU-factorization takes a matrix A and decomposes it into a lower triangular
330 /// matrix L and upper triangular matrix U such that A = LU. The
331 /// LUP-factorization also allows permutation of A, which makes the
332 /// decomposition always possible so long as A is not singular. In addition to
333 /// matrices L and U, LUP also finds permutation matrix P containing all zeros
334 /// except one 1 per row and column such that PA = LU.
335 ///
336 /// The result is done in place such that the lower triangular matrix, L, is
337 /// stored in the lower-left triangle of A including the diagonal. The upper
338 /// triangular matrix, U, is stored in the upper-right triangle of L not
339 /// including the diagonal. The diagonal of U in Crout's method is all 1's (and
340 /// therefore not explicitly stored).
341 ///
342 /// The permutation matrix P is represented by the permutation vector. If
343 /// permutation[i] = j then row j in the original matrix A has been moved to
344 /// row i in the resulting matrices. The permutation matrix P can be
345 /// represented by a matrix with p_i,j = 1 if permutation[i] = j and 0
346 /// otherwise. If using LUP-factorization to compute a determinant, you also
347 /// need to know the parity (whether there is an odd or even amount) of
348 /// inversions. An inversion is an instance of a smaller number appearing after
349 /// a larger number in permutation. Although you can compute the inversion
350 /// parity after the fact, this function keeps track of it with much less
351 /// compute resources. The parameter inversionParity is set to 1.0 for even
352 /// parity and -1.0 for odd parity.
353 ///
354 /// Not all matrices (specifically singular matrices) have an
355 /// LUP-factorization. If the LUP-factorization succeeds, valid is set to true.
356 /// Otherwise, valid is set to false and the result is indeterminant.
357 ///
358 template <typename T, vtkm::IdComponent Size>
MatrixLUPFactor(vtkm::Matrix<T,Size,Size> & A,vtkm::Vec<vtkm::IdComponent,Size> & permutation,T & inversionParity,bool & valid)359 VTKM_EXEC_CONT void MatrixLUPFactor(vtkm::Matrix<T, Size, Size>& A,
360 vtkm::Vec<vtkm::IdComponent, Size>& permutation,
361 T& inversionParity,
362 bool& valid)
363 {
364 // Initialize permutation.
365 for (vtkm::IdComponent index = 0; index < Size; index++)
366 {
367 permutation[index] = index;
368 }
369 inversionParity = T(1);
370 valid = true;
371
372 for (vtkm::IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
373 {
374 MatrixLUPFactorFindPivot(A, permutation, rowIndex, inversionParity, valid);
375 MatrixLUPFactorFindUpperTriangleElements(A, rowIndex);
376 }
377 }
378
379 /// Use a previous factorization done with MatrixLUPFactor to solve the
380 /// system Ax = b. Instead of A, this method takes in the LU and P
381 /// matrices calculated by MatrixLUPFactor from A. The x matrix is returned.
382 ///
383 template <typename T, vtkm::IdComponent Size>
MatrixLUPSolve(const vtkm::Matrix<T,Size,Size> & LU,const vtkm::Vec<vtkm::IdComponent,Size> & permutation,const vtkm::Vec<T,Size> & b)384 VTKM_EXEC_CONT vtkm::Vec<T, Size> MatrixLUPSolve(
385 const vtkm::Matrix<T, Size, Size>& LU,
386 const vtkm::Vec<vtkm::IdComponent, Size>& permutation,
387 const vtkm::Vec<T, Size>& b)
388 {
389 // The LUP-factorization gives us PA = LU or equivalently A = inv(P)LU.
390 // Substituting into Ax = b gives us inv(P)LUx = b or LUx = Pb.
391 // Now consider the intermediate vector y = Ux.
392 // Substituting in the previous two equations yields Ly = Pb.
393 // Solving Ly = Pb is easy because L is triangular and P is just a
394 // permutation.
395 vtkm::Vec<T, Size> y;
396 for (vtkm::IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
397 {
398 y[rowIndex] = b[permutation[rowIndex]];
399 // Recall that L is stored in the lower triangle of LU including diagonal.
400 for (vtkm::IdComponent colIndex = 0; colIndex < rowIndex; colIndex++)
401 {
402 y[rowIndex] -= LU(rowIndex, colIndex) * y[colIndex];
403 }
404 y[rowIndex] /= LU(rowIndex, rowIndex);
405 }
406
407 // Now that we have y, we can easily solve Ux = y for x.
408 vtkm::Vec<T, Size> x;
409 for (vtkm::IdComponent rowIndex = Size - 1; rowIndex >= 0; rowIndex--)
410 {
411 // Recall that U is stored in the upper triangle of LU with the diagonal
412 // implicitly all 1's.
413 x[rowIndex] = y[rowIndex];
414 for (vtkm::IdComponent colIndex = rowIndex + 1; colIndex < Size; colIndex++)
415 {
416 x[rowIndex] -= LU(rowIndex, colIndex) * x[colIndex];
417 }
418 }
419
420 return x;
421 }
422
423 } // namespace detail
424
425 /// Solve the linear system Ax = b for x. If a single solution is found, valid
426 /// is set to true, false otherwise.
427 ///
428 template <typename T, vtkm::IdComponent Size>
SolveLinearSystem(const vtkm::Matrix<T,Size,Size> & A,const vtkm::Vec<T,Size> & b,bool & valid)429 VTKM_EXEC_CONT vtkm::Vec<T, Size> SolveLinearSystem(const vtkm::Matrix<T, Size, Size>& A,
430 const vtkm::Vec<T, Size>& b,
431 bool& valid)
432 {
433 // First, we will make an LUP-factorization to help us.
434 vtkm::Matrix<T, Size, Size> LU = A;
435 vtkm::Vec<vtkm::IdComponent, Size> permutation;
436 T inversionParity; // Unused.
437 detail::MatrixLUPFactor(LU, permutation, inversionParity, valid);
438
439 // Next, use the decomposition to solve the system.
440 return detail::MatrixLUPSolve(LU, permutation, b);
441 }
442
443 /// Find and return the inverse of the given matrix. If the matrix is singular,
444 /// the inverse will not be correct and valid will be set to false.
445 ///
446 template <typename T, vtkm::IdComponent Size>
MatrixInverse(const vtkm::Matrix<T,Size,Size> & A,bool & valid)447 VTKM_EXEC_CONT vtkm::Matrix<T, Size, Size> MatrixInverse(const vtkm::Matrix<T, Size, Size>& A,
448 bool& valid)
449 {
450 // First, we will make an LUP-factorization to help us.
451 vtkm::Matrix<T, Size, Size> LU = A;
452 vtkm::Vec<vtkm::IdComponent, Size> permutation;
453 T inversionParity; // Unused
454 detail::MatrixLUPFactor(LU, permutation, inversionParity, valid);
455
456 // We will use the decomposition to solve AX = I for X where X is
457 // clearly the inverse of A. Our solve method only works for vectors,
458 // so we solve for one column of invA at a time.
459 vtkm::Matrix<T, Size, Size> invA;
460 vtkm::Vec<T, Size> ICol(T(0));
461 for (vtkm::IdComponent colIndex = 0; colIndex < Size; colIndex++)
462 {
463 ICol[colIndex] = 1;
464 vtkm::Vec<T, Size> invACol = detail::MatrixLUPSolve(LU, permutation, ICol);
465 ICol[colIndex] = 0;
466 vtkm::MatrixSetColumn(invA, colIndex, invACol);
467 }
468 return invA;
469 }
470
471 /// Compute the determinant of a matrix.
472 ///
473 template <typename T, vtkm::IdComponent Size>
MatrixDeterminant(const vtkm::Matrix<T,Size,Size> & A)474 VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, Size, Size>& A)
475 {
476 // First, we will make an LUP-factorization to help us.
477 vtkm::Matrix<T, Size, Size> LU = A;
478 vtkm::Vec<vtkm::IdComponent, Size> permutation;
479 T inversionParity;
480 bool valid;
481 detail::MatrixLUPFactor(LU, permutation, inversionParity, valid);
482
483 // If the matrix is singular, the factorization is invalid, but in that
484 // case we know that the determinant is 0.
485 if (!valid)
486 {
487 return 0;
488 }
489
490 // The determinant is equal to the product of the diagonal of the L matrix,
491 // possibly negated depending on the parity of the inversion. The
492 // inversionParity variable is set to 1.0 and -1.0 for even and odd parity,
493 // respectively. This sign determines whether the product should be negated.
494 T product = inversionParity;
495 for (vtkm::IdComponent index = 0; index < Size; index++)
496 {
497 product *= LU(index, index);
498 }
499 return product;
500 }
501
502 // Specializations for common small determinants.
503
504 template <typename T>
MatrixDeterminant(const vtkm::Matrix<T,1,1> & A)505 VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 1, 1>& A)
506 {
507 return A(0, 0);
508 }
509
510 template <typename T>
MatrixDeterminant(const vtkm::Matrix<T,2,2> & A)511 VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 2, 2>& A)
512 {
513 return A(0, 0) * A(1, 1) - A(1, 0) * A(0, 1);
514 }
515
516 template <typename T>
MatrixDeterminant(const vtkm::Matrix<T,3,3> & A)517 VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 3, 3>& A)
518 {
519 return A(0, 0) * A(1, 1) * A(2, 2) + A(1, 0) * A(2, 1) * A(0, 2) + A(2, 0) * A(0, 1) * A(1, 2) -
520 A(0, 0) * A(2, 1) * A(1, 2) - A(1, 0) * A(0, 1) * A(2, 2) - A(2, 0) * A(1, 1) * A(0, 2);
521 }
522
523 //---------------------------------------------------------------------------
524 // Implementations of traits for matrices.
525
526 /// Tag used to identify 2 dimensional types (matrices). A TypeTraits class
527 /// will typedef this class to DimensionalityTag.
528 ///
529 struct TypeTraitsMatrixTag
530 {
531 };
532
533 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
534 struct TypeTraits<vtkm::Matrix<T, NumRow, NumCol>>
535 {
536 using NumericTag = typename TypeTraits<T>::NumericTag;
537 using DimensionalityTag = vtkm::TypeTraitsMatrixTag;
538 };
539
540 /// A matrix has vector traits to implement component-wise operations.
541 ///
542 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
543 struct VecTraits<vtkm::Matrix<T, NumRow, NumCol>>
544 {
545 private:
546 using MatrixType = vtkm::Matrix<T, NumRow, NumCol>;
547
548 public:
549 using ComponentType = T;
550 static constexpr vtkm::IdComponent NUM_COMPONENTS = NumRow * NumCol;
551 using HasMultipleComponents = vtkm::VecTraitsTagMultipleComponents;
552 using IsSizeStatic = vtkm::VecTraitsTagSizeStatic;
553
554 VTKM_EXEC_CONT
555 static vtkm::IdComponent GetNumberOfComponents(const MatrixType&) { return NUM_COMPONENTS; }
556
557 VTKM_EXEC_CONT
558 static const ComponentType& GetComponent(const MatrixType& matrix, vtkm::IdComponent component)
559 {
560 vtkm::IdComponent colIndex = component % NumCol;
561 vtkm::IdComponent rowIndex = component / NumCol;
562 return matrix(rowIndex, colIndex);
563 }
564 VTKM_EXEC_CONT
565 static ComponentType& GetComponent(MatrixType& matrix, vtkm::IdComponent component)
566 {
567 vtkm::IdComponent colIndex = component % NumCol;
568 vtkm::IdComponent rowIndex = component / NumCol;
569 return matrix(rowIndex, colIndex);
570 }
571 VTKM_EXEC_CONT
572 static void SetComponent(MatrixType& matrix, vtkm::IdComponent component, T value)
573 {
574 GetComponent(matrix, component) = value;
575 }
576 };
577
578 } // namespace vtkm
579
580 //---------------------------------------------------------------------------
581 // Basic comparison operators.
582
583 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
584 VTKM_EXEC_CONT bool operator==(const vtkm::Matrix<T, NumRow, NumCol>& a,
585 const vtkm::Matrix<T, NumRow, NumCol>& b)
586 {
587 for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
588 {
589 for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
590 {
591 if (a(rowIndex, colIndex) != b(rowIndex, colIndex))
592 return false;
593 }
594 }
595 return true;
596 }
597 template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
598 VTKM_EXEC_CONT bool operator!=(const vtkm::Matrix<T, NumRow, NumCol>& a,
599 const vtkm::Matrix<T, NumRow, NumCol>& b)
600 {
601 return !(a == b);
602 }
603
604 #endif //vtk_m_Matrix_h
605