[Boost] Boost.SIMDを使ってみる(MKLとの比較ベンチもあるよ)

#概要
[Boost.SIMD](http://isocpp.org/blog/2014/02/nt2-boos)を使ってみた。
Takahashi様の[スライド](http://www.slideshare.net/faithandbrave/boostsimd)が参考になる。

まだ正式にはBoost入りしていない。
SPARCをサポートして下さい。なんでもしますから。

実際の所、これは神ライブラリである。

#チュートリアル
とりあえずベクトルの内積をBoost.SIMDで実装してみる。

```cpp
#include <boost/simd/sdk/simd/pack.hpp>
#include <boost/simd/include/functions/sum.hpp>
#include <boost/simd/include/functions/load.hpp>
#include <boost/simd/include/functions/plus.hpp>
#include <boost/simd/include/functions/multiplies.hpp>

template<typename Iterator, typename Iterator2, typename Value = typename remove_const_reference<(*(std::declval<Iterator>()))>::type>
Value dot(Iterator first1, Iterator last1, Iterator2 first2) //内積
{
    typedef boost::simd::pack<Value> pack_t; // SIMDレジスタのラッパークラス
    pack_t tmp(0); 
    while(first1 != last1)
    {
        pack_t x1 = boost::simd::load<pack_t>(&*first1); // memory -> SIMD register load
        pack_t x2 = boost::simd::load<pack_t>(&*first2);

        tmp = tmp + x1 * x2; //SIMDレジスタ間の演算

        first1 += pack_t::static_size; //SIMDレジスタ長分だけイテレータを回す
        first2 += pack_t::static_size;
    }
    return sum(tmp); // SIMDレジスタ内の和
}
```

意味は分かりやすいと思う。

既存のイテレータをアダプトして、SIMDレジスタを返してくれるイテレータを作る事もできる。

```cpp
#include <boost/simd/sdk/simd/pack.hpp>
#include <boost/simd/include/functions/sum.hpp>
#include <boost/simd/include/functions/load.hpp>
#include <boost/simd/include/functions/plus.hpp>
#include <boost/simd/include/functions/multiplies.hpp>
#include <boost/simd/memory/iterator.hpp>

template<typename Iterator, typename Iterator2, typename Value = typename remove_const_reference<(*(std::declval<Iterator>()))>::type>
Value dot(Iterator first1, Iterator last1, Iterator2 first2) //内積
{
    auto it1 = boost::simd::input_begin(first1); // イテレータアダプタ
    auto e1  = boost::simd::input_end(last1);
    auto it2 = boost::simd::input_begin(first2);

    for(;it1!=e1;++it1,++it2)
    {
        tmp = tmp + (*it1) * (*it2); //SIMDイテレータをデリファレンスすると自動的にSIMDレジスタ上にロードする。
    }
    return sum(tmp); // SIMDレジスタ内の和
}
```

ちなみにobjdumpしてアセンブリをみると上記のコードは複雑なコードパスを持ってしまう。
これは、アライメントの処理についてのコードがくっついてきてしまうからだ。

最初からアライメントを揃えておけば良いのだが、コードにアライメントが揃っている事を教えるには次の様にすれば良い。
デフォルトではアライメントチェックにBOOST_ASSERTを使っているので、リリース時にはNDEBUGをdefineしておく。

```cpp
    auto it1 = boost::simd::aligned_input_begin(first1); // イテレータアダプタ
    auto e1  = boost::simd::aligned_input_end(last1);
    auto it2 = boost::simd::aligned_input_begin(first2);
```

STLコンテナを使ってアライメントを揃えるために、Boost.SIMDはカスタムアロケータを提供している。

```cpp

#include <boost/simd/memory/allocator.hpp>
std::vector<double, boost::simd::allocator<double>> v1; //アライメントが揃えられたコンテナ
```

ちなみに、自作アロケータのアライメントを揃えるための、アロケータアダプタもある。詳しくはドキュメントを参照。

#ベンチマーク

`std::vector<double, boost::simd::allocator<double>>`に対して上記の`dot`を走らせてベンチマークを取った。
比較としてIntel MKL 14.0.3の`ddot`関数を用いて同じサイズのベクタの内積を計算させてベンチマークを取った。

OS: CentOS6.4
CPU: Intel Xeon E5-2667 v2 

Boost.SIMDを用いたコードはclang 3.6.0でコンパイルし
Intel MKLを用いたコードはicpc 14.0.3でコンパイルした。

コンパイルオプションは -O3 -m64 -mavx -mfma -std=c++11

それぞれのサイズのベクタに対して内積を16回計算し、その平均処理時間をプロットしたのが下図である。

![bench.png](https://qiita-image-store.s3.amazonaws.com/0/46436/b55e5106-5af3-56aa-43b9-8dec23c6b0ff.png)

この結果を見る限り、MKLのddotは1M次元以下でも非常に大きなサイズのブロックをキャッシュに転送してから計算している様だ。
大きいサイズではMKLの方が速い。
計算そのものに大して違いがあるとは思えないので、この違いはキャッシュ転送周りのコードからくると思われる。
ちなみにMKLのBLAS1,2は並列化されない。