CUDA Fortran入門

概要

CUDA Fortranの書き方・使い方を簡単に紹介します．

使用環境

• PGI Visual Fortran for Windows version 18.7
• CUDA 9.2 (Driver version 10.0)
• NVIDIA GTX Titan¹

CUDAとは

Compute Unified Device Architectureの略で，NVIDIA GPUで動作するプログラムの開発環境です．C/C++言語に独自の拡張が加えられています．CUDA Fortranは，簡単に言うとCUDAのFortran版です．CUDA Fortranと明確に区別するために，CUDAがCUDA Cと表現される場合もあります．CUDA FortranはPGI社のFortranコンパイラをインストールすると利用できるようになります．CUDA FortranはCUDA Cを利用しますが，新機能はPGIコンパイラが対応しないと利用できません．NVIDIA社がPGI社を買収したので，より使いやすくなってくれることを期待します．

CUDA Fortranの特徴は，なんといってもプログラムの記述がCUDA Cと比較してかなり簡単だということでしょう．プログラムの移植の際に，変数の宣言や動的割付け，GPUへのデータのコピーなど，ワクワクしない決まり切った作業を簡単に片付けることができます．

GPUとは

謎のAI半導体です．この謎の半導体を使うとAIが加速します．

GPUは画像処理装置で，主にCGの描画とそれに関係した計算を行うハードウェアです．画像処理を効率よく行うために並列処理に特化しており，2000年代からその並列処理能力を画像処理以外に用いる試みがさかんに行われてきました．GPUを画像処理以外の処理に用いることをGPGPU（General Purpose computing on GPU）とよびます．

GPUに仕事をさせるにはOpenGL等の画像処理用のライブラリしか手段がなかったので，一般的な処理を画像処理に置き換える必要がありました．CUDA登場後は，ハードウェアの特性を理解する必要はありますが，画像処理を意識することなくプログラムが書けるようになりました．

GPUのハードウェア構造

GPUの構造はCPUと大きく違います．GPUには，単純な処理が得意なCUDAコアが大量に搭載されています．高価なGPUで5000個以上，安価なGPUでも100個程度あります．多くても数十個のコアしか搭載されていないCPUと比べると桁違いに多い数ですが，CUDAコアは，CPUのコアのように複雑な処理はできません．どちらかというと演算器の方が近いでしょうか．

CUDAコアは，いくつかの数がまとめられて，SM（Streaming Multiprocessor）という単位で管理されます．どちらかというと，このSMの方がCPUのコアに近い存在だと思います．SMが集まってGPUのチップが構成されています．

このようにまとめると，CPUとGPUにあまり違いがないように思われます．CPUを使うときにプログラマが意識するのはコアであって，そのコアの中で演算器がどのように使われるかは，あまり意識しないと思います²．スレッド並列の場合には，各スレッドが各コアに割り当てられることを想定します．GPUでは，プログラマが主に意識するのはSMではなくCUDAコアです．非常に多くのスレッドが各CUDAコアに割り当てられ，CUDAコアにどのような処理をさせるかを意識してプログラムを書きます．そして，GPUを効率的に使うためには，すべてのCUDAコアが同じ処理を行うようにアルゴリズムを考えます（データは異なっていても構いません）．異なる処理をさせることもできますが，効率的ではありません．

つまり，ハードウェア構造の観点からはCPUのコアに近いのSMですが，プログラミングの観点では，CPUのコアとCUDAコアが近い存在だと見なせます．

CUDA Fortranを使ってみる

CUDA Fortranには無償のPGI Community Editionがあります．英語版しかありませんが，日本のPGIコンパイラの代理店であるソフテック社のWebページに日本語解説や様々な技術記事があるので，インストールやちょっと使ってみる分にはあまり困らないと思います．ただし，ソフテック社ではPGI Community Editionに関する技術サポートは行っていません．コミュニティの力を頼りましょう．

CUDA Fortranをインストールするには，PGIコンパイラのインストーラをダウンロードして，インストールウィザードに沿ってインストールしていくだけです．このとき，CUDAなど必要な環境を一緒にインストールできます．別途CUDAをインストール済みであれば，ここでのCUDAのインストールを省略できますが，リリースのタイミングによってはCUDA Fortranが，インストール済みのCUDAのバージョンに対応していないこともあります（特に最新版の時は）．

