Fast rotation/transformation matrix multiplications
.net, c#, matrix, matrix-multiplication, performance
Solution
This is what I use and it works surprisingly fast.
public struct Matrix3
{
public readonly double a11, a12, a13;
public readonly double a21, a22, a23;
public readonly double a31, a32, a33;
...
public vec3 Multiply(vec3 rhs)
{
// y= A*x
// fill vector by element
return new vec3(
(a11*rhs.X+a12*rhs.Y+a13*rhs.Z),
(a21*rhs.X+a22*rhs.Y+a23*rhs.Z),
(a31*rhs.X+a32*rhs.Y+a33*rhs.Z));
}
public mat3 Multiply(mat3 rhs)
{
// Y = A*X
// fill matrix by row
return new mat3(
(a11*rhs.a11+a12*rhs.a21+a13*rhs.a31),
(a11*rhs.a12+a12*rhs.a22+a13*rhs.a32),
(a11*rhs.a13+a12*rhs.a23+a13*rhs.a33),
(a21*rhs.a11+a22*rhs.a21+a23*rhs.a31),
(a21*rhs.a12+a22*rhs.a22+a23*rhs.a32),
(a21*rhs.a13+a22*rhs.a23+a23*rhs.a33),
(a31*rhs.a11+a32*rhs.a21+a33*rhs.a31),
(a31*rhs.a12+a32*rhs.a22+a33*rhs.a32),
(a31*rhs.a13+a32*rhs.a23+a33*rhs.a33));
}
}
where `vec3` and `mat3` are aliases to my own `Vector3` and `Matrix3` structures that store the elements are fields. Similarly for 4 element structures. Also I have coded it the inverses like this:
public double Determinant()
{
return a11*(a22*a33-a23*a32)
+a12*(a23*a31-a21*a33)
+a13*(a21*a32-a22*a31);
}
/// <summary>
/// Solves the system of equations this*x=rhs for x
/// </summary>
public vec3 Solve(vec3 rhs)
{
double D=Determinant();
double ID=1/D;
return new vec3(
(((a22*a33-a23*a32)*rhs.X+(a13*a32-a12*a33)*rhs.Y+(a12*a23-a13*a22)*rhs.Z)*ID),
-(((a21*a33-a23*a31)*rhs.X+(a13*a31-a11*a33)*rhs.Y+(a11*a23-a13*a21)*rhs.Z)*ID),
(((a21*a32-a22*a31)*rhs.X+(a12*a31-a11*a32)*rhs.Y+(a11*a22-a12*a21)*rhs.Z)*ID));
}
/// <summary>
/// Solves the system of equations this*X = rhs for X
/// </summary>
public mat3 Solve(mat3 rhs)
{
double D=Determinant();
double ID=1/D;
return new mat3(
(((a22*a33-a23*a32)*rhs.a11+(a13*a32-a12*a33)*rhs.a21+(a12*a23-a13*a22)*rhs.a31)*ID),
(((a22*a33-a23*a32)*rhs.a12+(a13*a32-a12*a33)*rhs.a22+(a12*a23-a13*a22)*rhs.a32)*ID),
(((a22*a33-a23*a32)*rhs.a13+(a13*a32-a12*a33)*rhs.a23+(a12*a23-a13*a22)*rhs.a33)*ID),
-(((a21*a33-a23*a31)*rhs.a11+(a13*a31-a11*a33)*rhs.a21+(a11*a23-a13*a21)*rhs.a31)*ID),
-(((a21*a33-a23*a31)*rhs.a12+(a13*a31-a11*a33)*rhs.a22+(a11*a23-a13*a21)*rhs.a32)*ID),
-(((a21*a33-a23*a31)*rhs.a13+(a13*a31-a11*a33)*rhs.a23+(a11*a23-a13*a21)*rhs.a33)*ID),
(((a21*a32-a22*a31)*rhs.a11+(a12*a31-a11*a32)*rhs.a21+(a11*a22-a12*a21)*rhs.a31)*ID),
(((a21*a32-a22*a31)*rhs.a12+(a12*a31-a11*a32)*rhs.a22+(a11*a22-a12*a21)*rhs.a32)*ID),
(((a21*a32-a22*a31)*rhs.a13+(a12*a31-a11*a32)*rhs.a23+(a11*a22-a12*a21)*rhs.a33)*ID));
}
Problem
I am looking for the most efficient way to do `matrix * matrix` and `matrix * vector` operations for 3x3 rotation and 4x4 transformation matrices in C#. I currently store my matrices in multidimensional arrays (`new double[3,3]`, `new double[4,4]`). I am not totally adverse to changing that but if possible I would like to keep the syntax. My current multiplication using the 3 standard nested for loops works fine but can be a bottleneck. My thoughts so far: - Optimized algorithms like Strassen are not practical for these sizes - Parallelisation does not make much sense either at the level of a single 4x4 multiplication; better done at a higher level. - multidimensional arrays are (were?) slower in c# due to less efficient boundary checks, however this can be overcome with unsafe pointer arithmetic. (I am not sure how current this information is) - rotation matrices are symmetric, there might be a way to exploit that? - the biggest gains can probably be achieved by using cache-locality, making sure that values that are close together in memory are accessed together; but I am unsure how to do this. So before I hack together my own solution using unsafe, fixed and 3 for loops, is there already a tested and optimized solution for this standard problem out there? Or are there other optimizations that I have overlooked?