内容に関するご質問はida＠cc.u-tokyo.ac.jp
まで、お願いします。

第62回 お試しアカウント付き
並列プログラミング講習会「ライブラリ利用：科学技術計算の効率化入門」

階層型行列法とℋACApKライブラリ
東京大学情報基盤センター 特任准教授 伊田 明弘

1

チュートリアルの流れ
1.
2.
3.

ppOpen-APPL/BEM
境界要素法とBEM-BBフレームワーク
プログラム実習Ⅰ
・BEM-BB4.
5.

階層型行列法とℋACApKライブラリ
プログラム実習Ⅱ
・BEM-BB+ ℋACApK

2

お試しアカウント付き講習会

多次元空間でシミューレートする・Transcribed from VINAS co.LTD
・Transcribed from Ohtani at al (2011)

どうやって？

3

科学技術計算の定式化例■どんどん近似していき，どこかで観念して計算する
物理現象
数学モデル化(非線形)時間発展方程式
時間微分に（半）陰解法(非線形)境界値問題
離散化

離散化手法の選択が重要な分かれ目

非線形代数方程系
線形化、ニュートン法

線形方程式系

この辺りがスパコンの出番

低ランク行列近似

線形方程式系
4

Mesh used in FDM, FEM, BEM

FDM

BEM
Pictures are transcribed from

FEM

FDM: de Havilland Aircraft Company, VINAS Co. Ltd
FEM: Prof. Nakazima (Tokyo Univ.)
BEM: μ-TEC Co. Ltd
56

FDM, FEM, BEM for Elliptic PDE
continuous

・Bi-linear form

discrete

For 𝑥𝑥 ∈ Ω ⊂ ℝ3 , find 𝜙𝜙(𝑥𝑥) s.t.

Az = 𝑏𝑏∆𝜙𝜙(𝑥𝑥) = 𝑓𝑓(𝑥𝑥)∫Ω 𝜓𝜓(𝑥𝑥)𝛥𝛥𝜙𝜙(𝑥𝑥)d𝑥𝑥 = ∫Ω 𝜓𝜓(𝑥𝑥)𝜌𝜌 d𝑥𝑥

・Green's identity∫Ω 𝜓𝜓(𝑥𝑥)𝛥𝛥𝜙𝜙(𝑥𝑥)d𝑥𝑥 =∫𝜕𝜕Ω 𝜓𝜓 𝑥𝑥 𝜙𝜙𝑛𝑛 (𝑥𝑥) d𝑥𝑥 − ∫Ω 𝛻𝛻𝜓𝜓(𝑥𝑥) � 𝛻𝛻𝜙𝜙(𝑥𝑥)d𝑥𝑥

� 𝛻𝛻𝜓𝜓(𝑥𝑥) � 𝛻𝛻𝛻𝛻(𝑥𝑥)d𝑥𝑥 = 𝜓𝜓, 𝑓𝑓 ̅Ω・Fredholm integral equation∫ 𝜓𝜓� 𝑥𝑥, 𝑦𝑦 𝜙𝜙𝑛𝑛 (𝑥𝑥) d𝑥𝑥 = 𝑓𝑓 ̿𝜕𝜕Ω

� 𝜑𝜑(𝑦𝑦) � 𝜓𝜓� 𝑥𝑥, 𝑦𝑦 𝜙𝜙𝑛𝑛 (𝑥𝑥) d𝑥𝑥dy = 𝜑𝜑, 𝑓𝑓̃
𝜕𝜕Ω

𝜕𝜕Ω

FEM
・Base function

BEM
・Base function■6概要：ppOpen-HPC■JST-CREST
・ “ポストペタスケール高性能計算に資するシステムソフトウェア技術の創出”領域
・ “自動チューニング機構を有するアプリケーション開発・実行環境ppOpen-HPC”研究代表者: 中島 研吾(東京大学情報基盤センター)■Download site: http://ppopenhpc.cc.u-tokyo.ac.jp/

このチュートリアルで使用

User’s Program

ppOpen-HPC covers …
ppOpen-APPL

FEM

FDMii

FVM

BEM

ppOpen-MATH

MG

ii
GRAPH

VIS

MP

ppOpen-AT

STATICii
DYNAMIC

ppOpen-SYS

COMM

ii
FT

DEM

ppOpen-HPC

Optimized Application with
Optimized ppOpen-APPL, ppOpen-MATH

7

ppOpen-APPL/BEMを用いた解析例■静電場解析
・Potential operator：
・半無限解析領域において物体の表面電化を計算する.

𝐸𝐸0
Ground

解析条件

解析結果
8

ppOpen-APPL/BEMを用いた解析例■静電場解析
・Potential operator：

Air

1V

0.25m

Conductor

0.5m

Ground

解析条件
・半径0.25mの導体球に1Vを与える
・表面に現れる電荷を求める（半無限領域内の電位が導出される）

解析結果
導体球の下半分に大きな電荷密度が生じる

