Matrix.h 21 KB
Newer Older
Kenneth Moreland's avatar
Kenneth Moreland committed
1 2 3 4 5 6 7 8 9 10
//=============================================================================
//
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt for details.
//
//  This software is distributed WITHOUT ANY WARRANTY; without even
//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
//  PURPOSE.  See the above copyright notice for more information.
//
11
//  Copyright 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
Kenneth Moreland's avatar
Kenneth Moreland committed
12 13 14
//  Copyright 2015 UT-Battelle, LLC.
//  Copyright 2015 Los Alamos National Security.
//
15
//  Under the terms of Contract DE-NA0003525 with NTESS,
Kenneth Moreland's avatar
Kenneth Moreland committed
16 17 18 19 20 21 22 23 24
//  the U.S. Government retains certain rights in this software.
//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
//  Laboratory (LANL), the U.S. Government retains certain rights in
//  this software.
//
//=============================================================================
#ifndef vtk_m_Matrix_h
#define vtk_m_Matrix_h

Kenneth Moreland's avatar
Kenneth Moreland committed
25
#include <vtkm/Assert.h>
Kenneth Moreland's avatar
Kenneth Moreland committed
26 27
#include <vtkm/Math.h>
#include <vtkm/TypeTraits.h>
28
#include <vtkm/Types.h>
Kenneth Moreland's avatar
Kenneth Moreland committed
29 30
#include <vtkm/VecTraits.h>

31 32
namespace vtkm
{
Kenneth Moreland's avatar
Kenneth Moreland committed
33 34 35 36 37 38 39 40 41 42 43

/// \brief Basic Matrix type.
///
/// The Matrix class holds a small two dimensional array for simple linear
/// algebra and vector operations. VTK-m provides several Matrix-based
/// operations to assist in visualization computations.
///
/// A Matrix is not intended to hold very large arrays. Rather, they are a
/// per-thread data structure to hold information like geometric transforms and
/// tensors.
///
44 45 46
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
class Matrix
{
Kenneth Moreland's avatar
Kenneth Moreland committed
47
public:
48
  using ComponentType = T;
49 50
  static constexpr vtkm::IdComponent NUM_ROWS = NumRow;
  static constexpr vtkm::IdComponent NUM_COLUMNS = NumCol;
Kenneth Moreland's avatar
Kenneth Moreland committed
51

52
  VTKM_EXEC_CONT
53
  Matrix() {}
Kenneth Moreland's avatar
Kenneth Moreland committed
54

55
  VTKM_EXEC_CONT
56 57 58 59
  explicit Matrix(const ComponentType& value)
    : Components(vtkm::Vec<ComponentType, NUM_COLUMNS>(value))
  {
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
60 61 62 63

  /// Brackets are used to reference a matrix like a 2D array (i.e.
  /// matrix[row][column]).
  ///
64
  VTKM_EXEC_CONT
65 66
  const vtkm::Vec<ComponentType, NUM_COLUMNS>& operator[](vtkm::IdComponent rowIndex) const
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
67
    VTKM_ASSERT(rowIndex >= 0);
68
    VTKM_ASSERT(rowIndex < NUM_ROWS);
Kenneth Moreland's avatar
Kenneth Moreland committed
69 70 71 72 73 74
    return this->Components[rowIndex];
  }

  /// Brackets are used to referens a matrix like a 2D array i.e.
  /// matrix[row][column].
  ///
75
  VTKM_EXEC_CONT
76 77
  vtkm::Vec<ComponentType, NUM_COLUMNS>& operator[](vtkm::IdComponent rowIndex)
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
78
    VTKM_ASSERT(rowIndex >= 0);
79
    VTKM_ASSERT(rowIndex < NUM_ROWS);
Kenneth Moreland's avatar
Kenneth Moreland committed
80 81 82 83 84 85
    return this->Components[rowIndex];
  }

