PGIアクセラレータを使用した行列積の計算とGPU性能

キーワード GPGPU、CUDA、MATMUL、PGIアクセラレータ、PGI CUDA Fortran

 ここからのコラムのシリーズの方針を少々変更して、PGIアクセラレータの基本的な使用法に関しての説明は後回しにすることにしました。まずは、いろいろなプログラムへの実施例を見ていただき、実践的な知見を先に得ていただきたいと思います。そこで、GPGPUの計算処理の例として、一般的な行列積(matmul)のプログラムを取り上げます。どこの教科書でもこの MATMUL は、GPGPU の計算適用例として採用されていますので、この PGI アクセラレータコンパイラを使って得た性能を、そのリファレンスとして呈示しておきます。これ以降説明する GPGPU 並列用のプログラミング方法は、PGI アクセラレータ用ディレクティブを挿入して GPU 並列コードに翻訳する ①PGIアクセラレータ・プログラミングモデルを使用する方法と、②PGI CUDA Fortran でコーディングする方法、そして、③NVIDIA の CUBLAS ライブラリを使用する場合の例を示します。まず、最初に、①PGIアクセラレータ・プログラミングモデルを使用して実行した結果をご紹介します。
2010年6月26日 Copyright © 株式会社ソフテック 加藤

行列積(matmul) プログラムへの適用

 この行列積の計算は、プログラム構造的には単純な並列ループで構成され、さらにリダクション処理等もないため、非常に多くのスレッドを同時実行できれば高い性能を得ることができるものです。また、配列の計算部分を小ブロックに分けて(Tiling と言う)、shared memory をキャッシュとして使用する最適化を行うことと、CUDA の Half Warp size(16)を意識した SIMD 長をプログラム上で使うことにより、メモリアクセスの最適化が図れると言った理想的なプログラムです。しかし、一般的なアプリケーションのアルゴリズムは、このような簡単な構造ではないため、究極の性能を得るまでは多くのソースプログラム上の最適化作業が必要です。ここでの内容は、NVIDIA GPU の構造に沿った形で最適化すると、実効性能としてどの程度の能力を出し切ることができるかと言う知見を得るためのものとして捉えて下さい。これから3回に渡って、行列積のプログラム化とその性能について説明します。

 行列積のプログラムに関するコラム(全3回)

  1. PGIアクセラレータを使用した行列積の計算とGPU性能
  2. PGI CUDA Fortran を使用した行列積の計算とGPU性能
  3. PGI Fortran と CUDA CUBLASを使用した行列積の計算とGPU性能

Performance Summary

単精度 性能加速 29.7x

倍精度 性能加速 37.9x

 最初に、今回 PGIアクセラレータコンパイラを用いて得られた性能を示します。

表-1 MATMUL 計算 x64 Nehalem 1CPU 性能とGPU 性能の比較(4096サイズ)
使用プログラミングモデル Nehalem 1CPU GTX 480 性能比
PGI Accelerator 単精度版ソースプログラム 5.3 GFLOPS 157.8 GFLOPS 29.7倍
PGI Accelerator 倍精度版ソースプログラム 2.5 GFLOPS 94.7 GFLOPS 37.9倍

性能実測するシステム環境について

 まず始めに、性能実測に使用する私のテストマシン環境について説明します。

ホスト側プロセッサ Intel(R) Core(TM) i7 CPU 920 @ 2.67GHz x 1 プロセッサ
ホスト側メモリ PC3-10600(DDR3-1333)2GB x 2 枚
ホスト側 OS Cent OS 5.4、カーネル 2.6.18-164.el5
PGI コンパイラ PGI 10.6 (PGI 2010)
GPU ボード1 NVIDIA Geforce GTX480(Fermi)、1401.0 MHz clock, 1535.7 MB memory
GPU ボード2 NVIDIA Geforce GTX285、1476.0 MHz clock, 1023.3 MB memory
NVIDIA CUDA バージョン CUDA 3.0ドライバと3.0 Toolkit
表-1.1 NVIDIA GPU 理論性能の比較
ハードウェア仕様 Geforce GTX285 Geforce GTX480
Compute Capability Compute Capability 1.3 2.0
総CUDAコア(SP)数 CUDA cores
(Streaming Processors)
240 cores 480 cores
プロセッサクロック CUDA cores
(Streaming Processors)
1476 MHz 1401 MHz
単精度理論FP演算性能 Single Precision GFLOPS 708 GFLOPS 1344 GFLOPS
倍精度理論FP演算性能 Double Precision GFLOPS 88 GFLOPS 1344x1/8 = 168 GFLOPS

