vtk-m2/vtkm/Matrix.h

//=============================================================================
//
//  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.
//
//  Copyright 2015 Sandia Corporation.
//  Copyright 2015 UT-Battelle, LLC.
//  Copyright 2015 Los Alamos National Security.
//
//  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
//  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

#include <vtkm/Math.h>
#include <vtkm/Types.h>
#include <vtkm/TypeTraits.h>
#include <vtkm/VecTraits.h>

namespace vtkm {

/// \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.
///
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
class Matrix {
public:
  typedef T ComponentType;
  static const vtkm::IdComponent NUM_ROWS = NumRow;
  static const vtkm::IdComponent NUM_COLUMNS = NumCol;

  VTKM_EXEC_CONT_EXPORT
  Matrix() { }

  VTKM_EXEC_CONT_EXPORT
  explicit Matrix(const ComponentType &value)
    : Components(vtkm::Vec<ComponentType, NUM_COLUMNS>(value)) { }

  /// Brackets are used to reference a matrix like a 2D array (i.e.
  /// matrix[row][column]).
  ///
  VTKM_EXEC_CONT_EXPORT
  const vtkm::Vec<ComponentType, NUM_COLUMNS> &
  operator[](vtkm::IdComponent rowIndex) const {
    return this->Components[rowIndex];
  }

  /// Brackets are used to referens a matrix like a 2D array i.e.
  /// matrix[row][column].
  ///
  VTKM_EXEC_CONT_EXPORT
  vtkm::Vec<ComponentType, NUM_COLUMNS> &
  operator[](vtkm::IdComponent rowIndex) {
    return this->Components[rowIndex];
  }

  /// Parentheses are used to reference a matrix using mathematical tuple
  /// notation i.e. matrix(row,column).
  ///
  VTKM_EXEC_CONT_EXPORT
  const ComponentType &
  operator()(vtkm::IdComponent rowIndex, vtkm::IdComponent colIndex) const {
    return (*this)[rowIndex][colIndex];
  }