Windows環境では，PGI FortranとVisual Studioとの統合ができたのですが，Visual Studioとの統合はもう提供されなくなりました．残念なことに，Windowsではシンタックスハイライト付きのまともなCUDA Fortran開発環境がありません．Linuxでは新しめのVimがCUDA Fortranのシンタックスハイライトに対応しています．

Hello World

お決まりのHello Worldを実行して，とりあえず動くことを確認しましょう．

helloworld.cuf

program main
    implicit none

    print *,"Hello World"
end program main

CUDA Fortranのソースコードの拡張子は.cufです．コンパイルするには，pgf90あるいはpgfortranコマンドを使い，-Mcudaオプション付きでコンパイルします．

> pgf90 -Mcuda helloworld.cuf
> helloworld.exe
 Hello World

環境によっては，-Mcudaオプションがなくてもほどよく判断してコンパイルしてくれることもあります．
コンパイルが成功して実行ファイルが作られ，それを実行して画面にHello Worldが表示されれば環境構築は成功しています．

Hello Kernel

helloworld.cufにはCPUで処理する内容しかありませんでしたので，GPUで実行する部分を追加していきます．

hellokernel.cuf

module kernel
    implicit none

    contains

    attributes(global) subroutine hello() !!!!
        implicit none

        print *,"Hello Kernel"
    end subroutine hello
end module kernel

program main
    use kernel
    implicit none

    call hello<<<1,1>>>() !!!!
end program main

kernelという名前のモジュールを作り，そこにサブルーチンを追加しました．見慣れたサブルーチンとは異なり，subroutineの前にattributes(global)がついています．これが，GPUで実行するサブルーチンだという印で，attributes(global)がついたサブルーチンをカーネルとよびます．

メインルーチンを見てると，カーネルを呼び出しているところで見慣れない記号<<<1,1>>>がついています．この記号には正式名称はないようで，3重括弧などとよばれています．CUDA Fortran（PGIコンパイラ）はこの記号をCheveron（シェブロン）とよんでいます．中世ヨーロッパの楯に描かれている模様との類似性からきていると思われます．

このシェブロンは，GPUで実行するときの並列実行の度合いを表しています．hellokernel.cufでは1スレッドで実行する設定です．

コンパイルして実行してみると，環境によって，Hello Kernelが表示される場合と，何も表示されずにプログラムが終了する場合があります．

> pgf90 -Mcuda hellokernel.cuf
> hellokernel.exe

これはCUDAのプログラミングモデルに関係しています．CUDAでは，CPUをホスト，GPUをデバイスとよんでおり，ホストがデバイスを制御します．ホストがデバイスにカーネルを実行させたとき，ホストはデバイスでの実行終了を待たずに次の処理に移ります．そのため，hellokernel.cufではGPUでprint *,"Hello Kernel"が実行される前にプログラムが終了することがあります．

GPUで実行したカーネルが終了する前にプログラムが終了してしまっては困るので，同期を取る関数を呼んでカーネルの終了までCPUを待機させるようにします．

program main
    use kernel
    implicit none

    integer :: stat

    call hello<<<1,1>>>()
    stat = cudaDeviceSynchronize() ! ホストとデバイスの同期をとる
end program main

cudaDeviceSynchronize()は一つのGPUの全スレッドの同期をとる関数です．Fortran的な設計では，cudaDeviceSynchronizeをサブルーチンにしてcall cudaDeviceSynchronize(stat)のように呼びそうなものですが，CUDA Cに合わせているので関数になっています．C系統の言語出あれば関数の戻り値を受け取らないことも可能ですが，Fortranでは関数の戻り値は必ず受け取らなければ（変数に代入しなければ）なりません．このような関数の呼び方は，Fortranに慣れているとかなり違和感があるのですが，CUDA Fortranではよく出てきます．

CUDA関連の関数は，拡張子を.cufにしておけばどのモジュールで定義されているかを意識する必要はありません．実際には，拡張子が.f90でもコンパイルはできますが，その場合はuse cudaforが必要です．

hellokernel.f90

module kernel
    implicit none

    contains

    attributes(global) subroutine hello()
        implicit none

        print *,"Hello Kernel"
    end subroutine hello
