Fortranプログラムをgfortranを使ってLLVM経由でCUDA化する

出典: トータル・ディスクロージャ・サイト(事実をありのままに)

4.1以降のCUDAは、コンパイラがLLVMベースのものに変更されている。具体的には、nvccはCUDA CのソースコードをまずLLVM IRと呼ばれるLLVM仮想マシンが解釈する中間コードにコンパイルし、LLVM仮想マシンが改めてPTXアセンブリにコンパイルする。かつてはこの中間コードはLLVMからのforkとして、CUDA Toolkit付属のnvvmプログラムが処理するNVVM IRとして厳密には別物扱いであったが、現在はCUDA向けの変更はLLVM本体のプロジェクトにマージされ、最新のLLVM 3.2では、marchnvptx (32-bit)/nvptx64 (64-bit)、mcpuに各種Compute Capabilityを指定することでコンパイル結果がPTXファイルとなる。

この機能の最大のメリットは、バックエンドとしてLLVMを持つコンパイラ言語であれば、LLVMを使ってCUDAのGPGPUコードを出力できるということである。今回は、GCCのバックエンドにLLVMを使うプラグインであるDragonEggを使い、GCCのFortranコンパイラであるgfortranを使ってFortranコードをCUDA化してみる検証を行ってみた。

目次

準備

まずは以下のコンパイラ環境を揃える。

これらの環境を兼ね備えたOSは少ないと思われるので、適宜ソースコードからビルドしていただくことになると思われる。今回の実験ではDebian WheezyのGCC 4.6を使い、LLVMとDragonEggはこの環境にインストールしてPTX生成までを行い、最終的にCentOS 5.9上のCUDA 5.0で動作確認を行った。

Fortranプログラムを書く

テスト用に、次のプログラムを書いた。

      SUBROUTINE MATMULKERNEL(N,M,K,A,B,C)
      IMPLICIT NONE
      INTEGER BLOCK_SIZE
      PARAMETER(BLOCK_SIZE=16)
      EXTERNAL THREADIDXX,THREADIDXY,BLOCKIDXX,BLOCKIDXY
      INTEGER THREADIDXX,THREADIDXY,BLOCKIDXX,BLOCKIDXY
      INTEGER N,M,K
      REAL A(N,K),B(K,M),C(N,M)
      INTEGER BLOCKROW,BLOCKCOL,ROW,COL,IB,I
      REAL CI
      COMMON AS,BS
      REAL AS(BLOCK_SIZE,BLOCK_SIZE),BS(BLOCK_SIZE,BLOCK_SIZE)

      BLOCKROW=BLOCKIDXX()*BLOCK_SIZE
      BLOCKCOL=BLOCKIDXY()*BLOCK_SIZE
      ROW=THREADIDXX()+1
      COL=THREADIDXY()+1
      CI=0.0
      DO 20 IB=0,K-1,BLOCK_SIZE
         AS(ROW,COL)=A(BLOCKROW+ROW,IB+COL)
         BS(ROW,COL)=B(IB+ROW,BLOCKCOL+COL)
         CALL SYNCTHREADS()
         DO 10 I=1,BLOCK_SIZE
            CI=CI+AS(ROW,I)*BS(I,COL)
 10      CONTINUE
         CALL SYNCTHREADS()
 20   CONTINUE
      C(BLOCKROW+ROW,BLOCKCOL+COL)=CI
      RETURN
      END

CUDA C Programming Guideの演習の行列積プログラムをそのままFortranに置き換えたようなプログラムであるが、以下のポイントが後の工程で大事となる。

  • Thread IDやBlock IDは、INTEGER型の値を返す関数(ここではTHREADIDXX()やBLOCKIDXY()など)から得る、というスタイルとなっている。
  • CUDA Cにおける__syncthreads()は、適当なサブルーチン(ここではSYNCTHREADS)の呼び出し、というスタイルとなっている。
  • Shared MemoryはCOMMON領域に確保する。

なお、THREADIDXXやSYNCTHREADSといった関数/サブルーチン名は、後の工程でLLVM IRファイルを書き換える際にLLVM IRのイントリンシックに置き換えるので、他の識別子と区別が付けばどういう名前でも構わない。