9ppOpen-APPL/BEMの構成■BEM－BB・境界要素法・積分方程式法用ソフトウェアフレームワーク（ライブラリ、API集ではない）
・使用法：境界要素積分関数をユーザが用意しBEM-BBに埋め込む（ユーザプログラムから必要な関数を呼ぶのではない）
・密行列演算利用版、階層型行列（H-matrix）利用版
・Fortran90, MPI+OpenMP, BLAS
■HACApK・階層型行列法ライブラリ
・使用法：ユーザプログラムから必要な関数を呼ぶ
HACApKから呼ばれる境界要素積分関数を用意しておく・Fortran90, MPI+OpenMP

10

ppOpen-APPL/BEMの構成
ppOpen-APPL/BEM

BEM-BB
密行列演算を用いた大規模並列境界要素解析フレームワーク

境界要素解析
テンプレート群

テンプレートは特定の問題に対してフレームワークと共に利用

H行列法を用いた大規模並列境界要素解析フレームワークℋACApK

分散並列H行列法ライブラリ

フレームワーク内で使用他のアプリケーションでも利用可能
11

チュートリアルの流れ
1.
2.
3.

ppOpen-APPL/BEM
境界要素法とBEM-BBフレームワーク
プログラム実習Ⅰ
・BEM-BB4.
5.

階層型行列法とℋACApKライブラリ
プログラム実習Ⅱ
・BEM-BB+ ℋACApK

12

お試しアカウント付き講習会

境界要素法（積分方程式法）• 偏微分方程式の境界値問題に対する代表的な離散化解法の一つ• モデルの境界面での離散化のみで解析を実行– 少ない未知変数で大規模なモデルを扱うことが可能– 開領域問題，形状最適化，移動物体問題に有効• 一般的には密行列計算が必要– 実応用解析では，FMMやH行列法等の高速行列演算技術の利用が不可欠• ターゲットアプリケーション– 電磁場解析– 地震サイクルシミュレーション

13境界要素積分方程式の離散化■密行列を係数行列に持つ線形方程式系が得られる.

ℊ 𝑢𝑢 = 𝑓𝑓, where ℊ 𝑢𝑢 𝑥𝑥 : = ∫Ω 𝜅𝜅 𝑥𝑥, 𝑦𝑦 𝑢𝑢(𝑦𝑦)d𝑦𝑦
singular kernel： 𝜅𝜅 𝑥𝑥, 𝑦𝑦 ∈ span𝑥𝑥 − 𝑦𝑦 −𝑝𝑝 , 𝑝𝑝 > 0e.g. 𝜅𝜅 𝑥𝑥, 𝑦𝑦 ∝ 𝑥𝑥 − 𝑦𝑦 −1

Discretization
by W.R.M., e.g.
Ritz-Galerkin

𝐴𝐴𝜙𝜙 = 𝑏𝑏,

𝐴𝐴𝑖𝑖𝑖𝑖 = � 𝜑𝜑𝑖𝑖 𝑥𝑥 � 𝑔𝑔 𝑥𝑥, 𝑦𝑦 𝜑𝜑𝑗𝑗 𝑦𝑦 d𝑦𝑦d𝑥𝑥,
Ω𝐴𝐴 : 密行列

𝜑𝜑𝑖𝑖

𝑢𝑢 ≅ 𝑢𝑢ℎ = � 𝜙𝜙𝑖𝑖 𝜑𝜑𝑖𝑖
𝑖𝑖∈ℶ

Ω

Ω

𝑏𝑏𝑖𝑖 = � 𝑓𝑓 𝑥𝑥 𝜑𝜑𝑖𝑖 (𝑥𝑥)d𝑥𝑥
Ω14境界要素解析の手順

1. モデルデータの入力
2. 係数行列の作成

BEM-BB
フレームワーク– 要素間（数値／解析）積分の実行

3. 右辺ベクトルの作成– 境界条件の適用

4. 連立一次方程式の求解– 密な係数行列を使用

5. 解析結果の出力15

BEM-BBの利用目的• 想定利用ユーザ– 主：境界要素法を使用し，大規模並列計算環境で効率的に動作するシミュレーションプログラムを開発する研究者・技術者– 副：特定のアプリケーション領域の問題を境界要素法
により解く研究者・技術者• 商用のアプリケーションソフトウェアと同様にモデルに関する情報を入力するだけで必要な解を得たい。
• PC等の逐次計算環境や小規模なクラスタ、マルチコアプロ
セッサでも動作することが望ましい– 大規模計算環境へのシームレスな移行

16

BEM-BBの設計思想• 計算科学者が自らのアイデアを反映できる余地が残されていなければならない– サポートの範囲を限定する• 要素間積分には様々な方法があり，新しい提案も今後あり得る．

– 並列処理はやはり面倒• 必要なプロセス間通信は提供されるソフトウェアが正
しく，ブラックボックス的にやってくれればうれしい

17

BEM-BBの機能とユーザの準備■Framework for parallel BEM analyses＜user＞＜Frame work＞
Model data input
Making matrix
・Call user function

Liner solver
on distributed memory
parallel computer systems

Input file・mesh, physical parameters and so on

Boundary element integral
・Function to calculate entry 𝐴𝐴𝑖𝑖𝑖𝑖

e.g. : real*8 function entry_ij(i, j, param)
integer :: i, j
type(phys) :: paramj

OpenMP/MPI ハイブリッド
並列処理の活用

