効率的な 4x4 行列乗算 (C とアセンブリ)

C で 2 つの 4x4 行列を乗算するためのより高速でトリッキーな方法を探しています。現在の研究は、SIMD 拡張を使用した x86-64 アセンブリに焦点を当てています。これまでのところ、単純な C 実装よりも約 6 倍高速な関数 witch を作成しました。これは、パフォーマンスの向上に対する私の期待を上回りました。残念ながら、これはコンパイルに最適化フラグが使用されていない場合にのみ当てはまります (GCC 4.7)。 -O2 で、C が速くなり、私の努力は無意味になります。

最新のコンパイラは、複雑な最適化手法を利用してほぼ完璧なコードを作成し、通常は手作りの精巧なアセンブリよりも高速であることを私は知っています。しかし、パフォーマンスが重要な少数のケースでは、人間がコンパイラでクロックサイクルを争おうとすることがあります。特に、最新の ISA に裏打ちされたいくつかの数学を調査できる場合 (私の場合のように)。

私の関数は次のようになります (AT&T 構文、GNU アセンブラー):

    .text
     .globl matrixMultiplyASM
     .type matrixMultiplyASM, @function
 matrixMultiplyASM:
     movaps   (%rdi), %xmm0    # fetch the first matrix (use four registers)
     movaps 16(%rdi), %xmm1
     movaps 32(%rdi), %xmm2
     movaps 48(%rdi), %xmm3
     xorq %rcx, %rcx           # reset (forward) loop iterator
 .ROW:
     movss (%rsi), %xmm4       # Compute four values (one row) in parallel:
     shufps $0x0, %xmm4, %xmm4 # 4x 4FP mul's, 3x 4FP add's 6x mov's per row,
     mulps %xmm0, %xmm4        # expressed in four sequences of 5 instructions,
     movaps %xmm4, %xmm5       # executed 4 times for 1 matrix multiplication.
     addq $0x4, %rsi
     movss (%rsi), %xmm4       # movss + shufps comprise _mm_set1_ps intrinsic
     shufps $0x0, %xmm4, %xmm4 #
     mulps %xmm1, %xmm4
     addps %xmm4, %xmm5
     addq $0x4, %rsi           # manual pointer arithmetic simplifies addressing
     movss (%rsi), %xmm4
     shufps $0x0, %xmm4, %xmm4
     mulps %xmm2, %xmm4        # actual computation happens here
     addps %xmm4, %xmm5        #
     addq $0x4, %rsi
     movss (%rsi), %xmm4       # one mulps operand fetched per sequence
     shufps $0x0, %xmm4, %xmm4 #  |
     mulps %xmm3, %xmm4        # the other is already waiting in %xmm[0-3]
     addps %xmm4, %xmm5
     addq $0x4, %rsi           # 5 preceding comments stride among the 4 blocks
     movaps %xmm5, (%rdx,%rcx) # store the resulting row, actually, a column
     addq $0x10, %rcx          # (matrices are stored in column-major order)
     cmpq $0x40, %rcx
     jne .ROW
     ret
 .size matrixMultiplyASM, .-matrixMultiplyASM

128 ビット SSE レジスタにパックされた 4 つの float を処理することにより、反復ごとに結果の行列の列全体を計算します。完全なベクトル化は、ちょっとした数学 (操作の並べ替えと集計) と mullps で可能です。 /addps 4xfloat パッケージの並列乗算/加算の手順。このコードは、パラメーターを渡すためのレジスターを再利用します (%rdi 、 %rsi 、 %rdx :GNU/Linux ABI) は、(内部) ループのアンローリングの恩恵を受け、1 つの行列全体を XMM レジスタに保持して、メモリの読み取りを減らします。ご覧のとおり、私はこのトピックを調査し、できる限り時間をかけて実装しました。

私のコードを征服する単純な C 計算は次のようになります:

