以下の内容はhttps://lpha-z.hatenablog.com/entry/2025/09/07/231500より取得しました。

AVX-512でSIMD化されたTwoSumの速度評価

TwoSumをAVX-512でSIMD化するという話題を聞いたので、速度調査をしてみました。

TwoSumとTwoSumCond

TwoSumは、二つの浮動小数点数aとbを受け取って、a+bおよびその丸め誤差を返す関数です。 （コーナーケースでの正しさの違いを除いて）以下の二種類の実装があります。

std::pair<double, double> TwoSum( double a, double b ) {
  double hi = a + b;
  double tmp = hi - a;
  double lo = a + (hi - tmp) + (b - tmp);
  return { hi, lo };
}

std::pair<double, double> TwoSumCond( double& a, double& b ) {
    double hi = a + b;

    double L = std::abs(a) < std::abs(b) ? b : a; // 絶対値の大きい方
    double S = std::abs(a) < std::abs(b) ? a : b; // 絶対値の小さい方

    double lo = L - hi + S; // FastTwoSum

    return { hi, lo };
}

実験コード

#include <iostream>
#include <chrono>
#include <array>
#include <random>

std::array<double, 4096> arr1;
std::array<double, 4096> arr2;

void TwoSum( double& a, double& b ) {
    double hi = a + b;
    double tmp = hi - a;
    double lo = a - (hi - tmp) + (b - tmp);

    a = hi;
    b = lo;
}

void TwoSumCond( double& a, double& b ) {
    double hi = a + b;

    double L = std::abs(a) < std::abs(b) ? b : a;
    double S = std::abs(a) < std::abs(b) ? a : b;

    double lo = L - hi + S;

    a = hi;
    b = lo;
}

void test1() {
    auto start = std::chrono::system_clock::now();
    for( int i = 0; i < 65536*32; ++i ) {
        for( int j = 0; j < 4096; ++j ) {
            TwoSum( arr1[j], arr2[j] );
        }
    }
    auto finish = std::chrono::system_clock::now();
    std::cout << "TwoSum, Compiler " << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() * 1e-9 << std::endl;
}

void test2() {
    auto start = std::chrono::system_clock::now();
    for( int i = 0; i < 65536*32; ++i ) {
        for( int j = 0; j < 4096; ++j ) {
            TwoSumCond( arr1[j], arr2[j] );
        }
    }
    auto finish = std::chrono::system_clock::now();
    std::cout << "TwoSumCond, Compiler " << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() * 1e-9 << std::endl;
}

#include <immintrin.h>

void test3() {
    auto start = std::chrono::system_clock::now();
    for( int i = 0; i < 65536*32; ++i ) {
        for( int j = 0; j < 4096; j += 8 ) {
            __m512d zmm1 = _mm512_load_pd( &arr1[j] );
            __m512d zmm2 = _mm512_load_pd( &arr2[j] );
            __m512d zmm_hi  = zmm1 + zmm2;
            __m512d zmm_tmp = zmm_hi - zmm1;
            __m512d zmm_lo = zmm1 - (zmm_hi - zmm_tmp) + (zmm2 - zmm_tmp);

            _mm512_store_pd( &arr1[j], zmm_hi );
            _mm512_store_pd( &arr2[j], zmm_lo );
        }
    }
    auto finish = std::chrono::system_clock::now();
    std::cout << "TwoSum, Intrinsic " << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() * 1e-9 << std::endl;
}

