非同期的な数学的
アルゴリズムのソフトウェアの
可能性
科研費基盤（B）
O(1億)コア環境における
スケーラブルな数値計算ソフトウェアの理論と応用
*Toshiyuki Imamura (今村俊幸),
RIKEN Advanced Institute for Computational Science
研究目的
数万から数億のコアプロセッサが搭載される計算システム環境下に
おいて、過去に蓄積された高性能な数値計算サービスや新しい数学
原理に基づく方法を早期に実現するには、様々な粒度での強or 弱ス
ケーリングを実現する基本カーネルの整備が必要である。様々な並
列階層で取扱いできる軽量数値カーネル構築と共に非同期的な数値
計算アルゴリズムを組み合わせることが重要である。
1.大規模非同期的数値計算アルゴリズムの理論
2.超メニイコアでのスケーラブルな軽量コード生成のための自動
チューニング・通信同期回避技術などの核基盤技術の研究
次世代数値計算ソフトウェアの新技術創出に繋がる新機軸探究を進
める。
25.Dec.2015 ATTA2015@山上会館 2サブテーマと体制
①．非同期的アルゴ リズムの構成方法と 数理的性質 ②．非同期的アルゴ リズムを計算機上で 表現する手法 ③．小規模・中規模 粒度での低オーバ ヘッド技術（通信・同 期削減技術など） ④．メニイコア環境で 実現するためのソフ トウェア周辺技術 ATTA2015@山上会館 3 25.Dec.2015きっかけとなった基礎技術
分析と可能性
CA = Communication Avoidance
AS = ASynchronous algorithm
CA: Communication Avoidance
•
There are many studies
• CAQR by L. Grigori, “Avoiding Communication in Linear Algebra”, J.
Demmel, “Introduction to Communication-Avoiding Algorithms, SC11 tutorial, 2011, et al.
• MPK by M. Hoemmen, “Communication Avoiding Krylov Subspace
Methods”, Ph.D dissertation, UCB, 2010.
• 2.5/3D matirix-matrix multiplication, and CALU and so on
• CA-CG by E. Carson et al. “An Efficient Deflation Technique For
the Communication-Avoiding Conjugate Gradient Method”, ETNA, 2014.
• CA-SBR by G. Ballard, “Avoiding Communication in Successive
Band Reduction”, UCB/EECS-2013-131.
• and so forth..
ATTA2015@山上会館 5
CAQR
•
Algorithm: TS-QR
ATTA2015@山上会館 6
In case of p=4 MPI_Allreduce or p-t-p comms are needed. If binary-tree merge is applied,
∝ log p(when with p processes)
Conventional Householder QR needs b log p allreduce’s (b is n. of vectors)
7
加算器の並列化（キャリー予測）
(2/2)
加算器数 (3/2) n，時間 (1/2) nt + α （t は1ビット全加算器のレイテンシ，α はセレクタのレイテンシ） 再帰的並列化も可能． S3 S4 B1 A1 S1 0 B2 A2 S2 B3 A3 A4 B4 0 B3 A3 A4 B4 1 C ：セレクタ8
並列加算器の解釈
y = f(x), z = g(y)
有限個のサンプル解 g(y(1)), g(y(2)), …, g(y(q)) と， x によって，サンプ
ル解以外で g を使わずに真解 z を表現．
【加算器】z = (1 - x) g(0) + x g(1)
f(x) , g(y(1)), g(y(2)), …, g(y(q)) は依存関係なし．
アルゴリズム
1. f(x), g(y(1)), g(y(2)), …, g(y(q)) を計算．（同期不要・演算量増大）
2. z = h(g(y(1)), g(y(2)), …, g(y(q)), y) を（低コストに）計算．
g が線型ならば， z は適切に選ばれた g(x(1)), g(x(2)), …, g(x(q)) の線型
9
Chow の非同期不完全 LU 分解
未知数を反復法によって計 算． 各未知数の計算を並列実 行．依存対象変数が更新さ れていればなければ，更新 前の値を使う． 古典的不完全LU 分解(S is a set of the indices of the nonzero elements.) for i = 2, …, n do
for k = 1, …, i - 1 and (i, k) is in S do
Compute ai,k= ai,k / ak,k
for j = k + 1, …, n and (i, j) が非ゼロ do Compute ai,j= ai,j - ai,kak,j
Chow の非同期不完全 LU 分解
(S is a set of the indices of the nonzero elements.) for p = 1, 2, … until convergence do
parallel for (i, j) is in S do if i > j then Compute if i > j then Compute (i, j) 要素の変数は計算結果 は，同じ行の左側，同じ列 の上側の要素の結果に依 存． 頻繁な同期が必要．
10
Chow の不完全 LU 分解の一反復の解釈
y = [y1, y2, ..., yq]T = f(x), z = g(y)
【Chow の ILU の一反復】依存対象の変数の更新後の値（yk ）が未だ
計算されていない場合，更新後の変数（yk ）を使う計算 g の代わり に，更新前の値を使う（変数 yk に依存しない）g の近似アルゴリズ ムによって z の近似値を求めている考える． あらかじめ g の近似 g(1), g(2), ... （g よりも少ない個数の変数に依存） が存在するように f, g を設計し，並列実行時の変数の更新状況によっ て，実際に使用するアルゴリズムを（自動的に）選択する． 類似アルゴリズム Jacobi 法，Gauss-Seidel 法（G-S 法）の中間的アルゴリズム 古典 Gram-Schmidt （CGS），修正 Gram-Schmidt （MGS）の中 間的アルゴリズム．
11
非同期アルゴリズムの二つの類型
サンプル解による真解の表現 （例：並列可算機等） サンプル解の数に比例する演算量の増加． 真解をサンプル解によって厳密に表現すれば，厳密に解を求める 非同期アルゴリズムが得られる． 複数の近似アルゴリズム候補からの選択 （例：Chow の ILU，Jacobi/G-S 法の中間，CGS/MGS の中間） 複数の候補アルゴリズムが存在． 並列実行時の依存対象変数の更新状況によってアルゴリズムが決 定．従来法（Householder QR）の問題点
• 列ベクトル全体をまとめて処理
⇒全プロセスでの同期（内積等のための集団通信）
• 一列ごとに上三角化
⇒同じ処理を
O(列数）だけ繰り返す
全体：（全プロセスの同期）×O(列数)
25.Dec.2015 ATTA2015@山上会館 12まとめると・・・
全プロセスでの同期 × O(列数）回 全プロセスでの同期 × O(列数/ブロック幅）回 一部のプロセスのみの同期 × O(列数/ブロック幅）回 一部のプロセスのみの同期 × O(列数）回 いわゆる タイルQR 25.Dec.2015 ATTA2015@山上会館 13その他、応用分野からの分析、
実装など
Grid-Occupancy（性能モデル）
理想的なスレッドブロック数 (SMs × Max-TBs-per-SMの倍数)と，起
動したスレッドブロック数 (Launched-TBs) の関係を，以下のように定
義：
20 GPU Launched-TBs = 104 Grid-Occupancy =104/ceil(104,104) = 100% SM SM SM SM … GPU Launched-TBs = 105 Grid-Occupancy = 105/ceil(105,104) = 50% … SM SM SM SMGrid-Occupancy =
Launched-TBs
ceil (Launched-TBs, SMs
×
Max-TBs-per-SM)
Launched-TBs Matrix M … 行数Mに比例して スレッドブロックが起動 GPU上の複数のSMが スレッドブロックを実行
※ceil(x, y): xを最も近いyの倍数に切り上げ
NVIDIA GPUにおけるメモリ律速なBLASカーネルのスレッド
数自動選択手法（2/9）
Warp, Block, Gridの3レベルのOccupancyの導入
[SWoPP2015]
Warp-Occupancy カーネルが達成可能な1マルチプロセッサにおけるワープレベルの 並列度
Block-Occupancy カーネルが達成可能な1マルチプロセッサにおけるスレッドブロッ クレベルの並列度
Grid-Occupancy グリッドが達成可能なデバイス上の全マルチプロセッサにおける スレッドブロックの占有率
3つのOccupancyに基づくスレッドブロックサイズ決定
1.
スレッドブロック：レイテンシ隠蔽に必要なWarp-OccupancyとBlock-Occupancyを達成可能なスレッド数2.
グリッド：Grid-Occupancyが100%になるスレッドブロック数
サンプリング実行＋自動チューニングによる最小Warp/Block-Occupancy
等の探索
（詳細は[SWoPP2015]） 21NVIDIA GPUにおけるメモリ律速なBLASカーネルのスレッド
数自動選択手法（5/9）
25.Dec.2015 ATTA2015@山上会館
これまでの研究（GPU BLASのスレッ
ド数自動選択）の発展
Maxwell以降はあまり価値がなくなりそう … 命令レイテンシ削減によりOccupancyを上 げる必要がなくなった スレッドブロックサイズ選択が性能に与え る影響が少なくなった＆適当な実装でもそ こそこ性能が出る
GEMM, SpMV等への応用 スレッドブロックサイズ以外のパラメータ選 択が重要になることが多い
マルチGPUへの拡張（進行中）
OpenCLへの適用？ 22今後の計画（1/3）
25.Dec.2015 ATTA2015@山上会館
新しいテーマ：マルチGPUの活用
マルチGPU化するメリットがない（コアが余ってしまうような）小行列の計算を高速化したい
2.5Dアルゴリズム風の実装が可能か？
分散並列における2.5D行列積
[Solomonik & Demmel 11]
3次元(ijk)のプロセスグリッドにおいて，k方向に行列 データをc個冗長に持って計算
k=0〜c-1の各2次元（ij）プロセスグリッド上で行列積 を1/cだけ計算し最後にk方向に結果を足し合わせ て計算完了
通信量・通信回数が減らせる
マルチGPUにおける2.5Dアルゴリズム
行列積でij方向の分割では並列数が足りない場合 にk方向をGPU数で分割？（現在検討中） 1次元スレッドブロックより2次元スレッドブロックの方が Occupancyを稼ぎやすいのと似た発想 23今後の計画（2/3）
0 (p/c)1/2-1 (p/c)1/2-1 c-1 k i 0 0 j 25.Dec.2015 ATTA2015@山上会館研究グループ活動まとめ
H27開始年度
まとめ
•非同期原理につながる基本原理の再考
• Communication Avoidance • Synchronous Reducing • Synchronous algorithms •注目の実アプリ積分スキーム
• Parallel-in-Space/Time •非同期ソフト構築のために
• GPGPU/MICでのLightweightカーネル • パラメタ選択戦略性能チューニング 25.Dec.2015 ATTA2015@山上会館 25今後の課題：AT x {CA, SR and AS}=?
•
What are principles of and relations among CA, SR and
AS ?
•
How do we organize algorithms?
• Mathematical features • Behavior analysis
•
How do we “map and describe” them on the
modern-computer architectures and “link with AT” ?
• Programability • Performance modeling • Debugging • Tracability • Reducibility ATTA2015@山上会館 26 25.Dec.2015 1st FY : 調査検討 2nd FY: 少し具体化（プロタイプ） 3rd FY: 具体化から実機上実験
