CPU(AVX/AVX2)とGPU(OpenCL/CUDA)のベンチマーク -- 意外と速いCPU 2014/1/29

GPGPUはCPUと比べてどのくらい速いのか

近年、GPUの性能は順調に伸び、OpenCLやCUDAなど、GPGPUを利用しやすい環境が整ってきた。また、CPU統合GPUもAMD、Intel双方から提供されている。一方、CPUの方もSandy Bridgeで256bit長のベクトル演算(主に浮動小数点演算)が可能になり、Haswellでさらに強化されて、256bit長の整数演算が可能になった。こうした状況の中、CPUとGPUの速度差がどの程度なのか、実際に測定してみた。

ベンチマーキングの環境

OpenCL 1.2をOSが標準でサポートしていること、CPU(Core i7-4850HQ(Crystal Well))、CPU統合GPU(Intel Iris Pro Graphics 5200)、Discrete GPU(NVIDIA GeForce GT 750M)の3種類の演算ユニット搭載していることから、MacBook Pro Late 2013(Retinaディスプレイ 2300/15.4 ME294J/A)上で、ベンチマークを行った。一部のベンチマークでは、Xeon E5-1650やNVIDIA GeForce GTX 660、を使用した。

ベンチマークプログラム

今回使用したベンチマークプログラムは、AppleがOpenCLのサンプルプログラムとして公開しているOpenCL N-Body Simulationと、一部、NVIDIAのCUDA SDKに含まれるnbodyサンプルプログラムを使用した。どちらのプログラムも、一回分のループを20 floating point operationsと見なして性能(Gflops)を算出している。body数nが10,000の計算を2秒に1回計算できるとすると20*(10,000*10,000)/2 = 1Gflopsとなる。また今回は、単精度浮動小数点数(float)でのベンチマークを行った。

但し、これらのサンプルプログラムは、あくまでOpenCL/CUDAのサンプルプログラムなので、比較対象としてのCPU用の実装は、恣意的なのか単なる手抜きなのか、非常に効率の悪い実装になっている(NVIDIAのnbodyサンプルプログラムには、CPUでの処理速度を表示する「-cpu」オプションがあるが、このオプションで表示されるCPUの処理速度に至っては、全く参考にならない(シングルスレッドの実装だし、CPUのSIMD命令が活かせるように書かれていない))。またデフォルトでは、コード生成のターゲットCPUのアーキテクチャーも古く、最適化ベルも低い。これではCPU本来の性能を出し切れない。GPUとのフェアな比較のために、CPU用の実装として、AVXのIntrinsicsを使用して私が手書きしたものも追加した(本当は、LLVM clangの自動ベクトル化に任せた方がフェアな気もするが、LLVM clangの自動ベクトル化が未成熟なので、CPU本来の性能を引き出すことの方を重視した)。下記が、そのコードのコアな部分である(最内ループのみ最適化している)。

namespace {
   inline float add_horizontally(__m256 packed_values) {
       packed_values = _mm256_hadd_ps(packed_values, packed_values);
       packed_values = _mm256_hadd_ps(packed_values, packed_values);
       __m128 sum2 = _mm256_extractf128_ps(packed_values, 1);
       return _mm_cvtss_f32(_mm_add_ss(_mm256_castps256_ps128(packed_values),
                                       sum2));
   }

   inline void compute_one_chunk(unsigned int i, float *p_x,float *p_y, float *p_z,
                                 __m256 pos_x, __m256 pos_y, __m256 pos_z,
                                 __m256 &acc_x, __m256 &acc_y, __m256 &acc_z,
                                 float *host_mass,
                                 __m256 softeningSq) {
       __m256 dx = _mm256_load_ps(&p_x[i]) - pos_x;
       __m256 dy = _mm256_load_ps(&p_y[i]) - pos_y;
       __m256 dz = _mm256_load_ps(&p_z[i]) - pos_z;

       __m256 Mi = _mm256_load_ps(&host_mass[i]);

       __m256 distSqr = dx * dx + dy * dy + dz * dz;
       distSqr += softeningSq;

       __m256 invDist = _mm256_rsqrt_ps(distSqr);
       __m256 s = (Mi * invDist) * (invDist * invDist);

       acc_x += dx * s;
       acc_y += dy * s;
       acc_z += dz * s;
   }
}

