京都大学 白眉センター
Since 2010
Kyoto Univ.

学問の新たな潮流を拓く、広い視野、柔軟な発想を求む！

• どんな分野でもOK• 5年間、好きな研究をしてOK• 毎年２０人採用

第24回理論懇親会：理論天文学・宇宙物理学の革新
2011/11/07

プログラミングという作業の革新
京都大学白眉センター
村主崇行
check!

☞

http://paraiso-lang.org/wiki/

quick start guideInstall Haskell Platform and git, then type> git clone git@github.com:nushio3/Paraiso.git> cd Paraiso/> cabal install> cd examples/Life/#Conway’s game of life example> make lib> ls output/OM.txt
output/OM.txt#this is analysis result for dataflow graph> ls dist/Life.cpp Life.hpp#an OpenMP implementation> ls dist-cuda/Life.cu Life.hpp#a CUDA implementation> cd ../Hydro/#hydrodynamics simulator example> make lib#this takes half a minute or so> ls output/; ls dist/; ls dist-cuda/#same as above2009年8月@うちの教室
8個のGPUを搭載した1台のパソコン
単精度演算性能7.3Tflops

GPGPU: General-Purpose Computation on GPUs
M. Harris et al (2002) who coined the name

演算器数 1’892’352
言語：CUDA
Register
Shared Memory

演算器の数＝115,200
プログラミング言語： OpenCL

L1 CacheL2 CacheVRAM
HOST MEMORY

SSD
Distributed FS (Hard Disc)

Green500 4位 958MFlops/W
Top500 5位 1192TFlops

Green500 3位 1376MFlops/W(47TFlops)

浮動小数演算器の数＝4,386,816
プログラミング言語： FORTRAN

Top500/Green500 as of June 2011

Green500 6位 825MFlops/W
Top500 1位 8162TFlops

HP Project Moonshot

2800枚ものARM(携帯電話向けCPU)を組み込むことで、従来のサーバーシステムと比べ、最大で89％のエネルギー節約、94％の設置面積縮小、63％のコスト削減を目指す。

今までも、これからも• より並列度の高い計算機を使っていく• 新しいプログラミング言語を覚えさせられる• 新しい並列プログラミングのスタイルを覚えさせられる

The Problem• 僕らが書くコードは美しくて長い• 美しい繰り返しパターン• 結晶質シリケイトのような美• このような美しいものをたくさん読み書きできるのは一種の幸福、時代性的快感でさえある• それはプログラマーの求めるべき幸福ではない8check!

☞

http://paraiso-lang.org/wiki/

Paraiso計画
偏微分方程式の陽解法のためのコード生成＆自動チューニングライブラリどのような計算が(Paraisoの)対象か？

一般相対性理論＇時空・座標（流体力学＇気体・液体（磁気流体力学＇プラズマ（輻射輸送＇光の放出・反射・吸収・伝搬（• 宇宙物理の問題のうち「双曲型・楕円型偏微分方程式」に分類され、「陽解法」に属するアルゴリズムで解くことができるもの。
• まあ、空間をマス目に切って、隣同士と情報をやりとりする手法。
• 計算機の並列度があがっても(weak scalingは)やりやすい

ある分野の計算に特化した自動コード生成や自動チューニング
Problem

高速フーリエ変換

行列計算
信号処理
偏微分方程式の陽的解法

Code Generator & Automated Tuning

FFTW
ATLAS
SPIRAL
Paraiso偏微分方程式の陽解法＇の一様メッシュ部分（は少数の基本操作の組み合わせで書ける５種類、９つくらいImm即値をロードする
LoadIndex,LoadSize
配列の添字やサイズを取得
Load ＇グラフの始点（
名前付き変数からロードする
Store ＇グラフの終点（
名前付き変数に書きこむ
Reduce配列をスカラー値にBroadcast
スカラー値を配列にShift
配列全体を少しずらす
Arith
演算を行う