void matrixMultiplyNormal(mat4_t *mat_a, mat4_t *mat_b, mat4_t *mat_r) {
     for (unsigned int i = 0; i < 16; i += 4)
         for (unsigned int j = 0; j < 4; ++j)
             mat_r->m[i + j] = (mat_b->m[i + 0] * mat_a->m[j +  0])
                             + (mat_b->m[i + 1] * mat_a->m[j +  4])
                             + (mat_b->m[i + 2] * mat_a->m[j +  8])
                             + (mat_b->m[i + 3] * mat_a->m[j + 12]);
 }

上記の C コードの最適化されたアセンブリ出力を調査しましたが、XMM レジスタに float を格納している間、並列操作は含まれていません – スカラー計算、ポインター演算、および条件付きジャンプのみ。コンパイラのコードは意図的ではないように見えますが、ベクトル化されたバージョンが約 4 倍高速であると予想されるよりもわずかに効果的です。一般的な考えは正しいと確信しています。プログラマーは同様のことを行い、やりがいのある結果をもたらします。しかし、ここで何が問題なのですか？私が認識していないレジスタ割り当てまたは命令スケジューリングの問題はありますか?マシンとの戦いをサポートする x86-64 アセンブリツールまたはトリックを知っていますか?

答え:

コードを高速化し、コンパイラーを凌駕する方法があります。洗練されたパイプライン分析やコードの詳細なマイクロ最適化は必要ありません (これは、これらのメリットが得られないという意味ではありません)。最適化には 3 つの簡単なトリックが使用されます:

関数は 32 バイトにアラインされました (これによりパフォーマンスが大幅に向上しました)、

メインループは逆になり、(EFLAGS に基づく) ゼロテストとの比較が減少します。

命令レベルのアドレス演算は、「外部」ポインター計算よりも高速であることが証明されました (ただし、「3/4 ケースでは」2 倍の加算が必要ですが)。これにより、ループ本体が 4 命令短縮され、実行パス内のデータ依存性が減少しました。関連する質問を参照してください。

さらに、コードは、シンボルの再定義エラーを抑制する相対ジャンプ構文を使用します。このエラーは、GCC がインライン展開しようとしたときに発生します (asm 内に配置された後)。 -O3 でコンパイルされたステートメント ).

    .text
     .align 32                           # 1. function entry alignment
     .globl matrixMultiplyASM            #    (for a faster call)
     .type matrixMultiplyASM, @function
 matrixMultiplyASM:
     movaps   (%rdi), %xmm0
     movaps 16(%rdi), %xmm1
     movaps 32(%rdi), %xmm2
     movaps 48(%rdi), %xmm3
     movq $48, %rcx                      # 2. loop reversal
 1:                                      #    (for simpler exit condition)
     movss (%rsi, %rcx), %xmm4           # 3. extended address operands
     shufps $0, %xmm4, %xmm4             #    (faster than pointer calculation)
     mulps %xmm0, %xmm4
     movaps %xmm4, %xmm5
     movss 4(%rsi, %rcx), %xmm4
     shufps $0, %xmm4, %xmm4
     mulps %xmm1, %xmm4
     addps %xmm4, %xmm5
     movss 8(%rsi, %rcx), %xmm4
     shufps $0, %xmm4, %xmm4
     mulps %xmm2, %xmm4
     addps %xmm4, %xmm5
     movss 12(%rsi, %rcx), %xmm4
     shufps $0, %xmm4, %xmm4
     mulps %xmm3, %xmm4
     addps %xmm4, %xmm5
     movaps %xmm5, (%rdx, %rcx)
     subq $16, %rcx                      # one 'sub' (vs 'add' & 'cmp')
     jge 1b                              # SF=OF, idiom: jump if positive
     ret

これは、これまでに見た中で最速の x86-64 実装です。私は感謝し、投票し、その目的のためのより迅速なアセンブリを提供する回答を受け入れます!