end module kernel

program main
    use kernel
    use cudafor !!!
    implicit none

    integer :: stat

    call hello<<<1,1>>>()
    stat = cudaDeviceSynchronize()
end program main

CUDA Fortranの基本スタイル

hellokernelは，いわゆるGPUプログラミングの最も基本的な形です．水泳でいうところの蹴伸びのような扱いです．

この節で触れた基本スタイルを言葉で書くと，次のようになるでしょう．

• デバイス（GPU）で実行するカーネルをモジュール内で定義する．
• カーネルにはattributes(global)を付ける．
• カーネルは常にサブルーチンであり，戻り値は利用できない．
• ホスト（CPU）からカーネルを呼ぶときには，カーネル名と引数の間にシェブロンを付ける．
• シェブロンで並列実行の度合いを指定する．
• ホストがカーネルの終了を待つために，同期をとる．

GPUで変数の値を更新してみる

値を2倍するサブルーチンのCPU実装

前の節では，GPUでカーネルを実行してみました．しかし，GPUを使って画面に文字が表示できれば満足だという人は居ても少数でしょう．CPU側で持っているデータをGPUで高速に処理したいと考えているはずです．この節では，ホスト側のデータをどのようにデバイスで処理するか，それをどうやってホストに反映するかという一連のやりとりを確認します．

引数の変数を2倍するサブルーチンdoublify()を定義してみました．特に説明はいらないでしょう．これをGPUに移植してみます．

doublify.f90

module kernel
    implicit none

    contains

    subroutine doublify(a)
        implicit none

        integer,intent(inout) :: a

        a = a*2
    end subroutine doublify
end module kernel

program main
    use kernel
    implicit none

    integer :: a

    a = 1

    call doublify(a)

    print *,a ! 2
end program main

値を2倍するサブルーチンのGPU移植

移植をする前に，GPUのメモリについて言及します．

CPUで処理するデータはメインメモリに置かれます．一方で，PCI-Expressカード型のGPUにもメモリ（GPUメモリとよぶことにします）が実装されており，GPUで処理するデータはそこに置かれます．CPUはGPUメモリを直接読めませんし，GPUもメインメモリを直接読めません．そこでCUDAの出番です．CUDAの機能を利用して，メインメモリとGPUメモリ間のデータ転送を行うことで，メインメモリのデータをGPUで処理できるようになります．CUDAは，メインメモリをホストメモリ，GPUメモリをデバイスメモリとよんでいます．

CUDA Fortranでどのようにホストメモリーデバイスメモリ間のやりとりをするのか，doublify.f90のGPU移植結果doublify.cufで確認します．

doublify.cuf

module kernel
    implicit none

    contains

    attributes(global) subroutine doublify(a) ! attributes(global)を付ける
        implicit none

        integer,intent(inout) :: a

        a = a*2
    end subroutine doublify
end module kernel

program main
    use kernel
    use cudafor
    implicit none

    integer :: a
    integer,device :: dev_a        ! デバイスメモリの割付（デバイス変数の宣言）

    a = 1
    dev_a = a                      ! GPUメモリへ値をコピー
    call doublify<<<1,1>>>(dev_a)  ! カーネル呼び出し
    a = dev_a                      ! GPUメモリから結果をコピー
    print *,a
end program main

変更はコメントで書いた箇所だけです．前節の基本スタイルに沿って，カーネルにattributes(global)を付け，カーネル名と引数の間にシェブロンを付けました．今はまだ1スレッド実行です．

ホストメモリのデータをデバイスメモリで受け取るには，デバイス側で変数（デバイス変数）を用意する必要があります．CUDA Fortranでは，デバイスメモリの宣言は，型宣言の際にdevice属性を付与するだけです．

ホストメモリからデバイスメモリへのデータ転送と，デバイスメモリからホストメモリへのデータ転送は，=で代入するだけで行われます．

前節で言及した基本スタイルには，

• ホストがカーネルの終了を待つために，同期をとる．

とありましたが，上のプログラムにはcudaDeviceSynchronize()が見当たりません．この同期は，=によるデータ転送で代用しています．=によるデータ転送は，転送が終わるまで処理が先に進まない同期通信です．そのためa = dev_aでデバイスメモリからホストメモリへの転送が終わるまで，print *,aは実行されません．

