PGI Fortran と CUDA CUBLASを使用した行列積の計算とGPU性能

　NVIDIAが提供する CUDA CUBLAS ライブラリを使用して、行列積を計算するプログラムを書いてみましょう。CUBLASは、NVIDIAの GPU に最適化された BLAS ライブラリですので、これを使えばある程度最高と言われる性能を享受できるはずです。ここでは、前2回のコラムでご紹介した PGI アクセラレータ・プログラミングモデルでのプログラムによる性能と PGI CUDA Fortran でコーディングした場合の性能を比べるために、CUBLAS の SGEMM / DGEMM ルーチンの性能を示すことにします。結論から言うと、行列サイズが 4096位までは CUBLASの性能と PGI CUDA Fortran によるプログラムによる性能は、似たような傾向の性能を出すことができることが分かりました。ただ、CUBLAS の方は、行列サイズがさらに大きくなると、Fermi の場合はまだ性能が伸びていくようです。

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

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

MATMUL Performance Summary using NVIDIA GPUs

　以下の表-1 に、本稿の③PGI Fortran による CUBLAS の利用による性能、 ②PGI CUDA Fortran を使って測定した性能と、①PGIアクセラレータ・プログラミングモデルを用いて得られた性能を比較するために記載します。

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

　一般的な MATMUL は、BLAS3 ルーチンである SGEMM（単精度） あるいは DGEMM（倍精度） ルーチンで代用できます。本稿では、CUDA Toolkit で提供されている NVIDIA CUDA の BLAS ルーチン（CUBLAS）を Fortran プログラムからリンクする方法を説明します。また、これで得られた MATMUL の性能も呈示します。まず、最初に一般的に BLAS を使用する場合のコーディング上の注意を述べます。
　NVIDIA の CUBLAS ルーチンは、C言語で記述されているため、C言語によるライブラリと言うことになります。しかし、本来の BLAS ルーチン起源は Fortran 言語で記述されていたことにより、特に多次元配列のメモリ上の配置（順並び）は、Column-major Format (列方向に並ぶ順序）に則った形で行われます。この並びは、一般に Fortran 言語でのメモリ配置オーダーと同じであるため、Fortran-order とも言うこともあります。この並びを簡単に言えば、多次元配列の「最初のインデックス（添字）を順番にメモリ上に配置してゆく形を取ります。具体的には、A(i,j) と言う配列であれば、そのメモリ上のデータ配置は、j=1 固定で i が 1,2,3, ... と言う順番で配置されます。一方、一般的な C/C++ 言語プログラムの場合は、Row-major Format による配置と呼ばれ、Fortran とは全く逆の配置でセットされます。上記の例では、i=1 固定で j が 1,2,3, ... と言う順番で配置されます。BLAS ルーチンの場合のルールは、Fortran-order で BLAS 内部のロジックが構築されていますので、Fortran でプログラムで組む場合は、こうしたメモリ配置については何も意識しないでBLASルーチンへの配列の引き渡しが可能です。また、BLAS 関数のデフォルトは、配列の添字が 1 から始めるものとして扱います。これも Fortran 流です。Cプログラムの場合の配列はデフォルト 0 から始まるものとして考えますので、Cプログラムから BLAS を使用する場合は、注意が必要です。一般に、Cプログラムから BLAS を使う場合は、行と列のインデックスの入れ替え（転置）する形をとって、BLAS ルーチンに配列を渡す必要があります。また、配列の添字の開始が0か1かにより、BLASルーチンへ渡す際のインデックス変換を行わなければなりません。例えば、次の例のようにマクロを定義します。

【C言語プログラムから blas を使用する場合の配列のインデックス変換】
　(m x n)行列の場合で x を行インデックス、 y を列インデックスとする。 m は行のサイズとする。
　#define COLM(x, y, m) ( (y-1)*(m) + ((x)-1) )   //column-major
　double* A = (double *)malloc(m*n, sizeof(double));
　m=5 の行列で  A[4][2] を表す場合は、
　A[COLM(4, 2, 5)] = ....	

　本稿は、Fortran プログラムから CUBLAS を使うため、上記に述べた C 言語特有な問題はありませんので、素直に引数の引き渡しを行えます。さて、以下に、サンプルプログラムと Makefile を用意しました。ダウンロード可能です。PGI CUDA Fortran の API をプログラムの中で使用していますので、コンパイル時に -Mcuda を入れて下さい。また、NVIDIA の CUBLASライブラリをリンクするために -lcublas を必要とします。

• 単精度版 PGI CUDA Fortranソースコード　Sgemm_cublas.F90 と Makefile

• 倍精度版 PGI CUDA Fortranソースコード　Dgemm_cublas.F90 と Makefile

● Makefile の中の FFLAGS （PGIオプション）
# Please select by CUDA driver version
# GPU CUDA3.0+CC20,CC13
 FFLAGS = -fastsse -Mcuda:3.0 -D_CUDAFOR -Minfo -lcublas
# GPU CUDA2.3+CC13,CC20
#FFLAGS = -fastsse -Mcuda -D_CUDAFOR -Minfo -lcublas

　※ PGI CUDA Fortran 用コンパイル・オプションについて

単精度版のソース・プログラムと CUBLASルーチン（C言語ライブラリ）のバインドの方法

　 従来から、Fortran と C あるいは、C++ 間のサブルーチン間の連結には、変数型の定義や引数の引き渡し等のルールを守り、慎重にプログラム・コーディングを行う必要があります。こうした背景があったために、Fortran 2003 で言語間の相互運用性を高める機能として、組込モジュール ISO_C_BIND がサポートされました。今回、NVIDIA の CUBLAS（C言語ライブラリ）と Fortran ルーチンとの相互運用を行うために、ISO_C_BIND を使用します。とは言っても、ごく一般的な標準的な設定をするだけです。引数インタフェースの定義を明示的に行うことと、BLAS 標準関数名（一例、SGEMM)と CUBLAS のそれに相当する関数名（cublasSgemm）は異なっているため、この両者関数シンボル名を同一のものとしてバインドします。これは、以下の15行目で行われている bind(c,name='cublasSgemm') で行うことができます。ISO_C_BIND モジュールで一番便利な機能と言えば、この名前のバインド機能と言うこともできます。後は、17～20行目に定義しているような Fortran の型に相当する C の型の定義（種別パラメータ）を行って変数の明示的な相互運用が可能となります。このように interface ブロックを ISO_C_BIND で定義するだけで、Fortran BLAS ルーチン名と CUBLAS ルーチンの連結ができるようになりました。 （CUDA CUBLASのマニュアル）http://developer.download.nvidia.com/compute/cuda/3_1/toolkit/docs/CUBLAS_Library_3.1.pdf

