958 lines
42 KiB
C++
958 lines
42 KiB
C++
|
/*=========================================================================
|
||
|
|
|
||
|
|
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.
|
||
|
|
*
|
||
|
|
* vtkMathPrivate provides meta-classes helpers that are used in some vtkMath
|
||
|
|
* methods.
|
||
|
|
* @sa
|
||
|
|
* vtkMath
|
||
|
|
* vtkMatrixUtilities
|
||
|
|
*/
|
||
|
|
|
||
|
|
#ifndef vtkMathPrivate_hxx
|
||
|
|
#define vtkMathPrivate_hxx
|
||
|
|
|
||
|
|
#include "vtkCommonCoreModule.h" //required for correct implementation
|
||
|
|
#include "vtkMatrixUtilities.h"
|
||
|
|
|
||
|
|
#include <type_traits>
|
||
|
|
|
||
|
|
namespace vtkMathPrivate
|
||
|
|
{
|
||
|
|
static constexpr int VTK_MATH_PRIVATE_PACK_SIZE = 4;
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// This class computes the dot product between row RowT of matrix M1
|
||
|
|
// and column ColT of matrix M2.
|
||
|
|
// The template parameter IdxT is on index ahead from computation.
|
||
|
|
// Template parameters LayoutT1 and LayoutT2 respectively reindex
|
||
|
|
// input matrices M1 and M2 according to MatrixLayout enumeration
|
||
|
|
// M1 (or M1^T if LayoutT1 == vtkMatrixUtilities::Layout::Transpose) is a matrix of RowsT x
|
||
|
|
// MidDimT M2 (or M2^T if LayoutT2 == vtkMatrixUtilities::Layout::Transpose) is a matrix of MidDimT
|
||
|
|
// x ColsDimT
|
||
|
|
template <class ScalarT, int RowsT, int MidDimT, int ColsT, int RowT, int ColT,
|
||
|
|
class LayoutT1 = vtkMatrixUtilities::Layout::Identity,
|
||
|
|
class LayoutT2 = vtkMatrixUtilities::Layout::Identity, int IdxT = 0, int PackSizeT = MidDimT>
|
||
|
|
class ContractRowWithCol
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static ScalarT Compute(const MatrixT1& M1, const MatrixT2& M2)
|
||
|
|
{
|
||
|
|
return ContractRowWithCol<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT, LayoutT1, LayoutT2, IdxT,
|
||
|
|
Shift>::Compute(M1, M2) +
|
||
|
|
ContractRowWithCol<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT, LayoutT1, LayoutT2,
|
||
|
|
IdxT + Shift, PackSizeT - Shift>::Compute(M1, M2);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT >= 0 && RowT < RowsT;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT >= 0 && ColT < ColsT;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
|
||
|
|
static constexpr int Shift =
|
||
|
|
(MidDimT - IdxT) / VTK_MATH_PRIVATE_PACK_SIZE ? VTK_MATH_PRIVATE_PACK_SIZE : MidDimT - IdxT - 1;
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for when there are 4 components left to compute
|
||
|
|
template <class ScalarT, int RowsT, int MidDimT, int ColsT, int RowT, int ColT, class LayoutT1,
|
||
|
|
class LayoutT2, int IdxT>
|
||
|
|
class ContractRowWithCol<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT, LayoutT1, LayoutT2, IdxT, 4>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static ScalarT Compute(const MatrixT1& M1, const MatrixT2& M2)
|
||
|
|
{
|
||
|
|
using Wrap1 = vtkMatrixUtilities::Wrapper<RowsT, MidDimT, MatrixT1, LayoutT1>;
|
||
|
|
using Wrap2 = vtkMatrixUtilities::Wrapper<MidDimT, ColsT, MatrixT2, LayoutT2>;
|
||
|
|
|
||
|
|
return Wrap1::template Get<RowT, IdxT>(M1) * Wrap2::template Get<IdxT, ColT>(M2) +
|
||
|
|
Wrap1::template Get<RowT, IdxT + 1>(M1) * Wrap2::template Get<IdxT + 1, ColT>(M2) +
|
||
|
|
Wrap1::template Get<RowT, IdxT + 2>(M1) * Wrap2::template Get<IdxT + 2, ColT>(M2) +
|
||
|
|
Wrap1::template Get<RowT, IdxT + 3>(M1) * Wrap2::template Get<IdxT + 3, ColT>(M2);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT >= 0 && RowT < RowsT;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT >= 0 && ColT < ColsT;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for when there are 3 components left to compute
|
||
|
|
template <class ScalarT, int RowsT, int MidDimT, int ColsT, int RowT, int ColT, class LayoutT1,
|
||
|
|
class LayoutT2, int IdxT>
|
||
|
|
class ContractRowWithCol<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT, LayoutT1, LayoutT2, IdxT, 3>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static ScalarT Compute(const MatrixT1& M1, const MatrixT2& M2)
|
||
|
|
{
|
||
|
|
using Wrap1 = vtkMatrixUtilities::Wrapper<RowsT, MidDimT, MatrixT1, LayoutT1>;
|
||
|
|
using Wrap2 = vtkMatrixUtilities::Wrapper<MidDimT, ColsT, MatrixT2, LayoutT2>;
|
||
|
|
|
||
|
|
return Wrap1::template Get<RowT, IdxT>(M1) * Wrap2::template Get<IdxT, ColT>(M2) +
|
||
|
|
Wrap1::template Get<RowT, IdxT + 1>(M1) * Wrap2::template Get<IdxT + 1, ColT>(M2) +
|
||
|
|
Wrap1::template Get<RowT, IdxT + 2>(M1) * Wrap2::template Get<IdxT + 2, ColT>(M2);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT >= 0 && RowT < RowsT;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT >= 0 && ColT < ColsT;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for when there are 2 components left to compute
|
||
|
|
template <class ScalarT, int RowsT, int MidDimT, int ColsT, int RowT, int ColT, class LayoutT1,
|
||
|
|
class LayoutT2, int IdxT>
|
||
|
|
class ContractRowWithCol<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT, LayoutT1, LayoutT2, IdxT, 2>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static ScalarT Compute(const MatrixT1& M1, const MatrixT2& M2)
|
||
|
|
{
|
||
|
|
using Wrap1 = vtkMatrixUtilities::Wrapper<RowsT, MidDimT, MatrixT1, LayoutT1>;
|
||
|
|
using Wrap2 = vtkMatrixUtilities::Wrapper<MidDimT, ColsT, MatrixT2, LayoutT2>;
|
||
|
|
|
||
|
|
return Wrap1::template Get<RowT, IdxT>(M1) * Wrap2::template Get<IdxT, ColT>(M2) +
|
||
|
|
Wrap1::template Get<RowT, IdxT + 1>(M1) * Wrap2::template Get<IdxT + 1, ColT>(M2);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT >= 0 && RowT < RowsT;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT >= 0 && ColT < ColsT;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for when there is 1 component left to compute
|
||
|
|
template <class ScalarT, int RowsT, int MidDimT, int ColsT, int RowT, int ColT, class LayoutT1,
|
||
|
|
class LayoutT2, int IdxT>
|
||
|
|
class ContractRowWithCol<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT, LayoutT1, LayoutT2, IdxT, 1>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static ScalarT Compute(const MatrixT1& M1, const MatrixT2& M2)
|
||
|
|
{
|
||
|
|
using Wrap1 = vtkMatrixUtilities::Wrapper<RowsT, MidDimT, MatrixT1, LayoutT1>;
|
||
|
|
using Wrap2 = vtkMatrixUtilities::Wrapper<MidDimT, ColsT, MatrixT2, LayoutT2>;
|
||
|
|
|
||
|
|
return Wrap1::template Get<RowT, IdxT>(M1) * Wrap2::template Get<IdxT, ColT>(M2);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT >= 0 && RowT < RowsT;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT >= 0 && ColT < ColsT;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// This class handles row and col contraction when at least one of the 2 input
|
||
|
|
// matrices are diagonal
|
||
|
|
template <class ScalarT, int RowsT, int MidDimT, int ColsT, int RowT, int ColT, class LayoutT1,
|
||
|
|
class LayoutT2>
|
||
|
|
class DiagContractRowWithCol;
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for when M1 is diagonal
|
||
|
|
template <class ScalarT, int RowsT, int MidDimT, int ColsT, int RowT, int ColT, class LayoutT2>
|
||
|
|
class DiagContractRowWithCol<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT,
|
||
|
|
vtkMatrixUtilities::Layout::Diag, LayoutT2>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static ScalarT Compute(const MatrixT1& M1, const MatrixT2& M2)
|
||
|
|
{
|
||
|
|
return M1[RowT] * M2[Mapper2::template GetIndex<RowT, ColT>()];
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT >= 0 && RowT < RowsT;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT >= 0 && ColT < ColsT;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
|
||
|
|
using Mapper2 = vtkMatrixUtilities::Mapper<MidDimT, ColsT, LayoutT2>;
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for when M2 is diagonal
|
||
|
|
template <class ScalarT, int RowsT, int MidDimT, int ColsT, int RowT, int ColT, class LayoutT1>
|
||
|
|
class DiagContractRowWithCol<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT, LayoutT1,
|
||
|
|
vtkMatrixUtilities::Layout::Diag>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static ScalarT Compute(const MatrixT1& M1, const MatrixT2& M2)
|
||
|
|
{
|
||
|
|
return M1[Mapper1::template GetIndex<RowT, ColT>()] * M2[ColT];
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT >= 0 && RowT < RowsT;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT >= 0 && ColT < ColsT;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
|
||
|
|
using Mapper1 = vtkMatrixUtilities::Mapper<RowsT, MidDimT, LayoutT1>;
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// This class returns zero. Is is used for rectangular diagonal matrices, when
|
||
|
|
// lines / rows are filled with zeros.
|
||
|
|
template <class ScalarT>
|
||
|
|
class NullContractRowWithCol
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static ScalarT Compute(const MatrixT1&, const MatrixT2&)
|
||
|
|
{
|
||
|
|
return ScalarT(0);
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Helper class to choose between regular contraction class vs diagonal
|
||
|
|
// handling contraction class.
|
||
|
|
// By default, Type is an instance of ContractRowWithCol
|
||
|
|
template <class ScalarT, int RowsT, int MidDimT, int ColsT, int RowT, int ColT, class LayoutT1,
|
||
|
|
class LayoutT2, bool HasOneDiagonalMatrixT = false>
|
||
|
|
struct ContractRowWithColSwitch
|
||
|
|
{
|
||
|
|
typedef ContractRowWithCol<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT, LayoutT1, LayoutT2> Type;
|
||
|
|
|
||
|
|
static constexpr bool NO_INPUT_MATRIX_CAN_BE_DIAGONAL =
|
||
|
|
!std::is_same<LayoutT1, vtkMatrixUtilities::Layout::Diag>::value &&
|
||
|
|
!std::is_same<LayoutT2, vtkMatrixUtilities::Layout::Diag>::value;
|
||
|
|
static_assert(NO_INPUT_MATRIX_CAN_BE_DIAGONAL,
|
||
|
|
"LayoutT1 and LayoutT2 cannot equal vtkMatrixUtilities::Layout::Diag in this setup");
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for diagonal matrices.
|
||
|
|
// Rectangular diagonal matrices are handled using NullContractRowWithCol.
|
||
|
|
template <class ScalarT, int RowsT, int MidDimT, int ColsT, int RowT, int ColT, class LayoutT1,
|
||
|
|
class LayoutT2>
|
||
|
|
struct ContractRowWithColSwitch<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT, LayoutT1, LayoutT2,
|
||
|
|
true>
|
||
|
|
{
|
||
|
|
static constexpr bool IsDiagonal1 =
|
||
|
|
std::is_same<LayoutT1, vtkMatrixUtilities::Layout::Diag>::value;
|
||
|
|
static constexpr bool IsDiagonal2 =
|
||
|
|
std::is_same<LayoutT2, vtkMatrixUtilities::Layout::Diag>::value;
|
||
|
|
|
||
|
|
// If on of the diagonal matrix is rectangular and we are out of the diagonal
|
||
|
|
// length, all remaining values are null.
|
||
|
|
static constexpr bool UseNullContractRowWithCol =
|
||
|
|
(IsDiagonal1 && RowT >= MidDimT) || (IsDiagonal2 && ColT >= MidDimT);
|
||
|
|
|
||
|
|
typedef typename std::conditional<UseNullContractRowWithCol, NullContractRowWithCol<ScalarT>,
|
||
|
|
DiagContractRowWithCol<ScalarT, RowsT, MidDimT, ColsT, RowT, ColT, LayoutT1, LayoutT2>>::type
|
||
|
|
Type;
|
||
|
|
|
||
|
|
static constexpr bool NEEDS_AT_LEAST_ONE_DIAGONAL_INPUT_MATRIX = IsDiagonal1 || IsDiagonal2;
|
||
|
|
static_assert(NEEDS_AT_LEAST_ONE_DIAGONAL_INPUT_MATRIX,
|
||
|
|
"LayoutT1 or LayoutT2 must equal vtkMatrixUtilities::Layout::Diag in this setup");
|
||
|
|
};
|
||
|
|
|
||
|
|
namespace detail
|
||
|
|
{
|
||
|
|
//=============================================================================
|
||
|
|
// Class in charge for actually multiplying 2 matrices.
|
||
|
|
// This method is called by MultiplyMatrix::Compute for chunks of size
|
||
|
|
// VTK_MATH_PRIVATE_PACK_SIZE at most, specified in ColPackSizeT.
|
||
|
|
// This class mostly consists on explicitly onfold computation for those chunks.
|
||
|
|
template <int RowsT, int MidDimT, int ColsT, class LayoutT1, class LayoutT2, int RowT, int ColT,
|
||
|
|
int ColPackSizeT>
|
||
|
|
class MultiplyMatrix;
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for a chunk of size 4
|
||
|
|
template <int RowsT, int MidDimT, int ColsT, class LayoutT1, class LayoutT2, int RowT, int ColT>
|
||
|
|
class MultiplyMatrix<RowsT, MidDimT, ColsT, LayoutT1, LayoutT2, RowT, ColT, 4>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3>
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
|
||
|
|
{
|
||
|
|
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT3>::value_type;
|
||
|
|
using Wrap3 = vtkMatrixUtilities::Wrapper<RowsT, ColsT, MatrixT3>;
|
||
|
|
|
||
|
|
Wrap3::template Get<RowT, ColT>(M3) = ContractRowWithColSwitch<Scalar, RowsT, MidDimT, ColsT,
|
||
|
|
RowT, ColT, LayoutT1, LayoutT2, OneMatrixIsDiagonal>::Type::Compute(M1, M2);
|
||
|
|
Wrap3::template Get<RowT, ColT + 1>(M3) = ContractRowWithColSwitch<Scalar, RowsT, MidDimT,
|
||
|
|
ColsT, RowT, ColT + 1, LayoutT1, LayoutT2, OneMatrixIsDiagonal>::Type::Compute(M1, M2);
|
||
|
|
Wrap3::template Get<RowT, ColT + 2>(M3) = ContractRowWithColSwitch<Scalar, RowsT, MidDimT,
|
||
|
|
ColsT, RowT, ColT + 2, LayoutT1, LayoutT2, OneMatrixIsDiagonal>::Type::Compute(M1, M2);
|
||
|
|
Wrap3::template Get<RowT, ColT + 3>(M3) = ContractRowWithColSwitch<Scalar, RowsT, MidDimT,
|
||
|
|
ColsT, RowT, ColT + 3, LayoutT1, LayoutT2, OneMatrixIsDiagonal>::Type::Compute(M1, M2);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT < RowsT && RowT >= 0;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT + 3 < ColsT && ColT >= 0;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
|
||
|
|
static constexpr bool OneMatrixIsDiagonal =
|
||
|
|
std::is_same<LayoutT1, vtkMatrixUtilities::Layout::Diag>::value ||
|
||
|
|
std::is_same<LayoutT2, vtkMatrixUtilities::Layout::Diag>::value;
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for a chunk of size 3
|
||
|
|
template <int RowsT, int MidDimT, int ColsT, class LayoutT1, class LayoutT2, int RowT, int ColT>
|
||
|
|
class MultiplyMatrix<RowsT, MidDimT, ColsT, LayoutT1, LayoutT2, RowT, ColT, 3>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3>
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
|
||
|
|
{
|
||
|
|
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT3>::value_type;
|
||
|
|
using Wrap3 = vtkMatrixUtilities::Wrapper<RowsT, ColsT, MatrixT3>;
|
||
|
|
|
||
|
|
Wrap3::template Get<RowT, ColT>(M3) = ContractRowWithColSwitch<Scalar, RowsT, MidDimT, ColsT,
|
||
|
|
RowT, ColT, LayoutT1, LayoutT2, OneMatrixIsDiagonal>::Type::Compute(M1, M2);
|
||
|
|
Wrap3::template Get<RowT, ColT + 1>(M3) = ContractRowWithColSwitch<Scalar, RowsT, MidDimT,
|
||
|
|
ColsT, RowT, ColT + 1, LayoutT1, LayoutT2, OneMatrixIsDiagonal>::Type::Compute(M1, M2);
|
||
|
|
Wrap3::template Get<RowT, ColT + 2>(M3) = ContractRowWithColSwitch<Scalar, RowsT, MidDimT,
|
||
|
|
ColsT, RowT, ColT + 2, LayoutT1, LayoutT2, OneMatrixIsDiagonal>::Type::Compute(M1, M2);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT < RowsT && RowT >= 0;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT + 2 < ColsT && ColT >= 0;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
|
||
|
|
static constexpr bool OneMatrixIsDiagonal =
|
||
|
|
std::is_same<LayoutT1, vtkMatrixUtilities::Layout::Diag>::value ||
|
||
|
|
std::is_same<LayoutT2, vtkMatrixUtilities::Layout::Diag>::value;
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for a chunk of size 2
|
||
|
|
template <int RowsT, int MidDimT, int ColsT, class LayoutT1, class LayoutT2, int RowT, int ColT>
|
||
|
|
class MultiplyMatrix<RowsT, MidDimT, ColsT, LayoutT1, LayoutT2, RowT, ColT, 2>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3>
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
|
||
|
|
{
|
||
|
|
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT3>::value_type;
|
||
|
|
using Wrap3 = vtkMatrixUtilities::Wrapper<RowsT, ColsT, MatrixT3>;
|
||
|
|
|
||
|
|
Wrap3::template Get<RowT, ColT>(M3) = ContractRowWithColSwitch<Scalar, RowsT, MidDimT, ColsT,
|
||
|
|
RowT, ColT, LayoutT1, LayoutT2, OneMatrixIsDiagonal>::Type::Compute(M1, M2);
|
||
|
|
Wrap3::template Get<RowT, ColT + 1>(M3) = ContractRowWithColSwitch<Scalar, RowsT, MidDimT,
|
||
|
|
ColsT, RowT, ColT + 1, LayoutT1, LayoutT2, OneMatrixIsDiagonal>::Type::Compute(M1, M2);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT < RowsT && RowT >= 0;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT + 1 < ColsT && ColT >= 0;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
|
||
|
|
static constexpr bool OneMatrixIsDiagonal =
|
||
|
|
std::is_same<LayoutT1, vtkMatrixUtilities::Layout::Diag>::value ||
|
||
|
|
std::is_same<LayoutT2, vtkMatrixUtilities::Layout::Diag>::value;
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for a chunk of size 1
|
||
|
|
template <int RowsT, int MidDimT, int ColsT, class LayoutT1, class LayoutT2, int RowT, int ColT>
|
||
|
|
class MultiplyMatrix<RowsT, MidDimT, ColsT, LayoutT1, LayoutT2, RowT, ColT, 1>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3>
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
|
||
|
|
{
|
||
|
|
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT3>::value_type;
|
||
|
|
using Wrap3 = vtkMatrixUtilities::Wrapper<RowsT, ColsT, MatrixT3>;
|
||
|
|
|
||
|
|
Wrap3::template Get<RowT, ColT>(M3) = ContractRowWithColSwitch<Scalar, RowsT, MidDimT, ColsT,
|
||
|
|
RowT, ColT, LayoutT1, LayoutT2, OneMatrixIsDiagonal>::Type::Compute(M1, M2);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = RowT < RowsT && RowT >= 0;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = ColT < ColsT && ColT >= 0;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
|
||
|
|
static constexpr bool OneMatrixIsDiagonal =
|
||
|
|
std::is_same<LayoutT1, vtkMatrixUtilities::Layout::Diag>::value ||
|
||
|
|
std::is_same<LayoutT2, vtkMatrixUtilities::Layout::Diag>::value;
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization when both input matrices are diagonal for a chunk of size 4
|
||
|
|
template <int SizeT, int IdxT>
|
||
|
|
class MultiplyMatrix<SizeT, SizeT, SizeT, vtkMatrixUtilities::Layout::Diag,
|
||
|
|
vtkMatrixUtilities::Layout::Diag, IdxT, IdxT, 4>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3>
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
|
||
|
|
{
|
||
|
|
M3[IdxT] = M1[IdxT] * M2[IdxT];
|
||
|
|
M3[IdxT + 1] = M1[IdxT + 1] * M2[IdxT + 1];
|
||
|
|
M3[IdxT + 2] = M1[IdxT + 2] * M2[IdxT + 2];
|
||
|
|
M3[IdxT + 3] = M1[IdxT + 3] * M2[IdxT + 3];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization when both input matrices are diagonal for a chunk of size 3
|
||
|
|
template <int SizeT, int IdxT>
|
||
|
|
class MultiplyMatrix<SizeT, SizeT, SizeT, vtkMatrixUtilities::Layout::Diag,
|
||
|
|
vtkMatrixUtilities::Layout::Diag, IdxT, IdxT, 3>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3>
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
|
||
|
|
{
|
||
|
|
M3[IdxT] = M1[IdxT] * M2[IdxT];
|
||
|
|
M3[IdxT + 1] = M1[IdxT + 1] * M2[IdxT + 1];
|
||
|
|
M3[IdxT + 2] = M1[IdxT + 2] * M2[IdxT + 2];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization when both input matrices are diagonal for a chunk of size 2
|
||
|
|
template <int SizeT, int IdxT>
|
||
|
|
class MultiplyMatrix<SizeT, SizeT, SizeT, vtkMatrixUtilities::Layout::Diag,
|
||
|
|
vtkMatrixUtilities::Layout::Diag, IdxT, IdxT, 2>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3>
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
|
||
|
|
{
|
||
|
|
M3[IdxT] = M1[IdxT] * M2[IdxT];
|
||
|
|
M3[IdxT + 1] = M1[IdxT + 1] * M2[IdxT + 1];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization when both input matrices are diagonal for a chunk of size 1
|
||
|
|
template <int SizeT, int IdxT>
|
||
|
|
class MultiplyMatrix<SizeT, SizeT, SizeT, vtkMatrixUtilities::Layout::Diag,
|
||
|
|
vtkMatrixUtilities::Layout::Diag, IdxT, IdxT, 1>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3>
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
|
||
|
|
{
|
||
|
|
M3[IdxT] = M1[IdxT] * M2[IdxT];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
} // namespace detail
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Multiply matrices such that M3 = M1 * M2.
|
||
|
|
// Template parameters LayoutT1 and LayoutT2 respectively reindex
|
||
|
|
// input matrices M1 and M2 following MatrixLayout options
|
||
|
|
// Hence, if LayoutT1 == vtkMatrixUtilities::Layout::Transpose, then M3 = M1^T * M2, and so on.
|
||
|
|
// M1 (or M1^T if LayoutT1 == vtkMatrixUtilities::Layout::Transpose) is a matrix of RowsT x
|
||
|
|
// MidDimT M2 (or M2^T if LayoutT2 == vtkMatrixUtilities::Layout::Transpose) is a matrix of MidDimT
|
||
|
|
// x ColsDimT
|
||
|
|
//
|
||
|
|
// RemainingRowSizeT should be disregarded when first instantiating this class. It is
|
||
|
|
// the number of remaining elements in the current row to process.
|
||
|
|
//
|
||
|
|
// To compute the multiplication, each component of the output matrix
|
||
|
|
// is computed chunk by chunk (of size VTK_MATH_PRIVATE_PACK_SIZE),
|
||
|
|
// starting at the top left, sweeping the rows one by one.
|
||
|
|
template <int RowsT, int MidDimT, int ColsT, class LayoutT1 = vtkMatrixUtilities::Layout::Identity,
|
||
|
|
class LayoutT2 = vtkMatrixUtilities::Layout::Identity, int NextRowT = 1, int NextColT = 1,
|
||
|
|
int RemainingRowSizeT = ColsT>
|
||
|
|
class MultiplyMatrix
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3>
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
|
||
|
|
{
|
||
|
|
detail::MultiplyMatrix<RowsT, MidDimT, ColsT, LayoutT1, LayoutT2, Row, Col, Shift>::Compute(
|
||
|
|
M1, M2, M3);
|
||
|
|
MultiplyMatrix<RowsT, MidDimT, ColsT, LayoutT1, LayoutT2, ForwardRow, ForwardCol,
|
||
|
|
ForwardRemainingRowSize>::Compute(M1, M2, M3);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr int Row = NextRowT - 1;
|
||
|
|
static constexpr int Col = NextColT - 1;
|
||
|
|
|
||
|
|
// We go to a new row when the last chunk is lower or equal to
|
||
|
|
// VTK_MATH_PRIVATE_PACK_SIZE
|
||
|
|
static constexpr bool GoToNewRow = !((ColsT - Col - 1) / VTK_MATH_PRIVATE_PACK_SIZE);
|
||
|
|
|
||
|
|
// This is the chunk size, which is from what we should shift in recursive
|
||
|
|
// calls
|
||
|
|
static constexpr int Shift = GoToNewRow ? ColsT - Col : VTK_MATH_PRIVATE_PACK_SIZE;
|
||
|
|
|
||
|
|
// Row / col / pack size for the next recursive call, depending on if we
|
||
|
|
// change row or not.
|
||
|
|
static constexpr int ForwardRow = GoToNewRow ? NextRowT + 1 : NextRowT;
|
||
|
|
static constexpr int ForwardCol = GoToNewRow ? 1 : NextColT + Shift;
|
||
|
|
static constexpr int ForwardRemainingRowSize = GoToNewRow ? ColsT : RemainingRowSizeT - Shift;
|
||
|
|
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = Row < RowsT && Row >= 0;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = Col < ColsT && Col >= 0;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for when both input matrices are diagonal
|
||
|
|
// Warning: RowsT, MidDimT and ColsT MUST match
|
||
|
|
template <int RowsT, int MidDimT, int ColsT, int NextIdxT, int RemainingRowSizeT>
|
||
|
|
class MultiplyMatrix<RowsT, MidDimT, ColsT, vtkMatrixUtilities::Layout::Diag,
|
||
|
|
vtkMatrixUtilities::Layout::Diag, NextIdxT, NextIdxT, RemainingRowSizeT>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3>
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3)
|
||
|
|
{
|
||
|
|
detail::MultiplyMatrix<RowsT, MidDimT, ColsT, vtkMatrixUtilities::Layout::Diag,
|
||
|
|
vtkMatrixUtilities::Layout::Diag, NextIdxT - 1, NextIdxT - 1, Shift>::Compute(M1, M2, M3);
|
||
|
|
MultiplyMatrix<RowsT, MidDimT, ColsT, vtkMatrixUtilities::Layout::Diag,
|
||
|
|
vtkMatrixUtilities::Layout::Diag, NextIdxT + Shift, NextIdxT + Shift,
|
||
|
|
RemainingRowSizeT - Shift>::Compute(M1, M2, M3);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
static constexpr bool Shift = VTK_MATH_PRIVATE_PACK_SIZE;
|
||
|
|
|
||
|
|
static constexpr bool DIAGONAL_MATRICES_DIMENSIONS_DONT_MATCH =
|
||
|
|
RowsT == MidDimT && RowsT == ColsT;
|
||
|
|
static_assert(DIAGONAL_MATRICES_DIMENSIONS_DONT_MATCH, "There must be RowsT = MidDimT = colsT");
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = NextIdxT <= RowsT && NextIdxT > 0;
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds");
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = NextIdxT <= ColsT && NextIdxT > 0;
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds");
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Macro defining the specialization for when both input matrices are diagonal,
|
||
|
|
// implementing the last chunk to compute. They should be specialized
|
||
|
|
// for integers below VTK_MATH_PRIVATE_PACK_SIZE
|
||
|
|
#define vtkEndForBothDiagonalMultiplyMatrixSpecializationMacro(RemainingRowSize) \
|
||
|
|
template <int RowsT, int MidDimT, int ColsT, int NextIdxT> \
|
||
|
|
class MultiplyMatrix<RowsT, MidDimT, ColsT, vtkMatrixUtilities::Layout::Diag, \
|
||
|
|
vtkMatrixUtilities::Layout::Diag, NextIdxT, NextIdxT, RemainingRowSize> \
|
||
|
|
{ \
|
||
|
|
public: \
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3> \
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3) \
|
||
|
|
{ \
|
||
|
|
detail::MultiplyMatrix<RowsT, MidDimT, ColsT, vtkMatrixUtilities::Layout::Diag, \
|
||
|
|
vtkMatrixUtilities::Layout::Diag, NextIdxT - 1, NextIdxT - 1, \
|
||
|
|
RemainingRowSize>::Compute(M1, M2, M3); \
|
||
|
|
} \
|
||
|
|
\
|
||
|
|
private: \
|
||
|
|
static constexpr bool DIAGONAL_MATRICES_DIMENSIONS_DONT_MATCH = \
|
||
|
|
RowsT == MidDimT && RowsT == ColsT; \
|
||
|
|
static_assert( \
|
||
|
|
DIAGONAL_MATRICES_DIMENSIONS_DONT_MATCH, "There must be RowsT = MidDimT = colsT"); \
|
||
|
|
static constexpr bool ROW_OUT_OF_BOUNDS = NextIdxT <= RowsT && NextIdxT > 0; \
|
||
|
|
static_assert(ROW_OUT_OF_BOUNDS, "RowT is out of bounds"); \
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = NextIdxT <= ColsT && NextIdxT > 0; \
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds"); \
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
vtkEndForBothDiagonalMultiplyMatrixSpecializationMacro(1);
|
||
|
|
vtkEndForBothDiagonalMultiplyMatrixSpecializationMacro(2);
|
||
|
|
vtkEndForBothDiagonalMultiplyMatrixSpecializationMacro(3);
|
||
|
|
vtkEndForBothDiagonalMultiplyMatrixSpecializationMacro(4);
|
||
|
|
|
||
|
|
#undef vtkEndForBothDiagonalMultiplyMatrixSpecializationMacro
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Macro defining specialization for the last chunk to compute. They should be
|
||
|
|
// specialized for integers below VTK_MATH_PRIVATE_PACK_SIZE
|
||
|
|
#define vtkLastRowMultiplyMatrixSpecializationMacro(RemainingRowSize) \
|
||
|
|
template <int RowsT, int MidDimT, int ColsT, class LayoutT1, class LayoutT2, int NextColT> \
|
||
|
|
class MultiplyMatrix<RowsT, MidDimT, ColsT, LayoutT1, LayoutT2, RowsT, NextColT, \
|
||
|
|
RemainingRowSize> \
|
||
|
|
{ \
|
||
|
|
public: \
|
||
|
|
template <class MatrixT1, class MatrixT2, class MatrixT3> \
|
||
|
|
static void Compute(const MatrixT1& M1, const MatrixT2& M2, MatrixT3& M3) \
|
||
|
|
{ \
|
||
|
|
detail::MultiplyMatrix<RowsT, MidDimT, ColsT, LayoutT1, LayoutT2, RowsT - 1, NextColT - 1, \
|
||
|
|
RemainingRowSize>::Compute(M1, M2, M3); \
|
||
|
|
} \
|
||
|
|
\
|
||
|
|
private: \
|
||
|
|
static constexpr bool COL_OUT_OF_BOUNDS = NextColT - 1 < ColsT && NextColT >= 1; \
|
||
|
|
static_assert(COL_OUT_OF_BOUNDS, "ColT is out of bounds"); \
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
vtkLastRowMultiplyMatrixSpecializationMacro(1);
|
||
|
|
vtkLastRowMultiplyMatrixSpecializationMacro(2);
|
||
|
|
vtkLastRowMultiplyMatrixSpecializationMacro(3);
|
||
|
|
vtkLastRowMultiplyMatrixSpecializationMacro(4);
|
||
|
|
|
||
|
|
#undef vtkLastRowMultiplyMatrixSpecializationMacro
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Class computing the determinant of square matrices SizeT x SizeT.
|
||
|
|
// The template parameter LayoutT is a struct embedded in MatrixLayout
|
||
|
|
template <int SizeT, class LayoutT = vtkMatrixUtilities::Layout::Identity>
|
||
|
|
class Determinant;
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for diagonal 3x3 matrices of size
|
||
|
|
template <>
|
||
|
|
class Determinant<3, vtkMatrixUtilities::Layout::Diag>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT>
|
||
|
|
static typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type Compute(
|
||
|
|
const MatrixT& M)
|
||
|
|
{
|
||
|
|
return M[0] * M[1] * M[3];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for diagonal 2x2 matrices
|
||
|
|
template <>
|
||
|
|
class Determinant<2, vtkMatrixUtilities::Layout::Diag>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT>
|
||
|
|
static typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type Compute(
|
||
|
|
const MatrixT& M)
|
||
|
|
{
|
||
|
|
return M[0] * M[1];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 1x1 matrices
|
||
|
|
template <class LayoutT>
|
||
|
|
class Determinant<1, LayoutT>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT>
|
||
|
|
static typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type Compute(
|
||
|
|
const MatrixT& M)
|
||
|
|
{
|
||
|
|
return M[0];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 2x2 non-diagonal matrices
|
||
|
|
template <class LayoutT>
|
||
|
|
class Determinant<2, LayoutT>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT>
|
||
|
|
static typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type Compute(
|
||
|
|
const MatrixT& M)
|
||
|
|
{
|
||
|
|
using Wrap = vtkMatrixUtilities::Wrapper<2, 2, MatrixT, LayoutT>;
|
||
|
|
|
||
|
|
return Wrap::template Get<0, 0>(M) * Wrap::template Get<1, 1>(M) -
|
||
|
|
Wrap::template Get<1, 0>(M) * Wrap::template Get<0, 1>(M);
|
||
|
|
}
|
||
|
|
|
||
|
|
private:
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 3x3 non-diagonal matrices
|
||
|
|
template <class LayoutT>
|
||
|
|
class Determinant<3, LayoutT>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT>
|
||
|
|
static typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type Compute(
|
||
|
|
const MatrixT& M)
|
||
|
|
{
|
||
|
|
using Wrap = vtkMatrixUtilities::Wrapper<3, 3, MatrixT, LayoutT>;
|
||
|
|
|
||
|
|
return Wrap::template Get<0, 0>(M) * Wrap::template Get<1, 1>(M) * Wrap::template Get<2, 2>(M) +
|
||
|
|
Wrap::template Get<0, 1>(M) * Wrap::template Get<1, 2>(M) * Wrap::template Get<2, 0>(M) +
|
||
|
|
Wrap::template Get<0, 2>(M) * Wrap::template Get<1, 0>(M) * Wrap::template Get<2, 1>(M) -
|
||
|
|
Wrap::template Get<0, 0>(M) * Wrap::template Get<1, 2>(M) * Wrap::template Get<2, 1>(M) -
|
||
|
|
Wrap::template Get<0, 1>(M) * Wrap::template Get<1, 0>(M) * Wrap::template Get<2, 2>(M) -
|
||
|
|
Wrap::template Get<0, 2>(M) * Wrap::template Get<1, 1>(M) * Wrap::template Get<2, 0>(M);
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Class inverting square matrices SizeT x SizeT.
|
||
|
|
// The template parameter LayoutT is a struct embedded in MatrixLayout
|
||
|
|
template <int SizeT, class LayoutT = vtkMatrixUtilities::Layout::Identity>
|
||
|
|
class InvertMatrix;
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 2x2 non-diagonal matrices
|
||
|
|
template <class LayoutT>
|
||
|
|
class InvertMatrix<2, LayoutT>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static void Compute(const MatrixT1& M1, MatrixT2& M2)
|
||
|
|
{
|
||
|
|
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT2>::value_type;
|
||
|
|
using Wrap1 = vtkMatrixUtilities::Wrapper<2, 2, MatrixT1, LayoutT>;
|
||
|
|
using Wrap2 = vtkMatrixUtilities::Wrapper<2, 2, MatrixT2>;
|
||
|
|
|
||
|
|
Scalar detInv = 1.0 / Determinant<2, LayoutT>::Compute(M1);
|
||
|
|
Wrap2::template Get<0, 0>(M2) = detInv * Wrap1::template Get<1, 1>(M1);
|
||
|
|
Wrap2::template Get<1, 0>(M2) = -detInv * Wrap1::template Get<1, 0>(M1);
|
||
|
|
Wrap2::template Get<0, 1>(M2) = -detInv * Wrap1::template Get<0, 1>(M1);
|
||
|
|
Wrap2::template Get<1, 1>(M2) = detInv * Wrap1::template Get<0, 0>(M1);
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 3x3 non-diagonal matrices
|
||
|
|
template <class LayoutT>
|
||
|
|
class InvertMatrix<3, LayoutT>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static void Compute(const MatrixT1& M1, MatrixT2& M2)
|
||
|
|
{
|
||
|
|
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT2>::value_type;
|
||
|
|
using Wrap1 = vtkMatrixUtilities::Wrapper<3, 3, MatrixT1, LayoutT>;
|
||
|
|
using Wrap2 = vtkMatrixUtilities::Wrapper<3, 3, MatrixT2>;
|
||
|
|
|
||
|
|
Scalar detInv = 1.0 /
|
||
|
|
(Wrap1::template Get<0, 0>(M1) *
|
||
|
|
(Wrap1::template Get<1, 1>(M1) * Wrap1::template Get<2, 2>(M1) -
|
||
|
|
Wrap1::template Get<2, 1>(M1) * Wrap1::template Get<1, 2>(M1)) -
|
||
|
|
Wrap1::template Get<0, 1>(M1) *
|
||
|
|
(Wrap1::template Get<1, 0>(M1) * Wrap1::template Get<2, 2>(M1) -
|
||
|
|
Wrap1::template Get<2, 0>(M1) * Wrap1::template Get<1, 2>(M1)) +
|
||
|
|
Wrap1::template Get<0, 2>(M1) *
|
||
|
|
(Wrap1::template Get<1, 0>(M1) * Wrap1::template Get<2, 1>(M1) -
|
||
|
|
Wrap1::template Get<2, 0>(M1) * Wrap1::template Get<1, 1>(M1)));
|
||
|
|
|
||
|
|
Wrap2::template Get<0, 0>(M2) = detInv *
|
||
|
|
(Wrap1::template Get<1, 1>(M1) * Wrap1::template Get<2, 2>(M1) -
|
||
|
|
Wrap1::template Get<2, 1>(M1) * Wrap1::template Get<1, 2>(M1));
|
||
|
|
Wrap2::template Get<1, 0>(M2) = -detInv *
|
||
|
|
(Wrap1::template Get<1, 0>(M1) * Wrap1::template Get<2, 2>(M1) -
|
||
|
|
Wrap1::template Get<2, 0>(M1) * Wrap1::template Get<1, 2>(M1));
|
||
|
|
Wrap2::template Get<2, 0>(M2) = detInv *
|
||
|
|
(Wrap1::template Get<1, 0>(M1) * Wrap1::template Get<2, 1>(M1) -
|
||
|
|
Wrap1::template Get<2, 0>(M1) * Wrap1::template Get<1, 1>(M1));
|
||
|
|
Wrap2::template Get<0, 1>(M2) = -detInv *
|
||
|
|
(Wrap1::template Get<0, 1>(M1) * Wrap1::template Get<2, 2>(M1) -
|
||
|
|
Wrap1::template Get<2, 1>(M1) * Wrap1::template Get<0, 2>(M1));
|
||
|
|
Wrap2::template Get<1, 1>(M2) = detInv *
|
||
|
|
(Wrap1::template Get<0, 0>(M1) * Wrap1::template Get<2, 2>(M1) -
|
||
|
|
Wrap1::template Get<2, 0>(M1) * Wrap1::template Get<0, 2>(M1));
|
||
|
|
Wrap2::template Get<2, 1>(M2) = -detInv *
|
||
|
|
(Wrap1::template Get<0, 0>(M1) * Wrap1::template Get<2, 1>(M1) -
|
||
|
|
Wrap1::template Get<2, 0>(M1) * Wrap1::template Get<0, 1>(M1));
|
||
|
|
Wrap2::template Get<0, 2>(M2) = detInv *
|
||
|
|
(Wrap1::template Get<0, 1>(M1) * Wrap1::template Get<1, 2>(M1) -
|
||
|
|
Wrap1::template Get<1, 1>(M1) * Wrap1::template Get<0, 2>(M1));
|
||
|
|
Wrap2::template Get<1, 2>(M2) = -detInv *
|
||
|
|
(Wrap1::template Get<0, 0>(M1) * Wrap1::template Get<1, 2>(M1) -
|
||
|
|
Wrap1::template Get<1, 0>(M1) * Wrap1::template Get<0, 2>(M1));
|
||
|
|
Wrap2::template Get<2, 2>(M2) = detInv *
|
||
|
|
(Wrap1::template Get<0, 0>(M1) * Wrap1::template Get<1, 1>(M1) -
|
||
|
|
Wrap1::template Get<1, 0>(M1) * Wrap1::template Get<0, 1>(M1));
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 1x1 matrices
|
||
|
|
template <class LayoutT>
|
||
|
|
class InvertMatrix<1, LayoutT>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static void Compute(const MatrixT1& M1, MatrixT2& M2)
|
||
|
|
{
|
||
|
|
M2[0] = 1.0 / M1[0];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 2x2 diagonal matrices
|
||
|
|
template <>
|
||
|
|
class InvertMatrix<2, vtkMatrixUtilities::Layout::Diag>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static void Compute(const MatrixT1& M1, MatrixT2& M2)
|
||
|
|
{
|
||
|
|
M2[0] = 1.0 / M1[0];
|
||
|
|
M2[1] = 1.0 / M1[1];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 3x3 diagonal matrices
|
||
|
|
template <>
|
||
|
|
class InvertMatrix<3, vtkMatrixUtilities::Layout::Diag>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT1, class MatrixT2>
|
||
|
|
static void Compute(const MatrixT1& M1, MatrixT2& M2)
|
||
|
|
{
|
||
|
|
M2[0] = 1.0 / M1[0];
|
||
|
|
M2[1] = 1.0 / M1[1];
|
||
|
|
M2[2] = 1.0 / M1[2];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Class solving systems M*y = x for square matrices RowsT x ColsT.
|
||
|
|
// The template parameter LayoutT is a struct embedded in MatrixLayout
|
||
|
|
// This class is currently specialized for 1x1, 2x2 and 3x3 matrices
|
||
|
|
template <int RowsT, int ColsT, class LayoutT = vtkMatrixUtilities::Layout::Identity>
|
||
|
|
class LinearSolve;
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 1x1 matrices
|
||
|
|
template <class LayoutT>
|
||
|
|
class LinearSolve<1, 1, LayoutT>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT, class VectorT1, class VectorT2>
|
||
|
|
static void Compute(const MatrixT& M, const VectorT1& x, VectorT2& y)
|
||
|
|
{
|
||
|
|
y[0] = x[0] / M[0];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 2x2 matrices
|
||
|
|
template <class LayoutT>
|
||
|
|
class LinearSolve<2, 2, LayoutT>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT, class VectorT1, class VectorT2>
|
||
|
|
static void Compute(const MatrixT& M, const VectorT1& x, VectorT2& y)
|
||
|
|
{
|
||
|
|
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type;
|
||
|
|
using Wrap = vtkMatrixUtilities::Wrapper<2, 2, MatrixT, LayoutT>;
|
||
|
|
|
||
|
|
Scalar detInv = 1.0 / Determinant<2, LayoutT>::Compute(M);
|
||
|
|
y[0] = (x[0] * Wrap::template Get<1, 1>(M) - x[1] * Wrap::template Get<0, 1>(M)) * detInv;
|
||
|
|
y[1] = (-x[0] * Wrap::template Get<1, 0>(M) + x[1] * Wrap::template Get<0, 0>(M)) * detInv;
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 3x3 matrices
|
||
|
|
template <class LayoutT>
|
||
|
|
class LinearSolve<3, 3, LayoutT>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT, class VectorT1, class VectorT2>
|
||
|
|
static void Compute(const MatrixT& M, const VectorT1& x, VectorT2& y)
|
||
|
|
{
|
||
|
|
using Scalar = typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type;
|
||
|
|
using Wrap = vtkMatrixUtilities::Wrapper<3, 3, MatrixT, LayoutT>;
|
||
|
|
|
||
|
|
Scalar detInv = 1.0 /
|
||
|
|
(Wrap::template Get<0, 0>(M) *
|
||
|
|
(Wrap::template Get<1, 1>(M) * Wrap::template Get<2, 2>(M) -
|
||
|
|
Wrap::template Get<2, 1>(M) * Wrap::template Get<1, 2>(M)) -
|
||
|
|
Wrap::template Get<0, 1>(M) *
|
||
|
|
(Wrap::template Get<1, 0>(M) * Wrap::template Get<2, 2>(M) -
|
||
|
|
Wrap::template Get<2, 0>(M) * Wrap::template Get<1, 2>(M)) +
|
||
|
|
Wrap::template Get<0, 2>(M) *
|
||
|
|
(Wrap::template Get<1, 0>(M) * Wrap::template Get<2, 1>(M) -
|
||
|
|
Wrap::template Get<2, 0>(M) * Wrap::template Get<1, 1>(M)));
|
||
|
|
|
||
|
|
y[0] = detInv *
|
||
|
|
(x[0] *
|
||
|
|
(Wrap::template Get<1, 1>(M) * Wrap::template Get<2, 2>(M) -
|
||
|
|
Wrap::template Get<2, 1>(M) * Wrap::template Get<1, 2>(M)) -
|
||
|
|
x[1] *
|
||
|
|
(Wrap::template Get<0, 1>(M) * Wrap::template Get<2, 2>(M) -
|
||
|
|
Wrap::template Get<2, 1>(M) * Wrap::template Get<0, 2>(M)) +
|
||
|
|
x[2] *
|
||
|
|
(Wrap::template Get<0, 1>(M) * Wrap::template Get<1, 2>(M) -
|
||
|
|
Wrap::template Get<1, 1>(M) * Wrap::template Get<0, 2>(M)));
|
||
|
|
y[1] = detInv *
|
||
|
|
(-x[0] *
|
||
|
|
(Wrap::template Get<1, 0>(M) * Wrap::template Get<2, 2>(M) -
|
||
|
|
Wrap::template Get<2, 0>(M) * Wrap::template Get<1, 2>(M)) +
|
||
|
|
x[1] *
|
||
|
|
(Wrap::template Get<0, 0>(M) * Wrap::template Get<2, 2>(M) -
|
||
|
|
Wrap::template Get<2, 0>(M) * Wrap::template Get<0, 2>(M)) -
|
||
|
|
x[2] *
|
||
|
|
(Wrap::template Get<0, 0>(M) * Wrap::template Get<1, 2>(M) -
|
||
|
|
Wrap::template Get<1, 0>(M) * Wrap::template Get<0, 2>(M)));
|
||
|
|
y[2] = detInv *
|
||
|
|
(x[0] *
|
||
|
|
(Wrap::template Get<1, 0>(M) * Wrap::template Get<2, 1>(M) -
|
||
|
|
Wrap::template Get<2, 0>(M) * Wrap::template Get<1, 1>(M)) -
|
||
|
|
x[1] *
|
||
|
|
(Wrap::template Get<0, 0>(M) * Wrap::template Get<2, 1>(M) -
|
||
|
|
Wrap::template Get<2, 0>(M) * Wrap::template Get<0, 1>(M)) +
|
||
|
|
x[2] *
|
||
|
|
(Wrap::template Get<0, 0>(M) * Wrap::template Get<1, 1>(M) -
|
||
|
|
Wrap::template Get<1, 0>(M) * Wrap::template Get<0, 1>(M)));
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 2x2 diagonal matrices
|
||
|
|
template <>
|
||
|
|
class LinearSolve<2, 2, vtkMatrixUtilities::Layout::Diag>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT, class VectorT1, class VectorT2>
|
||
|
|
static void Compute(const MatrixT& M, const VectorT1& x, VectorT2& y)
|
||
|
|
{
|
||
|
|
y[0] = x[0] / M[0];
|
||
|
|
y[1] = x[1] / M[1];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
|
||
|
|
//=============================================================================
|
||
|
|
// Specialization for 3x3 diagonal matrices
|
||
|
|
template <>
|
||
|
|
class LinearSolve<3, 3, vtkMatrixUtilities::Layout::Diag>
|
||
|
|
{
|
||
|
|
public:
|
||
|
|
template <class MatrixT, class VectorT1, class VectorT2>
|
||
|
|
static void Compute(const MatrixT& M, const VectorT1& x, VectorT2& y)
|
||
|
|
{
|
||
|
|
y[0] = x[0] / M[0];
|
||
|
|
y[1] = x[1] / M[1];
|
||
|
|
y[2] = x[2] / M[2];
|
||
|
|
}
|
||
|
|
};
|
||
|
|
} // namespace vtkMathPrivate
|
||
|
|
#endif
|