Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add Matrix class.

Browse Source 
The matrix class is for thread-local matrix computations. It is intended
to store things like tensors or geometric transforms.
This commit is contained in:
Kenneth Moreland 2015-07-02 11:20:21 -06:00
parent c52c83bdd0
commit abcc6da52f
6 changed files with 1221 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
vtkm/CMakeLists.txt
View File

@ -25,6 +25,7 @@ set(headers
Extent.h
ListTag.h
Math.h
Matrix.h
Pair.h
RegularConnectivity.h
RegularStructure.h

12
vtkm/Math.h
View File

@ -1505,6 +1505,18 @@ VTKM_EXEC_CONT_EXPORT vtkm::Float64 NegativeInfinity64() {
return vtkm::NegativeInfinity<vtkm::Float64>();
}
/// Returns the difference between 1 and the least value greater than 1
/// that is representable.
///
VTKM_EXEC_CONT_EXPORT vtkm::Float32 Epsilon32()
{
return vtkm::Epsilon<vtkm::Float32>();
}
VTKM_EXEC_CONT_EXPORT vtkm::Float64 Epsilon64()
{
return vtkm::Epsilon<vtkm::Float64>();
}
//-----------------------------------------------------------------------------
/// Returns true if \p x is not a number.
///

12
vtkm/Math.h.in
View File

@ -631,6 +631,18 @@ VTKM_EXEC_CONT_EXPORT vtkm::Float64 NegativeInfinity64() {
return vtkm::NegativeInfinity<vtkm::Float64>();
}
/// Returns the difference between 1 and the least value greater than 1
/// that is representable.
///
VTKM_EXEC_CONT_EXPORT vtkm::Float32 Epsilon32()
{
return vtkm::Epsilon<vtkm::Float32>();
}
VTKM_EXEC_CONT_EXPORT vtkm::Float64 Epsilon64()
{
return vtkm::Epsilon<vtkm::Float64>();
}
//-----------------------------------------------------------------------------
/// Returns true if \p x is not a number.
///

602
vtkm/Matrix.h Normal file
View File

@ -0,0 +1,602 @@
//=============================================================================
//
// 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

1
vtkm/testing/CMakeLists.txt
View File

@ -31,6 +31,7 @@ set(unit_tests
UnitTestExtent.cxx
UnitTestListTag.cxx
UnitTestMath.cxx
UnitTestMatrix.cxx
UnitTestPair.cxx
UnitTestTesting.cxx
UnitTestTypeListTag.cxx

593
vtkm/testing/UnitTestMatrix.cxx Normal file
View File

@ -0,0 +1,593 @@
//=============================================================================
//
// 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.
//
//=============================================================================
#include <vtkm/Matrix.h>
#include <vtkm/VecTraits.h>
#include <vtkm/testing/Testing.h>
// If more tests need a value for Matrix, we can move this to Testing.h
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
vtkm::Matrix<T,NumRow,NumCol>
TestValue(vtkm::Id index, const vtkm::Matrix<T,NumRow,NumCol> &)
{
vtkm::Matrix<T,NumRow,NumCol> value;
for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++)
{
typedef vtkm::Vec<T,NumCol> RowType;
RowType row =
TestValue(index, RowType())
+ RowType(static_cast<typename RowType::ComponentType>(10*rowIndex));
vtkm::MatrixSetRow(value, rowIndex, row);
}
return value;
}
namespace {
#define FOR_ROW_COL(matrix) \
for(vtkm::IdComponent row=0; row < (matrix).NUM_ROWS; row++) \
for(vtkm::IdComponent col=0; col < (matrix).NUM_COLUMNS; col++)
template<typename T, vtkm::IdComponent NumRow, vtkm::IdComponent NumCol>
struct MatrixTest
{
static const vtkm::IdComponent NUM_ROWS = NumRow;
static const vtkm::IdComponent NUM_COLS = NumCol;
typedef vtkm::Matrix<T,NUM_ROWS,NUM_COLS> MatrixType;
static void BasicCreation()
{
std::cout << "Basic creation." << std::endl;
MatrixType matrix(5);
FOR_ROW_COL(matrix)
{
VTKM_TEST_ASSERT(matrix(row,col) == static_cast<T>(5),
"Constant set incorrect.");
}
}
static void BasicAccessors()
{
std::cout << "Basic accessors." << std::endl;
MatrixType matrix;
MatrixType value = TestValue(0, MatrixType());
FOR_ROW_COL(matrix)
{
matrix[row][col] = value(row,col)*2;
}
FOR_ROW_COL(matrix)
{
VTKM_TEST_ASSERT(matrix(row,col) == value(row,col)*2,
"Bad set or retreive.");
const MatrixType const_matrix = matrix;
VTKM_TEST_ASSERT(const_matrix(row,col) == value(row,col)*2,
"Bad set or retreive.");
}
FOR_ROW_COL(matrix)
{
matrix(row,col) = value(row,col);
}
const MatrixType const_matrix = matrix;
FOR_ROW_COL(matrix)
{
VTKM_TEST_ASSERT(matrix[row][col] == value(row,col),
"Bad set or retreive.");
VTKM_TEST_ASSERT(const_matrix[row][col] == value(row,col),
"Bad set or retreive.");
}
VTKM_TEST_ASSERT(matrix == const_matrix,
"Equal test operator not working.");
VTKM_TEST_ASSERT(!(matrix != const_matrix),
"Not-Equal test operator not working.");
VTKM_TEST_ASSERT(test_equal(matrix, const_matrix),
"Vector-based equal test not working.");
}
static void RowColAccessors()
{
typedef vtkm::Vec<T,NUM_ROWS> ColumnType;
typedef vtkm::Vec<T,NUM_COLS> RowType;
const MatrixType const_matrix = TestValue(0, MatrixType());
MatrixType matrix;
std::cout << "Access by row or column" << std::endl;
FOR_ROW_COL(matrix)
{
RowType rowvec = vtkm::MatrixGetRow(const_matrix, row);
VTKM_TEST_ASSERT(rowvec[col] == const_matrix(row,col), "Bad get row.");
ColumnType columnvec = vtkm::MatrixGetColumn(const_matrix, col);
VTKM_TEST_ASSERT(columnvec[row] == const_matrix(row,col), "Bad get col.");
}
std::cout << "Set a row." << std::endl;
for (vtkm::IdComponent row = 0; row < NUM_ROWS; row++)
{
RowType rowvec = vtkm::MatrixGetRow(const_matrix, NUM_ROWS-row-1);
vtkm::MatrixSetRow(matrix, row, rowvec);
}
FOR_ROW_COL(matrix)
{
VTKM_TEST_ASSERT(matrix(NUM_ROWS-row-1,col) == const_matrix(row,col),
"Rows not set right.");
}
std::cout << "Set a column." << std::endl;
for (vtkm::IdComponent col = 0; col < NUM_COLS; col++)
{
ColumnType colvec = vtkm::MatrixGetColumn(const_matrix, NUM_COLS-col-1);
vtkm::MatrixSetColumn(matrix, col, colvec);
}
FOR_ROW_COL(matrix)
{
VTKM_TEST_ASSERT(matrix(row,NUM_COLS-col-1) == const_matrix(row,col),
"Columns not set right.");
}
}
static void Multiply()
{
std::cout << "Try multiply." << std::endl;
const MatrixType leftFactor = TestValue(0, MatrixType());
vtkm::Matrix<T,NUM_COLS,4> rightFactor =
TestValue(1, vtkm::Matrix<T,NUM_COLS,4>());
vtkm::Matrix<T,NUM_ROWS,4> product
= vtkm::MatrixMultiply(leftFactor, rightFactor);
FOR_ROW_COL(product)
{
vtkm::Vec<T,NUM_COLS> leftVector = vtkm::MatrixGetRow(leftFactor, row);
vtkm::Vec<T,NUM_COLS> rightVector=vtkm::MatrixGetColumn(rightFactor, col);
VTKM_TEST_ASSERT(test_equal(product(row,col),
vtkm::dot(leftVector,rightVector)),
"Matrix multiple wrong.");
}
std::cout << "Vector multiply." << std::endl;
MatrixType matrixFactor;
vtkm::Vec<T,NUM_ROWS> leftVector(2);
vtkm::Vec<T,NUM_COLS> rightVector;
FOR_ROW_COL(matrixFactor)
{
matrixFactor(row,col) = T(row + 1);
rightVector[col] = T(col + 1);
}
vtkm::Vec<T,NUM_COLS> leftResult =
vtkm::MatrixMultiply(leftVector, matrixFactor);
for (vtkm::IdComponent index = 0; index < NUM_COLS; index++)
{
VTKM_TEST_ASSERT(test_equal(leftResult[index], T(NUM_ROWS*(NUM_ROWS+1))),
"Vector/matrix multiple wrong.");
}
vtkm::Vec<T,NUM_ROWS> rightResult =
vtkm::MatrixMultiply(matrixFactor, rightVector);
for (vtkm::IdComponent index = 0; index < NUM_ROWS; index++)
{
VTKM_TEST_ASSERT(test_equal(rightResult[index],
T(((index+1)*NUM_COLS*(NUM_COLS+1))/2)),
"Matrix/vector multiple wrong.");
}
}
static void Identity()
{
std::cout << "Check identity" << std::endl;
MatrixType originalMatrix = TestValue(0, MatrixType());
vtkm::Matrix<T,NUM_COLS,NUM_COLS> identityMatrix;
vtkm::MatrixIdentity(identityMatrix);
MatrixType multMatrix =
vtkm::MatrixMultiply(originalMatrix, identityMatrix);
VTKM_TEST_ASSERT(test_equal(originalMatrix, multMatrix),
"Identity is not really identity.");
}
static void Transpose()
{
std::cout << "Check transpose" << std::endl;
MatrixType originalMatrix = TestValue(0, MatrixType());
vtkm::Matrix<T,NUM_COLS,NUM_ROWS> transMatrix =
vtkm::MatrixTranspose(originalMatrix);
FOR_ROW_COL(originalMatrix)
{
VTKM_TEST_ASSERT(originalMatrix(row,col) == transMatrix(col,row),
"Transpose wrong.");
}
}
static void Run()
{
std::cout << "-- " << NUM_ROWS << " x " << NUM_COLS << std::endl;
BasicCreation();
BasicAccessors();
RowColAccessors();
Multiply();
Identity();
Transpose();
}
private:
MatrixTest(); // Not implemented
};
template<typename T, int NumRow>
void MatrixTest1()
{
MatrixTest<T,NumRow,1>::Run();
MatrixTest<T,NumRow,2>::Run();
MatrixTest<T,NumRow,3>::Run();
MatrixTest<T,NumRow,4>::Run();
MatrixTest<T,NumRow,5>::Run();
}
template<typename T>
void NonSingularMatrix(vtkm::Matrix<T,1,1> &matrix)
{
matrix(0,0) = 1;
}
template<typename T>
void NonSingularMatrix(vtkm::Matrix<T,2,2> &matrix)
{
matrix(0,0) = -5; matrix(0,1) = 6;
matrix(1,0) = -7; matrix(1,1) = -2;
}
template<typename T>
void NonSingularMatrix(vtkm::Matrix<T,3,3> &matrix)
{
matrix(0,0) = 1; matrix(0,1) = -2; matrix(0,2) = 3;
matrix(1,0) = 6; matrix(1,1) = 7; matrix(1,2) = -1;
matrix(2,0) = -3; matrix(2,1) = 1; matrix(2,2) = 4;
}
template<typename T>
void NonSingularMatrix(vtkm::Matrix<T,4,4> &matrix)
{
matrix(0,0) = 2; matrix(0,1) = 1; matrix(0,2) = 0; matrix(0,3) = 3;
matrix(1,0) = -1; matrix(1,1) = 0; matrix(1,2) = 2; matrix(1,3) = 4;
matrix(2,0) = 4; matrix(2,1) = -2; matrix(2,2) = 7; matrix(2,3) = 0;
matrix(3,0) = -4; matrix(3,1) = 3; matrix(3,2) = 5; matrix(3,3) = 1;
}
template<typename T>
void NonSingularMatrix(vtkm::Matrix<T,5,5> &mat)
{
mat(0,0) = 2; mat(0,1) = 1; mat(0,2) = 3; mat(0,3) = 7; mat(0,4) = 5;
mat(1,0) = 3; mat(1,1) = 8; mat(1,2) = 7; mat(1,3) = 9; mat(1,4) = 8;
mat(2,0) = 3; mat(2,1) = 4; mat(2,2) = 1; mat(2,3) = 6; mat(2,4) = 2;
mat(3,0) = 4; mat(3,1) = 0; mat(3,2) = 2; mat(3,3) = 2; mat(3,4) = 3;
mat(4,0) = 7; mat(4,1) = 9; mat(4,2) = 1; mat(4,3) = 5; mat(4,4) = 4;
}
template<typename T, int Size>
void SingularMatrix(vtkm::Matrix<T,Size,Size> &singularMatrix)
{
FOR_ROW_COL(singularMatrix)
{
singularMatrix(row,col) = row+col;
}
if (Size > 1)
{
vtkm::MatrixSetRow(singularMatrix,
0,
vtkm::MatrixGetRow(singularMatrix, (Size+1)/2));
}
}
// A simple but slow implementation of finding a determinant for comparison
// purposes.
template<typename T>
T RecursiveDeterminant(const vtkm::Matrix<T,1,1> &A)
{
return A(0,0);
}
template<typename T, vtkm::IdComponent Size>
T RecursiveDeterminant(const vtkm::Matrix<T,Size,Size> &A)
{
vtkm::Matrix<T,Size-1,Size-1> cofactorMatrix;
T sum = 0.0;
T sign = 1.0;
for (vtkm::IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
{
// Create the cofactor matrix for entry A(rowIndex,0)
for (vtkm::IdComponent cofactorRowIndex = 0;
cofactorRowIndex < rowIndex;
cofactorRowIndex++)
{
for (vtkm::IdComponent colIndex = 1; colIndex < Size; colIndex++)
{
cofactorMatrix(cofactorRowIndex,colIndex-1) =
A(cofactorRowIndex,colIndex);
}
}
for (vtkm::IdComponent cofactorRowIndex = rowIndex+1;
cofactorRowIndex < Size;
cofactorRowIndex++)
{
for (vtkm::IdComponent colIndex = 1; colIndex < Size; colIndex++)
{
cofactorMatrix(cofactorRowIndex-1,colIndex-1) =
A(cofactorRowIndex,colIndex);
}
}
sum += sign * A(rowIndex,0) * RecursiveDeterminant(cofactorMatrix);
sign = -sign;
}
return sum;
}
template<typename T, vtkm::IdComponent Size>
struct SquareMatrixTest {
static const vtkm::IdComponent SIZE = Size;
typedef vtkm::Matrix<T,Size,Size> MatrixType;
static void CheckMatrixSize()
{
std::cout << "Check reported matrix size." << std::endl;
VTKM_TEST_ASSERT(MatrixType::NUM_ROWS == SIZE, "Matrix has wrong size.");
VTKM_TEST_ASSERT(MatrixType::NUM_COLUMNS == SIZE, "Matrix has wrong size.");
}
static void LUPFactor()
{
std::cout << "Test LUP-factorization" << std::endl;
MatrixType A;
NonSingularMatrix(A);
const MatrixType originalMatrix = A;
vtkm::Vec<vtkm::IdComponent,SIZE> permutationVector;
T inversionParity;
bool valid;
vtkm::detail::MatrixLUPFactor(A,
permutationVector,
inversionParity,
valid);
VTKM_TEST_ASSERT(valid, "Matrix declared singular?");
// Reconstruct L and U matrices from A.
MatrixType L(0);
MatrixType U(0);
FOR_ROW_COL(A)
{
if (row < col)
{
U(row,col) = A(row,col);
}
else //(row >= col)
{
L(row,col) = A(row,col);
if (row == col) { U(row,col) = 1; }
}
}
// Check parity of permutation.
T computedParity = 1.0;
for (int i = 0; i < SIZE; i++)
{
for (int j = i+1; j < SIZE; j++)
{
if (permutationVector[i] > permutationVector[j])
{
computedParity = -computedParity;
}
}
}
VTKM_TEST_ASSERT(inversionParity == computedParity,
"Got bad inversion parity.");
// Reconstruct permutation matrix P.
MatrixType P(0);
for (vtkm::IdComponent index = 0; index < Size; index++)
{
P(index, permutationVector[index]) = 1;
}
// Check that PA = LU is actually correct.
MatrixType permutedMatrix = vtkm::MatrixMultiply(P,originalMatrix);
MatrixType productMatrix = vtkm::MatrixMultiply(L,U);
VTKM_TEST_ASSERT(test_equal(permutedMatrix, productMatrix),
"LUP-factorization gave inconsistent answer.");
// Check that a singular matrix is identified.
MatrixType singularMatrix;
SingularMatrix(singularMatrix);
vtkm::detail::MatrixLUPFactor(singularMatrix,
permutationVector,
inversionParity,
valid);
VTKM_TEST_ASSERT(!valid, "Expected matrix to be declared singular.");
}
static void SolveLinearSystem()
{
std::cout << "Solve a linear system" << std::endl;
MatrixType A;
vtkm::Vec<T,SIZE> b;
NonSingularMatrix(A);
for (vtkm::IdComponent index = 0; index < SIZE; index++)
{
b[index] = index+1;
}
bool valid;
vtkm::Vec<T,SIZE> x = vtkm::SolveLinearSystem(A, b, valid);
VTKM_TEST_ASSERT(valid, "Matrix declared singular?");
// Check result.
vtkm::Vec<T,SIZE> check = vtkm::MatrixMultiply(A,x);
VTKM_TEST_ASSERT(test_equal(b, check),
"Linear solution does not solve equation.");
// Check that a singular matrix is identified.
MatrixType singularMatrix;
SingularMatrix(singularMatrix);
vtkm::SolveLinearSystem(singularMatrix, b, valid);
VTKM_TEST_ASSERT(!valid, "Expected matrix to be declared singular.");
}
static void Invert()
{
std::cout << "Invert a matrix." << std::endl;
MatrixType A;
NonSingularMatrix(A);
bool valid;
vtkm::Matrix<T,SIZE,SIZE> inverse = vtkm::MatrixInverse(A, valid);
VTKM_TEST_ASSERT(valid, "Matrix declared singular?");
// Check result.
vtkm::Matrix<T,SIZE,SIZE> product = vtkm::MatrixMultiply(A, inverse);
VTKM_TEST_ASSERT(test_equal(product, vtkm::MatrixIdentity<T,SIZE>()),
"Matrix inverse did not give identity.");
// Check that a singular matrix is identified.
MatrixType singularMatrix;
SingularMatrix(singularMatrix);
vtkm::MatrixInverse(singularMatrix, valid);
VTKM_TEST_ASSERT(!valid, "Expected matrix to be declared singular.");
}
static void Determinant()
{
std::cout << "Compute a determinant." << std::endl;
MatrixType A;
NonSingularMatrix(A);
T determinant = vtkm::MatrixDeterminant(A);
// Check result.
T determinantCheck = RecursiveDeterminant(A);
VTKM_TEST_ASSERT(test_equal(determinant, determinantCheck),
"Determinant computations do not agree.");
// Check that a singular matrix has a zero determinant.
MatrixType singularMatrix;
SingularMatrix(singularMatrix);
determinant = vtkm::MatrixDeterminant(singularMatrix);
VTKM_TEST_ASSERT(test_equal(determinant, T(0.0)),
"Non-zero determinant for singular matrix.");
}
static void Run()
{
std::cout << "-- " << SIZE << " x " << SIZE << std::endl;
CheckMatrixSize();
LUPFactor();
SolveLinearSystem();
Invert();
Determinant();
}
private:
SquareMatrixTest(); // Not implemented
};
struct MatrixTestFunctor
{
template<typename T>
void operator()(const T&) const {
MatrixTest1<T,1>();
MatrixTest1<T,2>();
MatrixTest1<T,3>();
MatrixTest1<T,4>();
MatrixTest1<T,5>();
}
};
struct SquareMatrixTestFunctor
{
template<typename T>
void operator()(const T&) const {
SquareMatrixTest<T,1>::Run();
SquareMatrixTest<T,2>::Run();
SquareMatrixTest<T,3>::Run();
SquareMatrixTest<T,4>::Run();
SquareMatrixTest<T,5>::Run();
}
};
struct VectorMultFunctor
{
template<class VectorType>
void operator()(const VectorType &) const {
// This is mostly to make sure the compile can convert from Tuples
// to vectors.
const int SIZE = vtkm::VecTraits<VectorType>::NUM_COMPONENTS;
typedef typename vtkm::VecTraits<VectorType>::ComponentType ComponentType;
vtkm::Matrix<ComponentType,SIZE,SIZE> matrix(0);
VectorType inVec;
VectorType outVec;
for (vtkm::IdComponent index = 0; index < SIZE; index++)
{
matrix(index,index) = 1;
inVec[index] = ComponentType(index+1);
}
outVec = vtkm::MatrixMultiply(matrix,inVec);
VTKM_TEST_ASSERT(test_equal(inVec, outVec), "Bad identity multiply.");
outVec = vtkm::MatrixMultiply(inVec,matrix);
VTKM_TEST_ASSERT(test_equal(inVec, outVec), "Bad identity multiply.");
}
};
void TestMatrices()
{
std::cout << "****** Rectangle tests" << std::endl;
vtkm::testing::Testing::TryTypes(MatrixTestFunctor(),
vtkm::TypeListTagScalarAll());
std::cout << "****** Square tests" << std::endl;
vtkm::testing::Testing::TryTypes(SquareMatrixTestFunctor(),
vtkm::TypeListTagFieldScalar());
std::cout << "***** Vector multiply tests" << std::endl;
vtkm::testing::Testing::TryTypes(VectorMultFunctor(),
vtkm::TypeListTagVecAll());
}
} // anonymous namespace
int UnitTestMatrix(int, char *[])
{
return vtkm::testing::Testing::Run(TestMatrices);
}