  /// Parentheses are used to reference a matrix using mathematical tuple
  /// notation i.e. matrix(row,column).
  ///
86
  VTKM_EXEC_CONT
87 88
  const ComponentType& operator()(vtkm::IdComponent rowIndex, vtkm::IdComponent colIndex) const
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
89
    VTKM_ASSERT(rowIndex >= 0);
90
    VTKM_ASSERT(rowIndex < NUM_ROWS);
Kenneth Moreland's avatar
Kenneth Moreland committed
91
    VTKM_ASSERT(colIndex >= 0);
92
    VTKM_ASSERT(colIndex < NUM_COLUMNS);
Kenneth Moreland's avatar
Kenneth Moreland committed
93 94 95 96 97 98
    return (*this)[rowIndex][colIndex];
  }

  /// Parentheses are used to reference a matrix using mathematical tuple
  /// notation i.e. matrix(row,column).
  ///
99
  VTKM_EXEC_CONT
100 101
  ComponentType& operator()(vtkm::IdComponent rowIndex, vtkm::IdComponent colIndex)
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
102
    VTKM_ASSERT(rowIndex >= 0);
103
    VTKM_ASSERT(rowIndex < NUM_ROWS);
Kenneth Moreland's avatar
Kenneth Moreland committed
104
    VTKM_ASSERT(colIndex >= 0);
105
    VTKM_ASSERT(colIndex < NUM_COLUMNS);
Kenneth Moreland's avatar
Kenneth Moreland committed
106 107 108 109 110 111 112 113 114 115
    return (*this)[rowIndex][colIndex];
  }

private:
  vtkm::Vec<vtkm::Vec<ComponentType, NUM_COLUMNS>, NUM_ROWS> Components;
};

/// Returns a tuple containing the given row (indexed from 0) of the given
/// matrix.
///
116 117
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT const vtkm::Vec<T, NumCol>& MatrixGetRow(
118 119
  const vtkm::Matrix<T, NumRow, NumCol>& matrix,
  vtkm::IdComponent rowIndex)
Kenneth Moreland's avatar
Kenneth Moreland committed
120 121 122 123 124 125 126
{
  return matrix[rowIndex];
}

/// Returns a tuple containing the given column (indexed from 0) of the given
/// matrix.  Might not be as efficient as the \c MatrixGetRow function.
///
127 128 129
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT vtkm::Vec<T, NumRow> MatrixGetColumn(const vtkm::Matrix<T, NumRow, NumCol>& matrix,
                                                    vtkm::IdComponent columnIndex)
Kenneth Moreland's avatar
Kenneth Moreland committed
130 131 132 133 134 135 136 137 138 139 140
{
  vtkm::Vec<T, NumRow> columnValues;
  for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
  {
    columnValues[rowIndex] = matrix(rowIndex, columnIndex);
  }
  return columnValues;
}

/// Convenience function for setting a row of a matrix.
///
141 142
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT void MatrixSetRow(vtkm::Matrix<T, NumRow, NumCol>& matrix,
143 144
                                 vtkm::IdComponent rowIndex,
                                 const vtkm::Vec<T, NumCol>& rowValues)
Kenneth Moreland's avatar
Kenneth Moreland committed
145 146 147 148 149 150
{
  matrix[rowIndex] = rowValues;
}

/// Convenience function for setting a column of a matrix.
///
151 152 153 154
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT void MatrixSetColumn(vtkm::Matrix<T, NumRow, NumCol>& matrix,
                                    vtkm::IdComponent columnIndex,
                                    const vtkm::Vec<T, NumRow>& columnValues)
Kenneth Moreland's avatar
Kenneth Moreland committed
155 156 157 158 159 160 161 162 163
{
  for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
  {
    matrix(rowIndex, columnIndex) = columnValues[rowIndex];
  }
}

/// Standard matrix multiplication.
///
164 165 166
template <typename T,
          vtkm::IdComponent NumRow,
          vtkm::IdComponent NumCol,
167 168 169 170
          vtkm::IdComponent NumInternal>
VTKM_EXEC_CONT vtkm::Matrix<T, NumRow, NumCol> MatrixMultiply(
  const vtkm::Matrix<T, NumRow, NumInternal>& leftFactor,
  const vtkm::Matrix<T, NumInternal, NumCol>& rightFactor)