arith:配列の同じ位置の要素どうしの演算を行う
c <- arith add <- a,b
do j = 1, N-5
do i = 1, M-1c(i,j) = a(i,j)+b(i,j)
end doend dob <- arith sin <- ado j = 1, N-5
do i = 1, M-1b(i,j) =sin(a(i,j))
end doend dod <- arith select <- a,b,cdo j = 1, N-5
do i = 1, M-1
if(a(i,j))thend(i,j)=b(i,j)
else
d(i,j)=c(i,j)
end ifend doend doshift:配列に入っているデータの位置をずらす
b <- shift (1,5) <- ado j = 1, N-5
do i = 1, M-1b(i+1,j+5) = a(i,j)
end doend doreduce:配列をひとつの値に変える
broadcast:値を配列にばらまくb <- reduce MIN <- ab = a(0,0)
do j = 1, N
do i = 1, M
b = min(b,a(i,j))
end doend doc <- broadcast <- b

do j = 1, N
do i = 1, M
c(i,j) = b
end doend doimm:定数値を設定
loadIndex:配列の添え字を取得
loadSize :配列のサイズを取得
a <- imm 4.2
do j = 1, N
do i = 1, M
a(i,j) = 4.2
end doend dob <- loadIndexc <- loadSize

do j = 1, N
do i = 1, M
b_0(i,j) = ib_1(i,j) = j
end doend doc_0 = M
c_1 = Nload:名前つき変数から読み込み
save:名前つき変数へ書き込み
a <- load “density”

save “density” <- b

do j = 1, N
do i = 1, M
a(i,j) = density(i,j)
end doend dodo j = 1, N
do i = 1, M
density(i,j) = b(i,j)
end doend doOrthotope Machine• 多次元配列状にならんだレジスタを持つ巨大なベクトル計算機• 実数ステンシル計算に対応する仮想マシン• 演算命令は基本的に全セルに並列に作用• 隣のセルからロードする命令なども実行するための仮想マシンではなく、データフローグラフを構築するための仮想マシン

データフローグラフの例

NValue1024573689235

Shift(-1,0)689129global value＇scalar value（Add5111517
local value＇Array（Reduce(Min)2local value＇Array（22222222224
18

NInst

Load(“hoge”)

Broadcast

Mul

1022Store(“hoge”)

Paraisoの全貌

基礎方程式
さしあたり離散化方程式記述言語：
Discrete PDE Language

人手アルゴリズムの記述
OMビルダー

多次元配列仮想マシン：
Orthotope Machine (OM)

OM上のコード
OMコンパイラ

実マシン上のコード
既存コンパイラ結果

実マシン上の実行ファイル

• Paraisoプログラムはどうやって書くのか？

プログラミング言語Paraisoには通常の意味でのフロントエンドがない• Paraiso言語のソースは文字列じゃない• LexerやParserはない• Paraisoは関数型言語Haskellの中に埋め込まれたライブラリ• Builderモナドとそのコンビネーターでプログラムを記述する。

Builder同士の演算の定義とは？
Builder A

Builder B

Builder C

プログラムB

プログラムC

値B

値C

生成
プログラムA

計算
値A「Builder A と builder B の和であるBuilder 」とは、それが生成するプログラムが計算する値が、「Builder Aの生成するプログラムが計算する値」と「Builder Bの生成するプログラムが計算する値」の和になるようなものである。

typelevel-tensorNotation in Haskell
using typelevel-tensor

Einstein’s notation

Cik  Aij B jk

notation in standard
mathematics terminology

Implementation in C++

3

Cik   Aij B jkj 124実際に使っているところただの数式にみえるが、各項はBuilderモナドであり、全体がOMグラフのジェネレータになっている

「流体っぽいもの」型クラスを定義• 必要そうな物理量の定義を全部用意• あとでDead Code Eliminationが消すから大丈夫「流体っぽいもの」をApplicativeにする• 結構たくさんある流体変数全体に一つの演算を施せるように！

隣り合う4マスを補間して間の量を求める関数• これ１つで、無数の流体変数全体を一気に処理• 任意の次元、任意の方向に対応！
• 一発で書ける

