0. はじめに

もともと所有するFortranコードをGPGPU移植するにあたって、Cuda Fortranが最も移植コストが少ないと考えた。
Monte Carlo法を有するコードなのでcuRANDを使いたいと思って調べたところ、Cuda Fortranのライブラリから直接読み込めるようであった。
乱数長が事前にわからないので、curand_init/curand_uniformを使ってカーネル内から直接cuRANDを使用したかった。
しかしその解説記事が存在しておらず、ChatGPT (2024年5月15日時点のGPT-4o)を用いても、情報古いのか、CUDA Cで書いたコードをCUDA Fortranカーネルから呼び出すように言われた。
最終的に実装できたのでその方法とTipsを示す。

環境

OS： Ubuntu 22.4
コンパイラ： nvfortran 2024 / 24.3
GPU: RTX 1660 Super

1. 実装

はじめにchatGPTに聞いて出てきたコードが以下

質問：
CUDA Fortranで各スレッドが独自の乱数を生成するためには、cuRANDライブラリを使用します。以下に、10個の乱数を別々のスレッドで生成するためのテストコードを示します。

module curand_module
    use cudafor
    implicit none
    include 'curand_kernel.h'

contains

    attributes(global) subroutine setup_kernel(state, seed)
        type(curandState), device :: state(*)
        integer, value :: seed
        integer :: id

        id = threadIdx%x + (blockIdx%x - 1) * blockDim%x
        call curand_init(seed, id, 0, state(id))
    end subroutine setup_kernel

    attributes(global) subroutine generate_kernel(state, result)
        type(curandState), device :: state(*)
        real(kind=4), device :: result(*)
        integer :: id

        id = threadIdx%x + (blockIdx%x - 1) * blockDim%x
        result(id) = curand_uniform(state(id))
    end subroutine generate_kernel

end module curand_module

program main
    use cudafor
    use curand_module
    implicit none

    integer, parameter :: n = 10
    integer, parameter :: seed = 1234
    type(curandState), device, allocatable :: dev_states(:)
    real(kind=4), device, allocatable :: dev_random_numbers(:)
    real(kind=4), allocatable :: random_numbers(:)

    ! デバイスメモリのアロケート
    allocate(dev_states(n))
    allocate(dev_random_numbers(n))
    allocate(random_numbers(n))

    ! 乱数状態の初期化
    call setup_kernel<<<1, n>>>(dev_states, seed)

    ! 乱数生成
    call generate_kernel<<<1, n>>>(dev_states, dev_random_numbers)

    ! 結果をホストにコピー
    random_numbers = dev_random_numbers

    ! 結果を表示
    print *, 'Generated random numbers:'
    print *, random_numbers

    ! メモリの解放
    deallocate(dev_states)
    deallocate(dev_random_numbers)
    deallocate(random_numbers)

end program main

include 'curand_kernel.h'
となっているし、適切に出力できていない。
これでできないですと返したら、

1.CUDA Cコードで乱数生成関数を実装
2.FortranモジュールでCUDA C関数のインターフェースを定義
3.Fortranメインプログラムで乱数生成を呼び出す

を提案してきた。いくつか他の聞き方をしてもイマイチ理解していないので、自分で調べることにした。

コード改修