Kenneth Moreland's avatar
Kenneth Moreland committed
171
{
172
  vtkm::Matrix<T, NumRow, NumCol> result;
Kenneth Moreland's avatar
Kenneth Moreland committed
173 174 175 176
  for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
  {
    for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
    {
177
      T sum = T(leftFactor(rowIndex, 0) * rightFactor(0, colIndex));
178
      for (vtkm::IdComponent internalIndex = 1; internalIndex < NumInternal; internalIndex++)
Kenneth Moreland's avatar
Kenneth Moreland committed
179
      {
180
        sum = T(sum + (leftFactor(rowIndex, internalIndex) * rightFactor(internalIndex, colIndex)));
Kenneth Moreland's avatar
Kenneth Moreland committed
181 182 183 184 185 186 187 188 189
      }
      result(rowIndex, colIndex) = sum;
    }
  }
  return result;
}

/// Standard matrix-vector multiplication.
///
190 191
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT vtkm::Vec<T, NumRow> MatrixMultiply(
192 193
  const vtkm::Matrix<T, NumRow, NumCol>& leftFactor,
  const vtkm::Vec<T, NumCol>& rightFactor)
Kenneth Moreland's avatar
Kenneth Moreland committed
194
{
195
  vtkm::Vec<T, NumRow> product;
Kenneth Moreland's avatar
Kenneth Moreland committed
196 197
  for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
  {
198
    product[rowIndex] = vtkm::Dot(vtkm::MatrixGetRow(leftFactor, rowIndex), rightFactor);
Kenneth Moreland's avatar
Kenneth Moreland committed
199 200 201 202 203 204
  }
  return product;
}

/// Standard vector-matrix multiplication
///
205 206
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT vtkm::Vec<T, NumCol> MatrixMultiply(
207 208
  const vtkm::Vec<T, NumRow>& leftFactor,
  const vtkm::Matrix<T, NumRow, NumCol>& rightFactor)
Kenneth Moreland's avatar
Kenneth Moreland committed
209
{
210
  vtkm::Vec<T, NumCol> product;
Kenneth Moreland's avatar
Kenneth Moreland committed
211 212
  for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
  {
213
    product[colIndex] = vtkm::Dot(leftFactor, vtkm::MatrixGetColumn(rightFactor, colIndex));
Kenneth Moreland's avatar
Kenneth Moreland committed
214 215 216 217 218 219
  }
  return product;
}

/// Returns the identity matrix.
///
220 221
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT vtkm::Matrix<T, Size, Size> MatrixIdentity()
Kenneth Moreland's avatar
Kenneth Moreland committed
222
{
223
  vtkm::Matrix<T, Size, Size> result(T(0));
Kenneth Moreland's avatar
Kenneth Moreland committed
224 225
  for (vtkm::IdComponent index = 0; index < Size; index++)
  {
226
    result(index, index) = T(1.0);
Kenneth Moreland's avatar
Kenneth Moreland committed
227 228 229 230 231 232
  }
  return result;
}

/// Fills the given matrix with the identity matrix.
///
233 234
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT void MatrixIdentity(vtkm::Matrix<T, Size, Size>& matrix)
Kenneth Moreland's avatar
Kenneth Moreland committed
235
{
236
  matrix = vtkm::MatrixIdentity<T, Size>();
Kenneth Moreland's avatar
Kenneth Moreland committed
237 238 239 240
}

/// Returns the transpose of the given matrix.
///
241 242 243
template <typename T, vtkm::IdComponent NumRows, vtkm::IdComponent NumCols>
VTKM_EXEC_CONT vtkm::Matrix<T, NumCols, NumRows> MatrixTranspose(
  const vtkm::Matrix<T, NumRows, NumCols>& matrix)
Kenneth Moreland's avatar
Kenneth Moreland committed
244
{
245
  vtkm::Matrix<T, NumCols, NumRows> result;
Kenneth Moreland's avatar
Kenneth Moreland committed
246 247 248
  for (vtkm::IdComponent index = 0; index < NumRows; index++)
  {
    vtkm::MatrixSetColumn(result, index, vtkm::MatrixGetRow(matrix, index));
249 250 251 252 253 254 255 256 257
#ifdef VTKM_ICC
    // For reasons I do not really understand, the Intel compiler with with
    // optimization on is sometimes removing this for loop. It appears that the
    // optimizer sometimes does not recognize that the MatrixSetColumn function
    // has side effects. I cannot fathom any reason for this other than a bug in
    // the compiler, but unfortunately I do not know a reliable way to
    // demonstrate the problem.
    __asm__("");
#endif
Kenneth Moreland's avatar
Kenneth Moreland committed
258 259 260 261
  }
  return result;
}

