From abcc6da52f3e1dfd5ce7a3ddd0d0cd3bcd6217fa Mon Sep 17 00:00:00 2001 From: Kenneth Moreland Date: Thu, 2 Jul 2015 11:20:21 -0600 Subject: [PATCH] Add Matrix class. The matrix class is for thread-local matrix computations. It is intended to store things like tensors or geometric transforms. --- vtkm/CMakeLists.txt | 1 + vtkm/Math.h | 12 + vtkm/Math.h.in | 12 + vtkm/Matrix.h | 602 ++++++++++++++++++++++++++++++++ vtkm/testing/CMakeLists.txt | 1 + vtkm/testing/UnitTestMatrix.cxx | 593 +++++++++++++++++++++++++++++++ 6 files changed, 1221 insertions(+) create mode 100644 vtkm/Matrix.h create mode 100644 vtkm/testing/UnitTestMatrix.cxx diff --git a/vtkm/CMakeLists.txt b/vtkm/CMakeLists.txt index efa6da2ad..029a50026 100644 --- a/vtkm/CMakeLists.txt +++ b/vtkm/CMakeLists.txt @@ -25,6 +25,7 @@ set(headers Extent.h ListTag.h Math.h + Matrix.h Pair.h RegularConnectivity.h RegularStructure.h diff --git a/vtkm/Math.h b/vtkm/Math.h index 5c4b25c31..269677a00 100644 --- a/vtkm/Math.h +++ b/vtkm/Math.h @@ -1505,6 +1505,18 @@ VTKM_EXEC_CONT_EXPORT vtkm::Float64 NegativeInfinity64() { return vtkm::NegativeInfinity(); } +/// 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_EXEC_CONT_EXPORT vtkm::Float64 Epsilon64() +{ + return vtkm::Epsilon(); +} + //----------------------------------------------------------------------------- /// Returns true if \p x is not a number. /// diff --git a/vtkm/Math.h.in b/vtkm/Math.h.in index 8028c1d50..7c6eefcd6 100644 --- a/vtkm/Math.h.in +++ b/vtkm/Math.h.in @@ -631,6 +631,18 @@ VTKM_EXEC_CONT_EXPORT vtkm::Float64 NegativeInfinity64() { return vtkm::NegativeInfinity(); } +/// 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_EXEC_CONT_EXPORT vtkm::Float64 Epsilon64() +{ + return vtkm::Epsilon(); +} + //----------------------------------------------------------------------------- /// Returns true if \p x is not a number. /// diff --git a/vtkm/Matrix.h b/vtkm/Matrix.h new file mode 100644 index 000000000..67e2c9b3c --- /dev/null +++ b/vtkm/Matrix.h @@ -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 +#include +#include +#include + +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 +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(value)) { } + + /// Brackets are used to reference a matrix like a 2D array (i.e. + /// matrix[row][column]). + /// + VTKM_EXEC_CONT_EXPORT + const vtkm::Vec & + 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 & + 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, NUM_ROWS> Components; +}; + +/// Returns a tuple containing the given row (indexed from 0) of the given +/// matrix. +/// +template +VTKM_EXEC_CONT_EXPORT +const vtkm::Vec &MatrixGetRow( + const vtkm::Matrix &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 +VTKM_EXEC_CONT_EXPORT +vtkm::Vec MatrixGetColumn( + const vtkm::Matrix &matrix, vtkm::IdComponent columnIndex) +{ + vtkm::Vec 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 +VTKM_EXEC_CONT_EXPORT +void MatrixSetRow(vtkm::Matrix &matrix, + vtkm::IdComponent rowIndex, + const vtkm::Vec &rowValues) +{ + matrix[rowIndex] = rowValues; +} + +/// Convenience function for setting a column of a matrix. +/// +template +VTKM_EXEC_CONT_EXPORT +void MatrixSetColumn(vtkm::Matrix &matrix, + vtkm::IdComponent columnIndex, + const vtkm::Vec &columnValues) +{ + for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++) + { + matrix(rowIndex, columnIndex) = columnValues[rowIndex]; + } +} + +/// Standard matrix multiplication. +/// +template +VTKM_EXEC_CONT_EXPORT +vtkm::Matrix MatrixMultiply( + const vtkm::Matrix &leftFactor, + const vtkm::Matrix &rightFactor) +{ + vtkm::Matrix 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 +VTKM_EXEC_CONT_EXPORT +vtkm::Vec MatrixMultiply( + const vtkm::Matrix &leftFactor, + const vtkm::Vec &rightFactor) +{ + vtkm::Vec 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 +VTKM_EXEC_CONT_EXPORT +vtkm::Vec MatrixMultiply( + const vtkm::Vec &leftFactor, + const vtkm::Matrix &rightFactor) +{ + vtkm::Vec 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 +VTKM_EXEC_CONT_EXPORT +vtkm::Matrix MatrixIdentity() +{ + vtkm::Matrix 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 +VTKM_EXEC_CONT_EXPORT +void MatrixIdentity(vtkm::Matrix &matrix) +{ + matrix = vtkm::MatrixIdentity(); +} + +/// Returns the transpose of the given matrix. +/// +template +VTKM_EXEC_CONT_EXPORT +vtkm::Matrix MatrixTranspose( + const vtkm::Matrix &matrix) +{ + vtkm::Matrix 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 +VTKM_EXEC_CONT_EXPORT +void MatrixLUPFactorFindPivot(vtkm::Matrix &A, + vtkm::Vec &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()) { valid = false; } + + if (maxRowIndex != topCornerIndex) + { + // Swap rows in matrix. + vtkm::Vec 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 +VTKM_EXEC_CONT_EXPORT +void MatrixLUPFactorFindUpperTriangleElements(vtkm::Matrix &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 +VTKM_EXEC_CONT_EXPORT +void MatrixLUPFactor(vtkm::Matrix &A, + vtkm::Vec &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 +VTKM_EXEC_CONT_EXPORT +vtkm::Vec +MatrixLUPSolve(const vtkm::Matrix &LU, + const vtkm::Vec &permutation, + const vtkm::Vec &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 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 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 +VTKM_EXEC_CONT_EXPORT +vtkm::Vec SolveLinearSystem(const vtkm::Matrix &A, + const vtkm::Vec &b, + bool &valid) +{ + // First, we will make an LUP-factorization to help us. + vtkm::Matrix LU = A; + vtkm::Vec 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 +VTKM_EXEC_CONT_EXPORT +vtkm::Matrix MatrixInverse(const vtkm::Matrix &A, + bool &valid) +{ + // First, we will make an LUP-factorization to help us. + vtkm::Matrix LU = A; + vtkm::Vec 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 invA; + vtkm::Vec ICol(T(0)); + for (vtkm::IdComponent colIndex = 0; colIndex < Size; colIndex++) + { + ICol[colIndex] = 1; + vtkm::Vec invACol = detail::MatrixLUPSolve(LU, permutation, ICol); + ICol[colIndex] = 0; + vtkm::MatrixSetColumn(invA, colIndex, invACol); + } + return invA; +} + +/// Compute the determinant of a matrix. +/// +template +VTKM_EXEC_CONT_EXPORT +T MatrixDeterminant(const vtkm::Matrix &A) +{ + // First, we will make an LUP-factorization to help us. + vtkm::Matrix LU = A; + vtkm::Vec 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 +VTKM_EXEC_CONT_EXPORT +T MatrixDeterminant(const vtkm::Matrix &A) +{ + return A(0,0); +} + +template +VTKM_EXEC_CONT_EXPORT +T MatrixDeterminant(const vtkm::Matrix &A) +{ + return A(0,0)*A(1,1) - A(1,0)*A(0,1); +} + +template +VTKM_EXEC_CONT_EXPORT +T MatrixDeterminant(const vtkm::Matrix &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 +struct TypeTraits > { + typedef typename TypeTraits::NumericTag NumericTag; + typedef TypeTraitsMatrixTag DimensionalityTag; +}; + +/// A matrix has vector traits to implement component-wise operations. +/// +template +struct VecTraits > { +private: + typedef vtkm::Matrix 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 +VTKM_EXEC_CONT_EXPORT +bool operator==(const vtkm::Matrix &a, + const vtkm::Matrix &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 +VTKM_EXEC_CONT_EXPORT +bool operator!=(const vtkm::Matrix &a, + const vtkm::Matrix &b) +{ + return !(a == b); +} + +#endif //vtk_m_Matrix_h diff --git a/vtkm/testing/CMakeLists.txt b/vtkm/testing/CMakeLists.txt index 2390b97ec..64b37a44c 100644 --- a/vtkm/testing/CMakeLists.txt +++ b/vtkm/testing/CMakeLists.txt @@ -31,6 +31,7 @@ set(unit_tests UnitTestExtent.cxx UnitTestListTag.cxx UnitTestMath.cxx + UnitTestMatrix.cxx UnitTestPair.cxx UnitTestTesting.cxx UnitTestTypeListTag.cxx diff --git a/vtkm/testing/UnitTestMatrix.cxx b/vtkm/testing/UnitTestMatrix.cxx new file mode 100644 index 000000000..8fa17d334 --- /dev/null +++ b/vtkm/testing/UnitTestMatrix.cxx @@ -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 + +#include + +#include + +// If more tests need a value for Matrix, we can move this to Testing.h +template +vtkm::Matrix +TestValue(vtkm::Id index, const vtkm::Matrix &) +{ + vtkm::Matrix value; + for (vtkm::IdComponent rowIndex = 0; rowIndex < NumRow; rowIndex++) + { + typedef vtkm::Vec RowType; + RowType row = + TestValue(index, RowType()) + + RowType(static_cast(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 +struct MatrixTest +{ + static const vtkm::IdComponent NUM_ROWS = NumRow; + static const vtkm::IdComponent NUM_COLS = NumCol; + typedef vtkm::Matrix 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(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 ColumnType; + typedef vtkm::Vec 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 rightFactor = + TestValue(1, vtkm::Matrix()); + + vtkm::Matrix product + = vtkm::MatrixMultiply(leftFactor, rightFactor); + + FOR_ROW_COL(product) + { + vtkm::Vec leftVector = vtkm::MatrixGetRow(leftFactor, row); + vtkm::Vec 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 leftVector(2); + vtkm::Vec rightVector; + FOR_ROW_COL(matrixFactor) + { + matrixFactor(row,col) = T(row + 1); + rightVector[col] = T(col + 1); + } + + vtkm::Vec 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 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 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 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 +void MatrixTest1() +{ + MatrixTest::Run(); + MatrixTest::Run(); + MatrixTest::Run(); + MatrixTest::Run(); + MatrixTest::Run(); +} + +template +void NonSingularMatrix(vtkm::Matrix &matrix) +{ + matrix(0,0) = 1; +} + +template +void NonSingularMatrix(vtkm::Matrix &matrix) +{ + matrix(0,0) = -5; matrix(0,1) = 6; + matrix(1,0) = -7; matrix(1,1) = -2; +} + +template +void NonSingularMatrix(vtkm::Matrix &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 +void NonSingularMatrix(vtkm::Matrix &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 +void NonSingularMatrix(vtkm::Matrix &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 +void SingularMatrix(vtkm::Matrix &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 +T RecursiveDeterminant(const vtkm::Matrix &A) +{ + return A(0,0); +} + +template +T RecursiveDeterminant(const vtkm::Matrix &A) +{ + vtkm::Matrix 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 +struct SquareMatrixTest { + static const vtkm::IdComponent SIZE = Size; + typedef vtkm::Matrix 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 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 b; + NonSingularMatrix(A); + for (vtkm::IdComponent index = 0; index < SIZE; index++) + { + b[index] = index+1; + } + bool valid; + + vtkm::Vec x = vtkm::SolveLinearSystem(A, b, valid); + VTKM_TEST_ASSERT(valid, "Matrix declared singular?"); + + // Check result. + vtkm::Vec 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 inverse = vtkm::MatrixInverse(A, valid); + VTKM_TEST_ASSERT(valid, "Matrix declared singular?"); + + // Check result. + vtkm::Matrix product = vtkm::MatrixMultiply(A, inverse); + VTKM_TEST_ASSERT(test_equal(product, vtkm::MatrixIdentity()), + "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 + void operator()(const T&) const { + MatrixTest1(); + MatrixTest1(); + MatrixTest1(); + MatrixTest1(); + MatrixTest1(); + } +}; + +struct SquareMatrixTestFunctor +{ + template + void operator()(const T&) const { + SquareMatrixTest::Run(); + SquareMatrixTest::Run(); + SquareMatrixTest::Run(); + SquareMatrixTest::Run(); + SquareMatrixTest::Run(); + } +}; + +struct VectorMultFunctor +{ + template + void operator()(const VectorType &) const { + // This is mostly to make sure the compile can convert from Tuples + // to vectors. + const int SIZE = vtkm::VecTraits::NUM_COMPONENTS; + typedef typename vtkm::VecTraits::ComponentType ComponentType; + + vtkm::Matrix 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); +}