　Interface の定義が終われば、後は普通に BLAS ルーチンを Fortran で使用するだけです。配列の添字のスタートも上述したように Fortran の慣習である 1 から始まるものとして指定すれば良いです。また、SGEMM ルーチンの引数のコンベンションは、行番号90 で示すような形となり、今まで使用してきた BLAS の呼び出しと全く同じように使うことができます。その他に何か、足りないなとお気づきかもしれませんが、例えば、GPU 側のデバイスメモリ上の配列の割付とかホスト側とデバイス間のデータ転送はどのように記述するのでしょうか？　以下のプログラム例の中で説明していますが、PGI CUDA Fortran を使えば、Fortran の構文レベルで非常に簡単にこれらの動作を記述できます。CUDA C の場合は、一つひとつの動作を CUDA API 関数を使用して記述しなければ行けませんが、PGI CUDA Fortran は Fortranネイティブな構文記述でこうした動作を記述できます。（もちろん CUDA Fortran でも CUDA API関数を用意しており、CUDA C と同じような記述もできます）。

PGF90 (Version     10.6)          06/29/2010  21:00:50      page 1

(    1) # 1 "Sgemm_cublas.F90"
(    1) !
(    2) !      Copyright 2010, SofTek Systems, Inc.
(    3) !      All rights reserved.
(    4) !
(    5) !      Using CUDA BLAS sgemm routine (Single Precision)
(    6) !      Fortran 2003 iso_c_binding interface version
(    7) !      SGEMM : http://www.netlib.org/blas/sgemm.f
(    8) !
(    9) program Sgemm_cublas
(   10)
(   11) # 12 (Fortran Interface文を使用して、CUBLASのC言語ライブラリとの連結） 
(   12)   use cudafor
(   13)   ! ISO_C_BIND interface (symbol conversion)
(   14)   interface
(   15)    subroutine SGEMM(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc ) &
                      bind(c,name='cublasSgemm')