理論浮動小数点演算性能の求め方(GPGPU用語に関してはこちらのコラムへ

Geforce GTX285

  • 単精度: 30 SM x 8core x 2FP/clock x 1476MHz = 708 GFLOPS/s (2FP 積・和 = MAD命令)
  • 倍精度: 30 SM x 2FP/clock x 1476MHz = 88 GFLOPS/s (2FP 積・和 = FAD命令、倍精度演算器は 1個/SM )

Geforce GTX480

同じ Fermi アーキテクチャである 「Geforce GTX480」 と 「Tesla C2050」シリーズの違い

  • Fermi の倍精度演算能力は単精度に較べ 1/2 ですが、これが可能なのは Tesla C2000シリーズの GPU です。低価格帯の Geforce 470/480 では、1/8 に制限されていますので、ご注意下さい。

PGIアクセラレータ・プログラミングモデルを使用した実行例

 一般的な MATMUL の Fortran プログラムに PGI Accelerator 用のディレクティブを挿入してみましょう。以下に、そのサンプルプログラムを用意しましたので、これをダウンロードして実行して下さい。添付の Makefile は実行モジュールを作成するためのものですが、システムに実装している CUDA のドライバのバージョンによって、FFLAGS の切り替えを行って下さい。現時点(PGI 10.5 レベル)での PGI コンパイラが生成するデフォルトの CUDA ドライバは、CUDA 2.3 を利用することを想定した実行バイナリを作成します。もし、最新バージョン CUDA 3.0 ドライバを使用している場合は、明示的に -ta=nvidia:3.0 のオプションを指定して下さい。これに関しては、こちらのページで詳細に説明しています
PGI Workstation シリーズの Winodows 版をお使いの場合も以下の Makefile は使用できます

● Makefile の中の FFLAGS (PGIオプション)
# Please select by CUDA driver version
#CUDA 3.0 driver + CC1.0,1.3 + CC2.0
 FFLAGS = -fastsse -Minfo=accel -ta=nvidia:3.0,time 
#CUDA 2.3 driver + CC1.0,1.3 + CC2.0
#FFLAGS = -fastsse -Minfo=accel -ta=nvidia,time

 ※ PGIアクセラレータ用コンパイル・オプションについて

単精度版のソース・プログラム

 単精度版のプログラムを以下に示します。MATMULルーチンのコアな部分は、module matmul の中にコード化されています。このモジュール内の全体のループが NVIDIA GPU に処理させる並列領域となります。このループ構造は、i, j, k のインデックスによる繰り返し構造となっていますが、この中で、最外側 j ループと、最内側 i ループは、配列データ要素間のデータ依存性がありません。すなわち、並列化あるいはベクトル化(SIMID化)が可能です。しかし、中間のループである k ループでは、k の繰り返しにおいて独立にストアできる配列(左辺の配列)がありません。すなわち、k に関して並列(ベクトル)実行させることはできません。この三重ループの構成では、k ループの次元はシリアル実行させるような構成で、i と j のループ領域を使って並列実行主体であるブロック(16x16)を構成させることにより並列動作が可能となります。

 PGIアクセラレータ用のディレクティブの挿入の仕方は、まず、並列計算可能なループブロックの前に !$acc region (行番号7)を入れ、最後の end do 構文の後に !$acc end region (行番号20)を入れます。基本的には、これだけでPGIコンパイラがそのブロックに関してのベクトル・並列化分析、GPUへの適用性分析等を実施して、GPU用の並列コードを生成します。MATMULのような単純なアルゴリズムは、コンパイラ側で十分認識していますので、最適(高性能)なコードを生成できます。一般のプログラムでは、!$acc region の内部に、さらに最適化のためのアクセラレータ用のディレクティブを入れ、その性能を試行錯誤で検証してみる必要があります。以下の例では、!$acc do parallel,vector(16) (行番号8,14)の kernel tuning 用のディレクティブも入れて、SIMD長(ベクトル長)、ブロックサイズやグリッドサイズを調整しています。このサイズの調整で性能が大きく変化しますので、この調整が PGI アクセラレータ・プログラミングモデルにおけるチューニング方法の一つと言えましょう。

(    1)        module matmul
(    2)        contains
(    3)        subroutine mm( a, b, c, n, m, p)
(    4)        integer(4) :: m,p,n
(    5)        real(4), dimension(:,:) :: a, b, c
(    6)
(    7) !$acc region
(    8) !$acc do parallel,vector(16)
(    9)        do j = 1, m
(   10)           do i = 1, n
(   11)              a(i,j)= 0.0
(   12)           end do
(   13)           do k = 1, p
(   14) !$acc do parallel,vector(16)
(   15)              do i = 1, n
(   16)                 a(i,j) = a(i,j) + b(i,k)*c(k,j)
(   17)              enddo
(   18)           enddo
(   19)        enddo
(   20) !$acc end region
(   21)        end subroutine
(   22)        end module matmul
(   23)
(   24)        program GPUtest
(   25)        use matmul
(   26)
(   27) #if defined (_ACCEL)
(   28)         use accel_lib
(   29) #endif
(   30)
(   31)        integer(4) :: n,m,p
(   32)        integer(4) :: i,j,k
(   33)        character(10) :: arg
(   34)        integer(4) :: hz, clock0, clock1
(   35)        real(8) :: walltime, mflops
(   36)        real(4), dimension (:,:), allocatable :: a, b, c
(   37)
(   38)        n=0
(   39)        m=0
(   40)        p=0
(   41)
(   42)        args = command_argument_count()
(   43)        if( args .gt. 0 ) then
(   44)          call getarg( 1, arg )
(   45)          read (arg, '(I10)') n
(   46)          if( args .gt. 1 ) then
(   47)            call getarg( 2, arg )
(   48)            read (arg, '(I10)') m
(   49)            if( args .gt. 2 ) then
(   50)              call getarg( 3, arg )
(   51)              read (arg, '(I10)') p
(   52)            endif
(   53)          endif
(   54)        endif
(   55)
(   56)        if(n <= 0) n=100
(   57)        if(m <= 0) m=100
(   58)        if(p <= 0) p=100
(   59)        allocate( a(n,m) )
(   60)        allocate( b(n,p) )
(   61)        allocate( c(p,m) )
(   62)
(   63)        b= 3.0e0
(   64)        c= 4.0e0
(   65)
(   66)        call system_clock(count_rate=hz)
(   67)        print '(/1x,a,e12.6)', "system_clock resolution: ", real(1.d0/hz)
(   68)
(   69) #if defined (_ACCEL)
(   70)        ! GPU explicitly is initialized
(   71)        call system_clock(count=clock0)
(   72)        call acc_set_device( acc_device_nvidia )
(   73)        call acc_init( acc_device_nvidia )
(   74)        call system_clock(count=clock1)
(   75)        walltime= real((clock1-clock0)) / real (hz)
(   76)        print '(1x,a,F8.4,a)', 'initialize nvidia, initilize time = ',
(   77)      + walltime, " second"
(   78) #endif
(   79)
(   80)        call system_clock(count=clock0)
(   81)
(   82)          call mm(a,b,c,n,m,p)
(   83)
(   84)        call system_clock(count=clock1)
(   85)        walltime= real((clock1-clock0)) / real (hz)
(   86)
(   87)        mflops = (real(m)*real(p)*2*real(n)) / (walltime * 1.0d+06)
(   88)        print '(/1x,3(I5,1x,a,1x))', n,"x",m,"x",p,' matrix size'
(   89)        print '(3x,F10.6,a)', walltime, ' seconds on GPU'
(   90)        write (*,901) mflops*1.e-3," GFLOPS/s"
(   91) 901    format(5x,f8.3,a)
(   92)
(   93)        end program

コンパイラによるGPUアクセラレータ翻訳メッセージ

 PGIアクセラレータ・コンパイラは、GPIU並列化を行った際のアクセラレータに関するメッセージを出力できます。このメッセージの中には、ホスト~GPU間の並列・変数のデータ転送の状況、並列化,ベクトル化を行ったループ、ベクトル長(SIMD width),ブロックのサイズ・構成、GPU内のshared memory/local memoryの使用量、kernel 1スレッドのレジスタ量、GPU occupancy の割合 等が含まれます。以下はその例を示し、そのメッセージの意味を赤字で記しました。

● コンパイル時の GPU翻訳に関するメッセージ -Minfo=accel 
[kato@photon29 MAT]$ make
pgfortran -c -fastsse -Minfo=accel -ta=nvidia:3.0,time -Mlist matmul.F
mm:
      7, Generating copyout(a(1:n,1:m))
         Generating copyin(c(1:p,1:m))
         Generating copyin(b(1:n,1:p)) (Host~GPU memory間のデータ転送を行うコード生成)
         Generating compute capability 1.0 binary
         Generating compute capability 1.3 binary
         Generating compute capability 2.0 binary (Compute Capability 1.0 ~ 1.3、2.0対応コード生成)
      9, Loop is parallelizable
     10, Loop is parallelizable
         Accelerator kernel generated (行番号9-10 ループ:アクセラレータ CUDA kernel を生成する)
          9, !$acc do parallel, vector(16)
         10, !$acc do parallel, vector(16)
             CC 1.0 : 6 registers; 24 shared, 88 constant, 0 local memory bytes; 100 occupancy
             CC 1.3 : 6 registers; 24 shared, 88 constant, 0 local memory bytes; 100 occupancy
             CC 2.0 : 6 registers; 8 shared, 112 constant, 0 local memory bytes; 100 occupancy
     13, Loop carried reuse of 'a' prevents parallelization (kループでは "a配列" により並列化できない)
     15, Loop is parallelizable
         Accelerator kernel generated (行番号9-15 ループ:アクセラレータ CUDA kernel を生成する)
          9, !$acc do parallel, vector(16) (行番号9ループは、並列とベクトル用に分割される)
         13, !$acc do seq                  (行番号13ループは、sequetntialに実行される)
             Cached references to size [16x16] block of 'b' (16x16のキャッシュを shared memoryに生成)
             Cached references to size [16x16] block of 'c'
         15, !$acc do parallel, vector(16) (行番号9ループは、並列とベクトル用に分割される)
             Using register for 'a' (a配列は、高速レジスタに割り付けることができた)
             CC 1.0 : 17 registers; 2072 shared, 92 constant, 0 local memory bytes; 33 occupancy
             CC 1.3 : 17 registers; 2072 shared, 92 constant, 0 local memory bytes; 75 occupancy
             CC 2.0 : 22 registers; 2056 shared, 112 constant, 0 local memory bytes; 83 occupancy
	上記の occupancy に関しては、こちらに詳細に説明あり

実行、CUDAプロファイル情報

 ここでは、実行モジュール名を Linux 流の a.out と表示しています。Windows では、***.exe となります。このプログラム例では、実行モジュールの後に三つの引数を指定できます。行列のサイズ n, m, p を整数で指定します。以下の実行例は、各サイズ 4096 とした場合の例です。

[kato@photon29 MAT]$ ./a.out 4096 4096 4096

 system_clock resolution: 0.100000E-05
 initialize nvidia, initilize time =   0.5621 second

  4096 x  4096 x  4096  matrix size
     0.870584 seconds on GPU
      157.870 GFLOPS/s

(これ以降の情報は、-ta=time オプションを付すると表示される
 CUDAプロファイルによる各 kernel の実行、性能情報です。その意味に関しては、こちらへ

Accelerator Kernel Timing data
/home/kato/GPGPU/MAT/matmul.F
  mm
    7: region entered 1 time
        time(us): total=870576 init=1 region=870575
                  kernels=814264 data=56311
        w/o init: total=870575 max=870575 min=870575 avg=870575
        10: kernel launched 1 times
            grid: [256x256]  block: [16x16]
            time(us): total=1611 max=1611 min=1611 avg=1611
        15: kernel launched 1 times
            grid: [256x256]  block: [16x16] 16x16のブロックを 256x256のグリッドで並列実行
            time(us): total=812653 max=812653 min=812653 avg=812653
acc_init.c
  acc_init
    41: region entered 1 time
        time(us): init=562141

ちなみに、x64 CPU の Nehalem 1コアで実行してみると ...

 一般的な x64 CPU では、4096のサイズの行列積の性能はどの程度でしょうか?このサイズは、CPUの持つキャッシュのサイズを超えているため、キャッシュ効果が現れにくいものとなります。以下のように、最高の最適化を施しても 5 GFLOPS どまりです。

[[kato@photon29 MAT]$ pgf90 -fastsse -Minfo matmul.F
mm:
     10, Memory zero idiom, loop replaced by call to __c_mzero4
     15, Generated 3 alternate loops for the loop
         Generated vector sse code for the loop
         Generated 2 prefetch instructions for the loop
gputest:
     63, Memory set idiom, array assignment replaced by call to pgf90_mset4
     64, Memory set idiom, array assignment replaced by call to pgf90_mset4
[kato@photon29 MAT]$ ./a.out 4096 4096 4096

 system_clock resolution: 0.100000E-05

  4096 x  4096 x  4096  matrix size
    25.833588 seconds on GPU
        5.320 GFLOPS/s

x64 Nehalem 1CPU 性能と GPU 性能の比較(Fortran)

 ここで、一般的な x64 CPU の性能と GPU の性能の実効性能の違いはどの程度あるかを把握しておきましょう。CPU/GPU 間のこの絶対的な性能感覚は、今後の自分のアプリケーションに対して GPGPU 計算を適用するかどうかの判断において必要なものです。以下の表は、4096の行列サイズの時のCPUとGPUの性能を示したものです。

表-2 MATMUL 計算 x64 Nehalem 1CPU 性能とGPU 性能の比較(4096サイズ)
GPU/行列サイズ Nehalem 1CPU GTX 480 性能比
PGI Accelerator 単精度版 5.3 GFLOPS 157.8 GFLOPS 29.7倍
PGI Accelerator 倍精度版 2.5 GFLOPS 94.7 GFLOPS 37.9倍

PGI アクセラレータによるMATMUL性能(Fortran)

 単精度版と倍精度版のプログラムを使用して、いくつかの行列サイズで行列積の計算を行った結果を示します。ここでは、最新の Fermi アーキテクチャである GTX 480 (Compute Capability :2.0) と1世代古い GTX 285 (Compute Capability :1.3)を使用して実測を行いました。表-3 は、単精度版での GTX 480 と GTX285 の性能比較を表し、表-4 では、倍精度版での性能比較を表しています。
 また、各 GPU ボード上で単精度と倍精度の実行性能の比較を行ったのが 表-5 と 表-6 である。表-5 は、GTX480 の性能、表-6 は、GTX285 の性能です。

単精度版 MATMUL Fortran プログラム性能(GTX480とGTX285 GPUの性能比較)

 GTX 480 (Compute Capability :2.0) と1世代古い GTX 285 (Compute Capability :1.3)のアーキテクチャの大きな違いは、Fermi でその Warp命令の issue タイミングの改善やCUDAコア数、演算のパイプラインが 2 倍となっています。従って、理想的には最大 2 倍程度までの性能改善が期待されます。並列スレッドが多くなる行列サイズ(4096)では 1.5 倍程度までの性能差が現れています。。高級言語のソースプログラムからのPGIアクセラレータによる GPU コード化ですので、Fermiでの 157.8 GFLOPS (4096サイズ)の実効性能はなかなかのものです。とは言え、両方の GPU 共に、理論ピーク性能の約 10~14 % 程度の実効性能が得られていることになります。

表-3 単精度 MATMUL 計算 PGIアクセラレータコンパイラによる性能(GFLOPS)
GPU/行列サイズ 256 512 1024 2048 4096
GTX 480 (CC 2.0) 38.5 76.9 116.3 147.2 157.8
GTX 285 (CC 1.3) 31.2 60.7 86.0 102.7 108.0
性能Ratio 1.23 1.27 1.35 1.44 1.47

倍精度版 MATMUL Fortran プログラム性能(GTX480とGTX285 GPUの性能比較)

 倍精度演算ではどうでしょう。GTX 480 では 94.7 GFLOPS の実効性能が記録されました。これは、理論ピーク性能の約 56 % 程度の実効性能が得られていることになります。一方、GTX 285 も GTX480 も、倍精度演算では単精度演算器に較べて 1/8 の演算リソースしかありません。GTX 285 の場合は、実性能としては 52.6 GFLOPS の実効性能で、理論ピーク性能の約 59% の実効性能が得られています。これは、少々驚きの事実です。どちらの GPU も、倍精度の実効的な演算性能は高いと言うことが分かります。

表-4 倍精度 MATMUL 計算 PGIアクセラレータコンパイラによる性能(GFLOPS)
GPU/行列サイズ 256 512 1024 2048 4096
GTX 480 (CC 2.0) 23.4 50.0 74.2 87.6 94.7
GTX 285 (CC 1.3) 18.1 33.9 45.2 50.1 52.6
性能Ratio 1.29 1.47 1.64 1.75 1.80

GTX 480 (Fermi)による単精度/倍精度の実効性能比較

 GTX 480 の単精度演算と倍精度演算の性能比を示したものです。GTX 480 のアーキテクチャは、演算パイプラインが単純に 1/2 の物量ですが、政策的に そのピーク性能をさらに 1/4 に制限していますので、本来は、理論性能的に倍精度性能は単精度演算の 1/8 程度になります。しかし、以下の比較を見ますと、理論値ほどの劇的な低下はありません。むしろ、倍精度演算はスペックに較べてかなり高い通知を示しています。単精度演算/倍精度演算の peformance ratio としては 1.6~1.7 と言うところになっているようです。

表-5 PGIアクセラレータコンパイラによるMATMUL性能(GFLOPS)
精度/行列サイズ 256 512 1024 2048 4096
単精度計算 38.5 76.9 116.3 147.2 157.8
倍精度計算 23.4 50.0 74.2 87.6 94.7
性能Ratio 1.65 1.54 1.56 1.68 1.67

GTX 285 による単精度/倍精度の実効性能比較

 GTX 285 の単精度演算と倍精度演算の性能比を示したものです。GTX 285 のアーキテクチャは、その理論性能比が倍精度/単精度の演算器の物量比として表され 1/8 となっていますので、倍精度の性能は単精度に比べかなり低いものと想定されていました。しかし、実際は、単精度演算/倍精度演算の peformance ratio としては 1.7~2.0 と言うところになっているようです。GTX285による倍精度計算は、思った程性能が落ちないのではないかと言うことが予想できます。これは、単に高級言語レベルのソースを翻訳して CUDA バイナリに落とし込んだ as is 性能に近いものであり、GPU のアーキテクチャを有効に利用する究極の最適化を施した性能ではないことによるものと考えた方が良いでしょう。普通にユーザが CUDA プログラムを作っても、高度な最適化(Shared memory を使用するとか...)を施していない場合は、単精度と倍精度の性能差は、こんなものと思った方が良いでしょう。GTX 285(CC 1.3)でも、その倍精度性能はそれ程悲観するようなものでもないと言うことでしょうか。

表-6 PGIアクセラレータコンパイラによるMATMUL性能(GFLOPS)
精度/行列サイズ 256 512 1024 2048 4096
単精度計算 31.2 60.7 86.0 102.7 108.0
倍精度計算 18.1 33.9 45.2 50.1 52.6
性能Ratio 1.72 1.79 1.91 2.04 2.05

 行列積のプログラムに関するコラム(全3回)

  1. PGIアクセラレータを使用した行列積の計算とGPU性能
  2. PGI CUDA Fortran を使用した行列積の計算とGPU性能(本稿)
  3. PGI Fortran と CUDA CUBLASを使用した行列積の計算とGPU性能