CUDA Fortranの基本スタイル（追加1）

データ転送ができるようになりました．水泳でいうとバタ足ができるようになったくらいでしょうか．

この節で触れた基本スタイルを追加しておきます．

• デバイス変数を宣言するにはdevice属性を付与する．
• ホストのデータをデバイスで処理するときは，データを転送する．
• デバイスで処理した結果をホストで参照するときも転送する．

CUDA Cによる実装

CUDA Fortranの実装とCUDA Cの実装を比べてみましょう．

doublify.cu

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>

__global__ void doublify(int *a) // __global__を付ける
{
    *a = *a * 2;
}

int main()
{
    int a;
    int *dev_a;                                                 // デバイス変数

    a = 1;
    cudaMalloc((void **)&dev_a, sizeof(int));                   // デバイスメモリの割付

    cudaMemcpy(dev_a, &a, sizeof(int), cudaMemcpyHostToDevice); // GPUメモリへ値をコピー
    doublify<<<1,1>>>(dev_a);                                   // カーネル呼び出し
    cudaMemcpy(&a, dev_a, sizeof(int), cudaMemcpyDeviceToHost); // GPUメモリから結果をコピー
    printf("%d\n", a);

    cudaFree(dev_a);
    return 0;
}

CUDA Cでは，デバイス変数を利用するにはまずポインタ変数として宣言し，cudaMalloc()で割り付ける必要があります．割り付けたら解放しなければならないので，プログラムの終了前にcudaFree()を呼びます．これらはCUDA Fortranでは明示的に現れていません．ホストとデバイス間のメモリコピーにはcudaMemcpy()を利用します．

カーネルは，仮引数の名前と型の宣言がまとまっているので，CUDA Cの方がすっきりしているように思います．

CUDA FortranでもcudaMalloc(), cudaFree(), cudaMemcpy()を明示的に呼ぶことは可能ですが，ややこしくなるだけなので，非同期転送をしたいなど特段の理由がなければ必要はないと思います．

module kernel
    implicit none

    contains

    attributes(global) subroutine doublify(a)
        implicit none

        integer,intent(inout) :: a(:)

        a = a*2
    end subroutine doublify
end module kernel

program main
    use kernel
    use cudafor
    implicit none

    integer :: a
    integer,allocatable,device :: dev_a(:)
    integer :: stat

    a = 1
    stat = cudaMalloc(dev_a,1)

    stat = cudaMemcpy(dev_a, a, 1, cudaMemcpyHostToDevice)
    call doublify<<<1,1>>>(dev_a)
    stat = cudaMemcpy(a, dev_a, 1, cudaMemcpyDeviceToHost)
    print *,a
    
    stat = cudaFree(dev_a)
end program main

CUDA Cは基本的に割付・転送するメモリサイズを指定しますが，CUDA Fortranは割付・転送する要素数を指定します．

逐次プログラムを並列化してみる

GPUで処理ができるようになりましたが，GPUの特長である並列処理は全く行っていませんでした．この節では，ようやく並列処理を行います．

並列実行するカーネルの書き方

GPUには数百から数千個のCUDAコアが搭載されており，すべてのCUDAコアに同じ処理をさせます．そうすると，全スレッドに対して処理を書くのはあまり意味がないように思われます．

そこで，CUDAではカーネルには1スレッドが実行する処理を書き，カーネルを呼び出す際に並列実行の度合いを指定します．その度合いは，単純にスレッドの数ではありません．

並列実行の度合いの指定

カーネルの並列実行の度合いは，シェブロンの中に2個の数字を書きます．総スレッド数を書くだけなら数字は一つだけでよいはずです．2個の数字を書くのは，GPUのハードウェアの構造に合わせているからです．

GPUにはCUDAコアという演算器があり，CUDAコアはいくつかの数がまとめられて，SMという単位で管理されていると最初の方で書きました．つまり，階層的に計算資源が管理されているということです．CUDAもスレッドを階層的に管理しています．その対応を表に示します．

ハードウェア構成	並列化階層	CUDA
GPU	全スレッド	グリッド
SM	スレッドの集まり	スレッドブロック
CUDAコア	スレッド	スレッド

