Post Peta CRESTで開発した

行列ライブラリ ACApK（ライブラリと使用法について）
伊田 明弘 （東京大学）
PCCCワークショップin大阪2020
2020年 2月 21日
富士通株式会社 関西システムラボラトリ（大阪）

目次
1
2
3
4
5

2

ℋACApKライブラリの概要
低ランク構造行列法（格子）ℋℋ行列法の並列化ℋACApKライブラリの使用方法
今後の取り組み概要：ppOpen-HPC
概要：■JST-CREST
・ “ポストペタスケール高性能計算に資するシステムソフトウェア技術の創出”領域
・ “自動チューニング機構を有するアプリケーション開発・実行環境ppOpen-HPC”研究代表者: 中島 研吾(東京大学・情報基盤センター)■Download site: http://ppopenhpc.cc.u-tokyo.ac.jp/

この話の主題ppOpen-HPC covers …
ℋACApK3

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

4

FDM, FEM, BEM による離散化
continuous

For  ∈ Ω ⊂ ℝ , find () s.t.
・Bi-linear form

5

discrete

Az = .

∆() = () ()()d =  () d・Green's identity ()()d =    () d −  %() · %()d$ %() · %()d = ,  ̅・Fredholm integral equation  ,   () d =  ̿

FEM・Base function$ )() $  ,   () ddy = ), +




BEM
・Base function■
5

ℋACApK ライブラリ■境界要素法（BEM）シミュレーション向け高速計算ライブラリ■境界要素法のボトルネック►ℊ 1 = , where ℊ 1  : =  3 ,  1()d・Degenerate kernel：3 ,  ≅ 8<

39:  3;:  ,for far remote >, ?

:=9e.g. 3 ,  ∝  −  A9

・旧来手法では密行列を扱いが必要

計算量4(56 )■大規模計算に向けて► 密行列に対し近似手法の導入
・低ランク構造行列法：B(56 ) ⇒ B(5DEF5)ℋ行列,
行列, 格子ℋ
格子ℋ行列,
行列, BLR行列► 分散メモリ並列計算
・Hybrid MPI+X programming model

6

ℋACApK ライブラリ■対応プラットフォーム
プラットフォーム
・CPU clusters
・Multi- GPUs（整備中）█ プログラミング言語
・Fortran90+OpenMP・CUDA, CUDA Fortran, OpenACC (GPU対応)■必要な外部ライブラリ
・MPI█ 公開形態・Open source ・MIT license█ ダウンロード サイト

・https://github.com/Post-Peta-Crest/ppOpenHPC/tree/MATH/HACApK・http://ppopenhpc.cc.u-tokyo.ac.jp（旧Version）

7

ℋACApKの適用例■Static electric field analysis
・Potential operator：
・Induced surface charge is calculated in half-infinite domain.

HI

GroundAnalysis condition

Analysis result
8

ℋACApKの適用例（CDWM）■Calculate LLG equation with terms of spin torque effect
N
V
RM
RM
− J⋅L M−
M× J⋅L M= −TM × U −
M×
OP
RS
OP
RS
Spin torque term
・Effective field
b
1
`
2hi
2l
U=
$ L · M  de + ; M ⋅ k k + ; M + Umnn
4^_I ca
Oj
Oj

Demagnetizing field Uo

Other field components U

p■Discretization method
M|X , U|X
M|XY∆Z・Time : Explicit Runge–Kutta
・Space : Magnetic material is divided into [ cubic elements

Initial position of DW

moving direction of DW

9

ℋACApKの適用例（CDWM）■Current-induced Domain Wall Motion (CDWM)■CDWM is induced by a spin-polarized current in magnetic nanowires.
120nm
magnetization direction

5nm

initial position of Domain Wall
x=0

z
y

I

x

magnetization directione-

4000nm1.0
Snapshot of DW(final step)

1.0

qs
-1.0

qr

[Ida : IEEE Trans. on Mag. 2020 ] -1.0 10

ℋACApKの適用例[Ida : IEEE Trans. on Mag. 2017 ]■Micromagnetic simulation of spin torque oscillator

1・Potential operator of demagnetic field ： 4^ $ )9 ? %|> − ?|A9 ⋅ dt
c
10nm~800nm

・2-nd order singularities
・Vector quantities
5 layers
10nm

・Mesh size : 2～3nm

・Oscillating frequency spectrum is calculated in magnetic materials.
・H-matrices accuracy of u = 10 Aw
is required to reproduce oscillation
processes.

11

ℋACApKの適用例[Tominaga : IEEE Trans. on Applied Superconductivity 2018 ]■Superconductor Simulation
・Potential operator：
・Current density is calculated in coils wound with coated conductors.

Model (N=1,516,060)

Current distribution

Memory usage
12

ℋACApKの適用例[Mifune : Superconductor Science and Technology 2019]■Superconductor Simulation accelerated by AMG
・Governing equation： Faraday’s law and Biot–Savart’s
1
R _I Sy% × }~ { ~ × ` ⋅ } ~%×% × }{ ⋅ } +$d = 0
xy [{]
RS 4^ P a
Galerkin and backward finite difference methodsl[ A9 ] ⋅   + ⋅ (  −  A9 ) =0.
SparseDense

・Previous study- Apply ℋACApK to integral- Solve by Newton-Raphsonusing BiCG-STAB

・New approach- Set l as preconditioning matrix
in the BiCG-STAB- Apply AMG to the preconditioner
13

ℋACApKの適用例■Earthquake Cycle Simulation[Ohtani : Proc. Comp. Sci. 2014 ]Integral operator with  ,  ∝  −  A
N

G

Vi
➢eq. motion
2Vsi 1➢friction law
 i   neff   *  Ai ln(Vi / V* )  Bi ln(V*i / Li )

 i   K ij (u j  V pl t ) d i exp(Vi / Vc )  Vii / Li ln(Vii / Li )dt

Subdivide
fault surface

14

ℋACApKの適用例（地震発生物理学：動的破壊シミュレーション)
2016年熊本地震の例：• 熊本側の本震と，大分側の動的誘発地震の発生を動力学的に再現• 現実的な断層形状を考慮• 動的誘発地震のメカニズムを解明： 本震が放射した強い応力波が50km以上を伝播
断層を横から見たところ