i

𝑨𝑨𝒊𝒊𝒊𝒊18

BEM-BBの機能とユーザの準備■Framework for parallel BEM analyses＜user＞＜Frame work＞
Model data input

Input file・mesh, physical parameters and so on

Making matrix
・Call user function

Liner solver
on distributed memory
parallel computer systems

Templates
・specific application domains■Static electric filed

OpenMP/MPI ハイブリッド
並列処理の活用19

BEM-BB構成
bem_bb_dense_mpi
メインプログラム

Initialization_MPIMPI変数初期化

Configファイル

Read_bem_bb_config
制御データ入力

入力ファイル名ハードコーディング
デフォルト：sample.input.txt

Read_model_data
メッシュ、物理量入力
入力データ分配

入力ファイル

Make_equation_data
係数行列、右辺ベクトル生成

element_ij
行列要素計算
境界要素積分実行

element_ij_???
???用テンプレート

element_ij_SCM
表面電化法用テンプレート

ppohBEM_pbicgstab_dense
線形ソルバ
仕様を公開

ユーザが用意今後追加予定機能

matvec_direct
行列・ベクトル積のbodyPrint_result
可視化データ出力

H行列版はHACAｐKルーチン使用

出力ファイル

出力ファイル名ハードコーディング
デフォルト：sample.output.vtk

20

BEM-BBの入力データ形式■３次元空間上の多角形メッシュで表現される面を想定■入力データファイル（拡張子 .pbf）フォーマット
(i) Number of nodes (1integer*4 number): nond(ii) Coordinates of the nodes in three dimensions (3 real*8 numbers for each node)(ii) Number of faces (1integer*4 number): nofc(iv) Number of nodes on each face (1 integer*4 number): nond_on_fc(v) Number of integer parameters set on each face (1 integer*4 number): nint_para_fc(vi) Number of real parameters set on each face (1 integer*4 number): ndble_para_fc(vii) Node number constructing for each face (nond_on_fc x nofc integer*4 numbers)(viii) Integer parameters set on faces (nond_on_fc x nint_para_fc integer*4 numbers)(ix) Real parameters set on faces (nond_on_fc x ndble_para_fc real*8 numbers)
4
0.000000000000e+00 0.000000000000e+00 0.000000000000e+00
0.000000000000e+00 1.000000000000e+00 0.000000000000e+00
1.000000000000e+00 1.000000000000e+00 0.000000000000e+00
2.000000000000e+00 0.000000000000e+00 2.000000000000e+01
2
3
0
1
123
243
0.000000000000e+00
1.000000000000e+00
21

Configファイル■入力データに含まれていない制御パラメータを指定(i) number_element_dof (integer*4)(ii) linear_solver (character*16)(iii) tor (real*8)(iv) max_steps (integer*4)
1
BICGSTAB
1d-8
5000
付属configファイル

22ユーザ定義関数 element_ij■２つの要素番号と各要素上の情報を受け取り、行列要素値を返す
real(8) function matrix_element_ij(i, j, nond, nofc, nond_on_fc, np,int_para_fc,
nint_para_fc, dble_para_fc, ndble_para_fc, face2node)type :: coordinate
real(8) :: x ,y ,z
end type coordinate
integer,intent(in) :: i,j,nond,nofc,nond_on_fc,nint_para_fc,ndble_para_fc
type(coordinate), intent(in) :: np(*) ! 節点座標
integer,intent(in) :: face2node(nond_on_fc,*), ! 要素を構成する節点
int_para_fc(nint_para_fc,*) ! 要素上にintで与えられている情報
real(8),intent(in) :: dble_para_fc(dble_para_fc,*) !要素上にreal8で与えられている情報23

BEM-BBによる解析結果の可視化■オープンソース可視化ソフトウェアParaViewに対応
・VTK形式ファイルを出力

24

チュートリアルの流れ
1.
2.
3.

ppOpen-APPL/BEM
境界要素法とBEM-BBフレームワーク
プログラム実習Ⅰ
・BEM-BB4.
5.

階層型行列法とℋACApKライブラリ
プログラム実習Ⅱ
・BEM-BB+ ℋACApK

25

お試しアカウント付き講習会

サンプルプログラムで計算する問題■静電場解析
・Potential operator：

Air

1V

0.25m

Conductor

0.5m

Ground

解析条件
・半径0.25mの導体球に1Vを与える
・表面に現れる電荷を求める（半無限領域内の電位が導出される）

解析結果
導体球の下半分に大きな電荷密度
が生じる

26

実習を行うあたっての注意点(BEM-BB)
 ファイル名

ppohBEM_0.4.1.tar.gz (Oakleaf-fxのみ)
 ジョブスクリプトファイルbem-bb.bash中のキュー名をlecture から tutorial に変更してから
pjsubしてください。
 lecture : 実習時間外のキュー（同時実行数１）
 tutorial : 実習時間内のキュー（同時実行数4+）27

お試しアカウント付き講習会

BEM-BB用ジョブスクリプト(Oakleaf-fx)#!/bin/bash#PJM -L "rscgrp=lecture"#PJM -L "node=1"#PJM --mpi "proc=1"#PJM -L "elapse=15:00"#PJM -g gt00
export OMP_NUM_THREADS=1
mpirun ./bem-bb-SCM.out input_sample.pbf