(   16)     use iso_c_binding
(   17)     integer(c_int), value :: m, n, k, lda, ldb, ldc
(   18)     real(c_float), device, dimension(m,n) :: a, b, c
(   19)     real(c_float), value :: alpha, beta
(   20)     character(kind=c_char), value :: transa, transb
(   21)    end subroutine SGEMM
(   22)   end interface
(   23)
(   24)   ! GPU device arrays CUDA Fortran でデバイス割付配列の宣言 
(   25)   real(4), device, allocatable, dimension(:,:) :: dA, dB, dC
(   26)   ! device properties
(   27)   type(cudaDeviceProp) :: prop
(   28)   ! device number
(   29)   integer idev
(   30)   ! command argument
(   31)   character*20 arg
(   32)
(   33) # 34
(   34)   real(4), allocatable, dimension(:,:) :: a, b, c
(   35)   real(4) :: alpha = 1.0e0, beta  = 0.0e0
(   36)   ! for timing
(   37)   real(8) :: GFLOPS
(   38)   real(8) :: time
(   39)   integer(4) :: hz, clock0, clock1
(   40)   integer :: i, j, k
(   41)
(   42) ! Matrix Size (N x N x N)
(   43)   print *, "What size of matrix : ?"
(   44)   read(5,*) n
(   45)   print *, "MatSize = ", n, " x ", n, " x ",n
(   46)
(   47) # 48
(   48) ! command argument and GPU device set
(   49)   nargs = command_argument_count()
(   50)   idev = 0
(   51)   do i = 1, nargs
(   52)     call get_command_argument(i,arg)
(   53)     if ((arg(1:7) .eq. "-device") .and. (i.lt.nargs)) then
(   54)       call get_command_argument(i+1,arg)
(   55)       read(arg,'(i2)') idev
(   56)     end if
(   57)   end do
(   58)   ! set a GPU device and print the properties
(   59)   istat = cudaSetDevice(idev)
(   60)   istat = cudaGetDeviceProperties(prop,idev)
(   61)   ilen = verify(prop%name, ' ', .true.)
(   62)   write (*,900) prop%name(1:ilen), &
(   63)                 real(prop%clockRate)/1000.0, &
(   64)                 real(prop%totalGlobalMem)/1024.0/1024.0
(   65)
(   66) # 67
(   67)   ! arrays for HOST
(   68)   allocate(a(n,n), b(n,n), c(n,n))
(   69) # 70 !  CUDA Fortran でデバイス割付配列の割付け 
(   70)   ! arrays for GPU
(   71)   allocate(dA(n,n), dB(n,n), dC(n,n))
(   72)
(   73) # 74
(   74)   a = 1.0e0
(   75)   b = 5.5e0
(   76)
(   77)   ! Initialize time function
(   78)   call system_clock(count_rate=hz)
(   79)   print *, "system_clock resolution: ", real(1.d0/hz)
(   80)   ! Start a clock
(   81)   call system_clock(count=clock0)
(   82)
(   83) # 84 
(   84)   ! for GPU  CUDA Fortran でHOSTからデバイス配列へデータコピー 
(   85)   dA = a
(   86)   dB = b
(   87)   if (beta .ne. 0.0d0) then
(   88)     dC = c
(   89)   endif
(   90)   call SGEMM('n', 'n', n, n, n, alpha, dA, n, dB, n, beta, dC, n)
(   91)   c = dC  ! デバイス配列からホスト側へデータの戻し 
(   92) # 96
(   96)   ! Stop a clock
(   97)   call system_clock(count=clock1)
(   98)   time= real((clock1-clock0)) / real (hz)
(   99)
(  100) print *, "*** Verify results ***"
(  101)
(  102) do j = 1, n
(  103)   do i = 1, n
(  104)     if (c(i,j) .ne. (5.50*real(n))) then
(  105)       print *, "error:  ",i,j,c(i,j)
(  106)     endif
(  107)   enddo
(  108) enddo
(  109)
(  110)   GFLOPS = (real(n) * real(n) * real(n) * 2.0)
(  111)   print *, "Total Time (sec): ",time
(  112)   print *, "Total SGEMM GFLOPS: ",GFLOPS/(time*1.0d9)
(  113)   print *, ""
(  114)
(  115) 900 format('\nDevice:',a,', ',f6.1,' MHz clock, ',f6.1,' MB memory.\n')
(  116) end

