Matrix.h 22.2 KB
Newer Older
nadavi's avatar
nadavi committed
1
//============================================================================
Kenneth Moreland's avatar
Kenneth Moreland committed
2
3
4
5
6
7
8
//  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.
nadavi's avatar
nadavi committed
9
//============================================================================
Kenneth Moreland's avatar
Kenneth Moreland committed
10
11
12
#ifndef vtk_m_Matrix_h
#define vtk_m_Matrix_h

Kenneth Moreland's avatar
Kenneth Moreland committed
13
#include <vtkm/Assert.h>
Kenneth Moreland's avatar
Kenneth Moreland committed
14
15
#include <vtkm/Math.h>
#include <vtkm/TypeTraits.h>
16
#include <vtkm/Types.h>
Kenneth Moreland's avatar
Kenneth Moreland committed
17
18
#include <vtkm/VecTraits.h>

19
20
namespace vtkm
{
Kenneth Moreland's avatar
Kenneth Moreland committed
21
22
23
24
25
26
27
28
29
30
31

/// \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.
///
32
33
34
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
class Matrix
{
Kenneth Moreland's avatar
Kenneth Moreland committed
35
public:
36
  using ComponentType = T;
37
38
  static constexpr vtkm::IdComponent NUM_ROWS = NumRow;
  static constexpr vtkm::IdComponent NUM_COLUMNS = NumCol;
Kenneth Moreland's avatar
Kenneth Moreland committed
39

40
  VTKM_EXEC_CONT
41
  Matrix() {}
Kenneth Moreland's avatar
Kenneth Moreland committed
42

43
  VTKM_EXEC_CONT
44
45
46
47
  explicit Matrix(const ComponentType& value)
    : Components(vtkm::Vec<ComponentType, NUM_COLUMNS>(value))
  {
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
48
49
50
51

  /// Brackets are used to reference a matrix like a 2D array (i.e.
  /// matrix[row][column]).
  ///
52
  VTKM_EXEC_CONT
53
54
  const vtkm::Vec<ComponentType, NUM_COLUMNS>& operator[](vtkm::IdComponent rowIndex) const
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
55
    VTKM_ASSERT(rowIndex >= 0);
56
    VTKM_ASSERT(rowIndex < NUM_ROWS);
Kenneth Moreland's avatar
Kenneth Moreland committed
57
58
59
60
61
62
    return this->Components[rowIndex];
  }

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

  /// Parentheses are used to reference a matrix using mathematical tuple
  /// notation i.e. matrix(row,column).
  ///
74
  VTKM_EXEC_CONT
75
76
  const ComponentType& operator()(vtkm::IdComponent rowIndex, vtkm::IdComponent colIndex) const
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
77
    VTKM_ASSERT(rowIndex >= 0);
78
    VTKM_ASSERT(rowIndex < NUM_ROWS);
Kenneth Moreland's avatar
Kenneth Moreland committed
79
    VTKM_ASSERT(colIndex >= 0);
80
    VTKM_ASSERT(colIndex < NUM_COLUMNS);
Kenneth Moreland's avatar
Kenneth Moreland committed
81
82
83
84
85
86
    return (*this)[rowIndex][colIndex];
  }

