ScaleGraph

Size: px

Start display at page:

Download "ScaleGraph"

くうしょういちぬの
5 years ago
Views:

1 超大規模半正定値計画問題に対する高性能汎用ソルバの開発と評価

2 数理計画問題 ( 最適化問題 ) と 2015 年予想 ( 目標 ) 非常に応用が広範 ( 企業社会公共政策 ) 高性能なソルバーを作ること自体が最適化問題センサーデータによる最適化問題の複雑 & 巨大化半正定計画問題 (SDP) と混合整数計画問題 (MIP) が 2 大注目数理計画問題汎用ソルバーの必要性 ( 個別の問題に対する仮定やチューニングは効果が低い ) 入力は疎データと密データの混合疎データ : 多数の CPU コアによる処理 ( 浮動小数点演算性能に非依存 : 密データへ変換 ) GPU 系による処理密データ :GPU 系による処理理論的性能限界等からボトルネック箇所を特定数値演算能力とメモリバンド等のトレードオフ関係を把握計算量とデータ移動量の正確な推定疎性やサイズなどのデータ特性と性能値の見極め計算量とポストペタスパコンでの想定 : 計算量 O(n 3 ) 数百万程度, O(n log n) 100 億以上, O(n) 100 兆以上

3 大計算量小上層 :NP 困難な最適化アルゴリズムの実行 1. MIP( 混合整数計画問題 ) の場合 : 0-1 整数変数の数 = n で計算量はおよそ O(2 n ) 2. 前処理による変数の削除と並列計算の適用 (CPU コア中心の大規模スレッド並列 : 分枝カット法の適用 ) 3. データ量は 10 5 以下 ( 整数変数の数 ) 4. 施設配置問題集合被覆 ( 分割 ) 問題スケジューリング問題などの最適化問題中層 : 多項式時間最適化アルゴリズムの実行 1. SDP( 半正定値計画問題 ) の場合 : n = 行列の大きさ, m = 制約式の数で計算量はおよそ O(n 3 + m 3 ) 2. 疎性の追求と前処理さらに並列計算の適用 (CPU コア中心の大規模スレッド並列が中心だが今後は CPU + GPU による高速化 ) 3. データ量は 10 8 以下 ( 非負変数の数 ) : 10 6 以下 ( 制約条件数 ) 4. グラフ分割センサーネットワークサポートベクターマシンなどの最適化問題下層 : グラフ解析アルゴリズムの実行を想定 1. ダイクストラ法 (1 対全最短路問題 : ヒープ付き ) の場合 : n = 点数, m = 枝数で計算量はおよそ O((n + m)log n) 2. グラフ探索の局所的な評価では優先キュー ( ヒープ木 ) を用いる実行時間メモリ消費量が安定的 (CPU コア中心の大規模スレッド並列が中心だが高速ストレージ技術による超大規模グラフ処理 ) 3. データ量は 10 12~14 以下 ( グラフの点数と枝数 ) 4. 最短路計算ネットワーク内での各点の重要度を推定各点の周辺及び広域内における影響 ( 情報の伝播力 ) を計算小データ量大

4 応用分野交通ネットワークソーシャルネットワークサイバーセキュリティバイオインフォマティクス脳科学防災計画策定 Twitter 事象同士の関係 (Relationships) ソーシャルネットワーク 6,160 万点 & 14 億 7 千万枝全米道路ネットワーク Step1 グラフ Step2 グラフ解析 Step3 分析と理解解析結果並列グラフ探索 ( 幅優先探索 ) 最適化 ( 最短路, 最大フロー, 最小費用フロー ) クラスタリング ( グラフ分割, コミュニティ抽出 ) ニューラル Human Brain Project サイバーセキュリティ 2,400 万点 & 5800 万枝 150 億 / 日のアクセスログ 890 億点 & 100 兆枝

5 l o g 2 ( m ) Human Brain Project Symbolic Network Graph500 (Medium) 京スパコン : 65536ノード Graph500: 5524GTEPS Graph500 (Large) Graph500 (Huge) Graph500 (Mini) Graph500 (Toy) Graph500 (Small) USA-road-d.USA.gr USA-road-d.LKS.gr USA-road-d.NY.gr log 2 (n) Android タブレット Tegra3 1.7GHz 4 コア : 1GB RAM 0.15GTEPS: 64.12MTEPS/W

6 A: 0.57, B: 0.19 C: 0.19, D: 0.05 TEPS(traversed edges per second) concurrent search, optimization (single source shortest path) edge-oriented (maximal independent set)

