PGI CUDA Fortran による CUBLAS 4.0 以降の使用方法

キーワード GPGPU、CUDA、SGEMM、MATMUL、PGI CUDA Fortran、CUBLAS 4.0

 NVIDIA® CUDA ランタイムの上位ライブラリである CUBLAS に、CUDA 4.0 から新しいインタフェース(API) が導入された。PGI 11.8 以降に提供する CUDA Fortran は、この新しいインタフェースも含めた CUBLAS 用の Fortran Module を提供した。この結果、NVIDIA C 言語の CUBLAS へのインタフェースを明示的に行う必要(ISO_C_BIND interface の必要)がなくなった。直接、Fortran 言語から BLAS 関数をコールすることができる。このコラムでは、その概要を説明する。
2011年10月21日 © 株式会社ソフテック 加藤

CUDA Fortran 上の CUBLAS 4.0 のインタフェース

 PGI 11.8 リリース以降において、PGI CUDA Fortran から CUBLAS ライブラリへのインタフェースを定義したモジュールを提供しました。これらのインタフェースは、CUDA Fortran のホストコードとなるプログラムの中に、次の use 文を置くことによりアクセスすることが可能です。なお、PGI 11.7 以前の CUDA Fortran では、ISO_C_BIND 機能を使用して明示的なインタフェースを作成する必要がありました。この方法に関しては、以前のコラムをご覧下さい。

use cublas

このインタフェースは、現在のところ、以下の 3 つのフォームで提供しています。

  1. 従来から使用されている(伝統的な) BLAS の呼び出しインタフェース。但し、配列引数は、ホスト側の配列ではなく、デバイス配列とする形態です。レガシーなプログラムをポーティングする際に使用すると便利です。但し、このレガシーな API は thread safety ではありません。すなわち、複数のホストスレッドとデバイスを使用する時には、スレッドセーフでないこうした関数は使用できません。thread safety な CUBLAS 関数は、以下の CUDA 4.0 以降の新しいインタフェースで使用可能となります。
    call saxpy(n, a, x, incx, y, incy)
        但し、x, y 引数はデバイス配列とする
  2. NVIDIA CUBLAS 4.0 より前のバージョンの Portable legacy CUBLAS インタフェースを直接使用する方法。但し、このレガシーな API は thread safety ではありません。
    call cublasSaxpy(n, a, x, incx, y, incy)
        但し、x, y 引数はデバイス配列とする
  3. 新しいライブラリの全ての機能にアクセス可能な、新 CUBLAS 4.0 のインタフェース。このインタフェースを使用する場合は、全て function call の形態をとります(call 文での呼び出しではありません)。最初の引数には、CUBLAS 関数のコンテキストへのハンドルとなる、派生型変数である handle (以下の例では h 変数) を置きます。このハンドルの使用は、マルチスレッド、マルチGPU上で使用する際に必須となります。この新しいインタフェースの使用法に関しては、CUDA Fortran Programming Guide and Reference の「Fortran Modules」の項をご参照下さい。また、NVIDIA社の CUBLAS Library User Guide も併せてご覧下さい。
    istat = cublasSaxpy_v2(h, n, a, x, incx, y, incy)

新 CUBLAS 4.0 のインタフェースについて(NVIDIA CUBLAS では、cublas_v2.h に含まれる I/F)

 いわゆる、CUBLAS version2 インタフェースで変化があった部分を説明します。例えば、上記の例の saxpy のケースの場合、引数 "a"(scalar used for multiplication)は、ホストあるいは、デバイスのどちらに置かれても良いような仕様となりました。こうした引数の詳細は、NVIDIA社の CUBLAS のマニュアルの各関数の引数の説明を見ると記載されています。また、sdot() や isamax() のようなスカラ値の結果を返す関数は、引数列の最後に新たな返値の引数を新設し、こうした引数形態でも使用できます。また、Fortran 上で伝統的に使用されてきた character*1 の属性を持つ引数、例えば、配列の転置形態を制御するパラメータは、整数値に変更されました。"n" は 0、 "t" は 1、"c" は 2 と言う形に変わりました。これらのパラメータは、cublas module の中に設定されています。最後に、ハンドル変数についてです。ハンドルに使われる派生型 cublasHandle は、cubals モジュールに定義されていますが、Fortran 上では、type(cublasHandle) :: h と言う形で宣言します。このハンドルは、一般に、cublasCreate(handle) 関数で初期化します。また、CUDA Fortran では、同様な機能として、h = cublasGetHandle()と言う関数も用意しています。

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

 性能実測に使用する私のテストマシン環境は以下の通りです。