int CPUSimulation::executeMultiSimd(unsigned int startDestIndex, unsigned int destCount)
{
   float const deltaTime = m_active_params.m_timestep;
   float const damping   = m_active_params.m_damping;
   float const softening = m_active_params.m_softening;

   float const softeningSq = softening * softening;

   dispatch_queue_t queue = dispatch_get_global_queue(DISPATCH_QUEUE_PRIORITY_DEFAULT, 0);
   dispatch_group_t group = dispatch_group_create();

   for (unsigned int l = startDestIndex; l < startDestIndex + destCount; l++)
   {
       dispatch_group_async(group, queue, ^{
           unsigned int i;

           float position_x = m_host_position_x[m_read_index][l];
           float position_y = m_host_position_y[m_read_index][l];
           float position_z = m_host_position_z[m_read_index][l];

           float accX = 0.0f;
           float accY = 0.0f;
           float accZ = 0.0f;
           {
               __m256 pos_x = _mm256_set1_ps(position_x);
               __m256 pos_y = _mm256_set1_ps(position_y);
               __m256 pos_z = _mm256_set1_ps(position_z);
               __m256 acc_x = _mm256_set1_ps(0.0f);
               __m256 acc_y = acc_x;
               __m256 acc_z = acc_x;
               __m256 softeningSq_ = _mm256_set1_ps(softeningSq);

               size_t end = (m_body_count / 32) * 32;
               for (i = 0; i < end; )
               {
                   compute_one_chunk(i, m_host_position_x[m_read_index],
                                     m_host_position_y[m_read_index],
                                     m_host_position_z[m_read_index],
                                      pos_x,  pos_y,  pos_z,
                                     acc_x, acc_y, acc_z,
                                     m_host_mass,
                                     softeningSq_);
                   i+=8;
                   compute_one_chunk(i, m_host_position_x[m_read_index],
                                     m_host_position_y[m_read_index],
                                     m_host_position_z[m_read_index],
                                     pos_x,  pos_y,  pos_z,
                                     acc_x, acc_y, acc_z,
                                     m_host_mass,
                                     softeningSq_);
                   i+=8;
                   compute_one_chunk(i, m_host_position_x[m_read_index],
                                     m_host_position_y[m_read_index],
                                     m_host_position_z[m_read_index],
                                     pos_x,  pos_y,  pos_z,
                                     acc_x, acc_y, acc_z,
                                     m_host_mass,
                                     softeningSq_);
                   i+=8;
                   compute_one_chunk(i, m_host_position_x[m_read_index],
                                     m_host_position_y[m_read_index],
                                     m_host_position_z[m_read_index],
                                     pos_x,  pos_y,  pos_z,
                                     acc_x, acc_y, acc_z,
                                     m_host_mass,
                                     softeningSq_);
                   i+=8;
               }
               accX = add_horizontally(acc_x);
               accY = add_horizontally(acc_y);
               accZ = add_horizontally(acc_z);
           }
           for (; i < m_body_count; i++)
           {

               float dx = m_host_position_x[m_read_index][i] - position_x;
               float dy = m_host_position_y[m_read_index][i] - position_y;
               float dz = m_host_position_z[m_read_index][i] - position_z;

               float Mi = m_host_mass[i];

               float distSqr = dx * dx + dy * dy + dz * dz;
               distSqr += softeningSq;

               float invDist = 1.0f / sqrtf(distSqr);
               float s = (Mi * invDist) * (invDist * invDist);

               accX += dx * s;
               accY += dy * s;
               accZ += dz * s;

           }

           float velocity_x = m_host_velocity_x[m_read_index][l];
           float velocity_y = m_host_velocity_y[m_read_index][l];
           float velocity_z = m_host_velocity_z[m_read_index][l];

           velocity_x += accX * deltaTime;
           velocity_y += accY * deltaTime;
           velocity_z += accZ * deltaTime;
           velocity_x *= damping;
           velocity_y *= damping;
           velocity_z *= damping;

           position_x += velocity_x * deltaTime;
           position_y += velocity_y * deltaTime;
           position_z += velocity_z * deltaTime;

           m_host_position_x[m_write_index][l] = position_x;
           m_host_position_y[m_write_index][l] = position_y;
           m_host_position_z[m_write_index][l] = position_z;

           m_host_velocity_x[m_write_index][l] = velocity_x;
           m_host_velocity_y[m_write_index][l] = velocity_y;
           m_host_velocity_z[m_write_index][l] = velocity_z;

           m_host_position[l].data[0] = position_x;
           m_host_position[l].data[1] = position_y;
           m_host_position[l].data[2] = position_z;
           m_host_position[l].data[3] = 1.0f;
       });
   }
   dispatch_group_wait(group, DISPATCH_TIME_FOREVER);
   dispatch_release(group);
   return 0;
}

ベンチマーク結果

ベンチマーク結果は、下記のグラフの通りだ。CPUやIris Proの性能は、CPUが熱くなる前のTurbo Boostが効いている時の数値である。温度が上がってくると1割前後性能が低下する。Discrete GPUのGeForceの性能は比較的安定していた。

