c# - 快速旋转/变换矩阵乘法

Question

我正在寻找在 C# 中为 3x3 旋转和 4x4 变换矩阵执行 matrix * matrix 和 matrix * vector 操作的最有效方法。

我目前将矩阵存储在多维数组中(new double[3,3]、new double[4,4])。我并不完全反对改变它，但如果可能的话我想保留语法。我当前使用 3 个标准嵌套 for 循环的乘法工作正常，但可能成为瓶颈。

到目前为止我的想法:

• 像 Strassen 这样的优化算法不适用于这些尺寸
• 在单个 4x4 乘法级别上，并行化没有多大意义；在更高的层次上做得更好。
• 由于边界检查效率较低，多维数组在 C# 中(曾经？)速度较慢，但​​这可以通过不安全的指针算法来克服。 (我不确定此信息的最新情况)
• 旋转矩阵是对称的，可能有办法利用它吗？
• 最大的 yield 可能是通过使用缓存局部性来实现的，确保内存中靠近的值被一起访问；但我不确定该怎么做。

因此，在我使用不安全、固定和 3 个 for 循环组合我自己的解决方案之前，是否已经有针对这个标准问题的经过测试和优化的解决方案？

还是有其他我忽略的优化？

Answer 1

这就是我使用的，它的运行速度出奇的快。

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));
    }
}

其中 vec3 和 mat3 是我自己的 Vector3 和 Matrix3 存储元素的结构的别名是字段.对于 4 元素结构也是如此。我也像这样对它进行逆编码:

    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));
    }