シェブロンの中には，スレッドブロック数と，スレッドブロックあたりのスレッド数を書きます．これらカーネル呼出し時の情報はカーネルの中で参照できるので，その情報を基にスレッドを振り分けます．

情報は派生型type(dim3)のビルトイン変数を介して参照します．派生型type(dim3)は，成分としてx,y,zを持っています．

階層	変数名	参照できる情報
グリッド	gridDim	グリッド内にあるスレッドブロックの数
スレッドブロック	blockIdx	各スレッドブロックに割り当てられた番号
	blockDim	スレッドブロック内にあるスレッドの数
スレッド	threadIdx	各スレッドに割り当てられた番号

カーネルhellothreadsの中で，ビルトイン変数（の成分x）gridDim%x, blockIdx%x, blockDim%x, threadIdx%xを表示して，どのような数値が表示されるのかを確認してみましょう．

hellothreads.cuf

module kernel
    implicit none

    contains

    attributes(global) subroutine hellothreads()
        implicit none

        print *,gridDim%x, blockIdx%x, blockDim%x, threadIdx%x
    end subroutine hellothreads
end module kernel

program main
    use kernel
    use cudafor
    implicit none

    integer :: stat

    call hellothreads<<<2,4>>>()
    stat = cudaDeviceSynchronize()
end program main

カーネルの中でビルトイン変数を表示しています．<<<2,4>>>としてカーネルを実行してみると，いくつか数字が表示されます．

2            1            4            1
2            1            4            2
2            1            4            3
2            1            4            4
2            2            4            1
2            2            4            2
2            2            4            3
2            2            4            4

グリッド内に2個のスレッドブロックがあり，そのブロックにはそれぞれ1, 2と番号が付けられています．各スレッドブロックの中には4個のスレッドが存在しており，1,2,3,4と番号が付けられています．すべての番号が重複しているスレッドはありませんので，これらの情報を利用すると，スレッドを識別できます．

ベクトル和のCPU実装

最も単純な処理の例として，ベクトル和を並列化してみます．

vectoradd.f90

module kernel
    implicit none

    contains

    subroutine add(a, b, c, N)
        implicit none

        real :: a(:)
        real :: b(:)
        real :: c(:)
        integer,value :: N

        integer :: i
        
        do i=1,N
            c(i) = a(i) + b(i)
        end do

    end subroutine add
end module kernel

program main
    use kernel
    implicit none

    integer,parameter :: N = 2**20

    real,allocatable :: a(:), b(:), c(:)

    allocate(a(N),source=1.0)
    allocate(b(N),source=2.0)
    allocate(c(N),source=0.0)

    call  add(a, b, c ,N)
    print *,sum(c)/N

    deallocate(a)
    deallocate(b)
    deallocate(c)
end program main

動的割付け配列a,b,cに$2^{20}$個の要素を割り付けて，$a=1.0$, $b=2.0$として，$c_i = a_i+b_i$を計算します．特段言及しなければならない処理はないでしょう．

ベクトル和のGPU移植

ベクトル和のプログラムをCUDA Fortranに移植します．ここで大きな変更が入るのは，サブルーチンをカーネルに変更するとき，特に配列の各要素の足し算を行っているところです．

do i=1,N
    c(i) = a(i) + b(i)
end do

CUDAでは，カーネルには1スレッドが実行する処理を書くのでした．全スレッドがそれぞれ異なるiの担当となってc(i) = a(i) + b(i)を実行できれば，並列実行の度合いが一番高くなります．そこで，処理を1スレッドが処理する内容

    c(i) = a(i) + b(i)

に書き直します．このようなある1要素に対する処理が全要素に展開されるような挙動は，普段から配列演算c(:) = a(:) + b(:)やelementalな関数に触れているFortranユーザにとって，なじみやすいと感じています．

カーネルを変更した後，どのスレッドがどのiを計算するのかを，threadIdxなどの情報を基に決定します．スレッド番号と配列添字の対応を取るかは，決まった式があります．その式は次のように書きます．

i = (blockIdx%x-1)*blockDim%x + threadIdx%x

このスレッド番号と配列添字を対応付けるのはCUDA C/C++であろうがCUDA Fortranであろうが同じで，CUDAを使う上で最も重要です．

本当にiが重複しないのかは，以下のようなプログラムで確認できます．

