CUDA BLAS
CUDAのツールキットには、CUDA BLASという最適化された演算ライブラリが付属している。BLASはBasic Linear Algebra Subprogramsの略称で、ベクトルと行列を扱う線形計算ルーチンを集めたものである。ベクトルと行列の演算は多くの科学技術計算で必要とされ、その性能が計算時間に大きく影響するので、各社とも自社のハードウェアにチューニングしたBLASを開発して性能をアピールしており、このCUDA BLASはCUDA用に最高性能を出すようにチューニングされたBLASである。
SDKのsimpleBLASプロジェクトを使って、matrixMulに換えて、このCUDA BLASに含まれているsgemm(単精度の行列積)を使ってみた。sgemmは、行列A、Bに対して、C=αAB+β*C(α、βはスカラ値)を計算する機能を持ち、かつ、行列A、Bそれぞれに対して、そのまま使うか、転置して使うかの指定ができるので、matrixMulよりも機能が多いが、その差には目をつぶって行列積の計算性能を測定してみた。
後述のように、GPUで動作する部分のコードは全てCUDA BLASに含まれており、ホスト側で動作するCプログラムからCUDA BLASを呼び出すだけでGPUを使った行列積の計算ができる。それは便利であるが、CプログラムからCUDAの高精度のタイマーを呼び出す方法が分からなかった(simpleBLASには実行時間の測定機能が無く、筆者が追加した)ので、通常のclock関数を使用して実行時間を計測した。そのため、時間の分解能が粗くなってしまい、2048×2048×2048のケースしか実用的な精度の測定ができなかったのであるが、その結果は、以下のようになった。
| ハードウェア | 行列積計算プログラム | 実行時間 [ms] | 実効性能 [GFlops] | ピーク比 |
|---|---|---|---|---|
| 1.18GHz 8600 GT GPU | matrixMul | 924.0 | 18.6 | 0.164 |
| CUDA sgemm | 420 | 40.9 | 0.361 | |
| 2.66GHz E6750 CPU | 単純C(最適化なし) | 296000 | 0.058 | -- |
| 単純C(最適化 /Ox) | 171000 | 0.101 | -- |
CUDAのsgemmを使うと、ホストからGPUへの行列A、Bの転送、行列積計算、結果の行列Cの転送を含めて420ms程度で実行が完了した。つまり、最適化したCUDA BLASのsgemmの使用により、matrixMulに比べて2.2倍に性能が向上している。
ついでに、例題の中で結果のチェックに使っているホストCPUで行列積を計算する関数(ソースコードを次に示す)の実行時間を測定してみると、CUDA BLASの700倍も時間が掛かってしまった。
static void simple_sgemm(int n, float alpha, const float *A, const float *B,
float beta, float *C)
{
int i;
int j;
int k;
for (i = 0; i < n; ++i) {
for (j = 0; j < n; ++j) {
float prod = 0;
for (k = 0; k < n; ++k) {
prod += A[k * n + i] * B[j * n + k];
}
C[j * n + i] = alpha * prod + beta * C[j * n + i];
}
}
}
また、コンパイルオプションに/Oxを指定して最適化した場合は、約1.7倍性能が改善したが、まだ、CUDA BLASの400倍程度の時間が掛かっている。これはGPUとCUDA BLASの高性能を示す証拠とも言えるが、このコードは、あまり効率の良いコードとは言えない。最内ループで行列Aをn跳びにアクセスしており、メモリアクセスに時間が掛かっていることが性能が上がらない原因と思われ、matrixMulのように16×16、あるいは2次キャッシュに格納できる上限の部分行列単位で処理すれば、演算に対するメモリアクセス回数が減り、大幅に性能が上がる筈である。しかし、上記のような単純なCコードを書いてCPUで走らせると悲惨な性能になるという一例ではある。
前記のmatrixMulのコードはGPUで走行させる拡張子が.cuのプログラムであるが、CUDA BLASを使う場合は、ホストで走らせるCプログラムを書けばよい。次に示すように、デバイスメモリの割り付けは、cublasAllocという関数を呼び出し、要素数、要素のサイズを与えると、&d_Aに割り当てられた領域のポインタが返される。
/* Allocate device memory for the matrices */
status = cublasAlloc(n2, sizeof(d_A[0]), (void**)&d_A);
if (status != CUBLAS_STATUS_SUCCESS) {
fprintf (stderr, "!!!! device memory allocation error (A)\n");
return EXIT_FAILURE;
}
そして、次に示すように、cublasSetVector関数の呼び出しでh_Aでポイントされるホストメモリの領域からd_Aでポイントされるデバイスメモリにデータを転送することができる。
/* Initialize the device matrices with the host matrices */
status = cublasSetVector(n2, sizeof(h_A[0]), h_A, 1, d_A, 1);
if (status != CUBLAS_STATUS_SUCCESS) {
fprintf (stderr, "!!!! device access error (write A)\n");
return EXIT_FAILURE;
}
そして、行列積を計算するcublasSgemmを呼び出す。
cublasSgemm('n', 'n', N, N, N, alpha, d_A, N, d_B, N, beta, d_C, N);
さらに、cublasGetVectorの呼び出しで、結果を、デバイスメモリのd_Cでポイントされる領域からホストメモリのh_CGでポイントされる領域に転送すればよい。
/* Read the result back */
status = cublasGetVector(n2, sizeof(h_C[0]), d_C, 1, h_CG, 1);
if (status != CUBLAS_STATUS_SUCCESS) {
fprintf (stderr, "!!!! device access error (read C)\n");
return EXIT_FAILURE;
}
なお、sgemmで2048元の行列の積を計算した場合、実効性能は、ピーク性能の36.1%の性能となっているが、行列積の計算では乗算と加算の出現比率が1対1であり、GeForce 8000シリーズGPUの持っているMAD(乗算＋加算）とMUL(乗算）の演算器を全部、有効に利用することはできない。従って、MADを使いMULを遊ばせるか、MADはADDだけを使うということになるので、毎サイクル2演算しか実行できず、行列積の場合に利用可能なピーク演算性能は、全演算器をフルに利用できる場合の2/3になってしまう。その意味では、sgemmの性能は、現実的に利用可能な最大性能に対しては54%程度の利用率を実現できているといえる。また、前に述べたように8600 GTは最上位のGPUの1/6程度の性能であるので、8800 Ultraや9800 GTXは単精度の行列積計算で240～270GFlops程度の性能を実現できると考えられる。
一方、ホストCPUで行列積を計算した場合の性能は、GPUに比べるとボロ負けであったが、CPUでもチューニングを頑張れば、もっと性能は出る。富士通のFORTRANコンパイラでSSEを使ったBLASライブラリであるblasp4.libのsgemmをリンクしてテストしてみると、次の表のように、2.66GHzのE6750 CPUの1コアで、2048元の行列積を1.406秒で計算することができた。
| ハードウェア | 行列積計算プログラム | 実行時間 [ms] | 実効性能 [GFlops] |
|---|---|---|---|
| 1.18GHz 8600 GT GPU | CUDA sgemm | 420 | 40.9 |
| 2.66GHz E6750 CPU | 富士通FORTRAN blasp4.lib sgemm | 1406 | 12.2 |
これは12.2GFlopsの性能であり、E6750 CPUのピーク性能の57%にあたる。また、GotoBLASで有名なテキサス大の後藤氏は、3GHzクロックのClovertownの1コアで2048元の倍精度行列積(dgemm)でピークの83%にあたる10GFlopsを実現したと報告している。最新の3.2GHzクロックでクワッドコアのQX9770を使い、行列積を4コア並列で処理するプログラムを書けば、blasp4.libのsgemmのレベルで58.7GFlops、GotoBLASと同様の計算手法で単精度のsgemmでも同等のピーク性能比率を実現できるとすると85GFlops程度は出せてよい筈である。
しかし、それでも、240～270GFlopsの性能を持つ現状のトップエンドのGPUは、トップエンドのQX9770 CPUと比較しても3～4倍程度高い演算性能を持っており、性能がいくらでも欲しい科学技術計算にとってはその高性能は魅力的である。
現在のGPUは32ビットの単精度の浮動小数点演算しか処理できないので、一般的な科学技術計算用としては精度不足であるが、NVIDIAもGPUの汎用科学技術計算エンジンとしての用途拡大を目指しており、次世代GPUでは倍精度浮動小数点演算を行えるようにすると表明している。従って、今後、GPUの汎用計算エンジンとしての使用が拡大して行くと考えられる。また、CUDAはスレッドブロック、グリッドという階層的グループ化の概念を用いて、超並列のスレッド処理を記述できるプログラミング環境を提供しており、このプログラミングモデルはGPU以外の超並列コンピュータにも適用されていくのではないかと思われる。(了)