断層面の
変位速度

大分側

熊本側

大分

地表面の
変位速度

熊本

断層と地表面を俯瞰したところ

変位

断層を横から見たところ

せん断応力
動画

活断層地図
空間要素数：18,000
時間要素数： 1,200
計算機：Xeonクラスタ
並列数：24コア×２ノード
計算時間：約１０時間

15

ℋACApKの適用例（地震発生物理学：動的破壊シミュレーション)• オリジナルBEM
積分核： h • 高速領域分割法（FDPM）： Ando （2016, GJI)
積分核の領域ごとの時間依存性と積分領域形状の違いを利用h 
 Domain F→時空間依存だが高々面積分• 積分核： h   + h   Domain I →時間項を分離した体積積分→時間項のない面積分
h P
 Domain S
9h   = ⋮ [9


⋯

 ]

Σ Σ h   A
積分領域の分割

効率化：密行列→ベクトルの積• FDPMによる計算コストの削減• 4 O; [ → 4 O;
メモリ
• 4 O[  → 4 O[ ; ~4 O; 演算数（M: 要素数. N ：時間ステップ数）
まだ減らしたい→H行列法適用

16三角錐： 動弾性の影響範囲

16

ℋACApKと密行列のメモリ使用量と密行列のメモリ使用量■ ℋ-matrices with ACA reduce memory usage.

Memory[GB]1000

Dense matricesℋACApK

100(Static electric field)

10ℋACApK

1(Static electric field)

0.1

104

105
106
107
108
Number of unknowns 5

Memory usage (Log(Log-Log scale)ℋACApK(Earthquake cycle)
17

18

1

ℋACApKライブラリの概要

2 低ランク構造行列法
3
4
5

（格子）ℋ
ℋ行列法の並列化ℋACApKライブラリの使用方法
今後の取り組み低ランク行列近似とベクトルの直積
低ランク行列近似(LRA)と近似■ベクトルの直積は「ランク１の行列」

⊗. =

9
;


w


⊗ .9

.;

. =

9 .9 9 .; 9 .
; .9 ; .; ; .
 .9  .;  .
 .9  .;  .
w .9 w .; w .
 .9  .;  .

■任意の行列はベクトルの直積（ランク１行列）で近似できる（精度は度外視）
99 9; 9
;9 ;; ;
9 ; 9 ; w9 w w
9 ; 

≃

9
;


w


⊗ .9

.;

.

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

≃

・


19

低ランク構造構造行列法
構造█積分方程式法に現れる密行列に対する近似手法
適用対象例ℊ 1 () = $  ,  1()dDegenerate kernel：  ,  ≅ 8離散化¨9<  ;< <=9for far remote >, ?

低ランク構造行列法

LRA

20,000Permutation
Partition

Full rank 密行列

Full-Rank
Low-Rank

・行列全体ではfull rank

q








O (§}) ⇒O (k(§ + }))・行列要素数・行列ベクトル積演算数

20

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

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

ACAª 9 : 任意の列ベクトル (e.g. １列目)
« 9 : ±9 行ベクトル, ±9 ≔ argmax |ª 9 |¨ ≈ ©¨ = 8 ª : (« : )¬:=9・推定される近似誤差：

q







­¨ : =b||ª ¨ || ⋅ ||« ¨ ||∑¨:=9 ||ª : || ⋅ ||« : ||・経験則：­¨ < ­¨A9■選択するベクトルの本数により使用メモリ量と近似精度を制御可能.

21

Procedure of ACA with full pivoting (ACA_full)█Pivoting and updating matrix need computational cost B(56 ).
ªInput:  ∈ ℝj×X .
Output: ª ¨ ∈ ℝ|j| , « ¨ ∈ ℝ|X|
«
<
¨
¨
¬
µ
s. t. © = ∑¨=9 ª (« ) .

a = 0, ­< = 1

For a given error tolerance u < 1,
Do while »` > uSa =a+1( ∗ , ± ∗ ):= argmax( ,) | , |, δ ≔  ∗ , ∗Pivoting : B(56 )searching over whole matrix(ª < ) :=  , ∗ ,  ∈ µ,(« < ) :=  ∗ , /À, ± ∈ S,b<<­< : = ||ª ||; ||« ||; /∑<¨=9 ||ª ¨ ||;; ||« ¨ ||;;

 , :=  , − (ª < ) (« < ) ,  ∈ µ, ± ∈ S
End doEstimation of
approximation error
Updating matrix : B(56 )Subtracting rank 1 matrix
22Procedure of ACA with partial pivoting (ACA_part)■Computational cost : B(5)Partial pivoting : B(5)searching around a vector
Vector calculations : B(5)instead of updating matrix
Estimation of
approximation error
Omitting the selected row
from searching space23

低ランク構造行列法の近似原理
低ランク構造行列法の近似原理■ツリー法，FMM，パネルクラスタリング法などと近似原理は同じ．
線形積分作用素： ℊ 1  =   ,  1  d,b

Ritz-Galerkin法ℊ 1  ≅ l,

遠く離れた２点>,に対しては,
遠く離れた２点 ?に対してはに対してはk

Ä >, ? ≅ 8 ÄÆÅ > ÄÆ6 ?

1 ≅ 1Â = 8  )例： ,  =  −  A9
)

∈ℶ

Ω

ΩÂ

l ∈ ℝℶ×ℶ ,  ∈ ℝℶ

l  = $ )  $  ,  ) ()dd
¨

Æ=Å

変数分離した関数積の級数で近似可能
変数分離した関数積の級数で近似可能≅ 8 $ 9:  )  d
:=9



⋅ 

積分実行

$ ;:  ) ()d

遠く離れた2領域 ΩÂj ∋  , ΩÂX ∋  で考えると，低ランク行列近似 l|j×X ≅ eÊ ¬ が可能ΩÂjdiam(ΩÂj )・

dist(ΩÂj , ΩÂX )

ΩÂX
・

近似許容：
：min{diam(ΩÂj ), diam(ΩÂX )}Ç ηdist(ΩÂj , ΩÂX )

µ

≅
S

S




24

低ランク構造行列
低ランク構造行列の行列の生成プロセス■ℋℋ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
25

ℋACApKにおける行列分割のACApKにおける行列分割の方法（従来における行列分割の方法（従来ℋ
方法（従来ℋ行列の場合）
行列の場合）■幾何情報に基づくクラスタリングを利用

Admissible condition：：min{diam(ΩÂj ), diam(ΩÂX )}Ç ηdist(ΩÂj , ΩÂX )
３
２
１
４

I

１

２

３

４

１
２

I３
４

I

26

従来ℋ
従来ℋ行列法を
行列法を用いた行列分割の例
用いた行列分割の例

27

従来ℋ
従来ℋ行列法を
行列法を用いた行列分割の例
用いた行列分割の例

28

低ランク構造行列法と近縁手法の種類■行列分割構造と要素保存法

分割構造

行列要素保存法
数値を直接保存
基底ベクトルと係数基底関数と係数

Strongadmissibility

Block

Weak
admissibility

従来ℋ-行列 BLR-行列
格子ℋ-行列HODLRℋ ; -行列

HSS

FMM

Tree

・青字： ℋACApKで取り扱える手法
29

既存の低ランク構造
既存の低ランク構造行列法の比較
構造行列法の比較■分割構造の違いにより、一長一短
・従来Ë-行列
従来 行列

・BLR, TLR

・HSS, HODLR

・Complicated structure

・Simple structure

・Relatively simple structure

like block-divided matrix
密行列用アルゴリズムを流用可能

・Very high compressibility
for low dimensional problem

効率的な演算・通信
パターンの構築が困難

・High compressibility
Memory : B(5DEF5)

・Low compressibility
Memory : B(5Å.Ì )・Strong admissible cond.

・Strong admissible cond.

・Worse asymptotic complexity
for high dimensional problem
・Weak admissible cond.
(All off-diagonal : Low-rank)
30

格子Ë行列法
格子 行列法 [A.Ida : IPDPS’18]■低ランク構造行列法の新手法
低ランク構造行列法の新手法
・Hybrid of ℋ-matrices and BLR-matrices
・Advantages of lattice ℋ-matrices- High compressibility of ℋ-matrices- Convenience of matrix arithmetic with BLR-matrices
・Lattice Ë-matrix

・BLR-matrix

Lattice structure・Ë-matrixℋ-submatrices
in blocks on lattice

31

格子ℋ
格子ℋ行列のメモリ使用量■格子サイズの決め方に依存・ Í ∶ 格子サイズ
・Í=[

従来ℋ行列

・Í=1

密行列
9
Î

・ Í = V[ for a constant V ( Ç V Ç 1) : 格子数固定

要素数 :  [ = 2T5 DEF 5 + [(2T log V + 2/V − T)・格子サイズ Í はMPIプロセス数に応じて決められるべき
利用可能なMPIプロセス数は行列サイズ[に応じて制限
32

低ランク構造行列の行列分割方法
・既存Ë行列・BLR行列・Cluster tree of index set ℶ

格子 Ë行列
・格子

33

34

1
2

ℋACApKライブラリの概要
低ランク構造行列法

3 （格子） 行列法の並列化
45ℋACApKライブラリの使用方法
今後の取り組み境界要素法解析で必要なℋの２機能
境界要素法解析で必要なℋACApKの２機能
演算は部分行列ごとに実施■ 低ランク構造行列生成
・低ランク近似(LRA)は各ブロックで個別独立に実施■ 低ランク構造行列-ベクトル積
・部分行列-ベクトル積の集まり

・
・

・

Add

=
35

ℋACApKにおけるにおける低ランクにおける低ランク構造行列法
低ランク構造行列法の並列化
構造行列法の並列化■低ランク構造行列の生成において・step 3 （最も時間がかかる部分）のみ並列化.
(step1,2も並列化した版は公開準備中）

・MPI の通信は不要.
step1 クラスタ ツリー 作成
step2 ℋ行列構造 作成
step3 部分行列の計算(ACA)

全MPIプロセスで冗長計算
並列計算■低ランク構造行列-ベクトル積計算において・全MPIプロセスがベクトル全体を持つ（ベクトル分散版は公開準備中）
・MPIプロセス間の通信が必要.

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

36

ℋACApKにおける従来における従来ℋにおける従来ℋ行列の並列計算手法

基本的にℋ行列を行方向にスライスするが、部分行列は切断しない■ MPIプロセス割当て戦略
MPIプロセス割当て戦略① 各ブロックの高さÐ O¨ を可能な限り最小化⇒ 通信データサイズの最小化② MPIプロセス間の計算負荷不均衡を最小化
×>█ 従来ℋ
従来ℋ行列-ベクトル積
行列 ベクトル積 ? ≔ Ö・各MPIプロセスで部分行列-ベクトル積を行った後に隣接間通信
l+

P0~7


P0

P1

P2

P3

P4

P5

P6

P7

MPI

37

格子ℋ行列に対するプロセス割当て戦略
格子 行列に対するMPIプロセス割当て戦略
行列に対する■ 「プロセスグリッド」に基づく2次元サイクリック割当て
・各MPIプロセスへの割当て : 同色で塗られた格子-複雑な構造は逐次orスレッド並列では許容範囲-密行列向けの既存アルゴリズムを流用可能■ 例: 3×2 プロセスグリッド (6 MPI-processes)：Process(0,0)：Process(2,1)(0,0) (0,1)(1,0) (1,1)(2,0) (2,1)

Assignment global view

Lattice ℋ -matrix
Local view
38

×>
格子Ë行列
格子 行列-ベクトル積
行列 ベクトル積,
ベクトル積 ? ≔ Ö

アルゴリズム 1 : “LatticeHMVM_fullVector”■各MPIプロセスで部分行列-ベクトル積を実行し部分ベクトルを得る■全プロセスで全ベクトル  を得るために,39

× >（ベクトル分散版）
格子Ë行列（ベクトル分散版）
格子 行列-ベクトル積
行列 ベクトル積,
ベクトル積 ? ≔ Ö

アルゴリズム 2 : “LatticeHMVM _partVector”■正方形プロセスグリッド限定: [n< = [nØ = b [n・各MPIプロセスが全ベクトルを持つ必要がない
× >を計算するため必要な部分ベクトルのみを通信
■次の Ö

40

数値実験：Current induced Domain Wall Motion(CDWM)
■ (格子
格子)
行列の近似精度格子 Ë-行列の近似精度
・Tolerance ε =1.0e-4 for ||l − l+||Ù Ç u.
■計算メッシュ・ Thin films of magnetic material
magnetization direction

initial position of DW
x=0

80nm
5nm
3000nm

4000nmuniform mesh (4nm)
x: -2000～～2000

3000nmgradation mesh (4nm～～120nm)x: 2000～～5000

l：Dense matrix,l+: ℋ-matrix・Single film model
N =20,910
N =31,860
N =116,100
・Double film model
N =232,200
・Triple film model
N =348,300█計算機: Fujitsu PRIMERGY CX2550 (北海道大学)
・Processor : Intel Xeon™ Gold 6148 (40cores/node)・Memory : 384GB/node, ・Network : 12.5 GB/s, Fat Tree.
41

格子ℋ
行列の計算オーダー
格子ℋ-行列の計算オーダー■ 格子サイズ Í の決め方で 4 [ log [ ～ 4([ ; )・Í=[従来 ℋ-行列

・Í=1

密行列■ MPI数で格子サイズを決定数で格子サイズを決定 (Í = [⁄ b 9[n ) : B 5 DEF 5
10-1

105

10-2

104

10-3

Memorylattice H-mat-vec
Linear103
Square105
Matrix size

Memory Usage [MB]Time of lattice H-mat-Vec[s]・Memory usage of Lattice ℋ-matrices
・Calculation time of Lattice ℋMVM when using 400 MPI processes

106
42

MPI プロセス間の計算負荷均衡■計算負荷均衡効率H≔__ÜÝÞß [n_ : Memory usage(entire matrix)_ÜÝÞß : The largest memory(an MPI process)・従来Ë行列 : MPIプロセス数の増加により効率が急速に低下
・格子Ë行列 : 効率低下の度合いが穏やか
1
Load blancing efficiency

Good

Poor

Lattice H (rectangle)
Lattice H (square)
Normal H-materices

0.8

・” Single film model”
N =20,910

0.6
0.4
0.2
0

101

102
103
Number of MPI processes

43

(格子
格子)
行列-ベクトル積の並列計算性能
格子 Ë行列行列 ベクトル積の並列計算性能■ 以下のMPIプロセス数まで計算速度向上を観測
・従来ℋℋ行列 : 約100 MPI プロセスまで
・格子ℋℋ行列（全ベクトル保持） : 約1,600 MPIプロセスまで
・格子ℋℋ行列（ベクトル分散） : 少なくとも 3,600 MPIプロセスまで■ MPIプロセス数[n > 100では格子ℋℋ行列が従来ℋℋ行列より高速■ MPIプロセス数[n > 600では格子ℋℋ行列（ベクトル分散）が最速
Time of H-mat-Vec[s]

10-1

10-2

Normal H-matrices
Lattice H (rectangle)
Lattice H (square)
Linear speedup” Single film model”
N =20,910

10-3

10-4 0
10

101
102
103
Number of MPI processes

44

45

1
2
3

4
5

ℋACApKライブラリの概要
低ランク構造行列法（格子）ℋℋ行列法の並列化

ACApKライブラリの使用方法

今後の取り組みℋACApKでで仕様を公開している関数■関数一覧と機能

赤字：必ず使用される関数

関数名

機能、コメント

HACApK_init

環境設定；最初に呼ばれる必要あり

HACApK_generate

低ランク構造行列の生成

HACApK_solve

線形方程式の求解

HACApK_adot_pmt_lfmtx_p

行列・ベクトル積計算

HACApK_free_leafmtxp

低ランク構造行列のメモリ開放

HACApK_finalize

後始末・境界要素法シミュレーションに最低限必要な機能を供給■静的な問題を解く場合に使用される関数：HACApK_solve
例：静電場解析■動的な問題を解く場合に使用される関数：HACApK_adot_pmt_lfmtx_p
例：マイクロマグネティクス計算、超電導解析、地震サイクル解析、動的破壊解析

46

main
メイン

ℋACApK利用コードの構成例
利用コードの構成例
MPI_Init
HACApK_init
初期化

supermatrix_
supermatrix_
construction_cog_
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

47

ユーザ定義関数 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

j

i

Öàá

■境界要素積分に必要な情報を格納する構造体を宣言しておく
type(st_HACApK_calc_entry
type(st_HACApK_calc_entry)
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
48

ℋACApKの使用方法
・モジュールと構造体を宣言し、決められた順番に関数を呼んでいくだけで使用可能ℋACApKモジュール呼出ℋACApK使用構造体宣言ℋACApKイニシャライズ
配列のメモリ確保ℋ行列の要求精度設定ℋ行列生成ℋ行列を係数に持つ
線形方程式系を解くℋ行列メモリ解放ℋACApKファイナライズ
49

ℋACApKへのへの入力情報への入力情報■行列サイズ：nd■ℋ行列への要求精度：ztol■線形方程式の右辺ベクトル：rhs■行列ベクトル積の元ベクトル：x
lrtrn=HACAkK_adot_pmt_lfmtx_p(st_leafmtxp,st_bemv,st_ctl,ax,x)50

ℋACApKへへの入力情報
入力情報：
情報：座標情報coord(:,:)■クラスタリングに必要な座標情報

・行列インデックスと関連付けられている三次元デカルト座標
・要素が関連付けられている場合には、要素重心の座標

51

HACApKのログ■ 制御パラメタst_ctl%param(1)の値が１以上で出力***************** 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 (%)（中略）

HACApK time = 98.4935837 [sec]

要素数
OpenMPスレッド数
MPIプロセス数
H行列の使用メモリ
密行列に対する
H行列の圧縮率
計算時間52

ℋACApKの一歩進んだ使い方の一歩進んだ使い方■制御パラメタ：st_ctl%param(1:100)
番号 初期値説明

1

1

ログ出力（0:エラーのみ, 1:標準出力, 2:デバグ用出力）

101

低ランク構造行列の種類（1:通常ℋℋ行列, 10:格子ℋℋ行列, 20:BLR行列）

110近似行列の精度チェック（0：チェックしない, １：チェックする）

21

15

最小クラスタサイズ

51

2.0

クラスタ間の距離判定パラメータ611

部分行列での低ランク近似手法選択（1:ACA法, 2:ACA+法）

63

1000

部分行列の最大ランク■制御パラメタの設定方法
・初期化関数HACApK_initと近似行列生成関数HACApK_generateの間で設定
lrtrn=HACApK_init(nd,st_ctl,st_bemv)
allocate(coord(nd,3),rhs(nd),sol(nd))
ztol=1.0e-5; st_ctl%param(10)=10; st_ctl%param(11)=1;
lrtrn=HACApK_generate(st_leafmtxp,st_bemv,st_ctl,coord,ztol)
53

制御パラメタ：st_ctl%paramの解説■低ランク構造行列法の近似原理
遠く離れた2領域 ΩÂj ∋  , ΩÂX ∋  で考えると，低ランク行列近似 l|j×X ≅ eÊ ¬ が可能ΩÂjdiam(ΩÂj )・

dist(ΩÂj , ΩÂX )

ΩÂX

⇒

µ

・

近似許容：
：min{diam(ΩÂj ), diam(ΩÂX )}Ç ηdist(ΩÂj , ΩÂX )

近似

S

S




最大値設定
param(62)

クラスタ間の距離判定
param(51)最小クラスタサイズ
param(21)

手法選択param(61)

54

55

1
2
3
4

ℋACApKライブラリの概要
低ランク構造行列法（格子）ℋℋ行列法の並列化ℋACApKライブラリの使用方法

5 今後の取り組み

低ランク構造行列分解への取り組み■取り組み中の行列分解
►LU分解

►QR分解

■想定する利点■想定する利点►密行列計算より軽い・メモリ使用量
メモリ使用量：B(5âãÄ5)< 4([ ; )メモリ使用量
・演算量：B(äDEF 6 5) ≪ 4([  )►疎行列の不完全分解より高精度►ランク数の調整により精度が制御可能？

■想定する使い道►低精度：クリロフ部分空間法の前処理、固有値の近似値計算►高精度：密行列による完全分解の代替

56

格子ℋ
格子ℋ行列に基づく数値線形代数の構築■データ科学計算および科学技術計算の大規模化・高速化►B(56 )以上の密行列演算をB(5DEFæ 5)以下に低減►密行列線形代数ライブラリ（ BLAS, LAPACK, ScaLAPACK）の置き換え■格子ℋℋ行列計算アルゴリズム群の開発:►行列演算種類の拡充

・開発済み：格子Ｈ行列生成，行列-ベクトル積，LU 分解
・開発中：QR 分解(BLRでは開発済)
・開発予定：行列-行列積，固有値計算，特異値計算►格子ℋ行列法の深化
・複素数化および長方行列への対応（現状は実数正方行列のみ対応可能）・テンソルへの拡張►データ科学分野の問題への適用

・最尤推測法に用いられる各種確率分布に基づく共分散行列への適用
・３次元ユークリッド空間と高次元データ空間の違いを克服
・ベイズ推論および深層学習へと適用範囲拡大

57

科学技術計算での密行列演算回避
微分作用素
離散化

積分作用素
逆作用素

離散化

逆行列B(5ç )密行列B(56 )

疎行列
格子ℋ行列
表現

近似離散化

ここを
回避

近似

格子ℋ
近似逆行列
4([ log ; [)格子ℋ行列4([ log [)58

階層的低ランク行列近似
階層的低ランク行列近似の近似の適用範囲① 線形積分作用素から得られる密行列の効率的な近似行列表現を得る② 上記①で得られた階層型行列に対し，行列計算（加法, 乗法, LU分解, QR分解など）を行う③ 線形微分作用素から得られる疎行列に対し，近似逆行列を計算する④ 線形微分作用素から得られる疎行列に対し，近似LU分解を計算する③の直感的理解(Bebendorf（2008）より） ・微分の逆作用素は積分のはず−∆1 =  in Ω ≔ , . , 1  = 1 . = 0
２次精度中心差分で離散化

線形積分作用素の⇒階層型行列
離散化行列に類似

59

データ科学計算への格子
データ科学計算への格子ℋ
格子ℋ行列適用■尤度関数，
フィッシャー情報行列■Dimensionality Reduction►主成分分析のホットスポット：共分散行列の固有値・固有ベクトル計算
・密行列： メモリB(56 )，演算量B(5ç )・格子ℋ行列 ： メモリB(5DEF5)，演算量B(5DEF 6 5)

■ディープニューラルネットワークのハイパーパラメータ予測►ヘッセ行列とヤコビ行列の関係式から最適な学習率やバッチサイズを見積り
・ McCandlishらの研究では対角近似 ⇒ 格子ℋ行列適用により予測精度向上

60

格子ℋ
格子ℋ行列の最新アーキテクチャへの高性能実装法■ハードウェアの性能を引き出せるアルゴリズム開発■最新計算機アーキテクチャ向け格子H行列計算コード開発:►次世代CPU 向けコード開発
・スパコン「京」の次世代機「富岳」で採用されるCPU「A64FX」
・メニーコア型CPU, ベクトル命令セット採用，小容量高速メモリ搭載►GPU向けコード開発

・マルチコアCPU 向けが既開発の演算（格子H 行列生成，行列-ベクトル積，LU 分解）
・テンソルコアの効率的な利用法►FPGA を用いた格子ℋℋ行列分割生成および行列分解

・格子ℋℋ 行列分割構造の作成ルーチン（条件分岐が多く、CPU・GPU に不向き）
・格子ℋℋ行列分解（深い演算パイプラインと細粒度での並列処理が可能なFPGA向き）►混合精度演算による計算の効率化

・部分行列のランク数に応じて数値の保存精度を決定
・誤差が蓄積する縮約演算の量に応じて演算精度を決定►動的手法を用いた高速化（ヘテロジニアスなシステムへの対応）
・タスク並列言語(Cilk Plus, Tascell)を用いた動的負荷分散
・ランタイムライブラリ（StarPU, PaRSECなど）の活用

61

Acceleration of ℋ-mat. generation using KNL[Hoshino : CLUSTER 2018 ]■ Test model of electrostatic field analysis►Perfect conducting sphere►Dielectric sphere•è[1]() ∶= Ω

D[1]() ∶= Ω

including branch divergence9é rAsrAs,(s)é rAs ë

１V►KNL : Intel Xeon Phi 7250, 68 core

-qopenmp -O3 -ipo -align array64byte-xAVX2 (BDW) –xMIC-AVX512 (KNL)■ Performance comparison►BDW : approximately 2x speedup

Ground16
14
Elapsed time [sec]►Compiler : Intel compiler 18.0.1

0.25m

Coefficient H-matrix generation
on BDW and KNL►BDW : Intel Xeon E5-2695 v4, 18 core

12
1.9x

10

►KNL : over 4x speedup•

1  d,  ∈ Ω

User-defined functions depend onthese integral equations■ Evaluation Environments•

1  d,  ∈ Ω

In the case of dense matrix generation,
new design achieved at most 6.6x speedup

8

4.1x

2.2x

4.3x

6
4
2
0
Perfect
BDW

Dierectric
BDW
Original

Perfect
KNL

Dierectric
KNL

SIMD design62

Acceleration of ℋ-mat.-vec. using GPU, KNL[Hoshino : ICCS 2018 ]0.3■Memory bandwidth bounded
Elapsed time / 1 BiCGSTAB step [sec]

0.25

0.2

略 称

P100

KNL

BDW

動作周波数
（GHz）

1.480

1.40

2.10

コア数（最大有効スレッド数）

3,584

68 (272)

18 (18)5,304

3,046.4

604.8

16

16

128

534

490

65.5

理論演算性能（GFLOPS）
主記憶容量（GB）
メモリバンド幅（GB/sec.,
Stream Triad）
）

0.15

13.5x
0.1

6.78x
0.05

7.48x
8.15x

0
P100
100ts

KNL
human1x1

BDW
216h

HACApK関連
HACApK関連の関連の査読付き論文リスト
1.

2.

3.
4.

5.

6.

7.

8.

9.

Rise Ooi, Takeshi Fukaya, Takeshi Iwashita, Akihiro Ida and Rio Yokota, “Accelerating Multithreaded Linear
Solver with Mixed Precision Hierarchical Matrix Computation and Data Structure”, International Conference
on High Performance Computing in Asia-Pacific Redion (HPC Asia 2020), pp.92-101.
Yu Pei, George Bosilca, Ichitaro Yamazaki, Akihiro Ida and Jack Dongarra, “Evaluation of Programming
Models to Address Load Imbalance on Distributed Multi-Core CPUs: A Case Study with Block Low-Rank
Factorization”, In: 2019 IEEE/ACM Parallel Applications Workshop, Alternatives To MPI (PAW-ATM). IEEE,
2019. (In press)
Akihiro Ida, Tadashi Ataka and Atsushi Furuya, “Lattice H-matrices for Massively Parallel Micromagnetic
Simulations of Current-induced Domain Wall Motion”, IEEE Transactions on Magnetics, (In press)
Zhengyang Bai, Tasuku Hiraishi, Hiroshi Nakashima, Akihiro Ida and Masahiro Yasugi. “Parallelization of
Matrix Partitioning in Construction of Hierarchical Matrices using Task Parallel Languages” , Journal of
Information Processing, Vol. 27 (2019) pp．840-851．
Akihiro Ida, Hiroshi Nakashima, Tasuku Hiraishi, Ichitaro Yamazaki, Rio Yokota and Takeshi Iwashita, “QR
Factorization of Block Low-Rank Matrices with Weak Admissibility Condition”, Journal of Information
Processing, Vol. 27 (2019) pp．831-839．
Satoshi Ohshima, Ichitaro Yamazaki, Akihiro Ida and Rio Yokota, “Optimization of Numerous Small DenseMatrix–Vector Multiplications in H-matrix Arithmetic on GPU”, In: Auto-Tuning for Multicore and GPU (ATMG)
In conjunction with the IEEE MCSoC-19, Singapore (Oct 1-4, 2019)
Ichitaro Yamazaki, Akihiro Ida, Rio Yokota and Jack Dongarra, “Distributed Memory Lattice H-matrix
Factorization”, The International Journal of High Performance Computing Applications (IHPCA) (2019)1094342019861139.
Takeshi Mifune, Naoki Tominaga, Yusuke Sogabe, Yudai Mizobata, Masahiro Yasunaga, Akihiro Ida, Takeshi
Iwashita and Naoyuki Amemiya, “Large-scale electromagnetic field analyses of coils wound with coatedconductors using a current-vector-potential formula-tion with a thin-strip approximation”, Superconductor
Science and Technology(SUST) .
Tetsuya Hoshino, Akihiro Ida, Toshihiro Hanawa and Kengo Nakajima,“Load-balancing-aware Parallel
Algorithms of H-matrices with Adaptive Cross Approximation for GPUs”, In Proceedings of IEEE International
Conference on Cluster Computing (CLUSTER 2018)64

HACApK関連
HACApK関連の査読付き
関連の査読付き論文リストの査読付き論文リスト
10. Tetsuya Hoshino, Akihiro Ida, Toshihiro Hanawa and Kengo Nakajima,“Design of Paral-lel BEM Analyses
Framework for SIMD Processors”, In: Shi Y. et al. (eds) Computational Science ? ICCS 2018. ICCS 2018.
Lecture Notes in Computer Science, vol 10860. Spring-er, Cham
11. Akihiro Ida, “Lattice H-Matrices on Distributed-Memory Systems”, 32nd IEEE Interna-tional Parallel &Distributed Processing Symposium (IPDPS 2018), Vancouver, Canada (May 21-25, 2018), pp.389-898.
12. Ichitaro Yamazaki, Ahmad Abdelfattah, Akihiro Ida, Satoshi Ohshima, Stanimire Tomov, Rio Yokota, and Jack
Dongarra “Analyzing Performance of BiCGStab with Hierarchical Matrix on GPU clusters”, 32nd IEEE
International Parallel & Distributed Processing Sym-posium (IPDPS 2018), Vancouver, Canada (May 21-25,
2018) ,pp930-939.
13. Naoki Tominaga, Takeshi Mifune, Akihiro Ida, Yusuke Sogabe, Takeshi Iwashita, Naoyuki Amemiya,“Application of hierarchical matrices to large-scale electromagnetic field analyses of coils wound with coatedconductors”, IEEE Transactions on Applied Su-perconductivity, Volume 28, Issue 3 (April 2018), pp.1-5.
14. Satoshi Ohshima, Ichitaro Yamazaki, Akihiro Ida and Rio Yokota, “Optimization of Hierarchical matrix
computation on GPU”, In: Asian Conference on Supercomputing Frontiers. Springer, Cham, 2018. p. 274292.
15. Akihiro Ida, Hiroshi Nakashima and Masatoshi Kawai, “Parallel Hierarchical Matrices with Block Low-rank
Representation on Distributed Memory Computer Systems”, International Conference on High Performance
Computing in Asia-Pacific Redion (HPC Asia 2018), Tokyo, Japan (Jan. 28-31, 2018), pp.232-240.
16. Akihiro Ida, Tadashi Ataka, Takeshi Mifune, Yasuhito Takahashi, Takeshi Iwashita and Tadashi Furuya,“Application of Improved H-matrices in Micromagnetic Simulations of Spin Torque Oscillator”, IEEE
Transactions on Magnetics, Vol. 54, Issue 3
17. Takeshi Iwashita, Akihiro Ida, Takeshi Mifune and Yasuhito Takahashi, "Software Framework for Parallel BEM
Analyses with H-matrices Using MPI and OpenMP", Procedia Computer Science 108C (2017) pp.2200–2209.
18. Akihiro Ida, Takeshi Iwashita, Takeshi Mifune and Yasuhito Takahashi, “Variable Preconditioning of Krylov
Subspace Methods for Hierarchical Matrices with Adaptive Cross Approximation,” IEEE Transactions on
Magnetics, Vol. 52, Issue 3
65

HACApK関連
HACApK関連の査読付き
関連の査読付き論文リストの査読付き論文リスト
19. Takeshi Iwashita, Naokazu Takemura, Akihiro Ida and Hiroshi Nakashima, “A new fill-in strategy for IC
factorization preconditioning considering SIMD instructions”, The 13-th IEEE International Symposium on
Parallel and Distributed Processing with Applications (IEEE ISPA-15), Helsinki Finland (2015)
20. Akihiro Ida, Takeshi Iwashita, Makiko Ohtani and Kazuro Hirahara, “Improvement of Hierarchical Matrices
with Adaptive Cross Approximation for Large-scale Simulation”, Journal of Information Processing, Vol. 23,
No. 3 (2015) pp．366-372．
21. Akihiro Ida, Takeshi Iwashita, Takeshi Mifune and Yasuhito Takahashi, “Parallel Hierarchical Matrices with
Adaptive Cross Approximation on Symmetric Multiprocessing Clusters”, Journal of Information Processing,Vol. 22, No. 4 (2014) pp．642-650．

66