  /// Parentheses are used to reference a matrix using mathematical tuple
  /// notation i.e. matrix(row,column).
  ///
87
  VTKM_EXEC_CONT
88
89
  ComponentType& operator()(vtkm::IdComponent rowIndex, vtkm::IdComponent colIndex)
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
90
    VTKM_ASSERT(rowIndex >= 0);
91
    VTKM_ASSERT(rowIndex < NUM_ROWS);
Kenneth Moreland's avatar
Kenneth Moreland committed
92
    VTKM_ASSERT(colIndex >= 0);
93
    VTKM_ASSERT(colIndex < NUM_COLUMNS);
Kenneth Moreland's avatar
Kenneth Moreland committed
94
95
96
97
98
99
100
101
102
103
    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.
///
104
105
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT const vtkm::Vec<T, NumCol>& MatrixGetRow(
106
107
  const vtkm::Matrix<T, NumRow, NumCol>& matrix,
  vtkm::IdComponent rowIndex)
Kenneth Moreland's avatar
Kenneth Moreland committed
108
109
110
111
112
113
114
{
  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.
///
115
116
117
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
118
119
120
121
122
123
124
125
126
127
128
{
  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.
///
129
130
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT void MatrixSetRow(vtkm::Matrix<T, NumRow, NumCol>& matrix,
131
132
                                 vtkm::IdComponent rowIndex,
                                 const vtkm::Vec<T, NumCol>& rowValues)
Kenneth Moreland's avatar
Kenneth Moreland committed
133
134
135
136
137
138
{
  matrix[rowIndex] = rowValues;
}

/// Convenience function for setting a column of a matrix.
///
139
140
141
142
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
143
144
145
146
147
148
149
150
151
{
  for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
  {
    matrix(rowIndex, columnIndex) = columnValues[rowIndex];
  }
}

/// Standard matrix multiplication.
///
152
153
154
template <typename T,
          vtkm::IdComponent NumRow,
          vtkm::IdComponent NumCol,
155
156
157
158
          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
159
{
160
  vtkm::Matrix<T, NumRow, NumCol> result;
Kenneth Moreland's avatar
Kenneth Moreland committed
161
162
163
164
  for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
  {
    for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
    {
165
      T sum = T(leftFactor(rowIndex, 0) * rightFactor(0, colIndex));
166
      for (vtkm::IdComponent internalIndex = 1; internalIndex < NumInternal; internalIndex++)
Kenneth Moreland's avatar
Kenneth Moreland committed
167
      {
168
        sum = T(sum + (leftFactor(rowIndex, internalIndex) * rightFactor(internalIndex, colIndex)));
Kenneth Moreland's avatar
Kenneth Moreland committed
169
170
171
172
173
174
175
176
177
      }
      result(rowIndex, colIndex) = sum;
    }
  }
  return result;
}

/// Standard matrix-vector multiplication.
///
178
179
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT vtkm::Vec<T, NumRow> MatrixMultiply(
180
181
  const vtkm::Matrix<T, NumRow, NumCol>& leftFactor,
  const vtkm::Vec<T, NumCol>& rightFactor)
Kenneth Moreland's avatar
Kenneth Moreland committed
182
{
183
  vtkm::Vec<T, NumRow> product;
Kenneth Moreland's avatar
Kenneth Moreland committed
184
185
  for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
  {
186
    product[rowIndex] = vtkm::Dot(vtkm::MatrixGetRow(leftFactor, rowIndex), rightFactor);
Kenneth Moreland's avatar
Kenneth Moreland committed
187
188
189
190
191
192
  }
  return product;
}

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

/// Returns the identity matrix.
///
208
209
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT vtkm::Matrix<T, Size, Size> MatrixIdentity()
Kenneth Moreland's avatar
Kenneth Moreland committed
210
{
211
  vtkm::Matrix<T, Size, Size> result(T(0));
Kenneth Moreland's avatar
Kenneth Moreland committed
212
213
  for (vtkm::IdComponent index = 0; index < Size; index++)
  {
214
    result(index, index) = T(1.0);
Kenneth Moreland's avatar
Kenneth Moreland committed
215
216
217
218
219
220
  }
  return result;
}

/// Fills the given matrix with the identity matrix.
///
221
222
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT void MatrixIdentity(vtkm::Matrix<T, Size, Size>& matrix)
Kenneth Moreland's avatar
Kenneth Moreland committed
223
{
224
  matrix = vtkm::MatrixIdentity<T, Size>();
Kenneth Moreland's avatar
Kenneth Moreland committed
225
226
227
228
}

/// Returns the transpose of the given matrix.
///
229
230
231
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
232
{
233
  vtkm::Matrix<T, NumCols, NumRows> result;
Kenneth Moreland's avatar
Kenneth Moreland committed
234
235
236
  for (vtkm::IdComponent index = 0; index < NumRows; index++)
  {
    vtkm::MatrixSetColumn(result, index, vtkm::MatrixGetRow(matrix, index));
237
238
239
240
241
242
243
244
245
#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
246
247
248
249
  }
  return result;
}

250
251
namespace detail
{
Kenneth Moreland's avatar
Kenneth Moreland committed
252
253

// Used with MatrixLUPFactor.
254
255
256
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT void MatrixLUPFactorFindPivot(vtkm::Matrix<T, Size, Size>& A,
                                             vtkm::Vec<vtkm::IdComponent, Size>& permutation,
257
258
                                             vtkm::IdComponent topCornerIndex,
                                             T& inversionParity,
259
                                             bool& valid)
Kenneth Moreland's avatar
Kenneth Moreland committed
260
261
262
{
  vtkm::IdComponent maxRowIndex = topCornerIndex;
  T maxValue = vtkm::Abs(A(maxRowIndex, topCornerIndex));
263
  for (vtkm::IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++)
Kenneth Moreland's avatar
Kenneth Moreland committed
264
265
266
267
268
269
270
271
272
  {
    T compareValue = vtkm::Abs(A(rowIndex, topCornerIndex));
    if (maxValue < compareValue)
    {
      maxValue = compareValue;
      maxRowIndex = rowIndex;
    }
  }

273
274
275
276
  if (maxValue < vtkm::Epsilon<T>())
  {
    valid = false;
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
277
278
279
280

  if (maxRowIndex != topCornerIndex)
  {
    // Swap rows in matrix.
281
282
    vtkm::Vec<T, Size> maxRow = vtkm::MatrixGetRow(A, maxRowIndex);
    vtkm::MatrixSetRow(A, maxRowIndex, vtkm::MatrixGetRow(A, topCornerIndex));
Kenneth Moreland's avatar
Kenneth Moreland committed
283
284
285
286
287
288
289
290
291
292
293
294
295
    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
296
297
298
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
299
300
{
  // Compute values for upper triangle on row topCornerIndex
301
  if (A(topCornerIndex, topCornerIndex) == 0)
Kenneth Moreland's avatar
Kenneth Moreland committed
302
  {
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
    for (vtkm::IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
    {
      A(topCornerIndex, colIndex) = std::numeric_limits<T>::quiet_NaN();
    }
  }
  else
  {
    // Let's make the reciprocal approximation here.
    // Doesn't make things much fast for small 'Size',
    // but definitely improves performance as 'Size' gets large.
    T rAdiag = 1 / A(topCornerIndex, topCornerIndex);
    for (vtkm::IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
    {
      A(topCornerIndex, colIndex) *= rAdiag;
    }
Kenneth Moreland's avatar
Kenneth Moreland committed
318
319
320
  }

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

/// Performs an LUP-factorization on the given matrix using Crout's method. The
331
/// LU-factorization takes a matrix A  and decomposes it into a lower triangular
Kenneth Moreland's avatar
Kenneth Moreland committed
332
333
/// 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
334
/// decomposition always possible so long as A is not singular. In addition to
Kenneth Moreland's avatar
Kenneth Moreland committed
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
/// 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.
///
360
361
362
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT void MatrixLUPFactor(vtkm::Matrix<T, Size, Size>& A,
                                    vtkm::Vec<vtkm::IdComponent, Size>& permutation,
363
364
                                    T& inversionParity,
                                    bool& valid)
Kenneth Moreland's avatar
Kenneth Moreland committed
365
366
367
368
369
370
{
  // Initialize permutation.
  for (vtkm::IdComponent index = 0; index < Size; index++)
  {
    permutation[index] = index;
  }
371
  inversionParity = T(1);
Kenneth Moreland's avatar
Kenneth Moreland committed
372
373
374
375
376
377
378
379
380
381
382
383
384
  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.
///
385
386
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT vtkm::Vec<T, Size> MatrixLUPSolve(
387
388
  const vtkm::Matrix<T, Size, Size>& LU,
  const vtkm::Vec<vtkm::IdComponent, Size>& permutation,
389
  const vtkm::Vec<T, Size>& b)
Kenneth Moreland's avatar
Kenneth Moreland committed
390
391
392
393
394
395
396
{
  // 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.
397
  vtkm::Vec<T, Size> y;
Kenneth Moreland's avatar
Kenneth Moreland committed
398
  for (vtkm::IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
399
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
400
401
402
    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++)
403
404
    {
      y[rowIndex] -= LU(rowIndex, colIndex) * y[colIndex];
Kenneth Moreland's avatar
Kenneth Moreland committed
405
    }
406
407
408
409
410
411
412
413
    if (LU(rowIndex, rowIndex) == 0)
    {
      y[rowIndex] = std::numeric_limits<T>::quiet_NaN();
    }
    else
    {
      y[rowIndex] /= LU(rowIndex, rowIndex);
    }
414
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
415
416

  // Now that we have y, we can easily solve Ux = y for x.
417
418
  vtkm::Vec<T, Size> x;
  for (vtkm::IdComponent rowIndex = Size - 1; rowIndex >= 0; rowIndex--)
Kenneth Moreland's avatar
Kenneth Moreland committed
419
420
421
422
  {
    // Recall that U is stored in the upper triangle of LU with the diagonal
    // implicitly all 1's.
    x[rowIndex] = y[rowIndex];
423
    for (vtkm::IdComponent colIndex = rowIndex + 1; colIndex < Size; colIndex++)
Kenneth Moreland's avatar
Kenneth Moreland committed
424
    {
425
      x[rowIndex] -= LU(rowIndex, colIndex) * x[colIndex];
Kenneth Moreland's avatar
Kenneth Moreland committed
426
427
428
429
430
431
432
433
434
435
436
    }
  }

  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.
///
437
438
template <typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT vtkm::Vec<T, Size> SolveLinearSystem(const vtkm::Matrix<T, Size, Size>& A,
439
440
                                                    const vtkm::Vec<T, Size>& b,
                                                    bool& valid)
Kenneth Moreland's avatar
Kenneth Moreland committed
441
442
{
  // First, we will make an LUP-factorization to help us.
443
444
445
  vtkm::Matrix<T, Size, Size> LU = A;
  vtkm::Vec<vtkm::IdComponent, Size> permutation;
  T inversionParity; // Unused.
Kenneth Moreland's avatar
Kenneth Moreland committed
446
447
448
449
450
451
452
453
454
  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.
///
455
456
457
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
458
459
{
  // First, we will make an LUP-factorization to help us.
460
461
462
  vtkm::Matrix<T, Size, Size> LU = A;
  vtkm::Vec<vtkm::IdComponent, Size> permutation;
  T inversionParity; // Unused
Kenneth Moreland's avatar
Kenneth Moreland committed
463
464
465
466
467
  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.
468
469
  vtkm::Matrix<T, Size, Size> invA;
  vtkm::Vec<T, Size> ICol(T(0));
Kenneth Moreland's avatar
Kenneth Moreland committed
470
  for (vtkm::IdComponent colIndex = 0; colIndex < Size; colIndex++)
471
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
472
    ICol[colIndex] = 1;
473
    vtkm::Vec<T, Size> invACol = detail::MatrixLUPSolve(LU, permutation, ICol);
Kenneth Moreland's avatar
Kenneth Moreland committed
474
475
    ICol[colIndex] = 0;
    vtkm::MatrixSetColumn(invA, colIndex, invACol);
476
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
477
478
479
480
481
  return invA;
}

/// Compute the determinant of a matrix.
///
482
483
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
484
485
{
  // First, we will make an LUP-factorization to help us.
486
487
  vtkm::Matrix<T, Size, Size> LU = A;
  vtkm::Vec<vtkm::IdComponent, Size> permutation;
Kenneth Moreland's avatar
Kenneth Moreland committed
488
489
490
491
492
493
  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.
494
495
496
497
  if (!valid)
  {
    return 0;
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
498
499
500
501
502
503
504

  // 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++)
505
506
507
  {
    product *= LU(index, index);
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
508
509
510
511
512
  return product;
}

// Specializations for common small determinants.

513
514
template <typename T>
VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 1, 1>& A)
Kenneth Moreland's avatar
Kenneth Moreland committed
515
{
516
  return A(0, 0);
Kenneth Moreland's avatar
Kenneth Moreland committed
517
518
}

519
520
template <typename T>
VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 2, 2>& A)
Kenneth Moreland's avatar
Kenneth Moreland committed
521
{
522
  return vtkm::DifferenceOfProducts(A(0, 0), A(1, 1), A(1, 0), A(0, 1));
Kenneth Moreland's avatar
Kenneth Moreland committed
523
524
}

525
526
template <typename T>
VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix<T, 3, 3>& A)
Kenneth Moreland's avatar
Kenneth Moreland committed
527
{
528
529
  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
530
531
532
533
534
535
536
537
}

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

/// Tag used to identify 2 dimensional types (matrices). A TypeTraits class
/// will typedef this class to DimensionalityTag.
///
538
539
540
struct TypeTraitsMatrixTag
{
};
Kenneth Moreland's avatar
Kenneth Moreland committed
541

542
543
544
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
struct TypeTraits<vtkm::Matrix<T, NumRow, NumCol>>
{
545
546
  using NumericTag = typename TypeTraits<T>::NumericTag;
  using DimensionalityTag = vtkm::TypeTraitsMatrixTag;
547
548
549
550
551
552

  VTKM_EXEC_CONT
  static vtkm::Matrix<T, NumRow, NumCol> ZeroInitialization()
  {
    return vtkm::Matrix<T, NumRow, NumCol>(vtkm::TypeTraits<T>::ZeroInitialization());
  }
Kenneth Moreland's avatar
Kenneth Moreland committed
553
554
555
556
};

/// A matrix has vector traits to implement component-wise operations.
///
557
558
559
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
struct VecTraits<vtkm::Matrix<T, NumRow, NumCol>>
{
Kenneth Moreland's avatar
Kenneth Moreland committed
560
private:
561
  using MatrixType = vtkm::Matrix<T, NumRow, NumCol>;
562

Kenneth Moreland's avatar
Kenneth Moreland committed
563
public:
564
  using ComponentType = T;
565
  using BaseComponentType = typename vtkm::VecTraits<T>::BaseComponentType;
566
  static constexpr vtkm::IdComponent NUM_COMPONENTS = NumRow * NumCol;
567
568
  using HasMultipleComponents = vtkm::VecTraitsTagMultipleComponents;
  using IsSizeStatic = vtkm::VecTraitsTagSizeStatic;
Kenneth Moreland's avatar
Kenneth Moreland committed
569

570
  VTKM_EXEC_CONT
571
  static vtkm::IdComponent GetNumberOfComponents(const MatrixType&) { return NUM_COMPONENTS; }
Kenneth Moreland's avatar
Kenneth Moreland committed
572

573
  VTKM_EXEC_CONT
574
575
  static const ComponentType& GetComponent(const MatrixType& matrix, vtkm::IdComponent component)
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
576
577
    vtkm::IdComponent colIndex = component % NumCol;
    vtkm::IdComponent rowIndex = component / NumCol;
578
    return matrix(rowIndex, colIndex);
Kenneth Moreland's avatar
Kenneth Moreland committed
579
  }
580
  VTKM_EXEC_CONT
581
582
  static ComponentType& GetComponent(MatrixType& matrix, vtkm::IdComponent component)
  {
Kenneth Moreland's avatar
Kenneth Moreland committed
583
584
    vtkm::IdComponent colIndex = component % NumCol;
    vtkm::IdComponent rowIndex = component / NumCol;
585
    return matrix(rowIndex, colIndex);
Kenneth Moreland's avatar
Kenneth Moreland committed
586
  }
587
  VTKM_EXEC_CONT
588
  static void SetComponent(MatrixType& matrix, vtkm::IdComponent component, T value)
Kenneth Moreland's avatar
Kenneth Moreland committed
589
590
591
  {
    GetComponent(matrix, component) = value;
  }
592
593
594
595
596
597
598
599
600

  template <typename NewComponentType>
  using ReplaceComponentType = vtkm::Matrix<NewComponentType, NumRow, NumCol>;

  template <typename NewComponentType>
  using ReplaceBaseComponentType =
    vtkm::Matrix<typename vtkm::VecTraits<T>::template ReplaceBaseComponentType<NewComponentType>,
                 NumRow,
                 NumCol>;
Kenneth Moreland's avatar
Kenneth Moreland committed
601
602
603
604
605
};

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

606
607
608
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
609
610
611
612
613
{
  for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
  {
    for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
    {
614
615
      if (a(rowIndex, colIndex) != b(rowIndex, colIndex))
        return false;
Kenneth Moreland's avatar
Kenneth Moreland committed
616
617
618
619
    }
  }
  return true;
}
620
621
622
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
623
624
625
626
{
  return !(a == b);
}

627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
/// Helper function for printing out matricies during testing
///
template <typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_CONT std::ostream& operator<<(std::ostream& stream, const vtkm::Matrix<T, NumRow, NumCol>& mat)
{
  stream << std::endl;
  for (vtkm::IdComponent row = 0; row < NumRow; ++row)
  {
    stream << "| ";
    for (vtkm::IdComponent col = 0; col < NumCol; ++col)
    {
      stream << mat(row, col) << "\t";
    }
    stream << "|" << std::endl;
  }
  return stream;
}
} // namespace vtkm

Kenneth Moreland's avatar
Kenneth Moreland committed
646
#endif //vtk_m_Matrix_h