是否有更好的方法来编写 SIMD 代码来反转变换矩阵？

Question

我正在为光线追踪器项目编写一个数学库，因此我试图使我的繁重操作（如矩阵求逆）更加优化。经过一些研究后，我发现了这个反转变换矩阵的技巧（如果有人感兴趣的话，我在我制作的这份文档中有更详细的描述：https://docs.google.com/document/d/1ok8dzMk7EZiZaVRB61zGDlxRSDoelX3z6ixaCRlg0yM/edit），其中我将我的变换矩阵存储为三个不同的矩阵（缩放、旋转和平移），然后我只需反转它们中的每一个，然后组合矩阵以获得最终结果，即逆矩阵。这是代码（也可以在 godbolt 上找到）： https://godbolt.org/z/cjsfGzW3c

#include <immintrin.h>

typedef union u_vec4s
{
    float       a[4];
    __m128      simd;
    struct
    {
        float   x;
        float   y;
        float   z;
        float   w;
    };
}__attribute((aligned(16))) t_vec4s;

typedef union u_mat4s
{
    float       a[4][4];
    __m128      simd[4];
    __m256      _ymm[2];
    struct
    {
        t_vec4s   r1;
        t_vec4s   r2;
        t_vec4s   r3;
        t_vec4s   r4;
    };
}__attribute((aligned(16))) t_mat4s;

t_mat4s lag_get_transform_matrix_inverse(const t_mat4s s, const t_mat4s r, const t_mat4s t)
{
    // const __m128 zeros = _mm_set1_ps(0);
    const __m128 mul00 = _mm_set_ps(1, -t.r3.w, -t.r2.w, -t.r1.w);
    t_mat4s t1w0;
    t_mat4s s1rr;

    __m128 tmp0 = _mm_unpacklo_ps(r.simd[0], r.simd[1]); // [r00, r10, r01, r11]
    __m128 tmp1 = _mm_unpackhi_ps(r.simd[0], r.simd[1]); // [r02, r12, r03, r13]
    __m128 tmp2 = _mm_unpacklo_ps(r.simd[2], r.simd[3]); // [r20, r30, r21, r31]
    __m128 tmp3 = _mm_unpackhi_ps(r.simd[2], r.simd[3]); // [r22, r32, r23, r33]

    s1rr.simd[0] = _mm_mul_ps(_mm_movelh_ps(tmp0, tmp2), _mm_set1_ps(1.f / s.r1.x)); // [r00/sx, r10/sx, r20/sx, r30/sx]
    s1rr.simd[1] = _mm_mul_ps(_mm_movehl_ps(tmp2, tmp0), _mm_set1_ps(1.f / s.r2.y)); // [r01/sy, r11/sy, r21/sy, r31/sy]
    s1rr.simd[2] = _mm_mul_ps(_mm_movelh_ps(tmp1, tmp3), _mm_set1_ps(1.f / s.r3.z)); // [r02/sz, r12/sz, r22/sz, r32/sz]
    s1rr.simd[3] = _mm_movehl_ps(tmp3, tmp1); // [0, 0, 0, 1]

    t1w0.simd[0] = /*_mm_sub_ps(zeros, */_mm_dp_ps(s1rr.simd[0], mul00, 0xF1)/*)*/;
    t1w0.simd[1] = /*_mm_sub_ps(zeros, */_mm_dp_ps(s1rr.simd[1], mul00, 0xF1)/*)*/;
    t1w0.simd[2] = /*_mm_sub_ps(zeros, */_mm_dp_ps(s1rr.simd[2], mul00, 0xF1)/*)*/;

    s1rr.r1.w = _mm_cvtss_f32(t1w0.simd[0]);
    s1rr.r2.w = _mm_cvtss_f32(t1w0.simd[1]);
    s1rr.r3.w = _mm_cvtss_f32(t1w0.simd[2]);

    return (s1rr);
}

