Paraisoでシミュレーションを始める
シミュレーションを始めるには、まず初期条件を設定できる必要があります。初期条件の場を計算する仕事は、Paraiso式づくりのとっかかりとしても優れています。やってみましょう。サンプルを見てください。
table :: Named (StaticValue TArray Double) table = "table" `isNameOf` StaticValue TArray undefined createBuilder :: Builder Vec2 Int Annotation () createBuilder = do x01 <- bind $ (cast $ loadIndex (Axis 0)) / (cast $ broadcast $ loadSize (Axis 0)) y01 <- bind $ (cast $ loadIndex (Axis 1)) / (cast $ broadcast $ loadSize (Axis 1)) x <- bind $ 4 * (x01-0.5) y <- bind $ 5 * (y01-0.5) z <- bind $ atan((1- x^2 - (y-(x^2)**(1/3))^2)*10) store table z
Generator.hs
では、少々煩雑ですが、cast
とbroadcast
を使って配列変数の添字とサイズから、x座標とy座標を0から1までの範囲の実数として求めています。cast
やbroadcast
の引数や返値は、うまいこと型推論してもらえてるみたいですね。実際、配列変数の添字とサイズはcreateBuilder
の型注釈からしてInt
型なのに対し、table
はDouble
型で、そこから明示的なcast
を挟まずに繋がっているx01
,y01
,x
,y
,z
などの型もみなDouble
(のBuilder)になります。経験のある読者は、1/3
が一見整数どうしの除算として0
に切り捨てられることなく、ちゃんとDouble
同士の除算として扱われていることにも気づいたかもしれません。
それでは、これがParaisoと筆者から読者への気持ちです!
人生っていいもんだ
Paraisoを実用的な問題で試してみましょう。そうですね、生命現象のシミュレーションなどが実用的で良いのではないでしょうか。さっそくサンプルプログラムを見ていきましょう。まずはインポートです。
#!/usr/bin/env runhaskell {-# LANGUAGE NoImplicitPrelude, OverloadedStrings #-} {-# OPTIONS -Wall #-} import Control.Monad import Data.Tensor.TypeLevel import Language.Paraiso.Annotation (Annotation) import qualified Language.Paraiso.Annotation.Boundary as Boundary import Language.Paraiso.Generator (generateIO) import qualified Language.Paraiso.Generator.Native as Native import Language.Paraiso.Name import Language.Paraiso.OM import Language.Paraiso.OM.Builder import Language.Paraiso.OM.Builder.Boolean (select,eq,ge,le) import Language.Paraiso.OM.DynValue as DVal import Language.Paraiso.OM.Realm import qualified Language.Paraiso.OM.Reduce as Reduce import Language.Paraiso.OM.Value (StaticValue(..)) import Language.Paraiso.Optimization import Language.Paraiso.Prelude import NumericPrelude hiding ((||),(&&)) -- the main program main :: IO () main = do _ <- generateIO mySetup myOM return ()
80マス×48マスの大きさで、周期境界条件が課された計算領域を用意します。OMの名前はずばり"Life"にしましょう。
-- the code generation setup mySetup :: Native.Setup Vec2 Int mySetup = (Native.defaultSetup $ Vec :~ 80 :~ 48) { Native.directory = "./dist/" , Native.boundary = compose $ const Boundary.Cyclic } -- the orthotope machine to be generated myOM :: OM Vec2 Int Annotation myOM = optimize O3 $ makeOM (mkName "Life") [] myVars myKernels
セルオートマトンの状態を表す配列変数cell
、人口を数えるためのスカラー変数population
、および経過時間を数えるためのスカラー変数generation
という三つの変数を用意します。
-- the variables we use cell :: Named (StaticValue TArray Int) cell = "cell" `isNameOf` StaticValue TArray undefined population :: Named (StaticValue TScalar Int) population = "population" `isNameOf` StaticValue TScalar undefined generation :: Named (StaticValue TScalar Int) generation = "generation" `isNameOf` StaticValue TScalar undefined myVars :: [Named DynValue] myVars = [f2d cell, f2d population, f2d generation]
カーネルは2つ作ります。1つは初期化のためのカーネルで、もうひとつはシミュレーションを1ステップ進めるためのカーネルです。
>|haskell|
-
- our kernel
myKernels :: [Named (Builder Vec2 Int Annotation ())]
myKernels = ["init" `isNameOf` initBuilder,
"proceed" `isNameOf` proceedBuilder]
|
Paraisoの境界条件
直胞機械のshift
命令は配列変数全体を平行移動する命令です。実際のコードが扱う配列変数は常に有限サイズですから、平行移動によりはみ出す部分の処理を境界条件により指定してやる必要があります。Paraisoは、データフローグラフの中のshift
命令を解析して、境界条件を実現するのに必要な処理(袖領域など)を追加してくれます。
Paraisoでは、今のところ周期境界と自由境界という二種類の境界条件をサポートしています。
コード生成時、次元ごとにいづれかの境界条件を選ぶことができます。この二種類とほかの命令を組み合わて作れる範囲で、より複雑な境界条件も扱えます。境界条件は、Setup
データ型のboundary
レコードで、Condition
型を要素にもつベクトルとして指定します。
二種類の境界条件のうち、周期境界は、単に右端からはみ出したアクセスは左端のものを見るだけです。いわゆるドラクエ地図です。
これに対し、自由境界の場合は少々複雑です。まず、袖領域は、どのParaisoカーネルを1ステップ実行しても、
Setupで指定した計算したい領域、添字が0以上localSize
未満の領域が埋めるのに必要十分なサイズをもって決まります。
次に、各カーネルの計算する範囲は、袖領域をふくめた最大の領域から開始し、なるべく大きな領域を埋めるように決まります。その領域の外側での値は未定義です。また、Reduce
演算はつねに袖領域をふくむ全体を畳み込みます。もし袖領域を除外したい場合は手動でマスクする必要があります。
このサンプルプログラムで境界条件の挙動を確かめてみましょう。このプログラムには3つのカーネルがあります。
myKernels :: [Named (Builder Vec1 Int Annotation ())] myKernels = ["init" `isNameOf` do store table $ loadIndex (Axis 0) ,"increment" `isNameOf` do store table $ 1 + load table ,"calculate" `isNameOf` do center <- bind $ load table right <- bind $ shift (Vec :~ (-1)) center left <- bind $ shift (Vec :~ ( 1)) center ret <- bind $ 10000 * left + 100 * center + right store table ret store total $ reduce Reduce.Sum ret ]
それぞれ、「配列の内容を配列添字で初期化する」「配列の内容に1を加える」「ひとつ右とひとつ左の内容を読んで値を計算する」カーネルです。
maker.init();
dump();
maker.increment();
dump();
maker.calculate();
dump();
cout << "total: " << maker.total() << endl;
cout << endl;
make
を実行すると、境界条件だけが異なる2種類のプログラムが生成されます。このように、周期境界条件のもとでは境界を飛び越えて反対側の値にアクセスできます。
$ ./main-cyclic.out index: 0 1 2 3 4 5 6 7 value: 0 1 2 3 4 5 6 7 index: 0 1 2 3 4 5 6 7 value: 1 2 3 4 5 6 7 8 index: 0 1 2 3 4 5 6 7 value: 80102 10203 20304 30405 40506 50607 60708 70801 total: 363636
一方、自由境界のもとでは袖領域が追加され、0未満、あるいはlocalSize
以上の添字の範囲にもアクセスできます。
$ ./main-open.out index: -1 0 1 2 3 4 5 6 7 8 value: -1 0 1 2 3 4 5 6 7 8 index: -1 0 1 2 3 4 5 6 7 8 value: 0 1 2 3 4 5 6 7 8 9 index: -1 0 1 2 3 4 5 6 7 8 value: 0 102 10203 20304 30405 40506 50607 60708 70809 0 total: 283644
Builderモナド
Paraisoのプログラムは、Builder
モナドを組み合わせることで記述します。Builder
モナドは、実のところ作りかけのデータフローグラフを状態にもつState
モナドです。それぞれのBuilder
モナドは、データフローグラフの頂点をいくつか引数に取って、新たな頂点を返します。
type Graph vector gauge anot = Gr (Node vector gauge anot) Edge
できあがるグラフの型は、vector
、gauge
、anot
という3つの型引数を取るもので、これがBuilderにも常についてまわります。vector
は配列の次元(例:3次元)、gauge
は添字の型(例:Int)、anot
は解析や最適化に使う注釈の型でして、vector gauge
という組み合わせは配列にアクセスするための添字のベクトル(Intの3次元ベクトル)になります。
-- | value type, with its realm and content type discriminated in type level data Value rea con = -- | data obtained from the data-flow graph. -- 'realm' carries a type-level realm information, -- 'content' carries only type information and its ingredient is insignificant -- and can be 'undefined'. FromNode {realm :: rea, content :: con, node :: G.Node} | -- | data obtained as an immediate value. -- 'realm' carries a type-level realm information, -- 'content' is the immediate value to be stored. FromImm {realm :: rea, content :: con} deriving (Eq, Show)
この頂点はValue
型であり、各頂点にたいして領域情報rea
と、中身の型の情報con
を運んでいます。rea
の位置に入るのはLanguage.Paraiso.OM.RealmにあるTArray
かTScalar
のいづれかであり、それぞれ配列変数およびスカラー変数を表します。con
の位置に入るのはInt
とかDouble
とかいった要素型ですが、多くの場合は型情報のみがグラフの構築に使われ値は単に無視されます。この多相な[Value
型]を接点にグラフを組み立てることで、できあがったグラフの型安全性が(かなり)保証されます。
Paraisoプログラムを組み立てる材料はLanguage.Paraiso.OM.Builderモジュールに揃っています。モジュールのページに行って、"Synopsis"タブを押してください。
bindが便利なわけ
これまでもいくつかParaisoのプログラムはでてきましたが、やたらと各行ごとにbind
があるなあと思われたかもしれません。
bind :: (Monad m, Functor m) => m a -> m (m a) bind = fmap return
こんなふうに。
x <- bind $ loadIndex (Axis 0) y <- bind $ loadIndex (Axis 1) z <- bind $ x*y w <- bind $ x+y
このプログラムは素直に書くと、以下のようになります。
x <- loadIndex (Axis 0) y <- loadIndex (Axis 1) z <- return x * return y w <- return x + return y
右辺の式loadIndex (Axis 0)
などはBuilder
モナドから成り立っている一方で、左辺のx,y,z,w
などはデータフローグラフの頂点ですから、Value
型の値です。このため、いったん左辺でx
などの変数名を束縛したら、その後右辺で使うたびにモナドに変換する必要があります。(このとき使うのはむろん、最小の文脈を付与するreturn
です!)ところがbind = fmap return
により束縛時点で一度return
を施せば、以降、右辺ではそのまま使えます。こちらの方がずっと便利でしょう?
というわけで、Paraisoのサンプルプログラムに出てくるx,y
といった何気ない変数名はすべてモナディックな値であることに注意して下さい。
直胞機械の命令
直胞機械の(実際コンパクトな)命令セットはLanguage.Paraiso.OM.Graph
モジュールのInst
型として定義されており、Language.Paraiso.OM.Builder
モジュールには、これにほぼ一対一対応する形で材料が用意してあります。以下に、それぞれのコンビネータを表形式で整理しておきます。B
を一般のモナド記号m
だと(実際の定義はtype B ret = forall v g a. Builder v g a ret
)思うと感覚がつかめるのではないでしょうか。
これが命令セットで、
data Inst vector gauge = Load StaticIdx | Store StaticIdx | Reduce R.Operator | Broadcast | LoadIndex (Axis vector) | LoadSize (Axis vector) | Shift (vector gauge) | Imm Dynamic | Arith A.Operator
-- options input nodes output nodes load :: Named (StaticValue r c) -> B (Value r c) store :: Named (StaticValue r c) -> B (Value r c) -> B () reduce :: Reduce.Operator -> B (Value TArray c) -> B (Value TScalar c) broadcast :: B (Value TScalar c) -> B (Value TArray c) loadIndex :: Axis v -> B (Value TArray g) loadSize :: Axis v -> B (Value TScalar g) shift :: v g -> B (Value TArray c) -> B (Value TArray c) imm :: c -> B (Value r c)
各命令の機能をひとことで言うと、
-
load
は名前つき静的変数から読み込む。 -
store
は名前つき静的変数へ書き込む。 -
reduce
は配列をスカラーに畳み込む。 -
broadcast
はスカラーを配列にする。 -
loadIndex
は配列の添字を取得。 -
loadSize
は配列のサイズを取得。 -
shift
はベクトルv g
を使って配列をずらす。
-
imm
は即値を読み込む。 -
Arith
はさまざまな算術演算を施す。
となります。各命令の詳細については、これから例で見ていきます。
演算子たち
そういえば、Arith
命令に対応するコンビネータが見当たりませんが、どうなってるのでしょうか。Paraisoでは、numeric-preludeというライブラリを使い、Buiderモナドをさまざまな代数構造のインスタンスにすることで、普通の数式を書くのと同じ感覚でBuiderモナドの数式を組み立てられるようになっています。
(+) :: B (Value r c) -> B (Value r c) -> B (Value r c) sin :: B (Value r c) -> B (Value r c)
ただし、HaskellのBoolを扱う演算子に関しては、(==) :: Eq a => a -> a -> Bool
などの固定された型を持っているため、BoolのBuilderを返させることができません。そのため、モジュールLanguage.Paraiso.OM.Builder.BooleanやLanguage.Paraiso.Prelude.で定義されている関数で代用してください。またif
も同じ理由でBuilderを組み立てる用途には使えませんので、代わりにselect
を使ってください。
eq :: B (Value r c) -> B (Value r c) -> B (Value r Bool) ne :: B (Value r c) -> B (Value r c) -> B (Value r Bool) select :: B (Value r Bool) -> B (Value r c) -> B (Value r c) -> B (Value r c)
最後に、型変換のための関数cast
があります。具体的にどういう型変換コードが生成されるかは対象言語しだいです。ここでc2
は変換先の型を指定するための変数で、値は使われません。
cast :: c2 -> B (Value r c1) -> B (Value r c2)
せきじいのサイト
先輩の関口さんがサイトを作ったらしいが検索にでてこないらしい。どうしたのだろう。
http://sites.google.com/site/yisekig/
などと、google juiceを与えてみるテスト。
はじめてのGPUプログラミング
それでは、いよいよGPUプログラミングをやってみましょう。C++コードの生成器をCUDA用のもの.に書き換えるのは、実際すごく簡単です。
$ diff HelloWorld/Generator.hs HelloGPU/Generator.hs 40c40,42 < { Native.directory = "./dist/" --- > { Native.directory = "./dist/" , > Native.language = Native.CUDA, > Native.cudaGridSize = (32,1)
このように、ネイティブコード生成の詳細を指定するSetup
型の値において、CUDAコードを生成したいということと、CUDAのグリッドサイズを指定するだけです。(以下の指定の場合、すべてのCUDAカーネルは<<<32,1>>>
のサイズをもって起動されます。現在のParaisoにはカーネルごとにグリッドサイズを変える機能はありません。)
mySetup :: Native.Setup Vec2 Int mySetup = (Native.defaultSetup $ Vec :~ 10 :~ 20) { Native.directory = "./dist/" , Native.language = Native.CUDA, Native.cudaGridSize = (32,1) }
生成されたコードは東京工業大学のTSUBAME2.0で、nvcc4.1(と内臓されたthrust)を用いて、C++版と同じ出力をすることを確かめました。
■
NGS現場の会 第二回研究会に行ってきた。