利用ノード数
ＭＰＩプロセス数
OpenMP
スレッド数

入力データ名

ノード数、MPIプロセス数、OpenMPスレッド数の設定はOakleaf-fxは１ノードあたり16コアであることに注意して行う
28

お試しアカウント付き講習会

BEM-BBの実行（Oafleaf-fx）
 Oakleaf-fxにログインする 以下のコマンドを実行する$ cp /home/z30107/ppohBEM_0.4.1.tar.gz ./$ tar xvzf ppohBEM_0.4.1.tar.gz$ cd ppohBEM_0.4.1$ cd bem-bb-framework_dense$ cd src$ cd framework_with_templates$ make$ pjsub bem-bb.bash
 実行が終了したら、以下を実行する$ cat bem-bb-fx.bash.oXXXXXX (XXXXXXは数値)
29

お試しアカウント付き講習会

BEM-BBの実行（Oafleaf-fx）■以下のような結果が見えれば成功
Number of processes 1
Number of unknowns set on each face element = 1
Selected linear solver = BICGSTABConvergence criterion = 1.000000000000000E-08
Upper limit of iteration counts = 5000
Load data file :input_sample.pbf
Number of total threads, Thread number 1 0
ext_ndim = 600 , ndim = 600
Original relative residual norm = 1.0000000E+00
Step 1 relative residual = 5.3560399E-03

MPIプロセス数
入力データ名OpenMPスレッド数
要素数

（中略）
Step 11 relative residual = 1.3362906E-08
Step 12 relative residual = 2.5294315E-09
Relative residual norm = 2.5294314E-09
before user_func time = 1.17245559E-02 [sec]user_func time = 0.704002202 [sec]bicgstab time = 8.27010907E-03 [sec]bicgstab time/steps = 6.89175737E-04 [sec/times]_____________all time = 0.723996818 [sec]

計算時間

30

演習課題１(BEM-BB)
1. 解析結果可視化用ファイル「 sample.output.vtk 」
を自分のPCにダウンロードし、ParaViewを用いて下図のような可視化を行ってください。

31

