CUDA FortranでcuRANDを使う

# 誰向けの記事?
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>).

```fortran
  istat = curandSetPseudoRandomGeneratorSeed(gen, iseed)
```
`iseed` は 4バイト整数.

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

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

https://qiita.com/osada-yum/items/8d2b6906678f985b56cb

# 参考
https://developer.nvidia.com/cuda-fortran

https://qiita.com/implicit_none/items/8229d1931cd236d62ca9