module kernel
    implicit none

    contains

    attributes(global) subroutine arrayindex()
        implicit none

        print *,(blockIdx%x-1)*blockDim%x + threadIdx%x
    end subroutine arrayindex
end module kernel

program main
    use kernel
    use cudafor
    implicit none

    integer :: stat

    call arrayindex<<<2,4>>>()
    stat = cudaDeviceSynchronize()
end program main

実行してみると，表示される順番はそろっていませんが，表示される数字に重複がないことを確認できます．

5
6
7
8
1
2
3
4

CUDA Fortranの基本スタイルに則って，vectoradd.f90をGPUに移植します．

module kernel
    implicit none

    contains

    attributes(global) subroutine add(a, b, c, N)   ! attributes(global)を付ける
        implicit none
        real :: a(N)
        real :: b(N)
        real :: c(N)
        integer,value :: N

        integer :: i
        
        i = (blockIdx%x-1)*blockDim%x + threadIdx%x ! スレッドの番号と配列添字を対応づける
        c(i) = a(i) + b(i)

    end subroutine add
end module kernel

program main
    use kernel
    use cudafor
    implicit none

    integer,parameter :: N=2**20

    real,allocatable,device :: a(:)  ! device属性を付与
    real,allocatable,device :: b(:)  ! device属性を付与
    real,allocatable,device :: c(:)  ! device属性を付与
    real,allocatable :: host_c(:)    ! 結果を受け取るためのホスト変数を追加

    allocate(a(N),source = 1.0)
    allocate(b(N),source = 2.0)
    allocate(c(N),source = 0.0)

    call add<<<N/256,256>>>(a,b,c,N) ! シェブロンに並列実行の度合いを追加

    ! 足し算の結果をコピーして表示
    allocate(host_c(N))
    host_c = c
    print *,sum(host_c)/N

    deallocate(a)
    deallocate(b)
    deallocate(c)
    deallocate(host_c)
end program main

シェブロンの中に書く数字のうち，制約が大きいのも性能に影響を及ぼすのも，1スレッドブロックあたりのスレッド数です．その上限は今のところ1024です．1スレッドが配列の1要素を処理するつもりなので，総スレッド数＝配列要素数です．総スレッド数は，スレッドブロック数×1スレッドブロックあたりのスレッド数なので，先に1スレッドブロックあたりのスレッド数を決めると，スレッドブロック数は，配列要素数/1スレッドブロックあたりのスレッド数から求められます．

上のプログラムでは，1スレッドブロックあたりのスレッド数を256としているので，スレッドブロック数はN/256です．

CUDA Fortranの基本スタイル（追加2）

並列処理ができるようになりました．水泳ならクロールができるようになって，スイスイと泳げるようになっています．

この節で触れた基本スタイルを追加しておきます．

• カーネルには，1スレッドが実行する内容を書く．
• カーネル内では，各スレッドに割り振られた情報と配列添字の対応を計算して，全スレッドが一つの式から異なる配列要素を参照することで並列に処理を行う．
• シェブロンには，スレッドブロック数と1スレッドブロック数あたりのスレッド数を指定する．

実行時間を測る

せっかくGPUに移植したので，1スレッドブロックあたりのスレッド数が実行速度にどのような影響があるかを調べてみましょう．

Fortranで時間の測定を行う方法は先日の記事に書きました．しかし，Hello Kernelで説明したように，CPUはカーネルを実行した後，カーネルの終了を待たずに次の処理を行います．つまり，先日の記事のように，カーネルの前後で時間を測定するサブルーチンを呼んだとしても，カーネルを呼び出した直後に終了時間測定用のサブルーチンが呼び出され，正しい実行時間が測定できません．

そこで，CUDAの機能を使って実行時間を測定することにします．測定の理屈は先日の記事と同じです．CUDAでは，GPUへの指示をイベントとよぶのですが，そのイベントが発生した時間の差から時間を測定します．

まず，イベント発生時間を記録するための変数を宣言します．イベントを取り扱う派生型type(cudaEvent)の変数と計算時間を取得する変数を宣言し，cudaEventCreate()関数で利用準備を整えます．

type(cudaEvent) :: begin, end ! イベント記録用
real(4) :: elapsed_time_ms    ! 実行時間記録用