演習課題２(BEM-BB)
2-1 サンプル以外のデータを用いた計算演習１で使用したデータはinput_sample.pbf(要素数：600）でした。
Input_10ts.pbf(要素数：約1000)およびInput_216h.pbf(要素数：約2,1000)のデータを用いた計算を行い、計算時間を測定してください。
・ヒント：ジョブスクリプトに書かれているデータ名を変更する
・ヒント：データ216hの計算時間は約１５分かかります。ジョブが流れたことを確認したら次の実習に進んでください

2-2 計算オーダーの計測
演習2-1で測定した計測時間を使って、密行列を用いた境界要素解析の計算時間が、要素数を𝑁𝑁として、O(𝑁𝑁 2 )であることを確かめて下さい。

32

演習課題３、４(BEM-BB)
3-1 コードの並列実行
いろいろなOpenMPスレッド数やMPIプロセス数で計算を実行し、計算時間を計測してください。
３-2 台数効果による性能向上の評価
逐次実行に比べ、どれだけ高速化されたか、性能評価をしてください。
４ 計算環境の移動
ppohBEM_0.4.1.tar.gzをReedBushへコピーし、BEM-BBを動作させてください。
・ヒント：Makeファイルの先頭数行で計算環境を設定しています
・ヒント：ジョブスクリプトは他の演習で使ったものを参考に自作
33

チュートリアルの流れ
1.
2.
3.

ppOpen-APPL/BEM
境界要素法とBEM-BBフレームワーク
プログラム実習Ⅰ
・BEM-BB4.
5.

階層型行列法とℋACApKライブラリ
プログラム実習Ⅱ
・BEM-BB+ ℋACApK

34

お試しアカウント付き講習会

階層型行列法ライブラリ ℋACApK■積分方程式を用いたシミュレーションのための高速計算ライブラリ► 行列近似手法
階層型行列法 ：𝑶𝑶(𝑵𝑵𝟐𝟐 ) ⇒ 𝑶𝑶(𝑵𝑵𝐥𝐥𝐥𝐥𝐥𝐥𝑵𝑵)(Hierarchical matrices, ℋ-matrices)► 並列計算MPI+OpenMPハイブリッド・プログラミングモデルを採用■Target applications・Electromagnetic fields

・Earthquake cycle

・Quantum mechanics

・Transcribed from NIST Image Gallery
・Transcribed from Ohtani at al (2011)35

階層型行列法■積分方程式法に現れる密行列に対する近似手法の一種
適用対象例ℊ 𝑢𝑢 (𝑥𝑥) = � 𝑔𝑔 𝑥𝑥, 𝑦𝑦 𝑢𝑢(𝑦𝑦)d𝑦𝑦
Ω

singular kernel： 𝑔𝑔 𝑥𝑥, 𝑦𝑦 ∈ span

離散化𝑥𝑥 − 𝑦𝑦 −𝑝𝑝 , 𝑝𝑝 > 0

階層型行列法20,000LRA
Permutation
Partition

Full rank 密行列

Full-Rank
Low-Rank

𝑚𝑚
𝑛𝑛

𝑘𝑘

𝑛𝑛

𝑘𝑘

O (𝒎𝒎𝒎𝒎) ⇒O (𝒌𝒌(𝒎𝒎 + 𝒏𝒏))・行列要素数・行列ベクトル積演算数

36

Memory[GB]メモリ使用量（密行列、階層型行列）

1000

Dense matrices

100ℋACApK(Static electric field)

10ℋACApK

1(Static electric field)

0.1

104

105

106

Matrix size 𝑁𝑁

107

108

Memory usage (Log-Log scale)ℋACApK(Earthquake cycle)37階層型行列法の近似原理■ツリー法，FMM，パネルクラスタリング法などと近似原理は同じ．
線形積分作用素： ℊ 𝑢𝑢 𝑥𝑥 = ∫Ω 𝑔𝑔 𝑥𝑥, 𝑦𝑦 𝑢𝑢 𝑦𝑦 d𝑦𝑦,
Ritz-Galerkin法ℊ 𝑢𝑢 𝑥𝑥 ≅ 𝐴𝐴𝜙𝜙,

遠く離れた２点𝒙𝒙, 𝒚𝒚に対しては,

𝒈𝒈 𝒙𝒙, 𝒚𝒚𝒌𝒌≅ � 𝒈𝒈𝝂𝝂𝟏𝟏 𝒙𝒙 𝒈𝒈𝝂𝝂𝟐𝟐 𝒚𝒚𝝂𝝂=𝟏𝟏𝑢𝑢 ≅ 𝑢𝑢ℎ = � 𝜙𝜙𝑖𝑖 𝜑𝜑𝑖𝑖𝑖𝑖∈ℶ

例：𝑔𝑔 𝑥𝑥, 𝑦𝑦 = 𝑥𝑥 − 𝑦𝑦 −1
𝜑𝜑𝑖𝑖
Ω

Ωℎ𝑖𝑖

𝐴𝐴 ∈ ℝℶ×ℶ , 𝜙𝜙 ∈ ℝℶ

𝐴𝐴𝑖𝑖𝑖𝑖 = � 𝜑𝜑𝑖𝑖 𝑥𝑥 � 𝑔𝑔 𝑥𝑥, 𝑦𝑦 𝜑𝜑𝑗𝑗 (𝑦𝑦)d𝑦𝑦d𝑥𝑥
Ω

変数分離した関数積の級数で近似可能

𝑘𝑘Ω≅ � � 𝑔𝑔1𝜈𝜈 𝑥𝑥 𝜑𝜑𝑖𝑖 𝑥𝑥 d𝑥𝑥
𝜈𝜈=1

Ω

𝑘𝑘

⋅ 𝑘𝑘

積分実行

� 𝑔𝑔2𝜈𝜈 𝑦𝑦 𝜑𝜑𝑗𝑗 (𝑦𝑦)d𝑦𝑦Ω

遠く離れた2領域 Ωℎ𝑠𝑠 ∋ 𝑥𝑥 , Ωℎ𝑡𝑡 ∋ 𝑦𝑦 で考えると，低ランク行列近似 𝐴𝐴|𝑠𝑠×𝑡𝑡 ≅ 𝑉𝑉𝑊𝑊 𝑇𝑇 が可能Ωℎ𝑠𝑠

diam(Ωℎ𝑠𝑠 )・𝑥𝑥

dist(Ωℎ𝑠𝑠 , Ωℎ𝑡𝑡 )

・𝑦𝑦

Ωℎ𝑡𝑡

近似許容：min{diam(Ωℎ𝑠𝑠 ), diam(Ωℎ𝑡𝑡 )}≤ ηdist(Ωℎ𝑠𝑠 , Ωℎ𝑡𝑡 )

𝑠𝑠

≅
𝑡𝑡

𝑘𝑘

𝑡𝑡

𝑘𝑘
38

ベクトルの直積と低ランク行列近似■ベクトルの直積は「ランク１行列」になる

𝑎𝑎 ⊗ 𝑏𝑏 =

𝑎𝑎1
𝑎𝑎2
𝑎𝑎3
𝑎𝑎4
𝑎𝑎5
𝑎𝑎6

⊗ 𝑏𝑏1

𝑏𝑏2

𝑏𝑏3 =

𝑎𝑎1 𝑏𝑏1 𝑎𝑎1 𝑏𝑏2 𝑎𝑎1 𝑏𝑏3
𝑎𝑎2 𝑏𝑏1 𝑎𝑎2 𝑏𝑏2 𝑎𝑎2 𝑏𝑏3
𝑎𝑎3 𝑏𝑏1 𝑎𝑎3 𝑏𝑏2 𝑎𝑎3 𝑏𝑏3𝑎𝑎4 𝑏𝑏1 𝑎𝑎4 𝑏𝑏2 𝑎𝑎4 𝑏𝑏3
𝑎𝑎5 𝑏𝑏1 𝑎𝑎5 𝑏𝑏2 𝑎𝑎5 𝑏𝑏3
𝑎𝑎6 𝑏𝑏1 𝑎𝑎6 𝑏𝑏2 𝑎𝑎6 𝑏𝑏3■任意の行列はベクトルの直積（ランク１行列）で近似できる（精度は度外視）
𝑐𝑐11 𝑐𝑐12 𝑐𝑐13
𝑐𝑐21 𝑐𝑐22 𝑐𝑐23
𝑐𝑐31 𝑐𝑐32 𝑐𝑐33
𝑐𝑐41 𝑐𝑐42 𝑐𝑐43
𝑐𝑐51 𝑐𝑐53 𝑐𝑐53
𝑐𝑐61 𝑐𝑐62 𝑐𝑐63

≃

𝑎𝑎1
𝑎𝑎2
𝑎𝑎3
𝑎𝑎4
𝑎𝑎5
𝑎𝑎6

⊗ 𝑏𝑏1

𝑏𝑏2

𝑏𝑏3

■基底の異なるランク１行列を𝑘𝑘個足すと「ランク𝑘𝑘行列」になる■任意の行列はランク𝑘𝑘行列で近似できる
・𝑘𝑘が十分小さければ低ランク近似

≃

・

𝑘𝑘

39

行列の特異値分解と𝜎𝜎𝑘𝑘 -近似

行列 𝐴𝐴 ∈ ℂ𝑚𝑚×𝑛𝑛 , 𝑚𝑚 ≥ 𝑛𝑛 に対して、■ SVD： 𝐴𝐴 = 𝑈𝑈Σ𝑉𝑉𝑈𝑈 ∈ ℂ𝑚𝑚×𝑛𝑛 , Σ ∈ ℝ𝑛𝑛×𝑛𝑛 , 𝑉𝑉 ∈ ℂ𝑛𝑛×𝑛𝑛

対角行列 Σ の成分 𝜎𝜎1 ≥ ⋯ ≥ 𝜎𝜎𝑛𝑛 ≥0■ 𝐴𝐴の𝜎𝜎𝑙𝑙 -近似： 𝐴𝐴𝑙𝑙 ≔ 𝑈𝑈Σ𝑙𝑙 𝑉𝑉,Σ𝑙𝑙 ≔ (𝜎𝜎1 , ⋯ , 𝜎𝜎𝑙𝑙 , 0, ⋯ , 0)・ 𝜎𝜎𝑙𝑙 -近似の誤差： 𝐴𝐴 − 𝐴𝐴𝑙𝑙 𝐹𝐹 =𝑛𝑛

� 𝜎𝜎𝑖𝑖𝑖𝑖=𝑙𝑙+1

2・ 許容誤差𝜀𝜀に対して 𝐴𝐴 − 𝐴𝐴𝑙𝑙 𝐹𝐹 < 𝜀𝜀 𝐴𝐴 𝐹𝐹 を満たす𝑙𝑙(𝜀𝜀)は、𝑙𝑙(𝜀𝜀) ≔ min{𝑙𝑙 ∈ ℕ:𝑛𝑛

� 𝜎𝜎𝑖𝑖𝑖𝑖=𝑙𝑙+12𝑛𝑛< 𝜀𝜀 � 𝜎𝜎𝑖𝑖
𝑖𝑖=1

2

}

■行列の特異値分解(SVD)は計算コストが高いので、HACApKでは近似手法を使う
40

ACAを用いた低ランク行列近似

Adaptive Cross Approximation
・ピボット列とピボット行を交互に選択

ACA

𝑣𝑣 1 : 任意の列ベクトル (e.g. １列目)
𝑤𝑤 1 : 𝑗𝑗1 行ベクトル, 𝑗𝑗1 ≔ argmax𝑗𝑗 |𝑣𝑣 1𝑗𝑗 |𝑘𝑘𝑎𝑎 ≈ 𝑎𝑎�𝑘𝑘 = � 𝑣𝑣 𝜈𝜈 (𝑤𝑤 𝜈𝜈 )𝑇𝑇
𝜈𝜈=1

𝑚𝑚
𝑛𝑛

𝑘𝑘

𝑛𝑛

・推定される近似誤差：𝑘𝑘

𝜖𝜖𝑘𝑘 : =||𝑣𝑣 𝑘𝑘 || ⋅ ||𝑤𝑤 𝑘𝑘 ||

𝑘𝑘
∑𝜈𝜈=1
||𝑣𝑣 𝜈𝜈 || ⋅ ||𝑤𝑤 𝜈𝜈 ||・経験則：

𝜖𝜖𝑘𝑘 < 𝜖𝜖𝑘𝑘−1■選択するベクトルの本数により使用メモリ量と近似精度を制御可能.

41

階層型行列の生成プロセス■ℋACApKでは青字の方法を採用している
Step 1 Clustering & Construction of a cluster tree
・Method based on geometric information
・Method based on nested dissection
・AMG-like method

Step 2 Construction of a partition■ There’re only blocks(frames of sub-matrices).

Step 3 Fill in sub-matrices■Low-rank approximation
・Cross approximation i) ACA ii) ACA+・Interpolation
・Hybrid method of interpolation and ACA
42