4つの解の候補のなかから場合分けに応じて正しいものを選ぶ• これ１つで、無数の流体変数全体を一気に処理• 任意の次元、任意の方向に対応！
• 一発で書ける各方向ごとの計算結果を足し合わせ全体の解を求める処理• これ１つで無数の流体変数全体を＇ｒｙ• 任意の次元、任意の方向に＇ｒｙ• モナド、Fold、演算子の部分適用などすごい
Haskellの楽しい機能を駆使• 自分で後からみても正直読めない• でもこんなに少ない行数で書ける！

Don’t Repeat Yourself• Paraisoには文字列フロントエンドがない• コード生成器Builder自体が言語の第一級の対象• 関数型言語の強力な利点！
• コード生成器を自由に操れる• DRY(同じことは２度書かない)原則をとことん追求できるベンチマーク結果
速度3.50E+06

ParaisoAthena↓

Paraiso GTX460 Paraiso Tesla 2050
(単精度)

Athena Core i7
920 x8core
2.67GHz

3.00E+06
2.50E+06

2.00E+06
1.50E+06
1.00E+06

5.00E+05
0.00E+00
Paraiso Core i7
920 x8core
2.67GHzベンチマーク結果
速度3500000

ParaisoAthena↓Paraiso GTX460 Paraiso Tesla 2050
(単精度)

Athena Core i7
920 x8core
2.67GHz

3000000
2500000

2000000
1500000
1000000

5000000
Paraiso Core i7
920 x8core
2.67GHz

Athena : 惑星業界をはじめ広く使われている既存のコード。
同じ条件の2次元、流体HLLC、2次精度に設定して比較。

残念
自動生成しただけではだめで自動チューニングが必要

解きたい式：

基礎方程式
さしあたり

流体計算記述言語：
Discrete PDE Language

人手アルゴリズムの記述
OMビルダー

多次元配列仮想マシン：
Orthotope Machine (OM)
からのNativeコード生成

OM上のコード
OMコンパイラ

実マシン上のコード
既存コンパイラ結果

実マシン上の実行ファイルコード生成部
OM グラフ解析と最適化
解析・最適化済みOM

OMTrans

PlanPlanTrans

ClarisClarisTrans
Native CodeAnalysis = 解析して注釈(Annotation)をグラフにつける
Optimization 注釈にもとづきグラフを変換
Plan = コード生成戦略の決定• どのくらいメモリを使うか？
• どこまでの計算を１つのサブルーチンに収めるか？
Claris• C++やCUDAをカバーする 抽
象構文Annotationの一例：計算結果をメモリに保持するか、捨てて再計算するか？
• メモリを食ってでも計算 • 計算を多くしてでもメモリを減らしたほうが速い？と帯域を節約した方がいい？
for(;;){
f[i] = calc_f(a[i], a[i+1]);}for (;;){
b[i] += f[i] – f[i-1];}

afb

演算量 ～ N
メモリ ～3N

for(;;){
f0 = calc_f(a[i-1], a[i]);f1 = calc_f(a[i], a[i+1]);b[i] += f1 – f0;}

afb

演算量 ～ 2N
メモリ ～ 2N

例：Paraisoで実装された流体計算• グラフのノード数＝3958個• レイアウトを自由に選べるValueは=1908個• 可能な実装は、それだけでも19082 通り=2318631474140359897594479094137816650163390396354617107978538972914676911296289889528
94988789846447793390988399384716551223336856806783982602912691606248364445770172335039
54535729241917880311363490383137914861274921255128950712734788397408670521950919714209
8322292697917713518111953435214333990623513447221563209222201346475070934362866728885394848451529803078779559205459073953255482226948670514566096452159327589352442445790848
16176470059329340736642337222850662358951938698298215645717772808920891115086440342006
478637177469672403326343875446350241918444483542305006944256通り

2011/09/08 00:08

ベンチマーク結果(改)速度25000000

Paraiso+Nushio

Athena↓20000000
15000000
10000000

5000000
0
Paraiso Core i7
920 x8core
2.67GHz

Paraiso GTX460 Paraiso Tesla 2050 Athena Core i7
920 x8core(単精度)
2.67GHzもっといろんなパターンを試したい！

全自動チューニング全体像
コンピュータマネージャー

ジョブ