コンパイラ翻訳メッセージの例

　PGI CUDA Fortran コンパイラを使用(-Mcudaオプションを入れる）してコンパイルし、-lcublas で CUBLAS ライブラリをリンクします。

●コンパイル・リンク メッセージ -Minfo 

pgfortran -c -fastsse -Mcuda:3.0 -D_CUDAFOR -Minfo -lcublas Sgemm_cublas.F90
sgemm_cublas:
     51, Loop not vectorized/parallelized: contains call
     74, Memory set idiom, array assignment replaced by call to pgf90_mset4
     75, Memory set idiom, array assignment replaced by call to pgf90_mset4
    103, Loop not vectorized/parallelized: contains call
pgfortran -o a.out Sgemm_cublas.o -fastsse -Mcuda:3.0 -D_CUDAFOR -Minfo -lcublas -Mlist

二つの GPU ボードで実行してみる

　ここでは、実行モジュール名を Linux 流の a.out と表示しています。Windows では、***.exe となります。このプログラム例では、実行モジュールの後に引数を指定できます。-device 0 と言った形で、GPU ボードの論理番号を指定します。私のマシンの場合は、二つのボードが搭載されていますので、以下のような 0, 1 の指定で、各 GPU の性能が測定できます。以下の例は、N=4096 の場合と N=5120 の場合の二つを測定しています。GTX480 は、N=5120 で 525 GLFLOPS となり、500 GFLOPSオーダーの性能を確保できるようです。一方、GTX285では、同様な条件で、378 GFLOPS となり、この辺りが飽和性能のようです。従って、Fermi のアーキテクチャの方が、GPU内の演算パイプラインを有効に使用できる余裕があるということになります。

GPU デバイス 0 (GTX480) を使う場合は、-device 0 引数を付ける

[kato@photon29 CUBLAS]$ ./a.out -device 0
 What size of matrix : ?
4096
 MatSize =          4096  x          4096  x          4096

Device:GeForce GTX 480, 1401.0 MHz clock, 1535.7 MB memory.

 system_clock resolution:    1.0000000E-06
 *** Verify results ***
 Total Time (sec):    0.3679809868335724
 Total SGEMM GFLOPS:     373.4947141009756
 
 [kato@photon29 CUBLAS]$ ./a.out -device 0
 What size of matrix : ?
5120
 MatSize =          5120  x          5120  x          5120

Device:GeForce GTX 480, 1401.0 MHz clock, 1535.7 MB memory.

 system_clock resolution:    1.0000000E-06
 *** Verify results ***
 Total Time (sec):    0.5108109712600708
 Total SGEMM GFLOPS:     525.5083995902089 (N=5120で500GFLOPS越え）
 
GPU デバイス 1 (GTX285) を使う場合は、-device 1 引数を付ける

[kato@photon29 CUBLAS]$ ./a.out -device 1
 What size of matrix : ?
5120
 MatSize =          5120  x          5120  x          5120

Device:GeForce GTX 285, 1476.0 MHz clock, 1023.3 MB memory.

 system_clock resolution:    1.0000000E-06
 *** Verify results ***
 Total Time (sec):    0.7090460062026978
 Total SGEMM GFLOPS:     378.5867964162276

[kato@photon29 CUBLAS]$ a.out -device 1
 What size of matrix : ?
5120
 MatSize =          5120  x          5120  x          5120

Device:GeForce GTX 285, 1476.0 MHz clock, 1023.3 MB memory.

 system_clock resolution:    1.0000000E-06
 *** Verify results ***
 Total Time (sec):    0.7088969945907593
 Total SGEMM GFLOPS:     378.6663761425110

　ここで、一般的な x64 CPU の性能と CUBLAS ライブラリによる GPU の性能の実効性能の違いはどの程度あるかを把握しておきましょう。以下の表は、4096 の行列サイズの時の CPU と GPU の性能を示したものです。なお、以下の表には載せていませんが、GTX480(Fermi)では N=5120の場合、525 GFLOPSの性能が実測されました。

　単精度版と倍精度版のプログラムを使用して、いくつかの行列サイズで CUBLAS の行列積の計算を行った結果を示します。ここでは、最新の 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 の性能です。

単精度版 SGEMM Fortran プログラム性能（GTX480とGTX285 GPUの性能比較）

　GTX 480 (Compute Capability :2.0) と1世代古い GTX 285での単精度版の性能は、サイズ N=4096 までは、大差ない性能結果となっております。しかし、以下の表には載せていませんが、N=5120 のサイズになると前者が 525 GFLOPS で後者が 378 GFLOPS となり、両者の性能差が現れてきます。 Active Warp が増えることにより、Fermi のアーキテクチャのマルチスレディングが活きてきていると言うことになります。なお、GTX 480の実効 525 GFLOPS は理論ピークの約 40% 弱、 GTCX 285の 378 GFLOPS は、53％となり、GTX 480（Fermi) では Active Warp が増えた場合、さらに性能が伸びるのではないかと思います。

倍精度版 DGEMM Fortran プログラム性能（GTX480とGTX285 GPUの性能比較）

　倍精度演算ではどうでしょう。GTX 480 では 147 GFLOPS の実効性能が記録されました。これは、理論ピーク性能の約 88 % の実効性能が得られていることになります。一方、GTX 285 も、倍精度演算では単精度演算器に較べて 1/8 の演算リソースしかありませんが、実性能としては 80 GFLOPS の実効性能であり、理論ピーク性能の約 91% の実効性能が得られています。両者とも倍精度の性能効率が良いと言うことです。特に、GTX 480 は Fermi でありながら、政策的にピーク性能を単精度の 1/8 に落としていますが、実効性能的にはそれ程の落ち込みは見られないと言うことでしょうか。Tesla 2000シリーズでどれだけの性能が出るのかは見物です（残念ながら弊社では測定できません）。

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

　GTX 480 の単精度演算と倍精度演算の性能比を示したものです。先のコラムの PGI CUDA Fortran プログラムでのコメントと全く同じです。CUDA BLAS によりネイティブなレベルで最大限最適化された状態では、単精度演算/倍精度演算の peformance ratio としては 2倍 以上 の性能差があると言うことです。ハードウェアのスペック（1/8 のピーク性能）から見たら、倍セ氏度はそれなりの高い性能を呈示しています。なお、PGI CUDA Fortran プログラムの性能と CUBLAS ライブラリ性能を比較した場合、行列のサイズが小さい場合の立ち上がり性能は、CUBLASも PGI CUDA Fortran もほぼ同じですが、サイズをさらに大きくすると CUBLAS による性能が伸びていくようです。

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

GTX 285 の単精度演算と倍精度演算の性能比を示したものです。先のコラムの PGI CUDA Fortran プログラムでのコメントと全く同じです。GTX 285 のアーキテクチャは、その理論性能比が倍精度/単精度の演算器の物量比として表され 1/8 となっていますが、実際は、単精度演算/倍精度演算の peformance ratio としては 3～5 と言うところになっているようです。これらの性能は、CUDA Fortran によるネイティブなレベルで最適化されたプログラム特性の場合の性能です。すなわち、GPUのアーキテクチャを有効に利用するような高度な最適化をソースレベルで施した場合の性能ですが、単精度/倍精度の理論的な性能比である 8 倍には届きません。プログラムを最適化すると両者の性能差は顕著になると言うことです（理論値に近づくと言うことです）