void test4() {
    auto start = std::chrono::system_clock::now();
    for( int i = 0; i < 65536*32; ++i ) {
        const __m512d sign_mask = _mm512_set1_pd(-0.0);
        for( int j = 0; j < 4096; j += 8 ) {
            __m512d a = _mm512_load_pd( &arr1[j] );
            __m512d b = _mm512_load_pd( &arr2[j] );

            __m512d absa = _mm512_andnot_pd(sign_mask, a);
            __m512d absb = _mm512_andnot_pd(sign_mask, b);

            __mmask8 m = _mm512_cmplt_pd_mask(absa, absb);

            __m512d zmm_hi = _mm512_add_pd(a, b);

            __m512d zmm_L = _mm512_mask_blend_pd(m, a, b);
            __m512d zmm_S = _mm512_mask_blend_pd(m, b, a);

            __m512d zmm_lo = _mm512_add_pd(_mm512_sub_pd(zmm_L, zmm_hi), zmm_S);

            _mm512_store_pd( &arr1[j], zmm_hi );
            _mm512_store_pd( &arr2[j], zmm_lo );
        }
    }
    auto finish = std::chrono::system_clock::now();
    std::cout << "TwoSumCond, Intrinsic (andnot_pd, 2blend) " << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() * 1e-9 << std::endl;
}

void test5() {
    auto start = std::chrono::system_clock::now();
    for( int i = 0; i < 65536*32; ++i ) {
        const __m512i sign_mask = _mm512_set1_epi64(0x8000000000000000ull);
        for( int j = 0; j < 4096; j += 8 ) {
            __m512d a = _mm512_load_pd( &arr1[j] );
            __m512d b = _mm512_load_pd( &arr2[j] );

            __m512d absa = _mm512_castsi512_pd(_mm512_andnot_epi64(sign_mask, _mm512_castpd_si512(a)));
            __m512d absb = _mm512_castsi512_pd(_mm512_andnot_epi64(sign_mask, _mm512_castpd_si512(b)));

            __mmask8 m = _mm512_cmplt_pd_mask(absa, absb);

            __m512d zmm_hi = _mm512_add_pd(a, b);

            __m512d zmm_L = _mm512_mask_blend_pd(m, a, b);
            __m512d zmm_S = _mm512_mask_blend_pd(m, b, a);

            __m512d zmm_lo = _mm512_add_pd(_mm512_sub_pd(zmm_L, zmm_hi), zmm_S);

            _mm512_store_pd( &arr1[j], zmm_hi );
            _mm512_store_pd( &arr2[j], zmm_lo );
        }
    }
    auto finish = std::chrono::system_clock::now();
    std::cout << "TwoSumCond, Intrinsic (andnot_epi64, 2blend) " << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() * 1e-9 << std::endl;
}

void test6() {
    auto start = std::chrono::system_clock::now();
    for( int i = 0; i < 65536*32; ++i ) {
        const __m512d sign_mask = _mm512_set1_pd(-0.0);
        for( int j = 0; j < 4096; j += 8 ) {
            __m512d a = _mm512_load_pd( &arr1[j] );
            __m512d b = _mm512_load_pd( &arr2[j] );

            __m512d absa = _mm512_andnot_pd(sign_mask, a);
            __m512d absb = _mm512_andnot_pd(sign_mask, b);

            __mmask8 m = _mm512_cmplt_pd_mask(absa, absb);

            __m512d zmm_hi = _mm512_add_pd(a, b);

            __m512d zmm_lo = _mm512_add_pd(_mm512_sub_pd(a, zmm_hi), b);
            zmm_lo = _mm512_mask_add_pd(zmm_lo, m, _mm512_sub_pd(b, zmm_hi), a);

            _mm512_store_pd( &arr1[j], zmm_hi );
            _mm512_store_pd( &arr2[j], zmm_lo );
        }
    }
    auto finish = std::chrono::system_clock::now();
    std::cout << "TwoSumCond, Intrinsic (andnot_pd, mask_add) " << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() * 1e-9 << std::endl;
}