HACApKにおける行列分割の方法■幾何情報に基づくクラスタリングを利用
Admissible condition：min{diam(Ωℎ𝑠𝑠 ), diam(Ωℎ𝑡𝑡 )}≤ ηdist(Ωℎ𝑠𝑠 , Ωℎ𝑡𝑡 )
３

２
１
４

I

１

２

３

４

１
２

I３
４

I

43

HACApKによる行列分割の例

44

HACApKによる行列分割の例

45

ℋACApKにおける階層型行列法の並列化手法■階層型行列の生成において・step 3 （最も時間がかかる部分）のみ並列化.
・MPI の通信は不要.
step1 Cluster tree 作成
step2 H-matrix構造 作成

全MPIプロセスで冗長計算

step3 部分行列の計算(ACA)

並列計算■階層型行列-ベクトル積計算において・全MPIプロセスがベクトル全体を持つ.
・MPIプロセス間の通信が必要.

■上記の両方の並列計算において・計算コアへのデータ割り当て方式は共通にする.
・演算は部分行列を単位として行う.
・計算コアへ割り当てられたデータは部分行列の集合となる.
46

データ割り当て戦略と意図■MPIプロセスについて① 可能な限り𝑅𝑅 𝑀𝑀𝑘𝑘 の最大値を最少にする⇒ 転送データサイズを少なくする② MPIプロセス間の計算負荷不均衡を最小化する■OpenMPスレッドについて
・ OpenMPスレッド間の計算負荷不均衡を最小化する
47