7 Follow-ship network 2009 User j User i (i, j)-edge

8 The Green Graph500 list on Nov Measures power-efficient using TEPS/W ratio Results on various system such as TSUBAME- KFC Cluster and Android mobiles Big Data category TSUBAME-KFC 6.72 MTEPS/W (44.01 GTEPS) Small Data category SONY Xperia-A-SO-04E 153 MTEPS/W (0.48 GTEPS)

内点法アルゴリズムの並列計算による超大規模半正定値計画問題の解決藤澤CREST + 遠藤 CREST

量子化学,制御分野など非常に幅広い応用を持っている高速かつ安定した反復解法である内点法アルゴリズムが存在している

内点法アルゴリズムには２つの大きなボトルネックが存在大規模な並列化 ELEMENTS 線形方程式系の行列 SCM)の生成

ノード内の NUMA 関係の情報を読み取り自動的に最適な processor affinity や memory

0 1360ノード(16,320 CPU cores) での高いスケーラビリティ CHOLESKY

10 内点法アルゴリズムの並列計算による超大規模半正定値計画問題の解決藤澤CREST + 遠藤 CREST 半正定値計画問題(SDP)は現在最も注目されている数理最適化問題の一つ組合せ最適化データマイニング量子化学,制御分野など非常に幅広い応用を持っている高速かつ安定した反復解法である内点法アルゴリズムが存在している SDPARA は現在開発公開を行なっている大規模な SDP に対する並列ソルバー内点法アルゴリズムには２つの大きなボトルネックが存在大規模な並列化 ELEMENTS 線形方程式系の行列 SCM)の生成 : メモリバンド幅依存 CHOLESKY SCMに対するコレスキー分解 : 倍精度浮動小数点演算依存 ELEMENTS ノード内の NUMA 関係の情報を読み取り自動的に最適な processor affinity や memory interleaving の設定を行う MPI + OpenMP による２段階並列計算東工大 TSUBAME ノード(16,320 CPU cores) での高いスケーラビリティ CHOLESKY 通信と計算のオーバーラップ多数 GPU による並列計算東工大 TSUBAME GPUs(NVIDIA K20X) での大規模分散並列化浮動小数点演算1.713 PFlops の達成と世界最大規模の SDP(233万制約超)を初めて解くことに成功した

11 SDPが注目される理由とは主双対内点法などのアルゴリズムによって多項式時間で最適解を求めることができる(つまり高速で安定したアルゴリズムが存在する) SDP は線形計画問題(LP), 凸二次計画問題や二次錐計画問題(SOCP) などを含んだより大きな凸計画問題の枠組である非凸最適化問題に対する強力な緩和値を導き出すことができる. そのため SDP を繰り返して解くことによって非凸最適化問題(例えば双線形行列方程式(BMI) など) を扱える可能性を持っている組合せ最適化問題, 整数計画問題, ノルムなどを用いた配置問題, システムと制御, ロバスト最適化, 量子化学など非常に多くのSDP の応用が存在する (つまり非常に多彩な応用分野を持っている). 多くのSDP に対するソフトウェアが開発され, インターネットより公開されている. つまり公開されているソフトウェアで実際に大きな問題を解くことができる２１世紀の線形計画問題として大きな期待を受けている

12 構造最適化問題量子化学システムと制御組合せ最適化問題データマイニング

13 SDP(SemiDefinite Programming) Primal Dual SDPA 1.x(1995) could solve SDPs with n = m = 20 ~30 (very small size!!)

14 Primal-Dual Interior-Point Methods

15 疎で大規模な SDP においては Schur complement equation がネック Bdz = g B ij = A i èx k A j (Y k ) Ä1 (i; j = 1;... ; m) Bは密行列と疎行列の両方に成り得る上記の連立線形方程式を非常に効率良く解く必要がある Bdz = g はコレスキー分解 Cholesky factorization や反復法 Conjugate gradient で解くのが一般的

16 Exploitation of Sparsity in Schur complement matrix The SDPA automatically selects a formula when computing each row of B by estimating computational costs of F1, F2 and F3 F1 F2 F3

The major bottleneck parts (80% - 90% of total execution time) ELEMENTS : Computation of the SCM Memory Access-intensive Time-complexity: or CHOLESKY : Cholesky factorizations of the SCM