void test7() {
    auto start = std::chrono::system_clock::now();
    for( int i = 0; i < 65536*32; ++i ) {
        const __m512i sign_mask = _mm512_set1_epi64(0x8000000000000000ull);
        for( int j = 0; j < 4096; j += 8 ) {
            __m512d a = _mm512_load_pd( &arr1[j] );
            __m512d b = _mm512_load_pd( &arr2[j] );

            __m512d absa = _mm512_castsi512_pd(_mm512_andnot_epi64(sign_mask, _mm512_castpd_si512(a)));
            __m512d absb = _mm512_castsi512_pd(_mm512_andnot_epi64(sign_mask, _mm512_castpd_si512(b)));

            __mmask8 m = _mm512_cmplt_pd_mask(absa, absb);

            __m512d zmm_hi = _mm512_add_pd(a, b);

            __m512d zmm_lo = _mm512_add_pd(_mm512_sub_pd(a, zmm_hi), b);
            zmm_lo = _mm512_mask_add_pd(zmm_lo, m, _mm512_sub_pd(b, zmm_hi), a);

            _mm512_store_pd( &arr1[j], zmm_hi );
            _mm512_store_pd( &arr2[j], zmm_lo );
        }
    }
    auto finish = std::chrono::system_clock::now();
    std::cout << "TwoSumCond, Intrinsic (andnot_epi64, mask_add) " << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() * 1e-9 << std::endl;
}

int main() {
    std::mt19937_64 mt;
    std::uniform_real_distribution<double> dist( -1.0, 1.0 );
    for( int i = 0; i < 4096; ++i ) { arr1[i] = dist( mt ); arr2[i] = dist( mt ); }

    test1();
    test2();
    test3();
    test4();
    test5();
    test6();
    test7();
}

test1はTwoSumを、test2はTwoSumCondを、それぞれコンパイラの自動ベクトル化でSIMD化するものです。 test3は、SIMD化されたTwoSumをイントリンシックを用いて記述したものです。 test4～test7は、ChatGPTにより生成された（ものに少し手を加えた）、SIMD化されたTwoSumCondをイントリンシックを用いて記述したコードです。 andnot_pd（AVX512-DQが必要らしい）を使うのがtest4・test6、それがなくても動くのがtest5・test7だそうです。 また(test4,test5)と(test6,test7)の違いは、loを計算する前に選択する（test2と同じ）か、a-hi+bとb-hi+aを両方計算してから1回だけ選択する*1かです。 前者はand×2, cmp×1, blend×2, add×3の合計8SIMD命令、後者はand×2, cmp×1, add×4, mask_add×1の合計8SIMD命令、になっています。

※subも全部addと書いています。以下同様に、加算と言ったらaddとsubを両方合わせたものを言います。

速度調査1：Alder Lake

Intel Core i9 12900Kの初期版（Eコアを無効化するとAVX-512が使える）を用いて評価しました。 使ったコンパイラはg++ 9.4.0で、最適化オプションは-O3 -march=nativeです。

実験結果は、test1とtest3（TwoSum）が0.7秒くらい、test2とtest4～7（TwoSumCond）が0.9秒くらいでした。test1はtest3より、test2はtest4～7より、それぞれ平均的に2%くらい遅かったです。 test1とtest2は一周当たり32要素（zmm使用で8要素、さらに4倍アンローリング）、test3～test7は一周当たり16要素（2倍アンローリング）のコードにコンパイルされていたので、その違いでしょうか。

Alder LakeのPコアであるGolden Coveの発行ポートは、uops.infoによれば

  • vaddpd(zmm)：p0/p5
  • vblendmpd(zmm)：p0/p5
  • vandnpd(zmm)：p0/p5
  • vcmpltpd(zmm)：p5

とのことなので、命令数の比がそのまま速度の比になるのが予想される結果です。 実験結果は、TwoSumの方がTwoSumCondよりも1.29倍高速、というものでした。 1.25倍が予期されるところなので、大体あっています。

速度調査2：Zen 4

AMD Ryzen 9 7950Xを用いて評価しました。 使ったコンパイラはg++ 13.3.0で、最適化オプションは -O3 -march=native -fno-loopinterchangeです。 -fno-loopinterchangeは、ループ交換（二重ループの入れ子の順番を入れ替える）最適化を抑制するために付けました。 ループ交換は、今回のベンチマークの意図を無視したものであるうえ、そもそも並列度が下がって却って速度低下するので、無効化しました。 今回のコードに対するループ交換は、レジスタに乗せたまま次のループを実行できてメモリ転送が減る利点があるので、gccはその点を重要視してループ交換最適化を適用したのかもしれません。