ホスト側プロセッサ Intel(R) Core(TM) i7 CPU 920 @ 2.67GHz x 1 プロセッサ
ホスト側メモリ PC3-10600(DDR3-1333)2GB x 2 枚
ホスト側 OS Scientific Linux 6.1 (Red Hat クローン)
PGI コンパイラ PGI 11.9 (PGI 2011)
GPU ボード1 NVIDIA GeForce GTX 580、1544 MHz clock, 1609 MB memory, 512 cores
NVIDIA CUDA バージョン CUDA 4.0ドライバと4.0 Toolkit

CUDA Fortranで、CUBLAS を使用するプログラム例

 以下に、PGI 11.8 以降の cublas モジュールを利用して CUBLAS 関数を使用する際のプログラム例を示しました。上記で述べた3つのインターフェースを使用して記述しています。コンパイルの方法は、以下のオプションを使用しますが、ここでは、CUDA 4.0 以上のツールキット(ライブラリ)を使用するようにしなければなりません。

● コンパイルオプション
 pgfortran -fastsse -Mcuda -Minfo -lcublas Sgemm_cublas.cuf

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

!
!      Copyright 2011, SofTek Systems, Inc.
!      All rights reserved.
!
program Sgemm_cublas

  use cudafor
  use cublas   ! for PGI 11.8 and later

  ! GPU device arrays
  real(4), device, allocatable, dimension(:,:) :: dA, dB, dC
  ! device properties
  type(cudaDeviceProp) :: prop
  ! device number
  integer idev
  ! command argument
  character*20 arg
  
  ! handle for CUBLAS
  type(cublasHandle) :: h1
  integer(4) :: iface = 0

  real(4), allocatable, dimension(:,:) :: a, b, c
  real(4) :: alpha = 1.0e0, beta  = 0.0e0
  ! for timing
  real(8) :: GFLOPS
  real(8) :: time
  integer(4) :: hz, clock0, clock1
  integer :: i, j, k

! BLAS interface
  print *, "Select BLAS interface? 0=traditional, 1=Portable legacy CUBLAS <4.0, 2=CUBLAS Version2"
  read (5,*) iface
  print *, "BLAS interface =", iface
  if (iface > 3) stop

! Matrix Size (N x N x N)
  print *, "What size of matrix : ?"
  read(5,*) n
  print *, "MatSize = ", n, " x ", n, " x ",n