戦略の違いによる割り当てデータ領域の違い

HACApKで各ＭＰＩプロセス
へ割り当てられた部分行列

計算負荷が最適になるような部分行列の割り当て

48

ℋACApKの並列計算性能テスト

次の２つの計算の並列計算性能を調べた
・ 階層型行列の生成
・ HMVM（階層型行列-ベクトル積）■Computer: 東京大学 Ｏａｋｌｅａｆ－ｆｘ

Processor : SPARC64TM Ixfx (16cores/node)
Memory : 32GB
Network : 5 GB/s, Tofu.

■要素数: Ncase 1 : N =1,000case 2 : N = 10,000
case 3 : N = 100,000
case 4 : N =1,000,000■Test model49

階層型行列を生成する計算におけるスケーラビリティ■データサイズが大きくなるにつれて,よりよいスケーラビリティがℋACApKでは得られた.

Speed-up

60

40

1,000,000unkowns
100,000
10,000
1,000

20

0
0

20
40
Number of Processors

6050

階層型行列ベクトル積計算におけるスケーラビリティ
Flat-MPI版ＨＡＣＡｐＫの１ノードによる計算時間を基準とした．

Speed-up vs. 16 core flat-MPI■ MPI+OpenMPハイブリッド版はFlat-MPI版より高い並列計算性能が得られた。ＭＰＩの通信の削減効果によるものと考えられる■Flat-MPIでは、計算速度向上が96コアまでしか得られなかった.
10
・1,000,000 要素

8
6
Flat-MPIMPI+OMP2threads
MPI+OMP4threads
MPI+OMP8threads
MPI+OMP16threads

4
2
0
0

50

100
150
200
Number of cores

250

Parallel scalability when performing an H-matrix vector multiplication

51

HACApK利用コードの構成例
main
メイン

MPI_Init
HACApK_init
初期化

supermatrix_
construction_cog_
leafmtx
H行列構造作成

create_cluster

create_ctree_ssgeom
クラスタツリー作成

RECURSIVE

set_bndbox_cog
バウンディングボックス作成

RECURSIVE
dist_2cluster

create_leafmtx
リーフ小行列作成

HACApK_generate
H行列生成

RECURSIVE

qsort_leafmtx
リーフ小行列ソート

fill_leafmtx
各小行列の近似
chk_leafmtx
setcutthreadacacomp_row

acapluscomp_colHACApK_element_ij
行列要素計算
境界要素積分実行

maxabsvalloc_daccuracy_leafmtx

HACApK_solve
線形ソルバ
仕様を公開

ユーザが用意HACApK_free_leafmtxp
行列メモリ解放bicgstab_lfmtx_p
Bicgstabソルバ

adot_lfmtx_p
H行列・ベクトル積

adot_body_lfmtx
H行列・ベクトル積のbody

dotp_d
内積

HACApK_finalize
制御用構造体メモリ解放
MPI_Finalize

52

HACApKの使用方法
・決められた順番に関数を呼んでいくだけで使用可能

HACApKイニシャライズ

H行列生成
H行列を係数に持つ
線形方程式系を解く
H行列メモリ解放
HACApKファイナライズ
53

ユーザ定義関数 HACApK_element_ij■２つの要素番号と各要素上の情報を受け取り、行列要素値を返す
real*8 function HACApK_entry_ij(i,j,st_bemv)
integer :: i,jtype(st_HACApK_calc_entry) :: st_bemv⋮
return HACApK_entry_ij
end function■境界要素積分に必要な情報を格納する構造体を宣言しておく
type(st_HACApK_calc_entry)
integer :: nd,lp61
real*8,pointer :: ao(:)
integer :: int_ex
real*8 :: dbl_ex
integer,pointer :: int_ex1(:),int_ex2(:,:)
real*8,pointer :: dbl_ex1(:),dbl_ex2(:,:)
end function
54座標情報coord■クラスタリングに必要な座標情報
・行列インデックスと関連付けられている三次元デカルト座標を配列として与える・行列インデックスと関連付けられているのが要素の場合には、重心の座標を与える