実験結果は、test1とtest3が1.218秒くらい、test2が0.996秒くらい、test4とtest5が0.985秒くらい、test6とtest7が1.059秒くらい、でした。 この結果からまず、コンパイラの自動SIMD化コードとイントリンシックを用いたコードの速度差は高々1%程度であることがわかります。 また、TwoSumCond(2blend)＜TwoSumCond(mask_add)＜＜TwoSumという速度差が見受けられます。

Zen 4の発行ポートは、uops.infoによれば

  • vaddpd(zmm)：FP2/FP3
  • vblendmpd(zmm)：FP0/FP1/FP2/FP3
  • vandnpd(zmm)：FP0/FP1/FP2/FP3
  • vcmpltpd(zmm)：FP0/FP1

なので、TwoSumが遅いのはvaddpdばかりやっているからだということが読み取れます。

各実装の最適なポート割り当てを見ていくと、

  • TwoSum：add×6：FP2/FP3に6命令
  • TwoSumCond(mask_add)：and×2, cmp×1, add×4, mask_add×1：FP0/FP1に3命令・FP2/FP3に5命令
  • TwoSumCond(2blend)：and×2, cmp×1, blend×2, add×3：FP0/FP1に4命令・FP2/FP3に4命令

となるので、この順番で速くなるのは納得です。

実際にかかったサイクル数と理論サイクル数を表にまとめてみました。 このプロセッサは5.44 GHzで動いていること、Zen 4には256ビット幅の演算器しかないため512ビット幅の演算は2回演算器を占有することに注意して計算しました。

実装 かかった時間（秒） サイクル数 理論サイクル数 比
TwoSum 1.218 6.62 G 6.44 G 1.03
TwoSumCond(mask_add) 1.059 5.76 G 5.37 G 1.07
TwoSumCond(2blend) 0.985 5.36 G 4.30 G 1.25

これを見ると、TwoSumとTwoSumCond(mask_add)は、理論サイクルとほぼ同じくらいで、理論による予測が正しそうだということがわかります。 TwoSumCond(2blend)は理論サイクルと比べてかなり遅いですが、これはポート振り分けが理論上最適な均等振り分けにならないことによる損が発生していることが原因と考えられます。

追加実験

TwoSumは命令数が少ないのにポートを偏って使用しているため低速でした。 そこで、加算のためにvaddpd命令（FP2/FP3を使う）を使うのではなく、贅沢にvfmaddpd命令（FP0/FP1を使う）を使ってしまえば、高速になる可能性があります。

fmaは破壊的3オペランド命令なので、破壊できるオペランドを持っている加算を置き換える場合は問題ありませんが、破壊できるオペランドがない場合はvmovapd命令が増えてしまうのが気になります（レジスタリネーミングの段階でmov eliminationされるからノーコストというのが定説ですが）。 なので、加算の少なくとも一方が以降でもう使わない変数になっているところ3か所を置き換えてみました。

void test8() {
    auto start = std::chrono::system_clock::now();
    for( int i = 0; i < 65536*32; ++i ) {
        const __m512d one = _mm512_set1_pd(1.0);
        for( int j = 0; j < 4096; j += 8 ) {
            __m512d zmm1 = _mm512_load_pd( &arr1[j] );
            __m512d zmm2 = _mm512_load_pd( &arr2[j] );
            __m512d zmm_hi  = zmm1 + zmm2;
            __m512d zmm_tmp = zmm_hi - zmm1;
            __m512d zmm_lo = _mm512_fmadd_pd(_mm512_fmsub_pd(zmm1, one, zmm_hi - zmm_tmp), one, _mm512_fmsub_pd(zmm2, one, zmm_tmp) );

            _mm512_store_pd( &arr1[j], zmm_hi );
            _mm512_store_pd( &arr2[j], zmm_lo );
        }
    }
    auto finish = std::chrono::system_clock::now();
    std::cout << "TwoSum, Intrinsic (3add3fma) " << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() * 1e-9 << std::endl;
}

実験してみると、0.728秒でした。