stat = cudaEventCreate(begin)
stat = cudaEventCreate(end)

利用準備が整ったら，実行時間を測定したいカーネルの前後でイベントを記録するためにcudaEventRecord()関数を呼びます．

stat = cudaEventRecord(begin,0)
call add<<<N/256,256>>>(a,b,c,N)
stat = cudaEventRecord(end, 0)

cudaEventRecord()の引数に先ほど宣言した派生型type(cudaEvent)の変数を渡します．二つ目の引数はstream IDとよばれる値です．複数のGPUを用いたり，GPUで複数のカーネルを並行に実行したりということがなければ，0でかまいません．

2回目のcudaEventRecord()を呼んだ後は，beginとendイベントの時間が正しく記録し終わっていることを保証するために，イベントの同期を取る関数cudaEventSynchronize()を呼びます．beginが先に呼ばれているので，endイベントの同期をとっておけば十分です．

stat = cudaEventSynchronize(end)

二つのイベントが発生した時間差を計算するのはcudaEventElapsedTime()関数です．結果は単精度浮動小数点型で，単位はミリ秒です．

stat = cudaEventElapsedTime(elapsed_time_ms, begin, end)

最後にcudaEventDestroy()関数の引数にイベントを渡し，後始末をします．

stat = cudaEventDestroy(begin)
stat = cudaEventDestroy(end)

この処理を前節のベクトル和のプログラムに追加して実行します．

module kernel
    implicit none

    contains

    attributes(global) subroutine add(a, b, c, N)
        implicit none

        real,device :: a(N)
        real :: b(N)
        real :: c(N)
        integer,value :: N

        integer :: i
        
        i = (blockIdx%x-1)*blockDim%x + threadIdx%x
        c(i) = a(i) + b(i)

    end subroutine add
end module kernel

program main
    use kernel
    use cudafor
    implicit none

    integer,parameter :: N=2**20

    real,allocatable,device :: a(:)
    real,allocatable,device :: b(:)
    real,allocatable,device :: c(:)

    integer,parameter :: numThread = 256
    integer,parameter :: numBlock  = N/numThread

    type(cudaEvent) :: begin, end
    real(4) :: elapsed_time_ms
    integer :: stat

    stat = cudaEventCreate(begin)
    stat = cudaEventCreate(end)

    allocate(a(N),source = 1.0)
    allocate(b(N),source = 2.0)
    allocate(c(N),source = 0.0)

    stat = cudaEventRecord(begin,0)
    call add<<<numBlock,numThread>>>(a,b,c,N)
    stat = cudaEventRecord(end, 0)

    stat = cudaEventSynchronize(end)
    stat = cudaEventElapsedTime(elapsed_time_ms, begin, end)
    print *,elapsed_time_ms,"ms"

    stat = cudaEventDestroy(begin)
    stat = cudaEventDestroy(end)

    deallocate(a)
    deallocate(b)
    deallocate(c)
end program main

実行時のパラメータを変更できるように，変数numThreadとnumBlockを導入しました．1スレッドブロックあたりのスレッド数を変更するだけで，配列の要素数とnumThreadから，スレッドブロック数numBlockを決定してくれます．

numThreadを変えて実行した結果を表にまとめました．ちなみに，CPUで実行した結果は，OpenMPの4並列で0.3438 msでした．

1ブロックあたりの スレッド数	実行時間 [ms]
32	0.1229
64	0.0765
128	0.0650
256	0.0657
512	0.0662
1024	0.0678
2048	0.0017

1ブロックあたりのスレッド数を大きくすると，段々と実行時間が短くなっていきますが，128を超えると遅くなっていきます．経験的には128か256が無難です．1ブロックあたりのスレッド数の最適な値を探すのが，CUDAの最も基本的なチューニングです．1ブロックあたりのスレッド数が2048のときに異様に高速化しているのは，設定できる値の上限を超えたためにカーネルが正しく実行されなかったことが理由です．カーネルのチューニングをしていると，あるときものすごく速くなって「おっ」と思うのですが，大体はカーネルが正しく実行されなかっただけというオチです．

こういったカーネルの起動や性能に影響するパラメータは他にもいくつかありますが，利用しているGPUよって設定できる上限値が変化します．利用しているGPUのパラメータは，PGIコンパイラと一緒にインストールされるpgaccelinfoで確認できます．

