Proposal: meta-programmed matrix operations?

Development
1

Hi,

I recently needed a few matrix operations for a filter where I am using ranges. In order to be able to feed those objects in vtkMath, there needs to be templated versions of each operation. It is pretty straightforward to do so with most methods, like vtkMath::Dot for instance.

I think that this would be the occasion to do loop-free meta-programmed matrix operations, like vtkMath::MultiplyMatrix for instance, templated on the matrices dimensions and on the parameters as well as the matrix / vector parameters.

Is anyone against such a feature?

2

Could you please provide examples that show how things are done now and how could be done with the new templates?

3

Right now, you can do matrix multiplication this way:

// A = 1 0
//     0 1
const double A[2][2] = { { 1, 0 }, { 1, 0} };
// B = 1 0 0
//     0 1 0
const double B[2][3] = { { 1, 0, 0 }, { 0, 1, 0 } };

double C[2][3]
vtkMath::MultiplyMatrix(A, B, 2 /* rowA*/, 2 /*colA*/, 2 /*rowB*/, 3 /*colB*/, C)

You need to put double** or float** as inputs and output. Now, using ranges, one could do that:

vtkNew<vtkDoubleArray> arrayOfMatricesA, arrayOfMatricesB, arrayOfMatricesC;
/**
 * Depending on the matrix dimensions, fill the arrays.
 * One could fill matrices row by row.
 */
// rangeA is a range of 2x2 matrices
auto rangeA = vtk::DataArrayTupleRange<4>(arrayOfMatricesA)
// rangeB is a range of 2x3 matrices
auto rangeB = vtk::DataArrayTupleRange<6>(arrayOfMatricesB);
// rangeC is a range of 2x3 matrices
auto rangeC = vtk::DataArrayTupleRange<6>(arrayOfMatricesC);
using Scalar = typename decltype<rangeC>::value_type;

for (vtkIdType id = 0; id < arrayOfMatricesA->GetNumberOfTuples(); ++id)
{
  // transposeA and transposeB would be false by default, but I put them for completeness
  vtkMath::MultiplyMatrix<Scalar, 2 /*rowA*/, 2 /*colA == rowB*/, 3 /*colB*/,
                          false /*transposedA*/, false /*transposedB*/>
                         (rangeA[id], rangeB[id], rangeC[id]);
}

You could use the first example as well by feeding *A, *B, and *C in this templated version.

4

I just realized that Scalar can be induced with meta-programming, so I think it would not be necessary.

5

I implemented matrix multiplications in vtkMath in master. Look for MultiplyMatrix, MultiplyMatrixWithVector, and Dot.

Matrices are 1D vector that get rearranged in row-column fashion depending on the template parameters.

You can change the layout of matrices at compile-time using some template flags. 3 layouts are supported:

  • vtkMath::MatrixLayout::Identity - matrix is read as is
  • vtkMath::MatrixLayout::Transpose - the transposed matrix is read when multiplying
  • vtkMath::MatrixLayout::Diag - the matrix is a diagonal matrix, i.e. values of the input array are values of the diagonal.