さらに、置き換え方を全通り探索しました。 具体的には加算6か所について、a+bを使う、fma(a,1.0,b)を使う、fma(b,1.0,a)を使う、の3通りを試す（合計729通りの）実験しました。 fmaを3回使うものが断然速く、その中でもhi_ = (hi-tmp)とb_ = (b-tmp)とa_ = a-hi_の三か所をfmaにするのが最速のようでした（a_ = a-hi_とb_ = (b-tmp)とlo = a_+b_の三か所をfmaにする上記の構成も次点で高速でした）。 以下のテンプレートの中でtest9<0,0,2,1,2,0>が最速で、0.708秒でした。

fmaを使う回数別に最速の実装の速度を見てみると、0回：1.18秒、1回：0.99秒、2回：0.79秒、3回：0.71秒、4回：0.81秒、5回：0.99秒、６回：1.18秒、という感じでした。

template<int> __m512d zmm_add_pd( __m512d a, __m512d b ) { return a + b; }
template<> __m512d zmm_add_pd<1>( __m512d a, __m512d b ) { return _mm512_fmadd_pd( a, _mm512_set1_pd(1.0), b ); }
template<> __m512d zmm_add_pd<2>( __m512d a, __m512d b ) { return _mm512_fmadd_pd( b, _mm512_set1_pd(1.0), a ); }
template<int> __m512d zmm_sub_pd( __m512d a, __m512d b ) { return a - b; }
template<> __m512d zmm_sub_pd<1>( __m512d a, __m512d b ) { return _mm512_fmsub_pd( a, _mm512_set1_pd(1.0), b ); }
template<> __m512d zmm_sub_pd<2>( __m512d a, __m512d b ) { return _mm512_fnmadd_pd( b, _mm512_set1_pd(1.0), a ); }

template<int A, int B, int C, int D, int E, int F>
void test9() {
    auto start = std::chrono::system_clock::now();
    for( int i = 0; i < 65536*32; ++i ) {
        const __m512d one = _mm512_set1_pd(1.0);
        for( int j = 0; j < 4096; j += 8 ) {
            __m512d zmm1 = _mm512_load_pd( &arr1[j] );
            __m512d zmm2 = _mm512_load_pd( &arr2[j] );
            __m512d zmm_hi  = zmm_add_pd<A>( zmm1, zmm2 );
            __m512d zmm_tmp = zmm_sub_pd<B>( zmm_hi, zmm1 );
            __m512d zmm_lo = zmm_add_pd<F>( zmm_sub_pd<D>( zmm1, zmm_sub_pd<C>( zmm_hi, zmm_tmp ) ), zmm_sub_pd<E>( zmm2, zmm_tmp ) );

            _mm512_store_pd( &arr1[j], zmm_hi );
            _mm512_store_pd( &arr2[j], zmm_lo );
        }
    }
    auto finish = std::chrono::system_clock::now();
    std::cout << "TwoSum, Intrinsic (3add3fma, " << A << B << C << D << E << F << ") " << std::chrono::duration_cast<std::chrono::nanoseconds>(finish - start).count() * 1e-9 << std::endl;
}

理論サイクル数を求めてみます。 上記実装はFP0/FP1に3命令・FP2/FP3に3命令になっていて、512ビット幅の演算は2回演算器を占有することに注意すると、理論的には3サイクルでいけそうです。 これを表にまとめると以下のようになります。

実装 かかった時間（秒） サイクル数 理論サイクル数 比
TwoSum 1.218 6.62 G 6.44 G 1.03
TwoSumCond(mask_add) 1.059 5.76 G 5.37 G 1.07
TwoSumCond(2blend) 0.985 5.36 G 4.30 G 1.25
TwoSum(3add3fma) 0.728 3.96 G 3.22 G 1.23
TwoSum(3add3fma, 002120) 0.708 3.85 G 3.22 G 1.20

また理論サイクル数とかなり乖離しました。 ポートの振り分けは自明なはずで、振り分けが難しいことによる損ではなさそうです。 念のため命令数を数えてみると、二倍アンローリングされているループ一周当たり、56 Byte・20 x86命令・20マイクロ命令？*2になっていたので、6サイクルあればZen 4の6マイクロ命令ディスパッチ/cycleでも十分に命令を供給できそうです。