17 The major bottleneck parts (80% - 90% of total execution time) ELEMENTS : Computation of the SCM Memory Access-intensive Time-complexity: or CHOLESKY : Cholesky factorizations of the SCM Compute-intensive Time-complexity: n : matrix size m: # of constraints ELEMENTS-bound SDP problems m < n (not m >> n), and Fully Dense SCM CHOLESKY-bound SDP problems m >> n, and Fully Dense SCM torus quantum Comb. Opt. Quad. assignment Fast computation of sparse SCM is future work. e.g.) The sensor network location problem and The polynomial optimization problem

18 ULIBC を用いた NUMA を考慮した高速化 ELEMENTS : Computation of the SCM NUMA アーキテクチャ Non-uniform Memory Access CPU ソケットとローカルメモリの対により構成されるためコアとメモリの距離が一定ではなくアクセスコストが不均一メモリアクセスの局所性の向上により高い性能を示す先行研究 NUMA の考慮にはこの組合せが重要だが存在しない計算機トポロジを検出するためのライブラリ群 CPUアフィニティ & ローカルメモリ確保のためのライブラリ群メモリアクセスの局所性を高めるためのスケジューリング方法提案高速化に重要な機能を ULIBC ライブラリで実現トポロジ検出 + CPUアフィニティ & ローカルメモリ確保 C 言語のみで実装関数で提供しており汎用的な実装が可能 4 種類の適用例で高速化を達成性能評価を行う

19 NUMA アーキテクチャの構成 4-way Intel Xeon E (Sandybridge-EP) 4 (CPU ソケット数 ) 8 (CPU ソケットあたりの物理コア数 ) 2 ( 物理コアあたりのスレッド数 ) 最大 4 x 8 x 2 = 64 スレッド ( ハイパースレッディング ) NUMA ノードローカルメモリへのアクセス ( アクセスコスト小 ) NUMA ノード CPU ソケット (16 論理コア ) ローカルメモリリモートメモリへのアクセス ( アクセスコスト大 ) データアクセスが不均一であるため並列時の偏りが生じるため並列効率を向上することが難しい本来の性能を予想することが難しい

NUMA の計算機トポロジの検出 4-way Intel Xeon E5-4640 (Sandybridge-EP) 4 (CPU ソケット数 ) 8 (CPU ソケットあたりの物理コア数 ) 2 ( 物理コアあたりのスレッド数 )

ローカルメモリリモートメモリへのアクセス ( アクセスコスト大 ) 計算機トポロジの検出 Processor ID : 各論理コアに対応した ID Package ID : 各 CPUソケット ( 一部例外あり ) に対応した ID