まず、CPUの処理速度を見てみよう。OpenCLを使った実装は、AVXのIntrinsicsを使った場合に比べて、半分強くらいの性能しか出ていない。OpenCLのclGetDeviceInfo APIで調べると、CPUはi7-4850HQと認識されているが、CL_DEVICE_PREFERRED_VECTOR_WIDTH_FLOAT, CL_DEVICE_NATIVE_VECTOR_WIDTH_FLOATが共に4となっていた。32bit(float) × 4 = 128bitsなので、OpenCLのCPU向けのカーネルコンパイラは、AVX(256bits)に対応しておらずSSEのコードを生成していると思われる。CPU本来の性能を知りたい場合は、グラフの黒色のバーに着目するとよいだろう。シングルスレッドとマルチスレッドの性能を比べてみると、AVX Intrinsicsの場合で、4倍強の差がある。Core i7-4850HQが4コア8スレッドのCPUであることを考慮すると、妥当な性能差である。参考のために、AVX Intrinsicsのプログラムを使って、別のマシンのXeon E5-1650の性能を確認したところ、210 Gflopsの性能が出た。動作周波数の差によるものか、(ターボ・ブースト利用時の最大周波数は大差ないので)コア数の差によるものか、Xeon E5-1650はCore i7-4850HQの1.5倍の性能となった。

各CPUの諸元
	Core i7-4850HQ	Xeon E5-1650
コア数	4	6
スレッド数	8	12
動作周波数	2.3 GHz	3.2 GHz
ターボ・ブースト利用時の最大周波数	3.5 GHz	3.8 GHz

次に、GeForce GT 750Mの結果を見てみよう。やはり、NVIDIAが自社のGPU向けに提供しているCUDAの方が、移植性のあるOpenCLよりも、35%程度高い性能を示した。ただし、今回のMacBook Proのように、NVIDIA以外のGPUやCPUも同時並行で利用できるヘテロジニアス(異種混合)な実行環境では、多少GeForceの性能が下がるとしても、OpenCLで記述してCPUとIris Proの計算能力も活用した方が、トータルの処理能力は高くなる。

Iris Pro Graphics 5200は統合GPUなのであまり期待していなかったのだが、OpenCL実装では、MacBook Proに搭載のデバイスの中で最高性能だった。CPUのIntrinsics実装と比べると約2倍の性能を発揮している。ただし、CPUに統合されているので、前述の通り温度変化に伴う性能の変動が大きい。ノートパソコンでは難しいが、2U以上のサーバーやデスクトップパソコンなら、冷却系を強化すれば高い性能を持続できるだろう。

比較のために、デスクトップ用のNVIDIAのGeForce GTX 660上でのCUDA実装の速度も測定した。1 Tflopsに迫る性能をたたき出し、ノート用のGPUの約3倍、Xeon E5-1650と比べても約4.4倍、Core i7-4850HQと比べると6.5倍となった。

善戦するCPU

GPGPUを使えば、CPUの何十倍もの性能を手に入れられるという印象があった。しかし、今回のベンチマークで、AVX(AVX2)SIMD命令とマルチコアをしっかりと活用するプログラムを使えば、CPUとGPUの性能差は一桁倍に収まることがわかった。もちろん、ハイエンドGPUを複数利用すれば、その差は広がるだろうし、コア数の多いCPUやマルチプロセッサ構成にすれば、差は縮まるだろう。いずれにせよ、OpenCLやCUDAのサンプルプログラムが示すCPUの性能は、あまり当てにならない。

ヘテロジニアス(異種混合)をサポートするOpenCL

今回のベンチマークで使用したMacのように、CPUとGPU 2種(IntelとNVIDIA)を搭載するようなヘテロジニアスな環境において、OpenCLは開発者にとって福音なのは間違いないだろう。OpenCLの使い方を覚えるだけで、様々な種類の演算装置(CPU(Intel, AMD)とGPU(NVIDIA, AMD, Intel))をフル活用できるのは魅力的だ。

一方で、OpenCLにも課題がある。CPUでは、OpenCL実装はIntrinsics実装に比べて半分程度の73 Gflopsだった。また、GeForce GT 750Mにおいても、CUDA実装より25%程度性能が低下した。この辺りは、まだ改善の余地がありそうだ。

確かに、OpenCLを使えば、様々な演算装置で動作するコードをポータブルに記述できる。しかし、単に「動作する」ことと、その演算装置の性能を引き出すことは別物である。OpenCLを使った場合でも、性能を引き出すには、メモリー上のデータレイアウト、キャッシュサイズ、SIMD命令の有無等の各演算装置の様々な違いに合わせてチューニングをする必要がある。絶対的な性能よりも、ポータブルなアプリケーションが作成できるということの方が重要な場合に、OpenCLは適していると言える。