この辺りをよく考えてみたところ、どうやら「FP0～3で合計で最大で8レジスタしか読み出せない」というZen 2の頃からの制約が関係しているようです。 Zen 2では、FP0がレジスタ読み出しポート0,1,2を、FP1がレジスタ読み出しポート3,4,5を、FP2がレジスタ読み出しポート6,7を、それぞれ使う一方、FP3はレジスタ読み出しポート2,3を使うため、FP0かFP1からfmaが発行されるときはFP3からaddが発行できません（fma+fma+add+addの組み合わせでの発行ができず、mul+mul+add+addやfma+fma+addの組み合わせまでしか発行できません）。 これだけだと4fma2addの実装が理論上最速となって3fma3addの実装の方が速かったこととの整合性が無いように見えますが、Zen 2やZen 3では第三オペランドがバイパス経由で受け取れる場合はFP3からの発行が阻害されないらしいです（Zen 2 - Microarchitectures - AMD - WikiChipおよびSoftware Optimization Guide for AMD EPYC 7003）。 また、Zen 4からは条件が変わってFP3とは競合しなくなったようです。 具体的には、FP0の第三オペランドはFP4用と、FP1の第三オペランドはFP5用と、それぞれ共有されていて、バイパス経由で受け取れるかFP4/FP5用のレジスタ読み出しポートから供給できるなら、発行できるようです（Software Optimization Guide for the AMD Zen4 MicroarchitectureおよびSoftware Optimization Guide for the AMD Zen5 Microarchitecture）。

バイパスで受け取るかに応じてスケジューリングを変えるというのは、スケジューラが複雑になるので普通はやりたくないところだと思いますが、まぁ不可能でもないのでやっている感じでしょうか。 実際には、バイパスで受け取れる場合ではなく、結果が生成されたサイクルだけ、競合しないらしいので、wakeup直後を特別扱いする実装にして追加コストを下げている気がします。 行列積を計算する時のことを考えると、ストア用のFP4/FP5よりもadd用のFP3と共有したほうが合理的のような気もしますが（ストアは1つしかレジスタを読み出さないので、fma+storeは問題にならないのかも）、Zen 4での変更はおそらく、AVX-512をサポートしたことによるものです。 AVX-512では、マスク付きであればソースオペランドを3つ取る*3ことがあり得るので、FP2/FP3で実行されるvaddpdなども第三オペランドが必要になることがあり、それをFP4/FP5との共有で確保しているということだと思います。

結局、FP4やFP5はストアで使うので常に空いているとは限らず、バイパス経由の方も常に使えるとは限らないので、限界までのスループットが出ない、ということだったようです。

まとめ

TwoSumをAVX-512でSIMD化する方法は複数ありますが、それらの速度を比較しました。 Alder Lakeでは、命令数の少ない、TwoSumを素直にSIMD化する方法が最適でした。 一方Zen 4では、TwoSumをSIMD化したものよりも、TwoSumCondをSIMD化したものの方が高速でした。 これは、TwoSumはFP2/FP3発行ポートしか使えない一方、TwoSumCondはFP0/FP1発行ポートも活用できることによるものでした。 したがって、FP0/FP1発行ポートを使うfmaを加算のために贅沢に使うことで、さらに高速化できました。

*1:実際には最終演算をマスク付き演算にすればいいので、選択命令は不要

*2:vmovapd zmm0 zmm2は消えるはず・ロード付き命令が2個・cmp-jneフュージョンが1個……であってる？

*3:レジスタリネームしているので、マスクされた箇所だけ書き込まないみたいな実装は無理です。



以上の内容はhttps://lpha-z.hatenablog.com/entry/2025/09/07/231500より取得しました。
このページはhttp://font.textar.tv/のウェブフォントを使用してます
不具合報告/要望等はこちらへお願いします。
モバイルやる夫Viewer Ver0.14