现在，我的问题是：有更好的方法吗？也许有更有效的方法可以达到相同的结果？或者也许我正在做一些不需要的操作。我知道我可以做的一件事是使其更快，是以列顺序而不是行顺序存储旋转矩阵，这样就不必转置。另一件事可能是避免使用

_mm_dp_ps

，因为该功能有点繁重，最好只使用水平添加。比如：

    t1w0.simd[0] = _mm_mul_ps(s1rr.simd[0], mul00);
    _mm_storeu_ps((float *)&t1w0.r1, _mm_hadd_ps(t1w0.simd[0], t1w0.simd[0]));
    s1rr.r1.w = t1w0.r1.a[0] + t1w0.r1.a[2];

这是值得做的事情吗？有什么想法或建议吗？请随意评论我之前提到的所有内容！任何建议将不胜感激。预先感谢<3

Answer 1

您的 3 个输入矩阵使用 192 字节内存或 12 个 SSE 向量。然而，它们只保留 15 个浮点数 = 60 字节的数据，其余元素为零。这是存储分解转换并将其传递给函数的非常低效的方法。

如果你确定你知道你在做什么，即你需要反转一大堆不同的矩阵，分析器显示反转是瓶颈，你已经使用 DirectXMath 的 XMMatrixInverse 测试了你的东西（链接的存储库使用复制粘贴友好的 MIT 许可证）并且对性能不满意，请重构您的函数以在单个向量（而不是 4x4 矩阵）中进行旋转和平移。

然后做类似的事情。以下代码是用 C++ 编写的，未经测试。

#include <immintrin.h>

// A row major 3x3 or 4x3 matrix stored in 3 SSE vectors
struct Mat3
{
    __m128 rows[ 3 ];
};

inline __m128 broadcastX( __m128 v ) { return _mm_permute_ps( v, _MM_SHUFFLE( 0, 0, 0, 0 ) ); }
inline __m128 broadcastY( __m128 v ) { return _mm_permute_ps( v, _MM_SHUFFLE( 1, 1, 1, 1 ) ); }
inline __m128 broadcastZ( __m128 v ) { return _mm_permute_ps( v, _MM_SHUFFLE( 2, 2, 2, 2 ) ); }

Mat3 invertDecomposed( Mat3 rot, __m128 scale, __m128 translation )
{
    // Invert the scale vector
    __m128 scaleInv = _mm_div_ps( _mm_set1_ps( 1 ), scale );
    // Negate the translation vector
    translation = _mm_sub_ps( _mm_setzero_ps(), translation );

    // Multiply rows of the rotation matrix by components of the inverted scale
    __m128 r0 = _mm_mul_ps( rot.rows[ 0 ], broadcastX( scaleInv ) );
    __m128 r1 = _mm_mul_ps( rot.rows[ 1 ], broadcastY( scaleInv ) );
    __m128 r2 = _mm_mul_ps( rot.rows[ 2 ], broadcastZ( scaleInv ) );

    // Compute last column of the output matrix
    // Note VDPPS can place output number into any lanes, we need it in W
    __m128 r03 = _mm_dp_ps( r0, translation, 0b01111000 );
    __m128 r13 = _mm_dp_ps( r1, translation, 0b01111000 );
    __m128 r23 = _mm_dp_ps( r2, translation, 0b01111000 );

    // Insert numbers into the last column
    r0 = _mm_blend_ps( r0, r03, 0b1000 );
    r1 = _mm_blend_ps( r1, r13, 0b1000 );
    r2 = _mm_blend_ps( r2, r23, 0b1000 );

    // Return the result
    Mat3 result;
    result.rows[ 0 ] = r0;
    result.rows[ 1 ] = r1;
    result.rows[ 2 ] = r2;
    return result;
}

顺便说一句，如果你能保证输入旋转矩阵的W分量为零，你可以用

_mm_blend_ps

替换

_mm_or_ps

，但是

_mm_blend_ps

非常快，我不认为这会带来可衡量的利润。