262 263
namespace detail
{
Kenneth Moreland's avatar
Kenneth Moreland committed
264 265

// Used with MatrixLUPFactor.
266 267 268
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT void MatrixLUPFactorFindPivot(vtkm::Matrix<T, Size, Size>& A,
                                             vtkm::Vec<vtkm::IdComponent, Size>& permutation,
269 270
                                             vtkm::IdComponent topCornerIndex,
                                             T& inversionParity,
271
                                             bool& valid)
Kenneth Moreland's avatar
Kenneth Moreland committed
272 273 274
{
  vtkm::IdComponent maxRowIndex = topCornerIndex;
  T maxValue = vtkm::Abs(A(maxRowIndex, topCornerIndex));
275
  for (vtkm::IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++)
Kenneth Moreland's avatar
Kenneth Moreland committed
276 277 278 279 280 281 282 283 284
  {
    T compareValue = vtkm::Abs(A(rowIndex, topCornerIndex));
    if (maxValue < compareValue)
    {
      maxValue = compareValue;
      maxRowIndex = rowIndex;
    }
  }

285 286 287 288
  if (maxValue < vtkm::Epsilon<T>())
  {
    valid = false;
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
289 290 291 292

  if (maxRowIndex != topCornerIndex)
  {
    // Swap rows in matrix.
293 294
    vtkm::Vec<T, Size> maxRow = vtkm::MatrixGetRow(A, maxRowIndex);
    vtkm::MatrixSetRow(A, maxRowIndex, vtkm::MatrixGetRow(A, topCornerIndex));
Kenneth Moreland's avatar
Kenneth Moreland committed
295 296 297 298 299 300 301 302 303 304 305 306 307
    vtkm::MatrixSetRow(A, topCornerIndex, maxRow);

    // Record change in permutation matrix.
    vtkm::IdComponent maxOriginalRowIndex = permutation[maxRowIndex];
    permutation[maxRowIndex] = permutation[topCornerIndex];
    permutation[topCornerIndex] = maxOriginalRowIndex;

    // Keep track of inversion parity.
    inversionParity = -inversionParity;
  }
}

// Used with MatrixLUPFactor
308 309 310
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT void MatrixLUPFactorFindUpperTriangleElements(vtkm::Matrix<T, Size, Size>& A,
                                                             vtkm::IdComponent topCornerIndex)
Kenneth Moreland's avatar
Kenneth Moreland committed
311 312
{
  // Compute values for upper triangle on row topCornerIndex
313
  for (vtkm::IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
Kenneth Moreland's avatar
Kenneth Moreland committed
314
  {
315
    A(topCornerIndex, colIndex) /= A(topCornerIndex, topCornerIndex);
Kenneth Moreland's avatar
Kenneth Moreland committed
316 317 318
  }

  // Update the rest of the matrix for calculations on subsequent rows
319
  for (vtkm::IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++)
Kenneth Moreland's avatar
Kenneth Moreland committed
320
  {
321
    for (vtkm::IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
Kenneth Moreland's avatar
Kenneth Moreland committed
322
    {
323
      A(rowIndex, colIndex) -= A(rowIndex, topCornerIndex) * A(topCornerIndex, colIndex);
Kenneth Moreland's avatar
Kenneth Moreland committed
324 325 326 327 328 329 330 331
    }
  }
}

/// Performs an LUP-factorization on the given matrix using Crout's method. The
/// LU-factorization takes a matrix A and decomposes it into a lower triangular
/// matrix L and upper triangular matrix U such that A = LU. The
/// LUP-factorization also allows permutation of A, which makes the
luz.paz's avatar
luz.paz committed
332
/// decomposition always possible so long as A is not singular. In addition to
Kenneth Moreland's avatar
Kenneth Moreland committed
333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357
/// matrices L and U, LUP also finds permutation matrix P containing all zeros
/// except one 1 per row and column such that PA = LU.
///
/// The result is done in place such that the lower triangular matrix, L, is
/// stored in the lower-left triangle of A including the diagonal. The upper
/// triangular matrix, U, is stored in the upper-right triangle of L not
/// including the diagonal. The diagonal of U in Crout's method is all 1's (and
/// therefore not explicitly stored).
///
/// The permutation matrix P is represented by the permutation vector. If
/// permutation[i] = j then row j in the original matrix A has been moved to
/// row i in the resulting matrices. The permutation matrix P can be
/// represented by a matrix with p_i,j = 1 if permutation[i] = j and 0
/// otherwise. If using LUP-factorization to compute a determinant, you also
/// need to know the parity (whether there is an odd or even amount) of
/// inversions. An inversion is an instance of a smaller number appearing after
/// a larger number in permutation. Although you can compute the inversion
/// parity after the fact, this function keeps track of it with much less
/// compute resources. The parameter inversionParity is set to 1.0 for even
/// parity and -1.0 for odd parity.
///
/// Not all matrices (specifically singular matrices) have an
/// LUP-factorization. If the LUP-factorization succeeds, valid is set to true.
/// Otherwise, valid is set to false and the result is indeterminant.
///
358 359 360
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT void MatrixLUPFactor(vtkm::Matrix<T, Size, Size>& A,
                                    vtkm::Vec<vtkm::IdComponent, Size>& permutation,
361 362
                                    T& inversionParity,
                                    bool& valid)
Kenneth Moreland's avatar
Kenneth Moreland committed
363 364 365 366 367 368
{
  // Initialize permutation.
  for (vtkm::IdComponent index = 0; index < Size; index++)
  {
    permutation[index] = index;
  }
369
  inversionParity = T(1);
Kenneth Moreland's avatar
Kenneth Moreland committed
370 371 372 373 374 375 376 377 378 379 380 381 382
  valid = true;

  for (vtkm::IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
  {
    MatrixLUPFactorFindPivot(A, permutation, rowIndex, inversionParity, valid);
    MatrixLUPFactorFindUpperTriangleElements(A, rowIndex);
  }
}

/// Use a previous factorization done with MatrixLUPFactor to solve the
/// system Ax = b.  Instead of A, this method takes in the LU and P
/// matrices calculated by MatrixLUPFactor from A. The x matrix is returned.
///
383 384
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT vtkm::Vec<T, Size> MatrixLUPSolve(
385 386
  const vtkm::Matrix<T, Size, Size>& LU,
  const vtkm::Vec<vtkm::IdComponent, Size>& permutation,
387
  const vtkm::Vec<T, Size>& b)
Kenneth Moreland's avatar
Kenneth Moreland committed
388 389 390 391 392 393 394
{
  // The LUP-factorization gives us PA = LU or equivalently A = inv(P)LU.
  // Substituting into Ax = b gives us inv(P)LUx = b or LUx = Pb.
  // Now consider the intermediate vector y = Ux.
  // Substituting in the previous two equations yields Ly = Pb.
  // Solving Ly = Pb is easy because L is triangular and P is just a
  // permutation.
395
  vtkm::Vec<T, Size> y;
Kenneth Moreland's avatar
Kenneth Moreland committed
396
  for (vtkm::IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
397
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
398 399 400
    y[rowIndex] = b[permutation[rowIndex]];
    // Recall that L is stored in the lower triangle of LU including diagonal.
    for (vtkm::IdComponent colIndex = 0; colIndex < rowIndex; colIndex++)
401 402
    {
      y[rowIndex] -= LU(rowIndex, colIndex) * y[colIndex];
Kenneth Moreland's avatar
Kenneth Moreland committed
403
    }
404 405
    y[rowIndex] /= LU(rowIndex, rowIndex);
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
406 407

  // Now that we have y, we can easily solve Ux = y for x.
408 409
  vtkm::Vec<T, Size> x;
  for (vtkm::IdComponent rowIndex = Size - 1; rowIndex >= 0; rowIndex--)
Kenneth Moreland's avatar
Kenneth Moreland committed
410 411 412 413
  {
    // Recall that U is stored in the upper triangle of LU with the diagonal
    // implicitly all 1's.
    x[rowIndex] = y[rowIndex];
414
    for (vtkm::IdComponent colIndex = rowIndex + 1; colIndex < Size; colIndex++)
Kenneth Moreland's avatar
Kenneth Moreland committed
415
    {
416
      x[rowIndex] -= LU(rowIndex, colIndex) * x[colIndex];
Kenneth Moreland's avatar
Kenneth Moreland committed
417 418 419 420 421 422 423 424 425 426 427
    }
  }

  return x;
}

} // namespace detail

/// Solve the linear system Ax = b for x. If a single solution is found, valid
/// is set to true, false otherwise.
///
428 429
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT vtkm::Vec<T, Size> SolveLinearSystem(const vtkm::Matrix<T, Size, Size>& A,
430 431
                                                    const vtkm::Vec<T, Size>& b,
                                                    bool& valid)
Kenneth Moreland's avatar
Kenneth Moreland committed
432 433
{
  // First, we will make an LUP-factorization to help us.
434 435 436
  vtkm::Matrix<T, Size, Size> LU = A;
  vtkm::Vec<vtkm::IdComponent, Size> permutation;
  T inversionParity; // Unused.
Kenneth Moreland's avatar
Kenneth Moreland committed
437 438 439 440 441 442 443 444 445
  detail::MatrixLUPFactor(LU, permutation, inversionParity, valid);

  // Next, use the decomposition to solve the system.
  return detail::MatrixLUPSolve(LU, permutation, b);
}

/// Find and return the inverse of the given matrix. If the matrix is singular,
/// the inverse will not be correct and valid will be set to false.
///
446 447 448
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT vtkm::Matrix<T, Size, Size> MatrixInverse(const vtkm::Matrix<T, Size, Size>& A,
                                                         bool& valid)
Kenneth Moreland's avatar
Kenneth Moreland committed
449 450
{
  // First, we will make an LUP-factorization to help us.
451 452 453
  vtkm::Matrix<T, Size, Size> LU = A;
  vtkm::Vec<vtkm::IdComponent, Size> permutation;
  T inversionParity; // Unused
Kenneth Moreland's avatar
Kenneth Moreland committed
454 455 456 457 458
  detail::MatrixLUPFactor(LU, permutation, inversionParity, valid);

  // We will use the decomposition to solve AX = I for X where X is
  // clearly the inverse of A.  Our solve method only works for vectors,
  // so we solve for one column of invA at a time.
459 460
  vtkm::Matrix<T, Size, Size> invA;
  vtkm::Vec<T, Size> ICol(T(0));
Kenneth Moreland's avatar
Kenneth Moreland committed
461
  for (vtkm::IdComponent colIndex = 0; colIndex < Size; colIndex++)
462
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
463
    ICol[colIndex] = 1;
464
    vtkm::Vec<T, Size> invACol = detail::MatrixLUPSolve(LU, permutation, ICol);
Kenneth Moreland's avatar
Kenneth Moreland committed
465 466
    ICol[colIndex] = 0;
    vtkm::MatrixSetColumn(invA, colIndex, invACol);
467
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
468 469 470 471 472
  return invA;
}

/// Compute the determinant of a matrix.
///
473 474
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, Size, Size>& A)
Kenneth Moreland's avatar
Kenneth Moreland committed
475 476
{
  // First, we will make an LUP-factorization to help us.
477 478
  vtkm::Matrix<T, Size, Size> LU = A;
  vtkm::Vec<vtkm::IdComponent, Size> permutation;
Kenneth Moreland's avatar
Kenneth Moreland committed
479 480 481 482 483 484
  T inversionParity;
  bool valid;
  detail::MatrixLUPFactor(LU, permutation, inversionParity, valid);

  // If the matrix is singular, the factorization is invalid, but in that
  // case we know that the determinant is 0.
485 486 487 488
  if (!valid)
  {
    return 0;
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
489 490 491 492 493 494 495

  // The determinant is equal to the product of the diagonal of the L matrix,
  // possibly negated depending on the parity of the inversion. The
  // inversionParity variable is set to 1.0 and -1.0 for even and odd parity,
  // respectively. This sign determines whether the product should be negated.
  T product = inversionParity;
  for (vtkm::IdComponent index = 0; index < Size; index++)
496 497 498
  {
    product *= LU(index, index);
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
499 500 501 502 503
  return product;
}

// Specializations for common small determinants.

504 505
template <typename T>
VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 1, 1>& A)
Kenneth Moreland's avatar
Kenneth Moreland committed
506
{
507
  return A(0, 0);
Kenneth Moreland's avatar
Kenneth Moreland committed
508 509
}

510 511
template <typename T>
VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 2, 2>& A)
Kenneth Moreland's avatar
Kenneth Moreland committed
512
{
513
  return A(0, 0) * A(1, 1) - A(1, 0) * A(0, 1);
Kenneth Moreland's avatar
Kenneth Moreland committed
514 515
}

516 517
template <typename T>
VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 3, 3>& A)
Kenneth Moreland's avatar
Kenneth Moreland committed
518
{
519 520
  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) -
    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);
Kenneth Moreland's avatar
Kenneth Moreland committed
521 522 523 524 525 526 527 528
}

//---------------------------------------------------------------------------
// Implementations of traits for matrices.

/// Tag used to identify 2 dimensional types (matrices). A TypeTraits class
/// will typedef this class to DimensionalityTag.
///
529 530 531
struct TypeTraitsMatrixTag
{
};
Kenneth Moreland's avatar
Kenneth Moreland committed
532

533 534 535
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
struct TypeTraits<vtkm::Matrix<T, NumRow, NumCol>>
{
536 537
  using NumericTag = typename TypeTraits<T>::NumericTag;
  using DimensionalityTag = vtkm::TypeTraitsMatrixTag;
Kenneth Moreland's avatar
Kenneth Moreland committed
538 539 540 541
};

/// A matrix has vector traits to implement component-wise operations.
///
542 543 544
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
struct VecTraits<vtkm::Matrix<T, NumRow, NumCol>>
{
Kenneth Moreland's avatar
Kenneth Moreland committed
545
private:
546
  using MatrixType = vtkm::Matrix<T, NumRow, NumCol>;
547

Kenneth Moreland's avatar
Kenneth Moreland committed
548
public:
549
  using ComponentType = T;
550
  static constexpr vtkm::IdComponent NUM_COMPONENTS = NumRow * NumCol;
551 552
  using HasMultipleComponents = vtkm::VecTraitsTagMultipleComponents;
  using IsSizeStatic = vtkm::VecTraitsTagSizeStatic;
Kenneth Moreland's avatar
Kenneth Moreland committed
553

554
  VTKM_EXEC_CONT
555
  static vtkm::IdComponent GetNumberOfComponents(const MatrixType&) { return NUM_COMPONENTS; }
Kenneth Moreland's avatar
Kenneth Moreland committed
556

557
  VTKM_EXEC_CONT
558 559
  static const ComponentType& GetComponent(const MatrixType& matrix, vtkm::IdComponent component)
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
560 561
    vtkm::IdComponent colIndex = component % NumCol;
    vtkm::IdComponent rowIndex = component / NumCol;
562
    return matrix(rowIndex, colIndex);
Kenneth Moreland's avatar
Kenneth Moreland committed
563
  }
564
  VTKM_EXEC_CONT
565 566
  static ComponentType& GetComponent(MatrixType& matrix, vtkm::IdComponent component)
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
567 568
    vtkm::IdComponent colIndex = component % NumCol;
    vtkm::IdComponent rowIndex = component / NumCol;
569
    return matrix(rowIndex, colIndex);
Kenneth Moreland's avatar
Kenneth Moreland committed
570
  }
571
  VTKM_EXEC_CONT
572
  static void SetComponent(MatrixType& matrix, vtkm::IdComponent component, T value)
Kenneth Moreland's avatar
Kenneth Moreland committed
573 574 575 576 577 578 579 580 581 582
  {
    GetComponent(matrix, component) = value;
  }
};

} // namespace vtkm

//---------------------------------------------------------------------------
// Basic comparison operators.

583 584 585
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT bool operator==(const vtkm::Matrix<T, NumRow, NumCol>& a,
                               const vtkm::Matrix<T, NumRow, NumCol>& b)
Kenneth Moreland's avatar
Kenneth Moreland committed
586 587 588 589 590
{
  for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
  {
    for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
    {
591 592
      if (a(rowIndex, colIndex) != b(rowIndex, colIndex))
        return false;
Kenneth Moreland's avatar
Kenneth Moreland committed
593 594 595 596
    }
  }
  return true;
}
597 598 599
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT bool operator!=(const vtkm::Matrix<T, NumRow, NumCol>& a,
                               const vtkm::Matrix<T, NumRow, NumCol>& b)
Kenneth Moreland's avatar
Kenneth Moreland committed
600 601 602 603 604
{
  return !(a == b);
}

#endif //vtk_m_Matrix_h