DragonEggがDRAGONEGG環境変数の示すディレクトリにインストールされ、このファイルをmatmulkernel.fというファイル名で保存すると仮定し、次のコマンドでgfortranを使ってLLVM IRのファイルとする。

$ gfortran -O3 matmulkernel.f -fplugin=$DRAGONEGG/dragonegg.so -fplugin-arg-dragonegg-emit-ir -S

これで、元のファイル名の.f.sとしたファイル名でLLVM IRファイルが保存される。

LLVMファイルの修正

LLVM IRファイルに行う修正は、以下の3点である。

  • 各種関数/サブルーチン呼び出しのイントリンシックへの置き換え。
  • matmulkernel_関数がCUDAにおける__global__空間に配置されるための指示。
  • COMMON領域がShared Memoryとなるための指示。

まずは、前工程で出力されたLLVM IRファイルの確認である。

; ModuleID = 'matmulkernel.f'
target datalayout = "e-p:64:64:64-S128-i1:8:8-i8:8:8-i16:16:16-i32:32:32-i64:64:64-f16:16:16-f32:32:32-f64:64:64-f128:128:128-v64:64:64-v128:128:128-a0:0:64-s0:64:64-f80:128:128-n8:16:32:64"
target triple = "x86_64--linux-gnu"

module asm "\09.ident\09\22GCC: (Debian 4.6.3-14) 4.6.3 LLVM: 3.2svn\22"

%0 = type { [256 x float], [256 x float] }

@__BLNK__ = common unnamed_addr global %0 zeroinitializer, align 32