! command argument and GPU device set
  nargs = command_argument_count()
  idev = 0
  do i = 1, nargs
    call get_command_argument(i,arg)
    if ((arg(1:7) .eq. "-device") .and. (i.lt.nargs)) then
      call get_command_argument(i+1,arg)
      read(arg,'(i2)') idev
    end if
  end do
  ! set a GPU device and print the properties
  istat = cudaSetDevice(idev)
  istat = cudaGetDeviceProperties(prop,idev)
  ilen = verify(prop%name, ' ', .true.)
  write (*,900) prop%name(1:ilen), &
                real(prop%clockRate)/1000.0, &
                real(prop%totalGlobalMem)/1024.0/1024.0

  ! arrays for HOST
  allocate(a(n,n), b(n,n), c(n,n))
  ! arrays for GPU
  allocate(dA(n,n), dB(n,n), dC(n,n))

  a = 1.0e0
  b = 5.5e0

  ! Initialize time function
  call system_clock(count_rate=hz)
  print *, "system_clock resolution: ", real(1.d0/hz)
  ! Start a clock
  call system_clock(count=clock0)

  ! for GPU
  dA = a
  dB = b
  if (beta .ne. 0.0d0) then
    dC = c
  endif

  if (iface == 0) then
    ! traditional BLAS interface
    call sgemm('n', 'n', n, n, n, alpha, dA, n, dB, n, beta, dC, n)

  else if (iface == 1) then
    ! Portable legacy CUBLAS interface
    call cublasSgemm('n', 'n', n, n, n, alpha, dA, n, dB, n, beta, dC, n)

  else if (iface ==2) then
    ! New CUBLAS 4.0 interface
    ! cublasOperation_t is defined as a integer namuber.
    ! "n" : CUBLAS_OP_N=0, "t" : CUBLAS_OP_T=1,  "c" : CUBLAS_OP_C=2
    h1 = cublasGetHandle()
    istat = cublasSgemm_v2(h1, 0, 0, n, n, n, alpha, dA, n, dB, n, beta, dC, n)

  end if

  ierr = cudathreadsynchronize()
  c = dC

  ! Stop a clock
  call system_clock(count=clock1)
  time= real((clock1-clock0)) / real (hz)

print *, "*** Verify results ***"

do j = 1, n
  do i = 1, n
    if (c(i,j) .ne. (5.50*real(n))) then
      print *, "error:  ",i,j,c(i,j)
    endif
  enddo
enddo

  GFLOPS = (real(n) * real(n) * real(n) * 2.0)
  print *, "Total Time (sec): ",time
  print *, "Total SGEMM GFLOPS: ",GFLOPS/(time*1.0d9)
  print *, ""

900 format('\nDevice:',a,', ',f6.1,' MHz clock, ',f6.1,' MB memory.\n')
end

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

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

●コンパイル・リンク メッセージ -Minfo 
[kato@photon29]$ pgfortran -o a.out Sgemm_cublas.cuf -fastsse -Mcuda -Minfo -lcublas
sgemm_cublas:
     48, Loop not vectorized/parallelized: contains call
     68, Memory set idiom, array assignment replaced by call to pgf90_mset4
     69, Memory set idiom, array assignment replaced by call to pgf90_mset4
    110, Loop not vectorized/parallelized: contains call

実行

 ここでは、実行モジュール名を Linux 流の a.out と表示しています。Windows では、***.exe となります。このプログラム例では、実行モジュールの後に引数を指定できます。-device 0 と言った形で、GPU ボードの論理番号を指定します。私のマシンの場合は、二つのボードが搭載されていますので、以下のような 0, 1 の指定で、各 GPU の性能が測定できます。以下の例は、device 0(GTX580) を使用して、N=4096 の場合と N=6000 の場合の性能を測定しています。N=6000 で 800GLFLOPS となり、800GFLOPSオーダーの性能を確保できるようです。

[kato@photon29]$ ./a.out
 Select BLAS interface? 0=traditional, 1=Portable legacy CUBLAS <4.0, 2=CUBLAS Version2
2
 BLAS interface =            2
 What size of matrix : ?
4096
 MatSize =          4096  x          4096  x          4096

Device:GeForce GTX 580, 1544.0 MHz clock, 1535.2 MB memory.

 system_clock resolution:    1.0000000E-06
 *** Verify results ***
 Total Time (sec):    0.1919399946928024
 Total SGEMM GFLOPS:     716.0516686059585 (716GFLOPS)
 
[kato@photon29 11.8]$ ./a.out
 Select BLAS interface? 0=traditional, 1=Portable legacy CUBLAS <4.0, 2=CUBLAS Version2
0
 BLAS interface =            0
 What size of matrix : ?
6000
 MatSize =          6000  x          6000  x          6000

Device:GeForce GTX 580, 1544.0 MHz clock, 1535.2 MB memory.

 system_clock resolution:    1.0000000E-06
 *** Verify results ***
 Total Time (sec):    0.5399159789085388
 Total SGEMM GFLOPS:     800.1245102345458 (800GFLOPS)