> pgaccelinfo

CUDA Driver Version:           10000

Device Number:                 0
Device Name:                   GeForce GTX TITAN
Device Revision Number:        3.5
Global Memory Size:            6442450944
Number of Multiprocessors:     14
Number of SP Cores:            2688
Number of DP Cores:            896
Concurrent Copy and Execution: Yes
Total Constant Memory:         65536
Total Shared Memory per Block: 49152
Registers per Block:           65536
Warp Size:                     32
Maximum Threads per Block:     1024
Maximum Block Dimensions:      1024, 1024, 64
Maximum Grid Dimensions:       2147483647 x 65535 x 65535
Maximum Memory Pitch:          2147483647B
Texture Alignment:             512B
Clock Rate:                    875 MHz
Execution Timeout:             Yes
Integrated Device:             No
Can Map Host Memory:           Yes
Compute Mode:                  default
Concurrent Kernels:            Yes
ECC Enabled:                   No
Memory Clock Rate:             3004 MHz
Memory Bus Width:              384 bits
L2 Cache Size:                 1572864 bytes
Max Threads Per SMP:           2048
Async Engines:                 1
Unified Addressing:            Yes
Managed Memory:                Yes
Concurrent Managed Memory:     No
PGI Default Target:            -ta=tesla:cc35

CUDA Fortranの基本スタイル（追加3）

• 実行時間はイベントを使って測定する．

まとめ

ね？簡単でしょ？

終わりに

Fortran Advent Calendar 2018に参加してくれた皆様のおかげで完走することができました．まずは皆様にお礼を申し上げます．

最終日の記事にこの内容を選んだのは，GPGPUがFortranに少なくない影響を与えたと感じているからです．

現在は機械学習の学習にGPUが利用されており，GPUが注目されています．機械学習に先駆けて，CAE業界では十数年前にGPUブームがありました．2010年頃までは，まともな開発環境はCUDA Cしかなかったので，Fortranユーザは頑張ってCUDA Cにコードを移植するか，見送るかという状況でした．

特に，小規模な計算機環境で高い性能を得るにはFortranでなくてもよいという状況になり，過去の資産とコンパイラ最適化でゴリ押しをしていたFortranという恐竜が，GPUという隕石による環境変化に対応できずに頭数を減らしたような印象を受けています．

FORTRANをCUDA Cに移植する人は，大抵は上から指示された若い人で，FORTRANは知らないけどCやCUDAはわかるという人たちです．今までFORTRANと縁遠かった人たちが濃縮されたFORTRANコードに触れ，実体験に基づくヘイトを溜める結果になっています．Fortranにとってこれはよい状況ではありません．

現在のCUDA Fortranの性能は，CUDA C/C++とあまり変わりません．FORTRANからCUDA C/C++に移植するよりは，FORTRAN→Fortran→CUDA Fortranと移植した方が，段階的に進められて移植の効率もよいはずです．ですが，FORTRANの資産の活用において，CUDA Fortranを使うという選択肢は現れてきません．そもそも使っている人が何人いるのかという状況でしょう．

Make Fortran great again!を目指している訳ではありませんが，ボチボチと情報を公開することで，自身の経験が誰かの役に立ったり，同じように情報を公開する人が増えて情報が入手しやすくなれば，みんな幸せになるんじゃないかと思っています．

Fortran Advent Calendarはよい勉強になりました．最初はノリでしたが，やって正解でした．
Fortranユーザの皆様，来年も豊かなFortranライフをお過ごしください．

参考資料

1. CUDA Programming Guide, https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html
2. PGI Community Edition（無償版） の導入方法, https://www.softek.co.jp/SPG/Pgi/comm_instruct.html
3. PGI Compiler Option 一覧, https://www.softek.co.jp/SPG/Pgi/TIPS/option1.html
4. Ruetsch, G. and Fatica, M., CUDA Fortran for Scientists and Engineers, https://www.amazon.co.jp/CUDA-Fortran-Scientists-Engineers-Programming/dp/0124169708

1. Pascal世代かVolta世代のGPU欲しいですね． ↩

2. 中には意識するという方もいらっしゃるでしょうが，アルゴリズムを考える段階では，そこまでは意識しないのではないでしょうか． ↩