一定時間ごとに、終了しているジョブの成績を回収
いままでの全ての遺伝子の成績をもとに新種を生成してジョブを作り投入する

１つのジョブは１つの遺伝子をもとにプログラムを生成してコンパイル、実行

１秒ごとに１個のジョブを昼夜休みなく投げる怪しい人

プログラムを実行
正しく計算できているかを検証したうえでベンチマークを記録する

新種を作りだす三つの方法• mutation
ATATATAAATTATATATATAAAAAAAAAAAAT↓ATATAGCAATTATATCTATAAAAAGTGAAAAT• cross
ATATATAAATTATATATATAAAAAAAAAAAATGGCCGCGCCCCGCGCGCCCGCGCGCCCGGCGG↓ATATGCGAATTATATATACGCGCGCCCGGCGT• triangulation
ATATATAAATTATATATATAAAAAAAAAAAATATATATAAATTATATATATAAAAAAGTTAAATATATAGCAATTATATCTATAAAAAAAAAAAAT↓ATATAGCAATTATATCTATAAAAAAGTTAAAT

実験１

速度

手動→
チューニングの最高記録

試したプログラムの数

速度

実験２

手動→
チューニングの最高記録

試したプログラムの数トータルで• Paraiso自体のソースコード：約4300行• Paraisoで書かれた流体計算：464行↓から↓• 15万0653通りのプログラム• 5億3353万3650行• 24.726464400ギガバイトのソースコードを生成した

一番速いコードの性能

M2050のピーク性能• 43.8GFlops• 128.2GB/s• 515.2GFlops(倍精度)• 144.8GB/s

帯域リミットであるベンチマーク結果(改改)速度
40000000
35000000

Paraiso+Nushio+Genome↓

Athena↓30000000
25000000

20000000
15000000
10000000

5000000
0

Paraiso Corei7 920 x8core
2.67GHz

ParaisoGTX460(単精度)

Paraiso Tesla Paraiso Tesla Athena Core i72050
2050
920 x8core
2.67GHz統計
Manifest

Hardwaresize of .cu file numberof CUDA
kernels

memory speed
consumpt (mesh/s)
ionAthena

Core i7 x8倍精度

noneCore i7 x8倍精度

13108 lines

52 x N

0.38×106

手動チュー Core i7 x8
ニング・２倍精度

2978 lines

68 x N

2.48×106

none

GTX 460

13108 lines

752 x N

3.03×106

手動チュー GTX 460
ニング・２単精度

2978 lines1168 x N 23.37×106

自動チュー Tesla
ニング
M2050
倍精度

3103 lines

1573 x N 35.30×106

2.90×106

Manifest

Hardwaresize of .cu file numberof CUDA
kernelsnone

GTX 460

13108 lines

752 x N

手動チュー GTX 460
ニング・１単精度

3417 lines

1584 x N 22.38×106

手動チュー GTX 460
ニング・２単精度

2978 lines1168 x N 23.37×106

手動チュー GTX 460
ニング・３単精度

17462 lines1268 x N

手動チュー Tesla
ニング・２M2050
倍精度

2978 lines1168 x N 16.97×106

自動チュー Tesla
ニングM2050
倍精度

3103 lines

1573 x N 35.30×106

AthenaCore i7 x8
倍精度

memory speed
consumpt (mesh/s)
ion

3.03×106

0.68×106

2.90×106Paraisoの現状まとめ• Paraisoソースから、OpenMPあるいはCUDAにより並列化されたコードを出力できる。
• OpenMP版は、広く使われているAthenaコードに匹敵する速度。
• CUDA版は一桁上の速度が出る。しかもCPU版からソースを全く書き換える必要がない。
• １，２行のAnnotationを追加するだけでメモリの使い方やサブルーチンの分割の仕方を、したがってコードの性能をがらりと変えられる。

Paraisoの現状まとめ• メモリの使い方やサブルーチンの分け方、同期のタイミングなどを、ベンチマークと遺伝的アルゴリズムで自動チューニングしてくれる• 人間がチューンしたコードをぶちこむとさらに2倍近く速くしてくれる• これからMPIを使った分散計算に対応していきたい• OpenCLやFortranのバックエンドも作りたい