JuliaGPU / CUDAnative でイジングモデルを計算する

# はじめに

`JuliaGPU / CUDAnative` を試してみて非常に良かったので、とりあえず概要を紹介する。
書ききれなかったことは別記事に仕立てるかもしれない。

`CUDAnative` は、`CuArrays` と比べて

- 低レベルのコーディングが必要
- 最適化されていて高速

開発者から最適化についての言及
https://discourse.julialang.org/t/cudanative-is-awesome/17861/2

`JuliaGPU / ArrayFire` でイジングモデルを計算する
https://qiita.com/Lirimy/items/3ecb550af79cc674e480
※ `ArrayFire` は `CuArrays` と同程度に高抽象度の外部ライブラリ

`CUDAnative` の詳しい使い方は見当たらないが
https://github.com/JuliaGPU/CUDAnative.jl/tree/master/examples
が参考になった

この記事のノートブック
https://nbviewer.jupyter.org/gist/Lirimy/c825ab9cc45c45a2173c52254d4802c9


# イジングモデルの並列処理

GPU 計算では、多数のコアを利用した並列化によってパフォーマンスを上げている。

一方、イジングモデルでは、スピンは隣のスピンの影響を受けるため、隣接するスピンを同時に評価することができない。

影響が直接及ばないスピンならば、並列に評価することができる。
スピン全体を格子模様で塗りわけて、色ごとに交互に評価すればよい。

![checkeredpattern.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202534/675b3c71-44e5-9f9f-b015-027fc77fabbb.png)

Simulation of Ising spin model using CUDA
http://www.vinkovic.org/Projects/MindExercises/radnje/2011_Jurisic.pdf

※ GPU のパフォーマンスは CPU と比較して 50-200 倍との記述あり


# CUDA プログラミング

## GPU での処理を記述する

GPU で呼び出される関数、いわゆるカーネルを書いていく。

並列処理以外の部分については、黒木玄による CPU 実装をほぼそのまま踏襲した。
https://nbviewer.jupyter.org/gist/genkuroki/4fa46c68c56ee0f3b1a6fc8ec628b9d7

```julia
using CUDAnative

function _update!(s, rnd, m, n, β, white=false)
    
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    j = (blockIdx().y-1) * blockDim().y + threadIdx().y
    
    if i ≤ m && j ≤ n && iseven(i+j) == white
        
        @inbounds ajs = (
            s[ifelse(i+1 ≤ m, i+1, 1), j] +
            s[ifelse(i-1 ≥ 1, i-1, m), j] +
            s[i, ifelse(j+1 ≤ n, j+1, 1)] +
            s[i, ifelse(j-1 ≥ 1, j-1, n)]
        )
        @inbounds prob = CUDAnative.exp_fast(-β * s[i, j] * ajs)
        @inbounds s[i, j] = ifelse(rnd[i, j] < prob, -s[i, j], s[i, j])
    end
    
    return
end
```

- CPU で計算するときの For ループの中身だけを抜き出したと考えればよい
- 関数内で、スレッド番号などから、自身が割り当てられたインデックスを求めている
- もしインデックスに対応するスピンがなければ処理をスキップする（スピンの数が 32 の倍数でない場合）
- 格子模様によるマスキングも、同じように処理のスキップで対応する（他の方法もある）

## GPU でスレッドを立ち上げる

`@cuda` マクロでカーネルを呼び出す。

```julia
using Random: rand!
using CuArrays.CURAND

m = 3072
n = 4096

rnd = curand(m, n)
s = @. Int32(2 * (rnd < 0.5f0) - 1)
βc = Float32(log(1+sqrt(2)))
threads = (32, 32)
blocks = @. ceil(Int, (m, n) / threads)

# Main
for _ in 1:100
	rand!(rnd)
	@cuda blocks=blocks threads=threads _update!(s, rnd, m, n, βc, true)
	@cuda blocks=blocks threads=threads _update!(s, rnd, m, n, βc, false)
end
```

- GPU では 32 スレッド (Warp) を単位として処理する
- スレッドのまとまりをブロックという
- ブロックあたりのスレッド数には制限がある（GTX 980 なら 1024）
- `CuArrays` で得た配列を `CUDAnative` のカーネルで利用できる
- 乱数をカーネル内で得るのは~~難しそう~~記述が長くなりそうだったので、引数で渡した


# ギャラリー

## 臨界温度での定常状態（静止画）

![ising.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202534/654b6c40-9695-4245-70c3-e0a405e78fd7.png)

全体をおよそ 100000 回評価後
$4096\times3072$ を 1/8 に縮小

## 臨界温度での定常状態（動画）

![ising.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202534/51f0e62e-7bdb-cec4-b23d-2b5ec9005b41.gif)

$4096\times3072$ を 1/16 に縮小
全体を 16 回評価／フレーム
30 フレーム

# ベンチマーク

![benchmark.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202534/63049f4b-6384-5c46-90d1-465a0f50f1e7.png)

- GPU は CPU と比較して、適した条件ではおよそ 80 倍高速
- モデルのサイズが小さいと CPU とパフォーマンスが変わらない
- 逆に大きすぎても（$32768\times32768$）メモリが不足する
- 最低でも $512\times512$ 、できれば $4096\times4096$ は欲しい