【c/c++】AVX2を使った円周率を求めるプログラム――修正

結果の加算周りが非効率なためこの記事に載せたプログラムは遅いです。
この部分を修正した場合、AVX2無しで動かしたほうが高速です。
wrongwrong163377.hatenablog.com

---

wrongwrong163377.hatenablog.com
上記記事で作ったモンテカルロ法で円周率を求めるプログラムの実装がマズかったので、マルチスレッドSIMD有り（MtWs）のみ作り直しました。

マズかったところ

前回記事を書いた後でkawa0810様による以下の2つの記事を見つけました。
kawa0810.hateblo.jp
kawa0810.hateblo.jp

SSE や AVX では条件分岐文の変わりに cmp 関数を用いて計算を行います.ただし,SSE 世代と AVX 世代では cmp 関数の使用方法が大きく変更されました.SSE 世代で a == b を計算したければ _mm_cmpeq_ps() 関数,a < b を計算したければ _mm_cmplt_ps() 関数など演算子ごとに専用の関数が用意されていましたが,AVX 世代では以下の関数で全ての条件演算子の計算を行います.

はい。この部分完全に勘違いしていて、AVX2ではcmpが無いもんだと思ってました。
この他にも一部気に入らない点があったので書き直しを行いました。

修正後

double MtWs_Kai() {
    using namespace std;

    mt19937 m = mt19937((int)time(NULL));
    constexpr int tnum = threadnum;
    mt19937_64 mts[tnum];
    for (mt19937_64 &mt : mts) { mt = mt19937_64(m()); }

    constexpr long L = (num/tnum) >> 2;
    double count = 0.0;
    const __m256d one = _mm256_set1_pd(1.0);

    #pragma omp parallel for 
    for (int i = 0; i < tnum; i++) {
        int t;
        double d[4];
        __m256d vx, vy, vnorm, vans;
        vans = _mm256_setzero_pd();
        double* ans = reinterpret_cast<double*>(&vans);
        uniform_real_distribution<> uniform = uniform_real_distribution<>(0.0, 1.0);

        for (long j = 0; j < L; j++) {
            //乱数で取って
            for (t = 0; t < 4; t++) d[t] = uniform(mts[i]);
            vx = _mm256_load_pd(d);
            for (t = 0; t < 4; t++) d[t] = uniform(mts[i]);
            vy = _mm256_load_pd(d);
            //ノルム
            vnorm = _mm256_add_pd(_mm256_mul_pd(vx, vx), _mm256_mul_pd(vy, vy));
            //比較してマスクを取り、これと1のandを取ることで外側を0にして、加算する
            vans = _mm256_add_pd(vans, _mm256_and_pd(_mm256_cmp_pd(vnorm, one, _CMP_LE_OS), one));
        }

        #pragma omp atomic
        count += ans[0] + ans[1] + ans[2] + ans[3];
    }
    return count / (double)(num >> 2);
}

改善した点

全体をループにするのではなく、外側と内側のループに分ける

このやり方では以下の3点の利点があります。

• 各スレッドで変数を宣言できるため、コピーが無い分高速。
• スレッド分割がどうなっているか把握しやすい。
• 変数のスコープが限定されるのでコードが見やすい。

点の位置判定と総和をAVX2命令で計算

はじめに書いたとおり、kawa0810様の記事を参考にcmp命令が使えるので判定部を置き直しました。

書き直した後の実行時間

前回記事でいうX=9で最速14.5秒程度でした。1.5秒程度の高速化がサックリ出来てしまいました。
勉強不足ェ……。