GPUに向き不向きの処理がある

今回のN-bodyプログラムは、たくさんの独立した処理を繰り返すタイプの処理で、GPUが得意とする処理である。実アプリケーションには、独立性が低く並列処理に向かない処理や、テーブルルックアップのように並列処理はできるがSIMDに不向きな処理、計算量に対しメモリーアクセスやI/Oが多い処理なども含まれる。今回の結果は、あくまでそうした様々なタイプの処理の中の一つの結果に過ぎないと思ってほしい。

今回のまとめ

GPU向きの処理タイプにおいても、CPUに最適化されたコードであれば、CPUの性能も思ったほど悪くない
OpenCLはポータブルであるが、性能的には改善の余地あり

Haswell(AVX2)対応(2014/5/26追記)

実は、上記の検証をおこったときは、Xcode等の開発ツールがAVX2に対応していなかった。なので上記のMacのCPUの速度は、OpenCL実装にしてもIntrinsics実装にしても、Haswellコアで実行した結果ではあるが、AVX2命令を活用していない結果だった。現在は、開発ツールもAVX2対応しているので、AVX2命令を使用することができる(但し、測定結果から推測すると、OpenCLはAVXにもAVX2にも未対応のようだ)。そこで、AVX2命令を使って、追加のベンチマークをしてみた。

Haswell(AVX2)の新機能

AVXは、256bit幅のSIMDとしては最初のバージョンであり、洗練されていない部分もあった。AVX2では、そうした部分が改良され、扱いやすく高い性能を出しやすくなっている。詳しくは、IDF 2012 - 次期Intel Core「Haswell」の内部構造を探る - 拡張命令(AVX2/TSX)編|マイナビニュースやHPC性能が大幅に強化されたHaswell|マイナビニュースを参照してほしい。今回のN-Bodyシミュレーションで出番があるのは、FMA命令だ。これは、X*Y+Zのような掛け算と足し算を1命令で(しかも演算精度も向上する)実行する命令だ。今回は出番がないが、Gather命令も高速化に役立ちそうだ。これまでは、変換テーブルのルックアップを行う部分はスカラー演算にせざるを得なかったが、Gather命令を使えば終止ベクトル演算で処理することが可能になる。Xeon PhiではScatter命令もサポートされているようだが、残念ながらHaswellにはScatter命令はないようだ。

AVX2(FMA)を利用した実装

AVX2で導入されたFMA(Fused Multiply Add)を利用して実装したコードは、下記の通りである。

#define ENABLE_FMA 1
    
    inline void compute_one_chunk(unsigned int i, float *p_x,float *p_y, float *p_z,
                                  __m256 pos_x, __m256 pos_y, __m256 pos_z,
                                  __m256 &acc_x, __m256 &acc_y, __m256 &acc_z,
                                  float *host_mass,
                                  __m256 softeningSq) {
        __m256 dx = _mm256_load_ps(&p_x[i]) - pos_x;
        __m256 dy = _mm256_load_ps(&p_y[i]) - pos_y;
        __m256 dz = _mm256_load_ps(&p_z[i]) - pos_z;
        
        __m256 Mi = _mm256_load_ps(&host_mass[i]);

#if ENABLE_FMA
        __m256 distSqr = _mm256_fmadd_ps(dz, dz, _mm256_fmadd_ps(dy, dy, dx * dx));
#else
        __m256 distSqr = dx * dx + dy * dy + dz * dz;
#endif
        distSqr += softeningSq;
        
        __m256 invDist = _mm256_rsqrt_ps(distSqr);
        __m256 s = (Mi * invDist) * (invDist * invDist);
        
#if ENABLE_FMA
        acc_x = _mm256_fmadd_ps(dx, s, acc_x);
        acc_y = _mm256_fmadd_ps(dy, s, acc_y);
        acc_z = _mm256_fmadd_ps(dz, s, acc_z);
#else
        acc_x += dx * s;
        acc_y += dy * s;
        acc_z += dz * s;
#endif
    }

下記が、測定結果である。「AVX2 Intrinsics」と「AVX Intrinsics」を比較するとFMAの効果が分かる。命令を変えただけで10〜20%高速になっているので、この命令を使わない手はない。

今回は使用しなかったが、三角関数等をあらかじめ計算してテーブルにに格納しておいて、Gather命令を使ってルックアップするような手法もある。Haswellで導入されたAVX2命令をうまく活用すれば、AVX命令のみを使った場合より、数割高速化することも可能ではないかと思われる。

他のコンテンツ:

C++でコードをパラメーター化する方法の比較