  /// Parentheses are used to reference a matrix using mathematical tuple
  /// notation i.e. matrix(row,column).
  ///
  VTKM_EXEC_CONT_EXPORT
  ComponentType &
  operator()(vtkm::IdComponent rowIndex, vtkm::IdComponent colIndex) {
    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.
///
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT_EXPORT
const vtkm::Vec<T, NumCol> &MatrixGetRow(
    const vtkm::Matrix<T,NumRow,NumCol> &matrix, vtkm::IdComponent rowIndex)
{
  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.
///
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<T, NumRow> MatrixGetColumn(
    const vtkm::Matrix<T,NumRow,NumCol> &matrix, vtkm::IdComponent columnIndex)
{
  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.
///
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT_EXPORT
void MatrixSetRow(vtkm::Matrix<T,NumRow,NumCol> &matrix,
                  vtkm::IdComponent rowIndex,
                  const vtkm::Vec<T,NumCol> &rowValues)
{
  matrix[rowIndex] = rowValues;
}

/// Convenience function for setting a column of a matrix.
///
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT_EXPORT
void MatrixSetColumn(vtkm::Matrix<T,NumRow,NumCol> &matrix,
                     vtkm::IdComponent columnIndex,
                     const vtkm::Vec<T,NumRow> &columnValues)
{
  for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
  {
    matrix(rowIndex, columnIndex) = columnValues[rowIndex];
  }
}

/// Standard matrix multiplication.
///
template<typename T,
         vtkm::IdComponent NumRow,
         vtkm::IdComponent NumCol,
         vtkm::IdComponent NumInternal>
VTKM_EXEC_CONT_EXPORT
vtkm::Matrix<T,NumRow,NumCol> MatrixMultiply(
    const vtkm::Matrix<T,NumRow,NumInternal> &leftFactor,
    const vtkm::Matrix<T,NumInternal,NumCol> &rightFactor)
{
  vtkm::Matrix<T,NumRow,NumCol> result;
  for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
  {
    for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
    {
      T sum = leftFactor(rowIndex, 0) * rightFactor(0, colIndex);
      for (vtkm::IdComponent internalIndex = 1;
           internalIndex < NumInternal;
           internalIndex++)
      {
        sum += leftFactor(rowIndex, internalIndex)
            * rightFactor(internalIndex, colIndex);
      }
      result(rowIndex, colIndex) = sum;
    }
  }
  return result;
}

/// Standard matrix-vector multiplication.
///
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<T,NumRow> MatrixMultiply(
    const vtkm::Matrix<T,NumRow,NumCol> &leftFactor,
    const vtkm::Vec<T,NumCol> &rightFactor)
{
  vtkm::Vec<T,NumRow> product;
  for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
  {
    product[rowIndex] =
        vtkm::dot(vtkm::MatrixGetRow(leftFactor,rowIndex), rightFactor);
  }
  return product;
}

/// Standard vector-matrix multiplication
///
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<T,NumCol> MatrixMultiply(
    const vtkm::Vec<T,NumRow> &leftFactor,
    const vtkm::Matrix<T,NumRow,NumCol> &rightFactor)
{
  vtkm::Vec<T,NumCol> product;
  for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
  {
    product[colIndex] =
        vtkm::dot(leftFactor,
                  vtkm::MatrixGetColumn(rightFactor, colIndex));
  }
  return product;
}

/// Returns the identity matrix.
///
template<typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT_EXPORT
vtkm::Matrix<T,Size,Size> MatrixIdentity()
{
  vtkm::Matrix<T,Size,Size> result(T(0));
  for (vtkm::IdComponent index = 0; index < Size; index++)
  {
    result(index,index) = T(1.0);
  }
  return result;
}

/// Fills the given matrix with the identity matrix.
///
template<typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT_EXPORT
void MatrixIdentity(vtkm::Matrix<T,Size,Size> &matrix)
{
  matrix = vtkm::MatrixIdentity<T,Size>();
}

/// Returns the transpose of the given matrix.
///
template<typename T, vtkm::IdComponent NumRows, vtkm::IdComponent NumCols>
VTKM_EXEC_CONT_EXPORT
vtkm::Matrix<T,NumCols,NumRows> MatrixTranspose(
    const vtkm::Matrix<T,NumRows,NumCols> &matrix)
{
  vtkm::Matrix<T,NumCols,NumRows> result;
  for (vtkm::IdComponent index = 0; index < NumRows; index++)
  {
    vtkm::MatrixSetColumn(result, index, vtkm::MatrixGetRow(matrix, index));
  }
  return result;
}

namespace detail {

// Used with MatrixLUPFactor.
template<typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT_EXPORT
void MatrixLUPFactorFindPivot(vtkm::Matrix<T,Size,Size> &A,
                              vtkm::Vec<vtkm::IdComponent,Size> &permutation,
                              vtkm::IdComponent topCornerIndex,
                              T &inversionParity,
                              bool &valid)
{
  vtkm::IdComponent maxRowIndex = topCornerIndex;
  T maxValue = vtkm::Abs(A(maxRowIndex, topCornerIndex));
  for (vtkm::IdComponent rowIndex = topCornerIndex + 1;
       rowIndex < Size;
       rowIndex++)
  {
    T compareValue = vtkm::Abs(A(rowIndex, topCornerIndex));
    if (maxValue < compareValue)
    {
      maxValue = compareValue;
      maxRowIndex = rowIndex;
    }
  }

  if (maxValue < vtkm::Epsilon<T>()) { valid = false; }

  if (maxRowIndex != topCornerIndex)
  {
    // Swap rows in matrix.
    vtkm::Vec<T,Size> maxRow = vtkm::MatrixGetRow(A, maxRowIndex);
    vtkm::MatrixSetRow(A,
                       maxRowIndex,
                       vtkm::MatrixGetRow(A,topCornerIndex));
    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
template<typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT_EXPORT
void MatrixLUPFactorFindUpperTriangleElements(vtkm::Matrix<T,Size,Size> &A,
                                              vtkm::IdComponent topCornerIndex)
{
  // Compute values for upper triangle on row topCornerIndex
  for (vtkm::IdComponent colIndex = topCornerIndex+1;
       colIndex < Size;
       colIndex++)
  {
    A(topCornerIndex,colIndex) /= A(topCornerIndex,topCornerIndex);
  }

  // Update the rest of the matrix for calculations on subsequent rows
  for (vtkm::IdComponent rowIndex = topCornerIndex+1;
       rowIndex < Size;
       rowIndex++)
  {
    for (vtkm::IdComponent colIndex = topCornerIndex+1;
         colIndex < Size;
         colIndex++)
    {
      A(rowIndex,colIndex) -=
          A(rowIndex,topCornerIndex)*A(topCornerIndex,colIndex);
    }
  }
}

/// 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
/// decomposition always posible so long as A is not singular. In addition to
/// 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.
///
template<typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT_EXPORT
void MatrixLUPFactor(vtkm::Matrix<T,Size,Size> &A,
                     vtkm::Vec<vtkm::IdComponent,Size> &permutation,
                     T &inversionParity,
                     bool &valid)
{
  // Initialize permutation.
  for (vtkm::IdComponent index = 0; index < Size; index++)
  {
    permutation[index] = index;
  }
  inversionParity = 1;
  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.
///
template<typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<T,Size>
MatrixLUPSolve(const vtkm::Matrix<T,Size,Size> &LU,
                    const vtkm::Vec<vtkm::IdComponent,Size> &permutation,
                    const vtkm::Vec<T,Size> &b)
{
  // 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.
  vtkm::Vec<T,Size> y;
  for (vtkm::IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
    {
    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++)
      {
      y[rowIndex] -= LU(rowIndex,colIndex)*y[colIndex];
      }
    y[rowIndex] /= LU(rowIndex,rowIndex);
    }

  // Now that we have y, we can easily solve Ux = y for x.
  vtkm::Vec<T,Size> x;
  for (vtkm::IdComponent rowIndex = Size-1; rowIndex >= 0; rowIndex--)
  {
    // Recall that U is stored in the upper triangle of LU with the diagonal
    // implicitly all 1's.
    x[rowIndex] = y[rowIndex];
    for (vtkm::IdComponent colIndex = rowIndex+1; colIndex < Size; colIndex++)
    {
      x[rowIndex] -= LU(rowIndex,colIndex)*x[colIndex];
    }
  }

  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.
///
template<typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<T,Size> SolveLinearSystem(const vtkm::Matrix<T,Size,Size> &A,
                                    const vtkm::Vec<T,Size> &b,
                                    bool &valid)
{
  // First, we will make an LUP-factorization to help us.
  vtkm::Matrix<T,Size,Size> LU = A;
  vtkm::Vec<vtkm::IdComponent,Size> permutation;
  T inversionParity;  // Unused.
  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.
///
template<typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT_EXPORT
vtkm::Matrix<T,Size,Size> MatrixInverse(const vtkm::Matrix<T,Size,Size> &A,
                                        bool &valid)
{
  // First, we will make an LUP-factorization to help us.
  vtkm::Matrix<T,Size,Size> LU = A;
  vtkm::Vec<vtkm::IdComponent,Size> permutation;
  T inversionParity;  // Unused
  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.
  vtkm::Matrix<T,Size,Size> invA;
  vtkm::Vec<T,Size> ICol(T(0));
  for (vtkm::IdComponent colIndex = 0; colIndex < Size; colIndex++)
    {
    ICol[colIndex] = 1;
    vtkm::Vec<T,Size> invACol = detail::MatrixLUPSolve(LU, permutation, ICol);
    ICol[colIndex] = 0;
    vtkm::MatrixSetColumn(invA, colIndex, invACol);
    }
  return invA;
}

/// Compute the determinant of a matrix.
///
template<typename T, vtkm::IdComponent Size>
VTKM_EXEC_CONT_EXPORT
T MatrixDeterminant(const vtkm::Matrix<T,Size,Size> &A)
{
  // First, we will make an LUP-factorization to help us.
  vtkm::Matrix<T,Size,Size> LU = A;
  vtkm::Vec<vtkm::IdComponent,Size> permutation;
  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.
  if (!valid) { return 0; }

  // 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++)
    {
    product *= LU(index,index);
    }
  return product;
}

// Specializations for common small determinants.

template<typename T>
VTKM_EXEC_CONT_EXPORT
T MatrixDeterminant(const vtkm::Matrix<T,1,1> &A)
{
  return A(0,0);
}

template<typename T>
VTKM_EXEC_CONT_EXPORT
T MatrixDeterminant(const vtkm::Matrix<T,2,2> &A)
{
  return A(0,0)*A(1,1) - A(1,0)*A(0,1);
}

template<typename T>
VTKM_EXEC_CONT_EXPORT
T MatrixDeterminant(const vtkm::Matrix<T,3,3> &A)
{
  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);
}

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

/// Tag used to identify 2 dimensional types (matrices). A TypeTraits class
/// will typedef this class to DimensionalityTag.
///
struct TypeTraitsMatrixTag {};

template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
struct TypeTraits<vtkm::Matrix<T, NumRow, NumCol> > {
  typedef typename TypeTraits<T>::NumericTag NumericTag;
  typedef TypeTraitsMatrixTag DimensionalityTag;
};

/// A matrix has vector traits to implement component-wise operations.
///
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
struct VecTraits<vtkm::Matrix<T, NumRow, NumCol> > {
private:
  typedef vtkm::Matrix<T, NumRow, NumCol> MatrixType;
public:
  typedef T ComponentType;
  static const vtkm::IdComponent NUM_COMPONENTS = NumRow*NumCol;
  typedef vtkm::VecTraitsTagMultipleComponents HasMultipleComponents;

  VTKM_EXEC_CONT_EXPORT
  static const ComponentType &GetComponent(const MatrixType &matrix,
                                           vtkm::IdComponent component) {
    vtkm::IdComponent colIndex = component % NumCol;
    vtkm::IdComponent rowIndex = component / NumCol;
    return matrix(rowIndex,colIndex);
  }
  VTKM_EXEC_CONT_EXPORT
  static ComponentType &GetComponent(MatrixType &matrix,
                                     vtkm::IdComponent component) {
    vtkm::IdComponent colIndex = component % NumCol;
    vtkm::IdComponent rowIndex = component / NumCol;
    return matrix(rowIndex,colIndex);
  }
  VTKM_EXEC_CONT_EXPORT
  static void SetComponent(MatrixType &matrix,
                           vtkm::IdComponent component,
                           T value)
  {
    GetComponent(matrix, component) = value;
  }
};

} // namespace vtkm

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

template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT_EXPORT
bool operator==(const vtkm::Matrix<T,NumRow,NumCol> &a,
                const vtkm::Matrix<T,NumRow,NumCol> &b)
{
  for (vtkm::IdComponent colIndex = 0; colIndex < NumCol; colIndex++)
  {
    for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
    {
      if (a(rowIndex, colIndex) != b(rowIndex, colIndex)) return false;
    }
  }
  return true;
}
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
VTKM_EXEC_CONT_EXPORT
bool operator!=(const vtkm::Matrix<T,NumRow,NumCol> &a,
                const vtkm::Matrix<T,NumRow,NumCol> &b)
{
  return !(a == b);
}

#endif //vtk_m_Matrix_h