上記のコードをベースに動かすために重要だったポイントを記す。
だいたいここ (https://docs.nvidia.com/hpc-sdk/compilers/fortran-cuda-interfaces/index.html#cfrand-runtime) のサンプルコードを見たらわかった。

type(curandState)は使えない

代わりに、type(curandStateXORWOW) が使用できる。

curand_deviceをインクルードする必要がある

use curand_deviceを最初に入れる。

コンパイルオプションに-gpu=nonvvmを入れる必要がある

旧情報だと-Mcuda=nonvvmだが、最新版コンパイラなら-gpu=nonvvmと記載する。

2. 最終版コード

上記以外にもいくつか修正を行いながら最終的に問題なく動くと思われる状態になったコードを以下に示す。

モジュール

curand.cuf

module curand_module
    use cudafor
    ! use curand
    use curand_device
    implicit none

contains

    attributes(global) subroutine setup_kernel(n, state, seed)
        integer, value :: n
        type(curandStateXORWOW), device :: state(*)
        integer, value :: seed
        integer :: id

        id = threadIdx%x + (blockIdx%x - 1)*blockDim%x
        if (id <= n) then
            call curand_init(seed, id, 0, state(id))
        end if
    end subroutine setup_kernel

    attributes(global) subroutine generate_kernel(n, state, result)
        integer, value :: n
        type(curandStateXORWOW), device :: state(*)
        double precision, device :: result(*)
        integer :: id

        id = threadIdx%x + (blockIdx%x - 1)*blockDim%x
        if (id <= n) then
            result(id) = curand_uniform(state(id))
        end if
    end subroutine generate_kernel

end module curand_module

メインプログラム

main.f90

program main
    use curand_module
    implicit none

    integer, parameter :: n = 1000000
    integer :: seed
    type(curandStateXORWOW), device, allocatable :: dev_states(:)
    double precision, device, allocatable :: dev_random_numbers(:)
    double precision, allocatable :: random_numbers(:)
    integer :: values(8) ! to hold date_and_time values
    double precision :: start_time, end_time, elapsed_time
    double precision :: mean, stddev, sumrand, sumsq
    integer :: i, n_thread, n_block
    integer(kind=c_size_t) :: mem_size

    ! 現在の時間を取得してシード値を生成
    call date_and_time(values=values)
    seed = values(5)*10000 + values(6)*100 + values(7) ! 時・分・秒を組み合わせてシード値を生成

    print *, 'Seed value: ', seed

    ! デバイスメモリのアロケート
    allocate (dev_states(n))
    allocate (dev_random_numbers(n))
    allocate (random_numbers(n))

    mem_size = size(dev_states) * c_sizeof(dev_states(1))
    print *, 'Memory size of dev_states: ', mem_size, ' bytes'

    mem_size = size(dev_random_numbers) * c_sizeof(dev_random_numbers(1))
    print *, 'Memory size of dev_random_numbers: ', mem_size, ' bytes'

    ! 乱数生成開始時刻の記録
    call cpu_time(start_time)

    n_thread = 256
    n_block = (n + n_thread - 1) / n_thread
    ! 乱数状態の初期化
    call setup_kernel <<<n_block, n_thread>>> (n, dev_states, seed)

    ! 乱数生成
    call generate_kernel <<<n_block, n_thread>>> (n, dev_states, dev_random_numbers)

    ! 乱数生成終了時刻の記録
    call cpu_time(end_time)

    ! 結果をホストにコピー
    random_numbers = dev_random_numbers

    ! 経過時間を計算
    elapsed_time = end_time - start_time

    ! 乱数の平均値と標準偏差を計算
    sumrand = 0.0d0
    sumsq = 0.0d0
    do i = 1, n
        sumrand = sumrand + random_numbers(i)
        sumsq = sumsq + random_numbers(i)**2
    end do
    mean = sumrand/n
    stddev = sqrt((sumsq/n) - mean**2)

    ! 結果を表示
    print *, 'Generated random numbers (first 10):'
    print *, random_numbers(1:10)
    print *, 'Elapsed time (seconds): ', elapsed_time
    print *, 'Mean: ', mean
    print *, 'Standard Deviation: ', stddev

    ! メモリの解放
    deallocate (dev_states)
    deallocate (dev_random_numbers)
    deallocate (random_numbers)
end program main

なお上記で配列のメモリを確認したり平均・標準偏差をとっているのは、乱数が適切に生成されているかたしかめるため。

コンパイルコマンド

nvfortran -cuda -gpu=nonvvm curand.cuf main.cuf -o main

3. 実行結果

 Seed value:        142021
 Memory size of dev_states:                  48000000  bytes
 Memory size of dev_random_numbers:                   8000000  bytes
 Generated random numbers (first 10):
   0.1877072751522064        0.2312117666006088        0.6245507597923279      
   0.6655896306037903        0.7520658969879150        0.4479118287563324      
   0.8056853413581848        0.8860529661178589        0.4329888820648193      
   0.4152245223522186     
 Elapsed time (seconds):    4.4918060302734375E-004
 Mean:    0.4998809951801293     
 Standard Deviation:    0.2885744465300747

4. 参考