55

チュートリアルの流れ
1. ppOpen-APPL/BEM
2. 境界要素法とBEM-BBフレームワーク
3. プログラム実習Ⅰ
・BEM-BB

4. 階層型行列法とℋACApKライブラリ
5. プログラム実習Ⅱ
・BEM-BB+ ℋACApK
56

お試しアカウント付き講習会

サンプルプログラムで計算する問題■静電場解析
・Potential operator：

Air

1V

0.25m

Conductor

0.5m

Ground

解析条件
・半径0.25mの導体球に1Vを与える
・表面に現れる電荷を求める（半無限領域内の電位が導出される）

解析結果
導体球の下半分に大きな電荷密度
が生じる

57

実習を行うあたっての注意点(HACApK) ファイル名

BEM-BBと同時にインストール済みです ジョブスクリプトファイルbem-bb.bash中のキュー名をlecture から tutorial に変更してから
pjsubしてください。
 lecture : 実習時間外のキュー（同時実行数１）
 tutorial : 実習時間内のキュー（同時実行数4+）58

お試しアカウント付き講習会

HACApK用ジョブスクリプト(Oakleaf-fx)(BEM-BB用と同一)#!/bin/bash#PJM -L "rscgrp=lecture"#PJM -L "node=1"#PJM --mpi "proc=1"#PJM -L "elapse=15:00"#PJM -g gt00
export OMP_NUM_THREADS=1
mpirun ./bem-bb-SCM.out input_sample.pbf

利用ノード数
ＭＰＩプロセス数
OpenMP
スレッド数

入力データ名

ノード数、MPIプロセス数、OpenMPスレッド数の設定はOakleaf-fxは１ノードあたり16コアであることに注意して行う
59

お試しアカウント付き講習会

BEM-BB＋HACApKの実行（Oafleaf-fx）
 Oakleaf-fxにログインする ppOpenBEM_***がインストールされていることを確認
 以下のコマンドを実行する$ cd ppohBEM_0.4.1$ cd HACApK_1.0.0$ cd src$ cd HACApK_with_BEM-BB-framework_1.0.0$ make$ pjsub bem-bb.bash
 実行が終了したら、以下を実行する$ cat bem-bb.bash.oXXXXXX (XXXXXXは数値)
60

お試しアカウント付き講習会

BEM-BB＋HACAｐKの実行（Oafleaf-fx）■以下のような結果が見えれば成功
Number of processes 1
Number of unknowns set on each face element = 1
Selected linear solver = BICGSTABConvergence criterion = 1.000000000000000E-06
Upper limit of iteration counts = 5000
Number of total threads, Thread number 1 0
ext_ndim = 21600 , ndim = 21600***************** HACApK start ********************
nd= 21600 nofc= 21600 nffc=
1
nrank=
1 nth=
1
No. of cluster=3903
No. of nsmtx 17002 33096 16699
1:Rk-matrix 2: dense-mat 3:H-matrix
ncpc= 24306128 ncpc/nrank= 24306128
time_supermatrix= 0.1259064860641956
time_fill_hmtx= 87.13198765483685
time_construction_Hmatrix = 87.25789824104868
Memory of the H-matrix= 327.1665954589844 (Mbyte)
Memory compression v.s. dense matrix=
9.191182270233195 (%)（中略）
before HACApK time = 0.266355455 [sec]HACApK time = 98.4935837 [sec]_____________all time = 98.7606354 [sec]MPIプロセス数
OpenMPスレッド数
要素数

H行列の使用メモリ
密行列に対する
H行列の圧縮率
計算時間61

演習課題１(HACAｐK)
１-1 サンプル以外のデータを用いた計算演習１で使用したデータはinput_sample.pbf(要素数：600）でした。
Input_10ts.pbf(要素数：約1000)およびInput_216h.pbf(要素数：約2,1000)のデータを用いた計算を行い、計算時間を測定してください。
・ヒント：ジョブスクリプトに書かれているデータ名を変更する

密行列を用いた場合の計算時間と比較してください。
１-２ 計算オーダーの計測
演習１-1で測定した計測時間を使って、階層型行列を用いた境界要素解析の計算時間が、要素数を𝑁𝑁としてO(𝑁𝑁log𝑁𝑁)であることを確かめて下さい。
62

演習課題２、３(HACAｐK)
3-1 コードの並列実行
いろいろなOpenMPスレッド数やMPIプロセス数で計算を実行し、計算時間を計測してください。
３-2 台数効果による性能向上の評価
逐次実行に比べ、どれだけ高速化されたか、性能評価をしてください。
３ 計算環境の移動
ppohBEM_0.4.1.tar.gzをReedBushへコピーし、BEM-BB＋HACAｐKを動作させてください。
・ヒント：Makeファイルの先頭数行で計算環境を設定しています
・ヒント：ジョブスクリプトは他の演習で使ったものを参考に自作
63お疲れさまでした64

お試しアカウント付き講習会