Haskell による並列・並行プログラミング #8 6章 AccelerateによるGPUプログラミング つづき
担当: maton
前回からのおさらい
Accelerate では Repa に近い配列や添字のデータ型が提供されている
Accelerate における主な処理の流れ
use で配列をGPUに飛ばし、A.なんとか でGPU内での計算を記述し、run で取り出す
頻出関数
unit :: Elt e => Exp e -> Acc (Scalar e)
単一要素の0次元配列（スカラ配列）を作る
the :: Elt e => Acc (Scalar a) -> Exp e
スカラ配列から要素を取り出す
(!) :: (Shape ix, Elt e) => Acc (Array ix e) -> Exp ix -> Exp e
Accelerate の世界で使える添字アクセス
Exp ix 型の添字を作るには、 index1 :: Exp Int -> Exp (Z :. Int) を使う
index2 や index3 も。 あと一応 index0 もある
p. 109 6.6 Accの中で配列を作る
配列の生成方法があるならば、リストから配列を作るよりも、Accの世界で直接作ったほうが早い
少なくとも配列のデータをGPUメモリに移動するコストは減らせる
すべて同じ値で埋める fill
fill :: (Shape sh, Elt e) => Exp sh -> Exp e -> Acc (Array sh e)
数列からの生成 enumFromN, enumFromStepN
code:haskell
enumFromN :: (Shepe sh, Elt e, IsNum e)
=> Exp sh -> Exp e -> Acc (Array sh e)
-- シェイプ 最初の要素
enumFromStepN :: (Shepe sh, Elt e, IsNum e)
=> Exp sh -> Exp e -> Exp e -> Acc (Array sh e)
-- シェイプ 最初の要素 増分
使用例
code:haskell
run $ enumFromStepN (index2 3 5) 15 (-1) :: Array DIM2 Int
関数を用いた生成 generate
code:haskell
generate :: (Shape sh, Elt a)
=> Exp ix -> (Exp ix -> Exp a) -> Acc (Array ix a)
-- シェイプ 添字から要素の値への関数
使用例
code:haskell
run $ generate (index2 3 5) (\ix -> let Z:.y:.x = unlift ix in x + y
Array (Z :. 3 :. 5) 0,1,2,3,4,1,2,3,4,5,2,3,4,5,6 unlift :: Exp (Z :. Int :. Int) -> Z :. Exp Int :. Exp Int ※定義ではない
Acc 世界のインデックスを分解して DIM2 を得る（ここでは）
これは Unlift 型クラスの DIM2 におけるインスタンス
様々な配列、要素、シェイプに対して Lift, Unlift のインスタンスがある
lift :: Z :. Exp Int :. Exp Int -> Exp (Z :. Int :. Int)
実はこれを使って index2 は作られている
p. 111 6.7 2つの配列を綴じ合わせる
２つの配列の各要素間で関数を適用し、新たな配列をつくる zipWith
code:haskell
zipWIth :: (Shape sh, Elt a, Elt b, Elt c)
=> (Exp a -> Exp b -> Exp c) -- 要素の組に適用する関数
-> Acc (Array ix a) -> Acc (Array ix b) -- 入力される2つの配列
-> Acc (Array ix c)
使用例
code:haskell
let a = enumFromN (index2 2 3) 1 :: Acc (Array DIM2 Int)
let b = enumFromStepN (index2 2 3) 6 (-1) :: Acc (Array DIM2 Int)
run $ A.zipWith (+) a b
使用例2
code:haskell
let a = enumFromN (index2 2 3) 1 :: Acc (Array DIM2 Int)
let b = enumFromStepN (index2 3 5) 10 10 :: Acc (Array DIM2 Int)
run $ A.zipWith (+) a b
-- 左上を始点とした、配列の共通部分だけが取り出されて計算される。
[ 1, 2, 3, [ 10, 20, 30, 40, 50 [ 11, 22, 33,
4, 5, 6 ] + 60, 70, 80, 90,100 = 64, 75, 86 ]
110,120,130,140,150 ]
p. 112 6.8 定数
リテラルではないInt型の値をExp Intにするには？ → constant を使おう
constant :: Elt t => t -> Exp t
実例編
p. 112 6.9 例：最短路問題
accelerate-cuda インストールバトル（飛ばす）
なんと depreceted → accelerate-llvm-ptx を入れる必要がある
ついでに accelerate-llvm, cuda (haskell package), nvvm (haskell package) も
libnvvm3 とか llvm-8-dev も必要
fwaccel.hsの方もImportとか古い関数の置き換えが必要
閑話休題
fwaccel.hs
テキストの説明を追っていく
step :: Acc (Scalar Int) -> Acc Graph -> Acc Graph
各頂点 i, j について、これまでの最短路と頂点 k を通る最短路のうち短い方でパス i→j のコストを更新
引数に k :: Acc (Acalar Int) として渡すのはなぜ？
Acc データ型はCUDAのコード辺（カーネル）になっており、Accelerateは一度カーネル化された関数は再利用できるように試みる
Acc を使った式に引数を適用すると、GHCi 上で可視化できる
配列を渡すことで、generate の外側で unit 2 が束縛される
Int を渡してしまうと、 与えた Int が generate の内側で束縛され、 generate に渡す式が異なる k 毎に変化するため、再利用が妨げられてしまう
コラムまとめ: GPU上で計算させたいものはスカラ値だろうとAcc配列で渡そう
shortestPathsAcc :: Int -> Acc Graph -> Acc Graph
k <- [0..n-1] について step のリストを作り、 foldl1 (>->) で畳み込む
(>->) :: (Arrays a, Arrays b, Arrays c) => (Acc a -> Acc b) -> (Acc b -> Acc c) -> Acc a -> Acc c
flip (.) みたいな関数
合成する2つの関数の中間で作られる配列を残す必要がないことをAcceralateに伝えられる→次の計算が始まるタイミングで中間配列をGCできる。
shortestPaths :: Graph -> Graph
引数のグラフを use して shortestPathsAcc を run でラップする
動かしてみよう
stack ghci fwaccel.hs
shortestPaths testGraph
p. 115 6.9.1 GPU上で実行する
Repa版
stack exec -- fwdense1 2000 +RTS -s -N12
Total time 223.874s ( 22.441s elapsed)
stack exec -- fwdense1 100 +RTS -s -N12
Total time 0.029s ( 0.011s elapsed)
Accelerate版
stack exec -- fwaccel-gpu 2000 +RTS -s
Total time 1.905s ( 2.040s elapsed)
stack exec -- fwaccel-gpu 100 +RTS -s
Total time 0.064s ( 0.170s elapsed)
サイズが小さいとオーバーヘッドが目立つ
p. 115 6.9.2 CUDAバックエンドのデバッグ
accelerate-cuda に -fdebug をつけてビルド
deprecated なので accelerate-llvm-ptx でできるのか調べてみる
そのようなオプションは無い模様…
p. 116 6.10 例：マンデルブロ集合の生成器
accelerate v1 以降では module の構成が変わっているのでコンパイルエラーになる
Data.Array.Accelerate.IO は各コンテナクラスに合わせてインポートが必要
実行方法はおいといて、一旦本文解説へ
マンデルブロ集合
複素平面上の各点 $ c に対して漸化式 $ z_{n+1}=c+z_n^2 \ (z_0=c) を繰り返し適用した時に絶対値が発散しないような点の集合
絶対値が2を超えると発散する
簡単のため、$ z = x + iy として $ |z|^2 = x^2 + y^2 > 4 を検査
よく目にするカラフルなマンデルブロ集合は、発散までの回数に色調を割り当てている
Accelerate の中で漸化式を計算する
目標
Accelerateの式に対する関数を定義できるようになる
Accelerateの式に対する分岐を扱えるようになる
登場する型シノニム
type F = Float Floatが頻出するので短く（DoubleはGPUの互換性のために使わない模様）
type Complex = (F, F) 複素数。自前で用意するので後で二項演算も用意する
実は、 accelerate-0.15.0.0以降ならば複素数をそのまま扱える…
type ComplexPlane = Array DIM2 Complex 複素平面。
点 $ c における $ z の次の値を計算する next
code:haskell
next :: Exp Complex -> Exp Complex -> Exp Complex
next c z = c plus (z times z)
自前で用意した複素数に対する二項演算 plus と times を実装する必要がある
ここでは plus のみを見ていく
Exp Complex はパターンマッチを直接行えないので、以下のようには書けない
plus (ax, ay) (bx, by) = (ax + bx, ay + by)
ベースライン: Exp (a, b) 型から値を取り出せる A.fst A.snd を利用する
code:haskell
plus :: Exp Complex -> Exp Complex -> Exp Complex
plus a b = lift (ax + bx, ay + by) -- lift して Exp Complex に
where
ax = A.fst a
ay = A.snd a
bx = A.fst b
by = A.snd b
改良1: unlift してパターンマッチする
code:haskell
plus :: Exp Complex -> Exp Complex -> Exp Complex
plus a b = lift (ax + bx, ay + by) -- lift して Exp Complex に
where
(ax, ay) = unlift a :: (Exp F, Exp F) -- 型注釈が必須
(bx, by) = unlift b :: (Exp F, Exp F)
改良2: 二項演算を lift2 する
lift2 :: (Unlift Exp e1, Unlift Exp e2, Lift Exp e3) => (e1 -> e2 -> e3) -> Exp (Plain e1) -> Exp (Plain e2) -> Exp (Plain e3)
Control.Applicative.liftA2 のような気持ちで使える（ただし型注釈は必須）
code:haskell
plus :: Exp Complex -> Exp Complex -> Exp Complex
plus a b = lift2 f
where
f :: (Exp F, Exp F) -> (Exp F, Exp F) -> (Exp F, Exp F)
f (ax, ay) (bx, by) = (ax + bx, ay + by)
times も似たような感じで（複素数上の乗算になるように）
Data.Array.Accelerate.Data.Complex でのNumのインスタンスはほぼこういう感じ
各点に対して、何回目で発散条件を満たすか計算する iter
Accelerate の条件演算子を使う
(?) :: Elt t => Exp Bool -> (Exp t, Exp t) -> Exp t
Boolが真なら組の第一要素、偽なら第二要素を返す
GPU コードで条件式を用いるのはSIMD発散の原因になるのでよくない。
Single Instruction, Multiple Data: GPU は１つの命令を複数のスレッドで同時に実行する
分岐が現れると、全スレッドで真の場合を実行し、さらに全スレッドで偽の場合を実行する
分岐が入れ子になると、場合分けの数だけ実行が増加する→SIMD発散
CUDAドキュメントによれば、Trueの場合の処理をまとめて実行し、Falseの場合の処理をまとめて実行することになっている模様
Trueの場合の実行中はFalseの場合のスレッドは停止しているので、GPUの使用率が低下する。ので、一応SIMD発散（Branch Divergenceというらしい）という問題自体は残っている模様
工夫：反復当たりの分岐を１つだけに制限し、固定回数の反復にすることで計算量を抑える
各要素について、組 (z, i) を保持する
各反復では以下を計算
z' = next c z を計算する
絶対値の二乗$ |z|^2 =dot z >* 4 が真なら (z, i) を返す
dot はドット積 dot (x,y) = x*x + y*y
(>*) :: (Elt t, IsScalar t) => Exp t -> Exp t -> Exp Bool
そうでなければ (z', i+1) を返す
初期複素平面の生成 genPlane
複素平面上で切り取りたい四隅を指定して配列を生成
配列のインデックスに複素平面上座標を対応させる関数を書き、generateに渡す
(z, i) の組を持つ初期配列の生成 mkinit
genPlaneで生成した配列に対して (, 0) を lift したものを map する
トップレベルの関数 mandelbrot
iterate :: (a -> a) -> a -> [a] に反復関数 go と初期配列 zs0 を渡して反復回数番目を取り出す。
プログラム実行
ビルドできませんでした…
代わりに（？）公式サンプル集を実行します（いいのかそれで）
次回
II部 並行Haskell
7章 並行制御の基本：スレッドとMVar
wadoさん