SIMD 演算で if 文の並列化 - その2 -

2012-10-20 - kawa0810の日記の続きです.今回は SIMD 演算の実用編(?)として,モンテカルロ法を用いて円周率を求めます.ただし,整数の AVX は Haswell から対応なので,円の内側の点数を数える際に整数ではなく浮動小数点を使用してます.また,乱数の生成にはメルセンヌ・ツイスタを使用してみます(もしかすると,一部の環境で動かないかもしれません...orz).

AVX でモンテカルロ法を実装する際の注意点

モンテカルロ法メルセンヌ・ツイスタについての詳細は wiki 等をご確認ください.モンテカルロ法で円周率を求める場合,円の内側にある離散点数を数える必要があります.今回は離散点が円の内側か外側かの判定を _mm256_cmp_pd() を用いてマスクを生成し,生成したマスクを元に円の内側の離散点の数を数える方法で実装しています.離散点の数を数える型には一般的に整数型を使用すると思いますが,AVX で整数の SIMD 演算は Haswell からの対応となるため,現在発売されているマイクロプロセッサでは使用することができません.よって,今回は代わりに浮動小数点型を用いて実装を行います.

AVX でモンテカルロ法のサンプルコード

倍精度で円周率を求めるサンプルコードです.

/**--------------------
   montecalro.cpp
   
   - Written by kawa0810
   --------------------*/

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <immintrin.h>//-mavx

#include <random>
#include <algorithm>
#include <functional>

std::size_t const mem_align = 32;//メモリアライメント(AVX なので 32byte 境界に揃える)

//逐次で計算した場合
double montecalro(std::size_t const n, double const* x, double const* y){
  std::size_t in = 0;

  for(std::size_t i=0; i<n; ++i){
    if(x[i] * x[i] + y[i] * y[i] <= 1.0) ++in;      
  }

  return static_cast<double>(4.0 * in) / n;
}

//AVX の SIMD 演算でモンテカルロ法を計算する場合
double simd_montecalro(std::size_t const n, double const* x, double const* y){
  std::size_t const para = mem_align / sizeof(double);
  std::size_t end = (n / para) * para;

  __m256d const inside = _mm256_set1_pd(1.0);//円の半径

  __m256d vin = _mm256_setzero_pd();//離散点を数える用
  double* in = reinterpret_cast<double*>( &vin );//離散点を最後に足し合わせる用

  for(std::size_t i=0; i<end; i+=para){
    __m256d vx = _mm256_load_pd(x+i);
    __m256d vy = _mm256_load_pd(y+i);

    //x[i ~ i+3] * x[i ~ i+3] + y[i ~ i+3] * y[i ~ i+3]
    __m256d vnorm = _mm256_add_pd(_mm256_mul_pd(vx, vx), _mm256_mul_pd(vy, vy));

    //norm <= 1.0
    __m256d mask = _mm256_cmp_pd(vnorm, inside, _CMP_LE_OS);//外側か内側か判定
    __m256d add_val = _mm256_and_pd(mask, inside);//外側の離散点に対応する部分を0にしておく
    vin = _mm256_add_pd(vin, add_val);//離散点を数える
  }
  
  //端数処理
  for(std::size_t i=end; i<n; ++i)
    if(x[i] * x[i] + y[i] * y[i] <= 1.0) in[0] += 1.0;

  //数えた結果をまとめる    
  in[0] += in[1] + in[2] + in[3];

  double pi = 4.0 * in[0] / n;

  return pi;
}

int main(void){
  std::size_t const n = 100000;
  double* x;
  double* y;
  
  //メルセンヌ・ツイスタで乱数を生成
  auto rnd = std::bind(std::uniform_real_distribution<double>(0.0, 1.0), 
		       std::mt19937(arc4random()));

  posix_memalign(reinterpret_cast<void**>(&x), mem_align, n * sizeof(double));
  posix_memalign(reinterpret_cast<void**>(&y), mem_align, n * sizeof(double));
  
  for(std::size_t i=0; i<n; ++i) x[i] = rnd();
  for(std::size_t i=0; i<n; ++i) y[i] = rnd();  

  std::cout << "pi = " << 2.0 * acos(0.0) << std::endl;
  std::cout << "Montecalro pi = " << montecalro(n, x, y) << std::endl;
  std::cout << "SIMD Montecalro pi = " << simd_montecalro(n, x, y) << std::endl;

  free(x);
  free(y);

  return 0;
}
SIMD 演算で if 文の並列化