Intel AVX を使用して SIMD 演算を試してみる - その5 - (最終回)

2012-03-06 - kawa0810の日記の続きです．
今回は Fused Multiply Add (FMA) 命令についてです．ただし，現在 (2012/3/7日) 発売されている Sandy Bridge では FMA 命令は未対応のようですので注意してください．
Why no FMA in AVX in Sandy Bridge?

FMA 命令とは

　プログラミングを行う際に積和演算 ( $D = A \times B + C$ ) といった計算をすることが多い人もいると思います．仮に， $\bf{z} = \alpha \bf{x} + \bf{y}$ といった計算を SIMD 演算で計算することを考えると以下のようなソースコードになると思います．

/*--------------------
  axpy.cpp
  --------------------*/

#include <iostream>
#include <immintrin.h>

void axpy(const std::size_t n, float* z,
	  const float alpha, float* x, float* y){
  static const std::size_t para = 8;
  const __m256 valpha = _mm256_set_ps(alpha, alpha, alpha, alpha, alpha, alpha, alpha, alpha);
  const __m256 *vx = reinterpret_cast<__m256*>( x );
  const __m256 *vy = reinterpret_cast<__m256*>( y );
  __m256 *vz = reinterpret_cast<__m256*>( z );  
  
  const std::size_t end = n / para;
  const std::size_t beg = end * para;
  
  for(std::size_t i=0; i<end; ++i)
    vz[i] = _mm256_add_ps( _mm256_mul_ps(valpha, vx[i]), vy[i]);
  for(std::size_t i=beg; i<n; ++i)
    z[i] = (alpha * x[i]) + y[i];
}

void axpy(const std::size_t n, double* z,
	  const double alpha, double* x, double* y){
  static const std::size_t para = 4;
  const __m256d valpha = _mm256_set_pd(alpha, alpha, alpha, alpha);
  const __m256d *vx = reinterpret_cast<__m256d*>( x );
  const __m256d *vy = reinterpret_cast<__m256d*>( y );
  __m256d *vz = reinterpret_cast<__m256d*>( z );  
  
  const std::size_t end = n / para;
  const std::size_t beg = end * para;
  
  for(std::size_t i=0; i<end; ++i)
    vz[i] = _mm256_add_pd( _mm256_mul_pd(valpha, vx[i]), vy[i]);
  for(std::size_t i=beg; i<n; ++i)
    z[i] = (alpha * x[i]) + y[i];
}

int main(void){
  const std::size_t n = 12;
  const float alpha = 2.0;
  float* x = static_cast<float*>( _mm_malloc(sizeof(float) * n, 32) );
  float* y = static_cast<float*>( _mm_malloc(sizeof(float) * n, 32) );
  float* z = static_cast<float*>( _mm_malloc(sizeof(float) * n, 32) );
  
  for(std::size_t i=0; i<n; ++i) x[i] = i;
  for(std::size_t i=0; i<n; ++i) y[i] = i+1;
  for(std::size_t i=0; i<n; ++i) z[i] = 0.0;

  axpy(n, z, alpha, x, y);

  for(std::size_t i=0; i<n; ++i)
    std::cout << z[i] << std::endl;

  _mm_free(x);
  _mm_free(y);
  _mm_free(z);
  
  return 0;
}

AVX の FMA 命令

　Intel AVX の FMA 命令には以下のような命令があるそうです．
Intel® Software Development Tools

上記の FMA 命令を用いることで現在は2回の SIMD 演算を行う関数を呼び出しするコードが以下のように書くことができるようになることが期待できます．

...
  //単精度
  for(std::size_t i=0; i<end; ++i)
    vz[i] = _mm256_fmadd_ps(valpha, vx[i], vy[i])
...
  //倍精度
  for(std::size_t i=0; i<end; ++i)
    vz[i] = _mm256_fmadd_pd(valpha, vx[i], vy[i])
...

上記のように FMA 命令を用いることで SIMD 演算を行う関数の呼び出し回数が減るだけでなく，_mm256_mul_ps(valpha, vx[i]) の計算結果を格納する際に発生する丸め誤差等も減ることが期待できます．
いずれにせよ Haswell に期待ですね!!

「Intel AVX を使用して SIMD 演算を試してみる」一旦終了

　「Intel AVX を使用して SIMD 演算を試してみる」を一旦終了します．理由は，
これ以上は結構マニアックになってくるため"試してみる"じゃなくなるよね？ && 引っ越しの準備等で忙しくなりそう...
だからですorz

ただ，総和計算とか書きたいと思っていることはあるので時間を見つけて書ければいいなと思っています．
サイドバーの「日記検索」や「カテゴリー」で探してみると良いかもです．