AVXでリダクションをしようと思ってみた

　Radeon HD7970!!!!!!!!!!!!
　僕の環境だと7970さんのPCIe周りがクソ遅いのでGTX580に完敗していてマジ涙目過ぎる今日この頃なのですが、コレはドライバの所為なのか、それとも僕のカードがはずれを引いたのかどっちだ！！

　――って感じな今日この頃の山田てるみです。
　そうそう、本日(ってももう昨日か) 開催された第3回x86/x64最適化勉強会にいってまいりました。
　いやーとても濃い話がいっぱいあって、すごく面白かったです。いつか僕もあんな感じに発表したいなとか思いましたが、さて、いつになるのかなと思っています

　それに触発されて、というわけではないですが、ブログをしたためているわけです。もともとはGCC 4.7 Cilk+ brancheの話を書くつもりでここ数日色々やっていたのですが、なんかうまく動いてくれないので、あれこれと手を尽くしていますがgccビルド時間かかりすぎﾜﾛﾀ
　なので、代わりに今、山田が地味に必要だなと思っていた部分のコードを書いていました。
　それがタイトルのAVXでのリダクションです

　もっと詳細に書くと、"インデックス付き最小値算出"というやつで、要はn個の配列の中に存在する最小値とそのインデックスを引っ張ってきたい、という話。
　で、配列への演算はAVX / SSEで行っているので、この際まとめてやってしまおうということです。
　もともと、出来るだろうということはわかっていたし、今のところnの値も決め打ちなので、あんまり考慮していない実装ではあるけども、あっさり書きました。あぁ、平和だ。

　環境はCorei7 2600Kで、gccは4.6.1.オプションは-mtune=native -O3でのテストです

　コードは以下の通りです。

#include <iostream>
#include <x86intrin.h>
#include <vector>
#include <cstdlib>
#include <sys/time.h>

typedef struct ReturnParams__ {
    float value;
    size_t idx;
} ReturnParams_t;

double gettimeofday_sec()
{
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return tv.tv_sec + tv.tv_usec * 1e-6;
}

void simd_reduction(ReturnParams_t& ret, const std::vector<float>& load_vec, const std::vector<float>& idx)
{
    const size_t num = load_vec.size();

    __m256 vec[num / 8];
    __m256 v_idx[num / 8];

    const size_t vec_num = num / 8;
    for(size_t i = 0;i < vec_num; ++i) {
        vec[i] = _mm256_loadu_ps(&load_vec[i*8]);
        v_idx[i] = _mm256_loadu_ps(&idx[i*8]);
    }

    // reduction
    __m256 v01;
    __m256 vidx_01;
    {
        __m256 v01_mask = _mm256_cmp_ps(vec[0],vec[1],1);
        v01 = _mm256_blendv_ps(vec[1],vec[0],v01_mask);
        vidx_01 = _mm256_blendv_ps(v_idx[1],v_idx[0],v01_mask);
    }

    __m256 v23;
    __m256 vidx_23;
    {
        __m256 v23_mask = _mm256_cmp_ps(vec[2],vec[3],1);
        v23 = _mm256_blendv_ps(vec[3],vec[2],v23_mask);
        vidx_23 = _mm256_blendv_ps(v_idx[3],v_idx[2],v23_mask);
    }

    __m256 v0123;
    __m256 vidx_0123;
    {
        __m256 v1234_mask = _mm256_cmp_ps(v01,v23,1);
        v0123 = _mm256_blendv_ps(v23,v01,v1234_mask);
        vidx_0123 = _mm256_blendv_ps(vidx_23,vidx_01,v1234_mask);
    }

    float stored_vec[8] __attribute__((aligned(32)));
    _mm256_store_ps(stored_vec,v0123);

    float stored_idx[8] __attribute__((aligned(32)));
    _mm256_store_ps(stored_idx,vidx_0123);

    float min_value = stored_vec[0];
    float min_idx = stored_idx[0];
    for(size_t i = 1;i < 8;++i) {
        if(min_value > stored_vec[i]) {
            min_value = stored_vec[i];
            min_idx = stored_idx[i];
        }
    }

    ret.value = min_value;
    ret.idx = static_cast<size_t>(min_idx);
}

void normal_reduction(ReturnParams_t& ret, const std::vector<float>& load_vec, const std::vector<float>& idx)
{
    const size_t num = load_vec.size();

    float min_value = load_vec[0];
    size_t min_idx = 0;
    for(size_t i = 1;i < num;++i) {
        if(min_value > load_vec[i]) {
            min_value = load_vec[i];
            min_idx = idx[i];
        }
    }

    ret.value = min_value;
    ret.idx = min_idx;
}

int main(int argc, char** argv)
{
    size_t num = 32;

    srand(1);

    std::vector<float> load_vec(num);
    std::vector<float> idx(num);

    for(size_t i = 0;i < num;++i) {
        load_vec[i] = static_cast<float>(rand() % 100);
        idx[i] = static_cast<float>(i);
    }

    const size_t loop_num = 10000;
    double s_0 = gettimeofday_sec();

    ReturnParams_t ret_simd;
    for(size_t i = 0;i < loop_num;++i) {
        simd_reduction(ret_simd, load_vec, idx);
    }

    double s_1 = gettimeofday_sec();

    double n_0 = gettimeofday_sec();

    ReturnParams_t ret_norm;
    for(size_t i = 0;i < loop_num;++i) {
        normal_reduction(ret_norm, load_vec, idx);
    }

    double n_1 = gettimeofday_sec();

    std::cout << "Simd: " << ret_simd.value << " " << ret_simd.idx << " needs : " << s_1 - s_0 << " [sec]" << std::endl;
    std::cout << "Norm: " << ret_norm.value << " " << ret_norm.idx << " needs : " << n_1 - n_0 << " [sec]" << std::endl;

    return 0;
}

　んで、結果が

Simd: 2 31 needs : 0.000283003 [sec]
Norm: 2 31 needs : 0.000405788 [sec]

　……うーん？
　まぁ、多少は速くなっているような…
　どうなの、この結果は……？

　とりあえず、Cilk+のGCCでの動作を早々に確認したいところなんですけど、なんで動かないんですかね？
　どなたか情報ご存じの方は教えてくださいオナシャス　

　ちなみに。
　Bulldozerでこれを実行すると

Simd: 2 31 needs : 0.000820875 [sec]
Norm: 2 31 needs : 0.00183201 [sec]

　何かの間違いだ