define void @matmulkernel_(i32* noalias nocapture %n, i32* noalias nocapture %m, i32* noalias nocapture %k, [0 x float]* noalias nocapture %a, [0 x float]* noalias nocapture %b, [0 x float]* noalias nocapture %c) unnamed_addr nounwind uwtable {
entry:
  %0 = load i32* %n, align 4, !tbaa !0
  %1 = sext i32 %0 to i64
  %2 = icmp slt i64 %1, 0
  %3 = select i1 %2, i64 0, i64 %1
  %4 = load i32* %k, align 4, !tbaa !0
  %5 = sext i32 %4 to i64
  %not = xor i64 %3, -1
  %6 = icmp slt i64 %5, 0
  %7 = select i1 %6, i64 0, i64 %5
  %8 = tail call i32 bitcast (i32 (...)* @blockidxx_ to i32 ()*)() nounwind
  %9 = shl nsw i32 %8, 4
  %10 = tail call i32 bitcast (i32 (...)* @blockidxy_ to i32 ()*)() nounwind
  %11 = shl nsw i32 %10, 4
  %12 = tail call i32 bitcast (i32 (...)* @threadidxx_ to i32 ()*)() nounwind
  %13 = add nsw i32 %12, 1
  %14 = tail call i32 bitcast (i32 (...)* @threadidxy_ to i32 ()*)() nounwind
  %15 = add nsw i32 %14, 1
  %16 = add nsw i32 %4, -1
  %17 = icmp slt i32 %16, 0
  br i1 %17, label %entry.8_crit_edge, label %"3"

entry.8_crit_edge:                                ; preds = %entry
  %.pre = add nsw i32 %13, %9
  %.pre5 = sext i32 %.pre to i64
  %.pre7 = add nsw i32 %15, %11
  %.pre9 = sext i32 %.pre7 to i64
  %.pre11 = add i64 %.pre5, %not
  br label %"8"

"3":                                              ; preds = %entry
  %not1 = xor i64 %7, -1
  %18 = lshr i32 %16, 4
  %19 = sext i32 %13 to i64
  %20 = sext i32 %15 to i64
  %21 = shl nsw i64 %20, 4
  %22 = add i64 %19, -17
  %23 = add i64 %22, %21
  %24 = add nsw i32 %13, %9
  %25 = sext i32 %24 to i64
  %26 = add i64 %25, %not
  %27 = getelementptr inbounds %0* @__BLNK__, i64 0, i32 0, i64 %23
  %28 = add nsw i32 %15, %11
  %29 = sext i32 %28 to i64
  %30 = mul nsw i64 %29, %7
  %31 = add i64 %30, %not1
  %32 = getelementptr inbounds %0* @__BLNK__, i64 0, i32 1, i64 %23
  %33 = add i64 %21, -17
  br label %"4"

"4":                                              ; preds = %"7", %"3"
  %indvars.iv3 = phi i64 [ %indvars.iv.next4, %"7" ], [ 0, %"3" ]
  %34 = phi float [ %54, %"7" ], [ 0.000000e+00, %"3" ]
  %35 = phi i32 [ %56, %"7" ], [ %18, %"3" ]
  %36 = add nsw i64 %indvars.iv3, %20
  %37 = mul nsw i64 %36, %3
  %38 = add i64 %26, %37
  %39 = getelementptr inbounds [0 x float]* %a, i64 0, i64 %38
  %40 = load float* %39, align 4, !tbaa !2
  store float %40, float* %27, align 4, !tbaa !2
  %41 = add nsw i64 %indvars.iv3, %19
  %42 = add i64 %31, %41
  %43 = getelementptr inbounds [0 x float]* %b, i64 0, i64 %42
  %44 = load float* %43, align 4, !tbaa !2
  store float %44, float* %32, align 4, !tbaa !2
  tail call void bitcast (void (...)* @syncthreads_ to void ()*)() nounwind
  br label %"5"

"5":                                              ; preds = %"4", %"5"
  %indvars.iv = phi i64 [ 1, %"4" ], [ %indvars.iv.next, %"5" ]
  %45 = phi float [ %34, %"4" ], [ %54, %"5" ]
  %46 = shl nsw i64 %indvars.iv, 4
  %47 = add i64 %22, %46
  %48 = getelementptr inbounds %0* @__BLNK__, i64 0, i32 0, i64 %47
  %49 = load float* %48, align 4, !tbaa !2
  %50 = add i64 %33, %indvars.iv
  %51 = getelementptr inbounds %0* @__BLNK__, i64 0, i32 1, i64 %50
  %52 = load float* %51, align 4, !tbaa !2
  %53 = fmul float %49, %52
  %54 = fadd float %45, %53
  %indvars.iv.next = add i64 %indvars.iv, 1
  %lftr.wideiv = trunc i64 %indvars.iv.next to i32
  %exitcond = icmp eq i32 %lftr.wideiv, 17
  br i1 %exitcond, label %"6", label %"5"

"6":                                              ; preds = %"5"
  tail call void bitcast (void (...)* @syncthreads_ to void ()*)() nounwind
  %55 = icmp eq i32 %35, 0
  br i1 %55, label %"8", label %"7"

"7":                                              ; preds = %"6"
  %indvars.iv.next4 = add i64 %indvars.iv3, 16
  %56 = add i32 %35, -1
  br label %"4"

"8":                                              ; preds = %"6", %entry.8_crit_edge
  %.pre-phi12 = phi i64 [ %.pre11, %entry.8_crit_edge ], [ %26, %"6" ]
  %.pre-phi10 = phi i64 [ %.pre9, %entry.8_crit_edge ], [ %29, %"6" ]
  %57 = phi float [ 0.000000e+00, %entry.8_crit_edge ], [ %54, %"6" ]
  %58 = mul nsw i64 %.pre-phi10, %3
  %59 = add i64 %.pre-phi12, %58
  %60 = getelementptr inbounds [0 x float]* %c, i64 0, i64 %59
  store float %57, float* %60, align 4, !tbaa !2
  ret void
}

declare i32 @blockidxx_(...)

declare i32 @blockidxy_(...)

declare i32 @threadidxx_(...)

declare i32 @threadidxy_(...)

declare void @syncthreads_(...)

!0 = metadata !{metadata !"alias set 7: integer(kind=4)", metadata !1}
!1 = metadata !{metadata !1}
!2 = metadata !{metadata !"alias set 15: real(kind=4)", metadata !1}

イントリンシックへの置き換え

threadIdx.x/y/zcall i32 @llvm.ptx.read.tid.x/y/z()blockDim.x/y/zcall i32 @llvm.ptx.read.ntid.x/y/z()blockIdx.x/y/zcall i32 @llvm.ptx.read.ctaid.x/y/z()gridDim.x/y/zcall i32 @llvm.ptx.read.nctaid.x/y/z()のイントリンシックで得ることができる。今回の出力で、具体的にそれらが必要な箇所は、22行目の

  %8 = tail call i32 bitcast (i32 (...)* @blockidxx_ to i32 ()*)() nounwind