20 NUMA の計算機トポロジの検出 4-way Intel Xeon E (Sandybridge-EP) 4 (CPU ソケット数 ) 8 (CPU ソケットあたりの物理コア数 ) 2 ( 物理コアあたりのスレッド数 ) Package ID Core ID SMT ID (0,1,2,3) (0,1,,7) (0,1) NUMA ノードローカルメモリへのアクセス ( アクセスコスト小 ) NUMA ノード CPU ソケット (16 論理コア ) ローカルメモリリモートメモリへのアクセス ( アクセスコスト大 ) 計算機トポロジの検出 Processor ID : 各論理コアに対応した ID Package ID : 各 CPUソケット ( 一部例外あり ) に対応した ID Core ID : CPU ソケット内の物理コアに対応した ID SMT ID : 物理コア内のスレッドに対応した ID Linux デバイスファイル /sys/devices/system/* /proc/cpuinfo APCID(CPUID 命令 )

21 メモリアクセスパターンを考慮したアフィニティ設定 Scatter アフィニティ使用するコアをできるだけ分散する 0, 4, 8,... 1, 5, 9,... 2, 6, 10,... 3, 7, 11,... Compact アフィニティ使用するコアをできるだけ集中する 0, 1, 2,... 16, 17, 18,... 32, 33, 33,... 48, 49, 50,... Compact+ アフィニティ Compact での優先度をソケット > 論理コア 0, 1, 2,... 8, 9, 10,... 16, 17, 18,... 32, 33, 34,...

interleaving の設定を行う MPI + OpenMP による 2 段階並列計算東工大

22 Parallel Computation for ELEMENTS ELEMENTS : ノード内の NUMA 関係の情報を読み取り自動的に最適な processor affinity や memory interleaving の設定を行う MPI + OpenMP による 2 段階並列計算東工大 TSUBAME ノード (16,320 CPU cores) での高いスケーラビリティ

Automate configuration of CPU Affinity and memory allocation policy Step1: Linux device files Step2: Processor mapping table Step3: thread Processor assignment table /sys/devices/system/* TSUBAME 2.

23 Automate configuration of CPU Affinity and memory allocation policy Step1: Linux device files Step2: Processor mapping table Step3: thread Processor assignment table /sys/devices/system/* TSUBAME 2.0 HP node Automate detection Processor ID Package ID Core ID Thread ID Scatter Compac t Local rank=0 Processor cores RAM0 0, 2, 4, 6, 8, 10 RAM1 Processor cores 1, 3, 5, 7, 9, 11 GPU1 Local rank=1 GPU2 Scatter-type Affinity distributes OpenMP threads as evenly as possible across the entire system. Compact-type Affinity binds the (n+1)-th OpenMP thread in a free-thread context as close as possible to the thread context in which the n-th OpenMP thread was bound.

24 Electronic4 : 量子化学系の問題 8604 n = 33,178 ( 行列の大きさ ) m = 76,554 ( 制約式の数 ) 14,116,348 ( 非零要素の総数 ) x 4 ブロック数 = 22 対角ブロック : 694

25 Electronic5 : 量子化学系の問題 4965 n = 19,640 ( 行列の大きさ ) m = 36,795 ( 制約式の数 ) 6,731,930 ( 非零要素の総数 ) x 4 ブロック数 = 22 対角ブロック : 490

26 CPU time for Electronic4 the memory interleaving is effective. the parameter set scatter and interleaving is fastest.

27 Scalability for Electronic4 and Electronic5 high efficiency (75.9% for Electronic4) higher efficiency when solving an SDP problem larger than Electronic4 because the efficiency for Electronic4 is higher than that for Electronic5. Strong Scaling based on 128 nodes Efficiency for large problems

0 Exa(10 18 ) Flops Our implementation is based on pdpotrf function (ScaLAPACK library) and HPL Library for TSUBAME 2.

28 CHOLESKY : Cholesky factorization of Schur complement matrix For symmetric and positive semidefinite matrix Time complexity is O(m 3 ) Ex. Matrix (1.48million x 1.48 million) 1.0 Exa(10 18 ) Flops Our implementation is based on pdpotrf function (ScaLAPACK library) and HPL Library for TSUBAME 2.0 ScaLAPACK and HPL library employ the 2-D block-cyclic distribution for Cholesky factorization In ELEMENTS, we directly computes all elements of the 2-D block-cyclic distribution of the SCM

29 Data Decomposition in Cholesky The dense matrix B(m x m) is uniformly distributed with 2D Block-Cyclic distribution with block size nb(= 2048) CHOLESKY algorithm consists of mb(= m/nb ) steps. MPI processes conceptually composes a 2D grid(ex. 6 (=2 x 3), 4080(=68 x 60) processes) Can use 4080GPUs in parallel!! m Matrix distribution on 6 (=2x3) processes nb=2048

30 ScaLAPACK の Cholesky 分解のブロック分割並列計算を行う

The k-th step proceeds as follows: We need mb(= m/nb ) steps. 1. Diagonal block factorization 2. Panel factorization 3. Panel broadcast and transposition 4.

31 The k-th step proceeds as follows: We need mb(= m/nb ) steps. 1. Diagonal block factorization 2. Panel factorization 3. Panel broadcast and transposition 4. Update: the most computation-intensive part B = B L x L We need call DGEMM function on GPU once. Matrix distribution on 6 (=2x3) processes Each process has a partial-matrix L N B L B B' B' L L' nb For Accelerating CHOLESKY by using massively parallel GPUs, we need optimization techniques to overlap computation, PCI-Express communication and MPI communication.

32 TSUBAME 2.0 & 2.5 Specification Specification CPU Intel Westmere EP (Xeon X5670, L2 Cache: 256 KB, L3: 12MB) 2.93 GHz processors, 12 CPU Cores (24 cores with Hyper Threading) x 2 sockets per 1 node (24 CPU Cores) RAM 54 GB OS SUSE Linux Enterprise 11 (Linux kernel: ) # of Total Nodes 1466 nodes (We used up to 1360 nodes) Network Topology Network GPGPU Full-Bisection Fat-Tree Topology Voltaire / Mellanox Dual-rail QDR Infiniband (40Gbps x2 = 80 Gbps) Three NVIDIA Fermi M2050 (or Kepler 20X) GPUs, CUDA C/C++ Compiler Intel icc 11.1 & gcc MPI MVAPICH 1.5.1

33 TSUBAME 2.5 System Configuration 1360 nodes, 2720 CPUs, 4080 GPUs

We have challenged to solve largest SDPs in 2013!

(matrix size) m = 2,339,331 (# of constraints) nnz = 28,078,376

computation time 1360 nodes, 2720 CPUs, 4080 GPUs (TSUBAME 2.

34 We have challenged to solve largest SDPs in 2013!! Quadratic Assignment Problem QAP(sko49) : 2013 n = 2,752,649 (matrix size) m = 2,339,331 (# of constraints) nnz = 28,078,376 (# of nonzero elements) CHOLESKY occupied over 95% of the entire computation time 1360 nodes, 2720 CPUs, 4080 GPUs (TSUBAME 2.5 in TITECH) Performance of Cholesky Fac. : 4.26 Exa flops / 2487 sec PFLOPS

transposition Update: Each process updates its own part of the rest matrix, taking

35 Parallel Algorithm of Cholesky Factorization Version 1: No overlapping Update: B = B L x Lt(Lt ) Diagonal block factorization Panel factorization Panel broadcast and transposition Update: Each process updates its own part of the rest matrix, taking the corresponding part of L and Lt(trace of L). Then B = B L x Lt(Lt ) is computed.

Parallel Algorithm of Cholesky Factorization Version 2: GPU computation and PCIe communication are overlapped Diagonal block factorization Panel factorization Panel

36 Parallel Algorithm of Cholesky Factorization Version 2: GPU computation and PCIe communication are overlapped Diagonal block factorization Panel factorization Panel broadcast and transposition Update: Each process updates its own part of the rest matrix, taking the corresponding part of L and Lt(trace of L). Then B = B L x Lt(Lt ) is computed.

Parallel Algorithm of Cholesky Factorization Version 3: GPU computation, PCIe communication, and MPI communication are overlapped Diagonal block factorization Panel factorization

37 Parallel Algorithm of Cholesky Factorization Version 3: GPU computation, PCIe communication, and MPI communication are overlapped Diagonal block factorization Panel factorization Panel broadcast and transposition Update: Each process updates its own part of the rest matrix, taking the corresponding part of L and Lt(trace of L). Then B = B L x Lt(Lt ) is computed.

factorization Panel factorization Panel broadcast and transposition Update: Each process updates its

38 Parallel Algorithm of Cholesky Factorization Version 4: GPU computation, PCIe communication, MPI communication, and panel factorization are overlapped Producer of panel L Diagonal block factorization Panel factorization Panel broadcast and transposition Update: Each process updates its own part of the rest matrix, taking the corresponding part of L and Lt(trace of L). Then B = B L x Lt(Lt ) is computed.

39 Performance of CHOLESKY on TSUBAME2.0 & 2.5 Quadratic assignment problems (QAP) org: Version 3 & new: Version 4

40 Performance of CHOLESKY on TSUBAME2.0 & 2.5 Quadratic assignment problems (QAP) 1.713PFLOPS(DP) with 4080GPUs!!

41 大規模グラフ解析 :CREST 終了時の目標と想定超大規模ネットワークに対する探索アルゴリズムとクラスタリングアルゴリズム最短路 ( 中心性 ), 最大フロー, PageRank, グラフ分割, 中心性, その他数理計画問題 (SDP, MIP) グラフクラスタリング, 高速グラフレイアウト, グラフ UI 数百万頂点数兆頂点数億枝数百兆枝からなる超大規模なグラフ解析 2 42 頂点数 PB 以上のワーキングデータセットグラフ解析 : 例 :890 億個のニューロンとその接続 100 兆枝ペタバイト級のストア領域 :(2023 年 : エクサスパコン ) 数百万人の被災者の避難経路の計算では数千万頂点のグラフ (1 スレッドあたりのメモリ要求量 1Gbytes) に対して同時に数百万スレッド単位で各被災者毎の最短路計算と各点の重要度判定 5 億頂点 250 億辺のグラフのインタラクティブな操作のための計算性能ベクトルスカラープロセッサから成る不均質でかつ分散共有メモリから成る大規模並列環境における大規模グラフ処理 : 数百万スレッドまでスケールさせる大規模グラフ処理を PGAS 言語 X10 を用いて一元的に記述メモリ多階層を考慮して高速性と省電力性を両立したアルゴリズムの確立世界最高性能の最適化ソルバーの開発 : 密データと疎データの分離と多階層並列 (MPI, OpenMP, OpenACC 等 ) 及びアクセラレータ (GPU, MIC 等 ) の活用

untitled

untitled c NUMA 1. 18 (Moore s law) 1Hz CPU 2. 1 (Register) (RAM) Level 1 (L1) L2 L3 L4 TLB (translation look-aside buffer) (OS) TLB TLB 3. NUMA NUMA (Non-uniform memory access) 819 0395 744 1 2014 10 Copyright