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)
はじめての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++版と同じ出力をすることを確かめました。
Paraisoが生成するC++コードのAPI
今回はParaisoが生成するコードの使い方を見ていきます。HelloWorldフォルダでGenerator.hsを実行すると、distというフォルダが作られてその下にC++のヘッダやソースが生成されているはずです。まずはヘッダを見てください。
class TableMaker{ private: std::vector<int> om_s0_table; private: int om_s1_total; private: std::vector<int> om_m0_5; public: TableMaker () : om_s0_table(om_memory_size()), om_m0_5(om_memory_size()) { }
このように、静的変数のための記憶領域やコンストラクタが定義されています。Paraisoの生成するクラスは、この他にもいくつかのアクセッサ関数を定義します。まずは、OMのサイズを取得するための関数がこちらです。全領域のサイズおよび、各次元ごとのサイズを取得できます。（これサイズがint型なのはまずいな・・・）
public: int om_size () { return 200; } public: int om_size_0 () { return 10; } public: int om_size_1 () { return 20; }
Paraisoが副次的に生成する関数の名前はすべてom_から始まるようになっています。ユーザーは静的変数やカーネルに、すべて異なり、かつどれもom_から始まらないような名前をつけてください。そうすればコード生成時に名前が衝突することはありません。
今のHelloWorldカーネルにはありませんが、配列変数で、近傍の値をとってくるshift命令を含むカーネルをコンパイルすると、Paraisoは通信のための余白領域を確保し、[0..size-1]までの領域が無事に計算できるようにします。この場合、実際に確保されたメモリの大きさはom_size()等が返す値よりも大きくなります。このメモリサイズを取得できる関数群がom_memory_size()などです。
public: int om_memory_size () { return 200; } public: int om_memory_size_0 () { return 10; } public: int om_memory_size_1 () { return 20; }
そしてこれが余白を取得する関数群です。
public: int om_lower_margin_0 () { return 0; } public: int om_lower_margin_1 () { return 0; } public: int om_upper_margin_0 () { return 0; } public: int om_upper_margin_1 () { return 0; }
Array変数に対しては要素へのアクセッサが、Scalar変数に対しては単なるアクセッサが定義されます。
public: int & table (int i0, int i1) { return (om_s0_table)[((om_lower_margin_0()) + (i0)) + ((om_memory_size_0()) * ((om_lower_margin_1()) + (i1)))]; } public: int & total () { return om_s1_total; }
また、配列の生データへのアクセッサも提供されますが、「生」データが具体的にどんなフォーマットなのかは実装の詳細に依存します。
public: std::vector<int> & table () { return om_s0_table; }
最後に、OMの定義で与えたカーネルに対応するメソッドが生成されます。
public: void create () ; };
はじめてのParaisoプログラミング
Paraisoは、一様メッシュ上での偏微分方程式の陽解法という分野に特化し、その分野の計算を、直胞機械(OrthotopeMachineの訳語、今決めた)と呼ばれる仮想機械のプログラムとして表します。直胞機械は、配列変数とその操作を表現する仮想機械であり、解きたい問題の内包する並列性をなるべく生かすように設計されています。たとえば、配列変数の要素を格納する順序や、アクセスする順序は自由。命令どうしの依存関係もはじめから明確で、依存性のない命令を実行する順序も、同時に実行することも自由です。Paraisoは、この直胞機械のプログラムを記述し、それを既存の言語のプログラムに翻訳し、さらに自動チューニングを行うための様々なメカニズムを備えたフレームワークです。より詳しい用語の解説などはParaisoの論文を見てください。
しかし、Paraisoにコードを生成させるためにユーザーがすべきことを、具体的に一言でいうなら、「OM型の値を作り、コードジェネレータに渡すこと」となります。
ハローワールドの中身
前回掲載したサンプルexamples/HelloWorld/Generator.hsには、Paraisoにコードを生成させるための全てが記されています。順にみていきましょう。
最初の十数行は#include みたいなもんで、必要なモジュールをインクルードしています。
#!/usr/bin/env runhaskell {-# LANGUAGE NoImplicitPrelude, OverloadedStrings #-} {-# OPTIONS -Wall #-} import Data.Tensor.TypeLevel import qualified Data.Text.IO as T import Language.Paraiso.Annotation (Annotation) 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.DynValue (f2d, DynValue) import Language.Paraiso.OM.Value (StaticValue(..)) import Language.Paraiso.OM.PrettyPrint (prettyPrintA1) import Language.Paraiso.OM.Realm import qualified Language.Paraiso.OM.Reduce as Reduce import Language.Paraiso.Optimization import NumericPrelude import System.Process (system)
まずはmainから見ていきましょう。一行目は"output"というフォルダを作っています。二行目では、myOMという名の直胞機械のデータフロー・グラフを出力しています。これは参考のために出力しているのであって、コード生成に必須の過程ではありません。三行目では、myOMの定義とmySetupというコード生成時の設定を渡してC++のコードを生成しています。
-- the main program main :: IO () main = do _ <- system "mkdir -p output" T.writeFile "output/OM.txt" $ prettyPrintA1 $ myOM _ <- generateIO mySetup myOM return ()
直胞機械を定義する時には使わないが、コード生成時にはじめて必要になる情報はLanguage.Paraiso.Generator.NativeモジュールにあるSetup型でまとめて指定します。ここでは具体的な計算のサイズや、ライブラリを生成するフォルダを指定しています。コード生成対象の言語はデフォルトではC++です。
-- the code generation setup mySetup :: Native.Setup Vec2 Int mySetup = (Native.defaultSetup $ Vec :~ 10 :~ 20) { Native.directory = "./dist/" }
次はいよいよmyOMの定義です。直胞機械を定義するには、名前(tableMaker)、マシン全体の注釈、静的変数のリスト、カーネルのリストの4つをmakeOM関数に渡します。
-- the orthotope machine to be generated myOM :: OM Vec2 Int Annotation myOM = optimize O3 $ makeOM (mkName "TableMaker") [] myVars myKernels
myOMには、九九の表を保持するためのtableと、その合計を計算するためのtotalという2つの変数があります。2つの変数のRealmは、それぞれArrayとScalarです。Array変数は解きたい問題の解像度とおなじサイズをもつ配列変数で、Scalar変数は普通に1つしかない値です。それぞれの変数はNameも持っています。
-- the variables we use table :: Named (StaticValue TArray Int) table = "table" `isNameOf` StaticValue TArray undefined total :: Named (StaticValue TScalar Int) total = "total" `isNameOf` StaticValue TScalar undefined myVars :: [Named DynValue] myVars = [f2d table, f2d total]
Paraisoでは、Language.Paraiso.Nameモジュールで定義されているName型を識別子の名前として使います。これは、ただの文字列と区別するために、Data.Textモジュールの提供するマルチバイト文字列型Textをnewtypeしたものです。
次にカーネルの定義です。OMのカーネルとは、一度に実行される計算のかたまりです。myOMにはcreateという名前のカーネルがひとつだけあります。
-- the only kernel our OM has myKernels :: [Named (Builder Vec2 Int Annotation ())] myKernels = ["create" `isNameOf` createBuilder]
カーネルを作るにはBuilder Monadを使います。ここでは、x軸とy軸の添字を読み込んでその積を計算するカーネルを作っています。ついでに九九の表の総計も求めてみましょう。
-- the kernel builder monad createBuilder :: Builder Vec2 Int Annotation () createBuilder = do x <- bind $ loadIndex (Axis 0) y <- bind $ loadIndex (Axis 1) z <- bind $ x*y store table z store total $ reduce Reduce.Sum z