24行目の

  %10 = tail call i32 bitcast (i32 (...)* @blockidxy_ to i32 ()*)() nounwind

26行目の

  %12 = tail call i32 bitcast (i32 (...)* @threadidxx_ to i32 ()*)() nounwind

28行目の

  %14 = tail call i32 bitcast (i32 (...)* @threadidxy_ to i32 ()*)() nounwind

である。これらをそれぞれ

  %8 = call i32 @llvm.ptx.read.ctaid.x()
  %10 = call i32 @llvm.ptx.read.ctaid.y()
  %12 = call i32 @llvm.ptx.read.tid.x()
  %14 = call i32 @llvm.ptx.read.tid.y()

に置き換える。

また、__syncthreads()

  call void @llvm.ptx.bar.sync(i32 0)

のイントリンシックで実行されるので、77行目と98行目の

  tail call void bitcast (void (...)* @syncthreads_ to void ()*)() nounwind

を上記の文に置き換えれば良い。

関数の__global__

11行目の

define void @matmulkernel_(i32* noalias nocapture %n, i32* noalias nocapture %m, i32* noalias nocapture %k, [0 x float]* noalias nocapture %a, [0 x float]* noalias nocapture %b, [0 x float]* noalias nocapture %c) unnamed_addr nounwind uwtable {

に対して、defineの直後にptx_kernelの一語を加えれば良い。

define ptx_kernel void @matmulkernel_(i32* noalias nocapture %n, i32* noalias nocapture %m, i32* noalias nocapture %k, [0 x float]* noalias nocapture %a, [0 x float]* noalias nocapture %b, [0 x float]* noalias nocapture %c) unnamed_addr nounwind uwtable {

という形になれば良い。

COMMON領域のShared Memory化

COMMON領域のデータは、9行目の

@__BLNK__ = common unnamed_addr global %0 zeroinitializer, align 32

の部分で宣言されている。これをShared Memoryとするには、common unnamed_addrの部分をaddrspace(3)に置き換えれば良い。結果この行は

@__BLNK__ = addrspace(3) global %0 zeroinitializer, align 32

という宣言になるのであるが、実際はプログラム中でこのアドレス空間を参照する部分で、float*の記述をfloat addrspace(3)*に置き換える作業が残っている。ここまでの変更を加えたらとにかくコンパイルしてみて、エラー箇所を次々に直していくのが、遠回りに見えて一番早い方法だろう。

今回は53、58、71、76、85、86、88、89行目に変更が加わることになる。

今回は最終的に、

; ModuleID = 'matmulkernel.f'
target datalayout = "e-p:64:64:64-S128-i1:8:8-i8:8:8-i16:16:16-i32:32:32-i64:64:64-f16:16:16-f32:32:32-f64:64:64-f128:128:128-v64:64:64-v128:128:128-a0:0:64-s0:64:64-f80:128:128-n8:16:32:64"
target triple = "x86_64--linux-gnu"

module asm "\09.ident\09\22GCC: (Debian 4.6.3-14) 4.6.3 LLVM: 3.2svn\22"

%0 = type { [256 x float], [256 x float] }

@__BLNK__ = addrspace(3) global %0 zeroinitializer, align 32

define ptx_kernel void @matmulkernel_(i32* noalias nocapture %n, i32* noalias nocapture %m, i32* noalias nocapture %k, [0 x float]* noalias nocapture %a, [0 x float]* noalias nocapture %b, [0 x float]* noalias nocapture %c) unnamed_addr nounwind uwtable {
entry:
  %0 = load i32* %n, align 4, !tbaa !0
  %1 = sext i32 %0 to i64
  %2 = icmp slt i64 %1, 0
  %3 = select i1 %2, i64 0, i64 %1
  %4 = load i32* %k, align 4, !tbaa !0
  %5 = sext i32 %4 to i64
  %not = xor i64 %3, -1
  %6 = icmp slt i64 %5, 0
  %7 = select i1 %6, i64 0, i64 %5
  %8 = call i32 @llvm.ptx.read.ctaid.x()
  %9 = shl nsw i32 %8, 4
  %10 = call i32 @llvm.ptx.read.ctaid.y()
  %11 = shl nsw i32 %10, 4
  %12 = call i32 @llvm.ptx.read.tid.x()
  %13 = add nsw i32 %12, 1
  %14 = call i32 @llvm.ptx.read.tid.y()
  %15 = add nsw i32 %14, 1
  %16 = add nsw i32 %4, -1
  %17 = icmp slt i32 %16, 0
  br i1 %17, label %entry.8_crit_edge, label %"3"

entry.8_crit_edge:                                ; preds = %entry
  %.pre = add nsw i32 %13, %9
  %.pre5 = sext i32 %.pre to i64
  %.pre7 = add nsw i32 %15, %11
  %.pre9 = sext i32 %.pre7 to i64
  %.pre11 = add i64 %.pre5, %not
  br label %"8"

"3":                                              ; preds = %entry
  %not1 = xor i64 %7, -1
  %18 = lshr i32 %16, 4
  %19 = sext i32 %13 to i64
  %20 = sext i32 %15 to i64
  %21 = shl nsw i64 %20, 4
  %22 = add i64 %19, -17
  %23 = add i64 %22, %21
  %24 = add nsw i32 %13, %9
  %25 = sext i32 %24 to i64
  %26 = add i64 %25, %not
  %27 = getelementptr inbounds %0 addrspace(3)* @__BLNK__, i64 0, i32 0, i64 %23
  %28 = add nsw i32 %15, %11
  %29 = sext i32 %28 to i64
  %30 = mul nsw i64 %29, %7
  %31 = add i64 %30, %not1
  %32 = getelementptr inbounds %0 addrspace(3)* @__BLNK__, i64 0, i32 1, i64 %23
  %33 = add i64 %21, -17
  br label %"4"

"4":                                              ; preds = %"7", %"3"
  %indvars.iv3 = phi i64 [ %indvars.iv.next4, %"7" ], [ 0, %"3" ]
  %34 = phi float [ %54, %"7" ], [ 0.000000e+00, %"3" ]
  %35 = phi i32 [ %56, %"7" ], [ %18, %"3" ]
  %36 = add nsw i64 %indvars.iv3, %20
  %37 = mul nsw i64 %36, %3
  %38 = add i64 %26, %37
  %39 = getelementptr inbounds [0 x float]* %a, i64 0, i64 %38
  %40 = load float* %39, align 4, !tbaa !2
  store float %40, float addrspace(3)* %27, align 4, !tbaa !2
  %41 = add nsw i64 %indvars.iv3, %19
  %42 = add i64 %31, %41
  %43 = getelementptr inbounds [0 x float]* %b, i64 0, i64 %42
  %44 = load float* %43, align 4, !tbaa !2
  store float %44, float addrspace(3)* %32, align 4, !tbaa !2
  call void @llvm.ptx.bar.sync(i32 0)
  br label %"5"

"5":                                              ; preds = %"4", %"5"
  %indvars.iv = phi i64 [ 1, %"4" ], [ %indvars.iv.next, %"5" ]
  %45 = phi float [ %34, %"4" ], [ %54, %"5" ]
  %46 = shl nsw i64 %indvars.iv, 4
  %47 = add i64 %22, %46
  %48 = getelementptr inbounds %0 addrspace(3)* @__BLNK__, i64 0, i32 0, i64 %47
  %49 = load float addrspace(3)* %48, align 4, !tbaa !2
  %50 = add i64 %33, %indvars.iv
  %51 = getelementptr inbounds %0 addrspace(3)* @__BLNK__, i64 0, i32 1, i64 %50
  %52 = load float addrspace(3)* %51, align 4, !tbaa !2
  %53 = fmul float %49, %52
  %54 = fadd float %45, %53
  %indvars.iv.next = add i64 %indvars.iv, 1
  %lftr.wideiv = trunc i64 %indvars.iv.next to i32
  %exitcond = icmp eq i32 %lftr.wideiv, 17
  br i1 %exitcond, label %"6", label %"5"

"6":                                              ; preds = %"5"
  call void @llvm.ptx.bar.sync(i32 0)
  %55 = icmp eq i32 %35, 0
  br i1 %55, label %"8", label %"7"

"7":                                              ; preds = %"6"
  %indvars.iv.next4 = add i64 %indvars.iv3, 16
  %56 = add i32 %35, -1
  br label %"4"

"8":                                              ; preds = %"6", %entry.8_crit_edge
  %.pre-phi12 = phi i64 [ %.pre11, %entry.8_crit_edge ], [ %26, %"6" ]
  %.pre-phi10 = phi i64 [ %.pre9, %entry.8_crit_edge ], [ %29, %"6" ]
  %57 = phi float [ 0.000000e+00, %entry.8_crit_edge ], [ %54, %"6" ]
  %58 = mul nsw i64 %.pre-phi10, %3
  %59 = add i64 %.pre-phi12, %58
  %60 = getelementptr inbounds [0 x float]* %c, i64 0, i64 %59
  store float %57, float* %60, align 4, !tbaa !2
  ret void
}

declare i32 @llvm.ptx.read.tid.x()
declare i32 @llvm.ptx.read.tid.y()
declare i32 @llvm.ptx.read.ctaid.x()
declare i32 @llvm.ptx.read.ctaid.y()
declare void @llvm.ptx.bar.sync(i32 %i)

!0 = metadata !{metadata !"alias set 7: integer(kind=4)", metadata !1}
!1 = metadata !{metadata !1}
!2 = metadata !{metadata !"alias set 15: real(kind=4)", metadata !1}

というLLVM IRファイルとなり、これをLLVMコンパイラを使ってコンパイルすることとなる。

LLVMコンパイラによるコンパイル

今回は64ビットのCompute Capability 2.0向けにコンパイルするので、以下のコマンドでコンパイルする。

$ llc -O3 -march=nvptx64 -mcpu=sm_20 < matmulkernel.s > matmulkernel.ptx

これで、今回は以下のようなPTXファイルが得られた。

//
// Generated by LLVM NVPTX Back-End
//

.version 3.1
.target sm_20, texmode_independent
.address_size 64

.visible .shared .align 32 .b8 __BLNK__[2048];

	// .globl	matmulkernel_
.entry matmulkernel_(
	.param .u64 .ptr .align 4 matmulkernel__param_0,
	.param .u64 .ptr .align 4 matmulkernel__param_1,
	.param .u64 .ptr .align 4 matmulkernel__param_2,
	.param .u64 .ptr .align 4 matmulkernel__param_3,
	.param .u64 .ptr .align 4 matmulkernel__param_4,
	.param .u64 .ptr .align 4 matmulkernel__param_5
)                                       // @matmulkernel_
{
	.reg .pred %p<396>;
	.reg .s16 %rc<396>;
	.reg .s16 %rs<396>;
	.reg .s32 %r<396>;
	.reg .s64 %rl<396>;
	.reg .f32 %f<396>;
	.reg .f64 %fl<396>;

// BB#0:                                // %entry
	ld.param.u64 	%rl0, [matmulkernel__param_0];
	ld.s32 	%rl0, [%rl0];
	setp.lt.s64 	%p0, %rl0, 0;
	selp.b64 	%rl0, 0, %rl0, %p0;
	ld.param.u64 	%rl1, [matmulkernel__param_2];
	ld.u32 	%r0, [%rl1];
	not.b64 	%rl2, %rl0;
	mov.u32	%r1, %ctaid.x;
	shl.b32 	%r4, %r1, 4;
	mov.u32	%r1, %ctaid.y;
	shl.b32 	%r2, %r1, 4;
	mov.u32	%r1, %tid.x;
	add.s32 	%r3, %r1, 1;
	mov.u32	%r1, %tid.y;
	add.s32 	%r1, %r1, 1;
	add.s32 	%r5, %r0, -1;
	setp.gt.s32 	%p0, %r5, -1;
	ld.param.u64 	%rl1, [matmulkernel__param_5];
	@%p0 bra 	BB0_3;
	bra.uni 	BB0_1;
BB0_3:                                  // %3
	cvt.s64.s32 	%rl3, %r0;
	setp.lt.s64 	%p0, %rl3, 0;
	selp.b64 	%rl8, 0, %rl3, %p0;
	ld.param.u64 	%rl4, [matmulkernel__param_4];
	ld.param.u64 	%rl5, [matmulkernel__param_3];
	not.b64 	%rl9, %rl8;
	shr.u32 	%r0, %r5, 4;
	cvt.s64.s32 	%rl6, %r3;
	cvt.s64.s32 	%rl7, %r1;
	mul.wide.s32 	%rl11, %r1, 16;
	add.s32 	%r4, %r3, %r4;
	cvt.s64.s32 	%rl3, %r4;
	add.s64 	%rl2, %rl3, %rl2;
	add.s32 	%r2, %r1, %r2;
	cvt.s64.s32 	%rl3, %r2;
	mad.lo.s64 	%rl8, %rl3, %rl8, %rl9;
	mul.wide.s32 	%rl9, %r3, 4;
	mov.u64 	%rl12, __BLNK__;
	add.s64 	%rl9, %rl9, %rl12;
	add.s64 	%rl9, %rl9, -4;
	mul.wide.s32 	%rl10, %r1, 64;
	add.s64 	%rl10, %rl10, %rl12;
	add.s64 	%rl10, %rl10, 960;
	add.s64 	%rl11, %rl6, %rl11;
	shl.b64 	%rl11, %rl11, 2;
	add.s64 	%rl12, %rl11, %rl12;
	add.s64 	%rl11, %rl12, 956;
	add.s64 	%rl12, %rl12, -68;
	mov.f32 	%f0, 0f00000000;
	mov.u64 	%rl13, 0;
	bra.uni 	BB0_4;
BB0_7:                                  // %7
                                        //   in Loop: Header=BB0_4 Depth=1
	add.s64 	%rl13, %rl13, 16;
	add.s32 	%r0, %r0, -1;
BB0_4:                                  // %4
                                        // =>This Loop Header: Depth=1
                                        //     Child Loop BB0_5 Depth 2
	add.s64 	%rl14, %rl13, %rl7;
	mad.lo.s64 	%rl14, %rl14, %rl0, %rl2;
	shl.b64 	%rl14, %rl14, 2;
	add.s64 	%rl14, %rl5, %rl14;
	add.s64 	%rl15, %rl13, %rl6;
	add.s64 	%rl15, %rl8, %rl15;
	shl.b64 	%rl15, %rl15, 2;
	add.s64 	%rl15, %rl4, %rl15;
	ld.f32 	%f1, [%rl14];
	st.shared.f32 	[%rl12], %f1;
	ld.f32 	%f1, [%rl15];
	st.shared.f32 	[%rl11], %f1;
	bar.sync	0;
	mov.u32 	%r1, 16;
	mov.u64 	%rl14, %rl9;
	mov.u64 	%rl15, %rl10;
BB0_5:                                  // %5
                                        //   Parent Loop BB0_4 Depth=1
                                        // =>  This Inner Loop Header: Depth=2
	ld.shared.f32 	%f1, [%rl15];
	ld.shared.f32 	%f2, [%rl14];
	fma.rn.f32 	%f0, %f2, %f1, %f0;
	add.s32 	%r1, %r1, -1;
	add.s64 	%rl15, %rl15, 4;
	add.s64 	%rl14, %rl14, 64;
	setp.ne.s32 	%p0, %r1, 0;
	@%p0 bra 	BB0_5;
// BB#6:                                // %6
                                        //   in Loop: Header=BB0_4 Depth=1
	bar.sync	0;
	setp.eq.s32 	%p0, %r0, 0;
	@%p0 bra 	BB0_2;
	bra.uni 	BB0_7;
BB0_1:                                  // %entry.8_crit_edge
	add.s32 	%r0, %r3, %r4;
	cvt.s64.s32 	%rl4, %r0;
	add.s32 	%r0, %r1, %r2;
	cvt.s64.s32 	%rl3, %r0;
	add.s64 	%rl2, %rl4, %rl2;
	mov.f32 	%f0, 0f00000000;
BB0_2:                                  // %8
	mad.lo.s64 	%rl0, %rl3, %rl0, %rl2;
	shl.b64 	%rl0, %rl0, 2;
	add.s64 	%rl0, %rl1, %rl0;
	st.f32 	[%rl0], %f0;
	ret;
}

最後に、このCUDAコードをテストする。

テスト

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <sys/time.h>

#include <cuda.h>

#define M 512
#define N 768
#define K 1024

extern "C" void
sgemm_(void *, void *, void *, void *, void *,
       void *, void *, void *, void *, void *, void *, void *, void *);

static float a[K * M], b[N * K], c[N * M], c_gold[N * M];

int
main()
{
  std::srand(0);
  for (ptrdiff_t i = 0; i < K * M; ++i) {
    a[i] = (float) std::rand() / (float) RAND_MAX;
  }
  for (ptrdiff_t i = 0; i < N * K; ++i) {
    b[i] = (float) std::rand() / (float) RAND_MAX;
  }
  int m = M;
  int n = N;
  int k = K;
  cuInit(0);
  CUdevice dev;
  cuDeviceGet(&dev, 0);
  CUcontext ctx;
  cuCtxCreate(&ctx, 0, dev);
  CUdeviceptr d_m, d_n, d_k, d_a, d_b, d_c;
  cuMemAlloc(&d_m, sizeof(int));
  cuMemAlloc(&d_n, sizeof(int));
  cuMemAlloc(&d_k, sizeof(int));
  cuMemAlloc(&d_a, K * M * sizeof(float));
  cuMemAlloc(&d_b, N * K * sizeof(float));
  cuMemAlloc(&d_c, N * M * sizeof(float));
  cuMemcpyHtoD(d_m, &m, sizeof(int));
  cuMemcpyHtoD(d_n, &n, sizeof(int));
  cuMemcpyHtoD(d_k, &k, sizeof(int));
  cuMemcpyHtoD(d_a, a, K * M * sizeof(float));
  cuMemcpyHtoD(d_b, b, N * K * sizeof(float));
  CUmodule module;
  cuModuleLoad(&module, "matmulkernel.ptx");
  CUfunction func;
  cuModuleGetFunction(&func, module, "matmulkernel_");
  void *kernelParams[] = {
    &d_m,
    &d_n,
    &d_k,
    &d_a,
    &d_b,
    &d_c,
  };
  timeval start;
  gettimeofday(&start, 0);
  cuLaunchKernel(func, M / 16, N / 16, 1, 16, 16, 1, 0, 0, kernelParams, 0);
  cuCtxSynchronize();
  timeval end;
  gettimeofday(&end, 0);
  double time =
    end.tv_sec - start.tv_sec + (end.tv_usec - start.tv_usec) * 1.0e-6;
  cuMemcpyDtoH(c, d_c, N * M * sizeof(float));
  char trans = 'N';
  float alpha = 1.0f;
  float beta = 0.0f;
  sgemm_(&trans, &trans, &m, &n, &k,
         &alpha, a, &m, b, &k, &beta, c_gold, &m);
  double error = 0.0;
  for (ptrdiff_t i = 0; i < N * M; ++i) {
    error +=
      std::fabs((double) c[i] - (double) c_gold[i]) /
      ((double) c[i] + (double) c_gold[i]);
  }
  error /= N * M;
  std::printf("LLVM CUDA Fortran Test\n");
  double flops = 2.0 * M * N * K / time;
  std::printf("%.6e seconds\t%.6e FLOPS\n", time ,flops);
  std::printf("Error: %e\n", error);
  return 0;
}

Tesla C2050を使ってテストした結果は、以下の通りであった。

LLVM CUDA Fortran Test
5.113000e-03 seconds	1.575017e+11 FLOPS
Error: 1.773923e-07

誤差は問題無さそうであるし、157 GFLOPSもこのプログラムにしては満足な数値。少々の手間を厭わなければ、実用に問題は無いと思われる。


この記事へのコメントをお寄せください

  • サイトへの書き込みに差し支えございましたら トータルディスクロージャーサイトサポート係へメールをお送りください
  • トータル・ディスクロージャ・サイトに投稿された文章と画像は、すべてその著作権がHPCシステムズ株式会社に帰属し、HPCシステムズ株式会社が著作権を所有することに同意してください。
  • あなたの文章が他人によって自由に編集、配布されることを望まない場合は、投稿を控えてください。
  • コメントを書き込む場合は名前にひらがなを織り交ぜてください。
  • あなたの投稿する文章と画像はあなた自身によって書かれたものであるか、パブリック・ドメインかそれに類する自由なリソースからの複製であることを約束してください。あなたが著作権を保持していない作品を許諾なしに投稿してはいけません!

<comments hideform="false" />


Comments

ノート:Fortranプログラムをgfortranを使ってLLVM経由でCUDA化する

個人用ツール