誰向けの記事?

FortranでGPUを使って乱数を使うようなシミュレーションをしたい人向け.

CUDA Fortranとは?

GPUで数値計算するためのもの.
いつのまにか, pgfortran から nvfortran になっている.
ダウンロードとかは https://developer.nvidia.com/hpc-sdk.

cuRANDとは?

NVIDIA公式のライブラリで, 乱数を生成できる.
- https://docs.nvidia.com/cuda/curand/index.html.

実行環境

Ubuntu 20.04.6 LTS
CUDA 12.2
NVIDIA GeForce RTX4090

やってみよう

test_curand_simple.f90

program test_curand
    use, intrinsic :: iso_fortran_env
    use cudafor
    use curand ! これだけ
    implicit none
    integer(int32), parameter :: n = 10
    integer(int32), device :: a_d(n) ! GPU側の変数.
    integer(int32) :: a_h(n)
    real(real64), device :: dr_d(n) ! GPU側の変数.
    real(real64) :: dr_h(n)
    type(curandGenerator) :: gen
    integer :: i, istat

    istat = curandCreateGenerator(gen, CURAND_RNG_PSEUDO_XORWOW)
    istat = curandGenerate(gen, a_d, n)  ! integer(int32)型の配列に乱数詰める.
    istat = curandGenerate(gen, dr_d, n) ! real(real64)型の配列に乱数詰める.

    ! cudaMemcpyDeviceToHost要らず.
    a_h(:)  = a_d(:)
    dr_h(:) = dr_d(:)

    write(output_unit, '(g0)') a_h(:)
    ! できない write(output_unit, '(g0)') a_d(:)
    write(output_unit, '(g0)') dr_h(:)
    ! できない write(output_unit, '(g0)') dr_d(:)

    stop "END PROGRAM"
end program test_curand

コンパイルと実行

$ nvfortran -cuda -cudalib=curand test_curand_simple.f90
$ ./a.out
-1115749450
-339329097
167591721
-133303299
-321518802
1917059131
-1428853312
472148682
2019573489
-2090817275
0.5084141922664625
0.6545496961753901
0.5126043582083046
0.2643004928742966
0.5198064948091491E-01
0.5789965062976243
0.3855563551535925
0.9082153289929680
0.6416203630205781
0.2833991036050670
Warning: ieee_inexact is signaling
END PROGRAM

たったこれだけ

use curand を implicit noneの前につける.
curandCreateGeneratorを使って, 変数 gen を初期化して, curandGenerate で配列に乱数を詰め込む.
CURAND_RNG_PSEUDO_XORWOW は乱数生成器のID.
- 他にもいっぱいある. https://docs.nvidia.com/cuda/curand/group__HOST.html
(2024/05/27追記) 64bit実数配列を引数に取る curandGenerate は curandGenerateUniformDouble を呼び出しているので, 乱数の範囲は (0, 1] らしい (https://tech.garilog.com/cuda-curand/)

平均値を取ってみる

test_curand_average.f90

program test_curand
    use, intrinsic :: iso_fortran_env
    use cudafor
    use curand
    implicit none
    integer(int64), parameter :: n = 100000000
    real(real64), parameter :: exact_average = 0.5d0
    integer(int64) :: n_times
    real(real64), device :: dr_d(n)
    real(real64) :: summ, average
    type(curandGenerator) :: gen
    integer(int32) :: i, istat
    istat = curandCreateGenerator(gen, CURAND_RNG_PSEUDO_XORWOW)

    write(error_unit, '(a)') "Input repeat number: "
    flush(error_unit)
    read(input_unit, *) n_times
    write(error_unit, '(a, i0, a)') "Average ", n * n_times, " random numers"

    summ = 0.0d0
    do i = 1, n_times
        istat = curandGenerate(gen, dr_d, n)

        summ = summ + sum(dr_d(:))
    end do

    average = summ / (n * n_times)
    write(output_unit, '(i0, 2(1x, g0))') n * n_times, average, abs(average - exact_average)

    stop "END PROGRAM"
end program test_curand

コンパイルと実行

$ nvfortran -Ofast -cuda -cudalib=curand test_curand_average.f90
$ time ./a.out <<< "1000"
Input repeat number: 
Average 100000000000 random numers
100000000000 0.5000002939475188 0.2939475187702101E-06
Warning: ieee_inexact is signaling
END PROGRAM

real	0m2.387s
user	0m2.301s
sys	    0m0.060s

結論

$10^{11}$ 個の乱数を2秒くらいで生成できた.

乱数の初期シードのせってい

数値計算では疑似乱数の初期シードを設定して, 同じ乱数列を生成できるようにする.
curand での初期シードの設定の仕方は curandSetPseudoRandomGeneratorSeed を使う(https://docs.nvidia.com/cuda/curand/host-api-overview.html).

  istat = curandSetPseudoRandomGeneratorSeed(gen, iseed)

iseed は 4バイト整数.

実際のシミュレーション！

CUDA Fortranを用いた 2次元Ising模型のシミュレーションは
Fortran Advent Calendar 2023 19日目へ...

参考