Commit e25e6e44 authored by Yohann Bearzi's avatar Yohann Bearzi

adding templated matrix inversion / det / solve

And adding 2D array handling for matrices
parent 208a56c6
Pipeline #181463 running with stage
......@@ -254,6 +254,7 @@ set(headers
vtkIOStreamFwd.h
vtkInformationInternals.h
vtkMathUtilities.h
vtkMatrixUtilities.h
vtkMeta.h
vtkNew.h
vtkRange.h
......
This diff is collapsed.
......@@ -41,6 +41,7 @@
#include "vtkCommonCoreModule.h" // For export macro
#include "vtkMathPrivate.hxx" // For Matrix meta-class helpers
#include "vtkMatrixUtilities.h" // For Matrix wrapping / mapping
#include "vtkObject.h"
#include "vtkSmartPointer.h" // For vtkSmartPointer.
#include "vtkTypeTraits.h" // For type traits
......@@ -708,24 +709,6 @@ public:
static void Multiply3x3(const double A[3][3], const double B[3][3], double C[3][3]);
//@}
/**
* Matrix reindexing options for templated methods MultiplyMatrix and
* MultiplyMatrixWithVector
*/
struct MatrixLayout
{
// Matrix is untouched.
using Identity = vtkMathPrivate::MatrixLayout::Identity;
// Matrix is transposed.
using Transpose = vtkMathPrivate::MatrixLayout::Transpose;
// Matrix is diagonal, i.e. value at index idx points to component of
// coordinates (idx, idx) in the diagonal matrix, and all other components
// are set to zero. The matrix does not need to be square.
using Diag = vtkMathPrivate::MatrixLayout::Diag;
};
/**
* Multiply matrices such that M3 = M1 x M2.
* M3 must already be allocated.
......@@ -736,9 +719,9 @@ public:
* Matrices are assumed to be a 1D array implementing `operator[]`. The matrices
* are indexed row by row. Let us develop the effect of LayoutT1 on M1 (the
* same is true for LayoutT2 on M2).
* If LayoutT1 == MatrixLayout::Identity (default), then M1 indexing is untouched.
* If LayoutT1 == MatrixLayout::Transpose, then M1 is read column-wise.
* If LayoutT1 == MatrixLayout::Diag, then M1 is considered to be composed
* If LayoutT1 == vtkMatrixUtilities::Layout::Identity (default), then M1 indexing is untouched.
* If LayoutT1 == vtkMatrixUtilities::Layout::Transpose, then M1 is read column-wise.
* If LayoutT1 == vtkMatrixUtilities::Layout::Diag, then M1 is considered to be composed
* of non zero elements of a diagonal matrix.
*
* @note M3 components indexing can be entirely controlled by swapping M1, M2,
......@@ -749,8 +732,10 @@ public:
* M3 will be diagonal as well (i.e. there won't be allocated memory for
* elements outsize of the diagonal).
*/
template <int RowsT, int MidDimT, int ColsT, class LayoutT1 = MatrixLayout::Identity,
class LayoutT2 = MatrixLayout::Identity, class MatrixT1, class MatrixT2, class MatrixT3>
template <int RowsT, int MidDimT, int ColsT,
class LayoutT1 = vtkMatrixUtilities::Layout::Identity,
class LayoutT2 = vtkMatrixUtilities::Layout::Identity, class MatrixT1, class MatrixT2,
class MatrixT3>
static void MultiplyMatrix(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
{
vtkMathPrivate::MultiplyMatrix<RowsT, MidDimT, ColsT, LayoutT1, LayoutT2>::Compute(M1, M2, M3);
......@@ -765,19 +750,19 @@ public:
*
* Matrix M is assumed to be a 1D array implementing `operator[]`. The matrix
* is indexed row by row. If the input array is indexed columns by colums.
* If LayoutT == MatrixLayout::Identity (default), then M indexing is untouched.
* If LayoutT == MatrixLayout::Transpose, then M is read column-wise.
* If LayoutT == MatrixLayout::Diag, then M is considered to be composed
* If LayoutT == vtkMatrixUtilities::Layout::Identity (default), then M indexing is untouched.
* If LayoutT == vtkMatrixUtilities::Layout::Transpose, then M is read column-wise.
* If LayoutT == vtkMatrixUtilities::Layout::Diag, then M is considered to be composed
* of non zero elements of a diagonal matrix.
*
* VectorT1 and VectorT2 are arrays of size RowsT, and must implement
* `operator[]`.
*
* @warning In the particular case where M1 and M2 BOTH have layout
* MatrixLayout::Diag, RowsT, MidDimT and ColsT MUST match.
* vtkMatrixUtilities::Layout::Diag, RowsT, MidDimT and ColsT MUST match.
*/
template <int RowsT, int ColsT, class LayoutT = MatrixLayout::Identity, class MatrixT,
class VectorT1, class VectorT2>
template <int RowsT, int ColsT, class LayoutT = vtkMatrixUtilities::Layout::Identity,
class MatrixT, class VectorT1, class VectorT2>
static void MultiplyMatrixWithVector(const MatrixT& M, const VectorT1& X, VectorT2& Y)
{
vtkMathPrivate::MultiplyMatrix<RowsT, ColsT, 1, LayoutT>::Compute(M, X, Y);
......@@ -791,8 +776,73 @@ public:
template <class ScalarT, int SizeT, class VectorT1, class VectorT2>
static ScalarT Dot(const VectorT1& x, const VectorT2& y)
{
return vtkMathPrivate::ContractRowWithCol<ScalarT, 1, SizeT, 1, 0, 0, MatrixLayout::Identity,
MatrixLayout::Transpose>::Compute(x, y);
return vtkMathPrivate::ContractRowWithCol<ScalarT, 1, SizeT, 1, 0, 0,
vtkMatrixUtilities::Layout::Identity, vtkMatrixUtilities::Layout::Transpose>::Compute(x, y);
}
/**
* Computes the determinant of input square SizeT x SizeT matrix M.
* The return type is the same as the underlying scalar type of MatrixT.
*
* LayoutT allow to perform basic matrix reindexing. It should be set to a component
* of MatrixLayout.
*
* Matrix M is assumed to be a 1D array implementing `operator[]`. The matrix
* is indexed row by row. If the input array is indexed columns by colums.
* If LayoutT == vtkMatrixUtilities::Layout::Identity (default), then M indexing is untouched.
* If LayoutT == vtkMatrixUtilities::Layout::Transpose, then M is read column-wise.
* If LayoutT == vtkMatrixUtilities::Layout::Diag, then M is considered to be composed
* of non zero elements of a diagonal matrix.
*
* This method is currently implemented for SizeT lower or equal to 3.
*/
template <int SizeT, class LayoutT = vtkMatrixUtilities::Layout::Identity, class MatrixT>
static typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type Determinant(
const MatrixT& M)
{
return vtkMathPrivate::Determinant<SizeT, LayoutT>::Compute(M);
}
/**
* Computes the inverse of input matrix M1 into M2.
*
* LayoutT allow to perform basic matrix reindexing of M1. It should be set to a component
* of MatrixLayout.
*
* Matrix M is assumed to be a 1D array implementing `operator[]`. The matrix
* is indexed row by row. If the input array is indexed columns by colums.
* If LayoutT == vtkMatrixUtilities::Layout::Identity (default), then M indexing is untouched.
* If LayoutT == vtkMatrixUtilities::Layout::Transpose, then M is read column-wise.
* If LayoutT == vtkMatrixUtilities::Layout::Diag, then M is considered to be composed
* of non zero elements of a diagonal matrix.
*
* This method is currently implemented for SizeT lower or equal to 3.
*/
template <int SizeT, class LayoutT = vtkMatrixUtilities::Layout::Identity, class MatrixT1,
class MatrixT2>
static void InvertMatrix(const MatrixT1& M1, MatrixT2& M2)
{
vtkMathPrivate::InvertMatrix<SizeT, LayoutT>::Compute(M1, M2);
}
/**
* This method solves linear systems M * x = y.
*
* LayoutT allow to perform basic matrix reindexing of M1. It should be set to a component
* of MatrixLayout.
*
* Matrix M is assumed to be a 1D array implementing `operator[]`. The matrix
* is indexed row by row. If the input array is indexed columns by colums.
* If LayoutT == vtkMatrixUtilities::Layout::Identity (default), then M indexing is untouched.
* If LayoutT == vtkMatrixUtilities::Layout::Transpose, then M is read column-wise.
* If LayoutT == vtkMatrixUtilities::Layout::Diag, then M is considered to be composed
* of non zero elements of a diagonal matrix.
*/
template <int RowsT, int ColsT, class LayoutT = vtkMatrixUtilities::Layout::Identity,
class MatrixT, class VectorT1, class VectorT2>
static void LinearSolve(const MatrixT& M, const VectorT1& x, VectorT2& y)
{
vtkMathPrivate::LinearSolve<RowsT, ColsT, LayoutT>::Compute(M, x, y);
}
/**
......@@ -804,18 +854,19 @@ public:
*
* Matrix M is assumed to be a 1D array implementing `operator[]`. The matrix
* is indexed row by row. If the input array is indexed columns by colums.
* If LayoutT == MatrixLayout::Identity (default), then M indexing is untouched.
* If LayoutT == MatrixLayout::Transpose, then M is read column-wise.
* If LayoutT == MatrixLayout::Diag, then M is considered to be composed
* If LayoutT == vtkMatrixUtilities::Layout::Identity (default), then M indexing is untouched.
* If LayoutT == vtkMatrixUtilities::Layout::Transpose, then M is read column-wise.
* If LayoutT == vtkMatrixUtilities::Layout::Diag, then M is considered to be composed
* of non zero elements of a diagonal matrix.
*/
template <class ScalarT, int SizeT, class LayoutT, class VectorT1, class MatrixT, class VectorT2>
template <class ScalarT, int SizeT, class LayoutT = vtkMatrixUtilities::Layout::Identity,
class VectorT1, class MatrixT, class VectorT2>
static ScalarT Dot(const VectorT1& x, const MatrixT& M, const VectorT2& y)
{
std::array<ScalarT, SizeT> tmp;
vtkMathPrivate::MultiplyMatrix<SizeT, SizeT, 1, LayoutT>::Compute(M, y, tmp);
return vtkMathPrivate::ContractRowWithCol<ScalarT, 1, SizeT, 1, 0, 0, MatrixLayout::Identity,
MatrixLayout::Transpose>::Compute(x, tmp);
return vtkMathPrivate::ContractRowWithCol<ScalarT, 1, SizeT, 1, 0, 0,
vtkMatrixUtilities::Layout::Identity, vtkMatrixUtilities::Layout::Transpose>::Compute(x, tmp);
}
/**
......
This diff is collapsed.
/*=========================================================================
Program: Visualization Toolkit
Module: vtkMathPrivate.hxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm 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.
=========================================================================*/
/**
* @class vtkMathPrivate
* @brief Internal toolkit used in some vtkMath methods.
*
* vtkMatrixUtilities provides matrix indexing / wrapping tools. One can use
* this utility to wrap a 1D array into a matrix shape, index it at compile
* time.
* @sa
* vtkMath
* vtkMathPrivate
*/
#ifndef vtkMatrixUtilities_hxx
#define vtkMatrixUtilities_hxx
#include <type_traits> // for type traits
namespace vtkMatrixUtilities
{
//=============================================================================
/**
* This struct determines a prior transform to input matrices, chaging the
* way they are indexed
*/
struct Layout
{
/**
* Input matrix is unchanged, i.e. sorted row-wise ordered
*/
struct Identity;
/*
* Input matrix is transposed, i.e. sorted in column-wise ordered.
*/
struct Transpose;
/**
* Input matrix is considered diagonal, and value at index idx points
* to component of coordinates (idx, idx) in the diagonal matrix.
*/
struct Diag;
};
namespace detail
{
// Extracting for STL-like containers
template <int ContainerTypeT, class ContainerT>
struct ScalarTypeExtractor
{
typedef typename ContainerT::value_type value_type;
};
// Extracting for C++ arrays
template <class ContainerT>
struct ScalarTypeExtractor<1, ContainerT>
{
typedef typename std::remove_all_extents<ContainerT>::type value_type;
};
// Extracting for C++ pointers
template <class ContainerT>
struct ScalarTypeExtractor<2, ContainerT>
{
typedef typename std::remove_pointer<ContainerT>::type value_type;
};
} // namespace detail
//=============================================================================
/**
* This class extract the underlying value type of containers. It works on
* multi-dimensional C++ arrays as well as with STL container like that
* have a value_type typedef.
*
* One can access the value type by fetching
* ScalarTypeExtractor<ContainerT>::value_type.
*/
template <class ContainerT>
struct ScalarTypeExtractor
{
typedef typename detail::ScalarTypeExtractor<
// This parameter equals 0, 1 or 2
std::is_array<ContainerT>::value + std::is_pointer<ContainerT>::value * 2,
ContainerT>::value_type value_type;
};
namespace detail
{
// Class actually implementing matrix mapping.
template <int RowsT, int ColsT, class LayoutT>
struct Mapper;
// Specialization of the matrix mapper for when the layout is the identity
template <int RowsT, int ColsT>
struct Mapper<RowsT, ColsT, Layout::Identity>
{
template <int RowT, int ColT>
static constexpr int GetIndex()
{
static_assert(RowT >= 0 && RowT < RowsT, "RowT out of bounds");
static_assert(ColT >= 0 && ColT < ColsT, "ColT out of bounds");
return ColsT * RowT + ColT;
}
};
template <int RowsT, int ColsT>
struct Mapper<RowsT, ColsT, Layout::Transpose>
{
template <int RowT, int ColT>
static constexpr int GetIndex()
{
static_assert(RowT >= 0 && RowT < RowsT, "RowT out of bounds");
static_assert(ColT >= 0 && ColT < ColsT, "ColT out of bounds");
return RowsT * ColT + RowT;
}
};
} // namespace detail
//=============================================================================
/**
* This class is a helper class to compute at compile time the index of a matrix
* stored as a 1D array from its 2D coordinates. This class maps matrices of
* dimension RowsT x ColsT. The LayoutT template parameter permits to switch to
* the indexing of the transpose of the matrix. LayoutT can be set to
* Layout::Identity for a row-wise ordering, or to Layout::Transpose for a
* column-wise ordering
*
* @warning This mapper does not work with matrices stored as 2D arrays, or with
* diagonal matrices.
*/
template <int RowsT, int ColsT, class LayoutT = Layout::Identity>
struct Mapper
{
template <int RowT, int ColT>
static constexpr int GetIndex()
{
return detail::Mapper<RowsT, ColsT, LayoutT>::template GetIndex<RowT, ColT>();
}
};
namespace detail
{
// Class implementing matrix wrapping.
template <int RowsT, int ColsT, class MatrixT, class LayoutT, bool MatrixIs2DArrayT>
struct Wrapper;
// Specializaion of matrix wrapping for matrices stored as 1D arrays
// in row-wise order
template <int RowsT, int ColsT, class MatrixT, class LayoutT>
class Wrapper<RowsT, ColsT, MatrixT, LayoutT, false>
{
private:
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type;
public:
template <int RowT, int ColT>
static const Scalar& Get(const MatrixT& M)
{
return M[Mapper<RowsT, ColsT, LayoutT>::template GetIndex<RowT, ColT>()];
}
template <int RowT, int ColT>
static Scalar& Get(MatrixT& M)
{
return M[Mapper<RowsT, ColsT, LayoutT>::template GetIndex<RowT, ColT>()];
}
};
// Specialization for matrices stored as 2D arrays with an unchanged layout
template <int RowsT, int ColsT, class MatrixT>
class Wrapper<RowsT, ColsT, MatrixT, Layout::Identity, true>
{
private:
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type;
public:
template <int RowT, int ColT>
static const Scalar& Get(const MatrixT& M)
{
return M[RowT][ColT];
}
template <int RowT, int ColT>
static Scalar& Get(MatrixT& M)
{
return M[RowT][ColT];
}
};
// Specialization for matrices stored as 2D arrays read as its transposed self.
template <int RowsT, int ColsT, class MatrixT>
class Wrapper<RowsT, ColsT, MatrixT, Layout::Transpose, true>
{
private:
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type;
public:
template <int RowT, int ColT>
static const Scalar& Get(const MatrixT& M)
{
return M[ColT][RowT];
}
template <int RowT, int ColT>
static Scalar& Get(MatrixT& M)
{
return M[ColT][RowT];
}
};
// Specialization for diagonal matrices.
// Note: a diagonal matrix has to be stored in a 1D array.
template <int RowsT, int ColsT, class MatrixT>
class Wrapper<RowsT, ColsT, MatrixT, Layout::Diag, false>
{
private:
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type;
template <int RowT, int ColT>
struct Helper
{
static constexpr Scalar ZERO = Scalar(0);
static Scalar& Get(const MatrixT& M) { return ZERO; }
};
template <int RowT>
struct Helper<RowT, RowT>
{
static Scalar& Get(MatrixT& M) { return M[RowT]; }
static const Scalar& Get(const MatrixT& M) { return M[RowT]; }
};
public:
template <int RowT, int ColT>
const Scalar& Get(const MatrixT& M)
{
return Helper<RowT, ColT>::Get(M);
}
template <int RowT, int ColT>
Scalar& Get(MatrixT& M)
{
return Helper<RowT, ColT>::Get(M);
}
};
} // namespace detail
//=============================================================================
/**
* Matrix wrapping class. This class implements a getter templated on the
* coordinates of the wanted element. A matrix can be a 2D C++ array, a 1D C++
* array row-wise ordered, or any STL-like container implementing operator[] and
* having a value_type typedef.
*
* This class wraps a RowsT x ColsT matrix stored in the container MatrixT. The
* LayoutT template parameter permits to reindex at compile-time the matrix. If
* it is set to Layout::Identity, the matrix is assumed to be row-wised ordered.
* If it is set to Layout::Transpose, the matrix is assumed to be column-wise ordered.
* One can also convert a 1D input array into a diagonal matrix by setting
* LayoutT to Layout::Diag. In ths particular case, method Get will return a
* read-only zero on elements outside of the diagonal.
*/
template <int RowsT, int ColsT, class MatrixT, class LayoutT = Layout::Identity>
class Wrapper
{
private:
using Scalar = typename ScalarTypeExtractor<MatrixT>::value_type;
static constexpr bool MatrixIs2DArray =
std::is_array<typename std::remove_extent<MatrixT>::type>::value;
static_assert(!MatrixIs2DArray || !std::is_same<LayoutT, Layout::Diag>::value,
"A diagonal matrix cannot be a 2D array");
public:
template <int RowT, int ColT>
static const Scalar& Get(const MatrixT& M)
{
return detail::Wrapper<RowsT, ColsT, MatrixT, LayoutT, MatrixIs2DArray>::template Get<RowT,
ColT>(M);
}
template <int RowT, int ColT>
static Scalar& Get(MatrixT& M)
{
return detail::Wrapper<RowsT, ColsT, MatrixT, LayoutT, MatrixIs2DArray>::template Get<RowT,
ColT>(M);
}
};
} // namespace vtkMatrixUtilities
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment