PCクラスタコンソーシアムHPCオープンソースソフトウェア普及部会ワークショップ「高性能数値計算ライブラリLexADV」, 2021年9月10日

階層型領域分割法ライブラリ
LexADV_IsDDMの概要と並列反復法の実装
荻野正雄(大同大学)

自己紹介
• 所属
• 大同大学情報学部 准教授• 研究分野• 計算工学，高性能計算• プログラミング言語• C言語は20年以上，Fortranは読み書き程度，Pythonはそれなり• MPIとOpenMPはそれなり, Pthreadは読み書き程度• ソフトウェア開発• ADVENTUREプロジェクト: AdvSolid，AdvThermal• HDDMPPSプロジェクト: LexADV_IsDDM, LexADV_TryDDM, LexADV_WOVis• PC歴• 1997年に大学の授業でPCでのFortran 77プログラミングを学んでから．
• 1998年から研究室でSun OSやLinux OSに触れる機会を得る．その後メイン機種で使ってきたLinux OSは，Kondara → Red Hat → SuSE, openSUSE →Ubuntu• 2000年代は研究室でよくPCクラスタ(数十台規模)を作っていた．

階層型領域分割法ライブラリの概要• LexADV_IsDDMの概要• 特徴• 固体解析向けAdvSolidモジュールの線形方程式ソルバ部分を置換• Schur補元方程式を陰的に構築• AdvMetisによる領域分割メッシュに対応
高効率なDDMを実• C言語, MPI+OpenMPによるハイブリッド並列装したい人向け参照• 履歴
コード• Ver.0.1b: 2016年8月25日【非公開】• LexADV_TryDDMの概要
• 特徴
•••
•

領域分割法(DDM)に基づいて線形方程式を求解Schur補元方程式を陽的に構築
AdvMetisによる領域分割メッシュに対応
C言語, MPIによるプロセス並列+並列BLAS• 履歴

線形ソルバとしてのDDMの評価やDDM前処理の開発向け• Ver.0.1b: 2014年3月23日公開• Ver.0.2a: BDD-DIAG前処理法追加【公開前】このスライドはざっくりとした説明なので細かいことは考えないでください

領域分割法による並列有限要素解析(1)• 有限要素法(FEM)• 汎関数𝐴 𝑥 の最小化問題と考える
min 𝐴 𝑥 s.t. 𝑥の境界条件• 領域分割法(DDM)• 𝑁個の部分領域に重なりがないよう分割する• 汎関数𝐴 𝑥,…,𝐴 𝑥の最小化問題と考える
s.t. 𝑥 の境界条件と
min 𝐴 𝑥𝑥 の領域間境界条件
領域間境界条件は不明なので，適当な初期値から開始し，反復法で修正していく

メッシュの領域分割

領域分割法による並列有限要素解析(2)• DDMの利点• 既存のFEMコードを内包する形で開発できる• (反復法の選び方によるが)部分領域ごとの計算は並列化できる• 例えば，OpenMP並列化されたFEMコードがあれば，部分領域ごとの計算をMPIで分散並列処理する，など• DDMの欠点• FEMコードが数万行規模であることも多く，実際の開発は大変• 元のFEMコードの(データ構造などの)影響を受けるため，FEMコードごとに開発する必要がある• 反復法は元のFEMを考慮して構築するため，反復法の変更はDDMアルゴリズムの変更になる
ADVENTURE界隈の話を聞いてDDMを試したくても気軽にできない
ADVENTUREの中の人でもコード改修は大変な労力である

DDMを数値線形代数に帰着させる• 解くべき線形方程式

一般的な反復法ライブラリはこの式をそのまま解く• 𝑁個の領域分割で未知数を領域内部(𝐼)と領域間境界(𝐵)で並び替え
𝐴0𝐴
⋮⋱0𝐴

⋯𝐴𝐴𝐴

∑𝐴𝐛𝐱⋮

⋮𝐛𝐱𝐱∑𝐛• 式の分割• 未知数𝐱 に関するSchur補元方程式
• 𝑆

∑

𝐴

𝐴

𝐴𝐴

これを反復法で解く∑,𝐠𝐛• 未知数𝐱 に関する領域内部方程式

𝐴 𝐱𝐛

𝐴 𝐱 ,𝑖1, … , 𝑁𝐴𝐴𝐛

に対する反復法の実装• の成分を陽に計算しない陰的構築 (LexADV_IsDDM)• 行列ベクトル積𝐪𝑆𝐩の計算手順• step 1. 行列ベクトル積𝐭𝐴 𝐩, 𝑖1, … , 𝑁𝐭 ,𝑖1, … , 𝑁• step 2. 線形方程式求解 𝐴 𝐮• step 3. 行列ベクトル積𝐪∑ 𝐴𝐩𝐴 𝐮• 部分領域単位で並列処理
行列𝑆の成分が不明のため, 評価できる前処理法が限定される• の成分を陽に計算する陽的構築 (LexADV_TryDDM)• 最初に𝑆 𝐴
𝐴 𝐴
𝐴 を計算する• 行列ベクトル積𝐪 𝑆𝐩はそのまま計算する• 行列分散で並列処理ほとんどの反復法・前処理法が評価可能である

LexADV_IsDDMの構成と利用方法• ライブラリに含まれる主なファイル
••
•
•
•
•
•
•
•
•
•

hddmpps_matrix.c ... Schur補元行列とのベクトル積計算
hddmpps_matrix_expansion ... Schur補元からの展開
hddmpps_matrix_reduction.c ... Schur補元への縮約
hddmpps_solver_cg.c ... Schur補元方程式をCG法で解く
hddmpps_solver_cr.c ... Schur補元方程式をCR法で解く
hddmpps_solver_minres.c ... MINRES法で解く
hddmpps_precon_bdd_1.c ... BDD前処理法ステップ1
hddmpps_precon_bdd_2.c ... BDD前処理法ステップ2
hddmpps_precon_bdd_3.c ... BDD前処理法ステップ3
hddmpps_precon_bdd_4.c ... BDD前処理法ステップ4
hddmpps_vector.c ... ベクトル演算• 利用方法• 自前アプリのソースファイルと一緒にビルドする方式• いわゆるライブラリ化すると，行列データのコピー作成が必要になり，時間・空間的にもコスト増になってしまうため• DDMソルバ向けのデータ構造が複雑で，API設計がうまくいっていないため・・・

LexADV_IsDDMの使い方
hddmpps_setup_options(&argc, &argv);

コマンドライン引数を渡す...
hddmpps_set_mesh(mynparts, pmesh); 領域分割メッシュ情報を渡す
hddmpps_allocate_matrix();...

係数行列のためのメモリ割り当て#pragma omp parallel for schedule(dynamic,1)
for (idom = 0; idom < pmesh->n_domain; idom++) {/* making coef. matrix of sub-domain in gk *//* making RHS vector of sub-domain in gf */hddmpps_set_matrix(idom, &gk);
hddmpps_set_vector(idom, gf);

部分領域ごとに係数行列𝐴 と右辺ベクトル𝐛 を作成して渡す}
hddmpps_solver_cg(nparts, optsw, cgoption,
nr_norm, pmesh, pfield);𝐴𝐱 𝐛をDDMに基づいてCG法で解く

CG法ソルバhddmpps_solver_cg.cの中身...
初期残差𝐫𝑔 𝑆𝐱 を計算
hddmpps_reduction(ap_h, r_h);
hddmpps_vec_part_gather(p_hddm_npart[0], p_hddm_partid[0],p_hddm_outfree_num[0], p_hddm_outfree_list[0], r_h);...

ベクトル𝐫は分割されているので，袖領域通信が必要/* CG loop */while (1) {for (i = 0; i < n_total; i++)ap_h[i] = 0.0;
行列ベクトル積の𝐪 𝑆𝐩計算
hddmpps_matvec(p_h, ap_h);
hddmpps_vec_part_gather(p_hddm_npart[0], p_hddm_partid[0],p_hddm_outfree_num[0], p_hddm_outfree_list[0], ap_h);
myprod = 0.0;ベクトル𝐪の袖領域通信
for (i = 0; i < n; i++)myprod += p[i] * ap[i];COM_GrpSumDbl(&myprod, &pap, 1);alpha = rr / pap;

LexADV_IsDDMを導入したアプリの実行• 基本的には元々のアプリの実行方法通り• IsDDMオプションとして以下を用意している
Usage:
mpiexec [mpi_options] app_name -iterate [solver] -precon [pc]Options:solver: cg, cr, minres
pc: jacobi, nn, bdd, bdd-diag

他に, Pipelined CG法，ドット積を倍々精度計算，なども実装済み

DDM専用の前処理
nn: Neumann-Neumann [De Roeck & Le Tallec, 1991]
bdd: Balancing Domain Decomposition [Mandel, 1993]
bdd-diag: BDD w/ diagonal-scaling [Ogino et al., 2008]

LexADV_IsDDMの実行例(1)
線形弾性問題の有限要素解析，未知数約100万
計算機: 仮想マシン@Intel i7-8560U, Ubuntu OS

計算結果のミーゼス応力分布

Jacobi前処理: 8,446反復, 1.18時間
ADVENTURE界隈は限られた反復法と前処理法の結果しか示せていなかった (根拠があって選択した手法だが新規案件やトラブルへの対応は遅い)

BDD-DIAG前処理: 63反復, 28秒
BDD前処理: 44反復, 33秒

LexADV_IsDDMの実行例(2)

反復法(CG, CR, MINRES)ごとの収束傾向比較(このケースはCG法で十分と確認)更新残差𝐫 と真の残差𝐛 𝐴𝐱 のノルム比較(一致しており偽収束でないことを確認)

LexADV_IsDDMによって，数万行規模で書かれたDDMコードであるAdvSolidで反復法の選択が容易になり，様々な評価や検討が可能になった．さらに，例えば係数が非対称になるなど新たな反復法が必要な場合の開発コストも低くなる．

LexADV_TryDDMディレクトリ構成• ライブラリに含まれるディレクトリ構成• doc ... マニュアル• lib ... ビルドされたライブラリがコピーされる• samples ... サンプルデータ• src• common ... global/local共通のプログラムソース• common-include ... global/local共通のヘッダーファイル• global-schur-ver ... グローバル構築版• local-schur-ver ... ローカル構築版

TryDDM利用時の処理の流れ
Domaindecomposition

FEM pre-processing

FEM main-processing

メッシュの領域分割

AdvMetis

LexADV_TryDDM

Input

領域分割メッシュ読み込み

Input DD

AdvIO

𝐴, 𝐛

Assembly𝐴, 𝐛

Schur補元方程式の構築

Reduction

MUMPS

𝑆, 𝐠Solver𝐱

SolverOutput𝐱𝐱

Schur補元方程式の求解

Lis
領域内部方程式の求解

Expansion

𝐱

MUMPS

FEM post-processing

グローバル 構築版とローカル 構築版• グローバルSchur補元行列構築版(global-schur-ver)• 𝐴 → MUMPS → 𝑆 𝐴
𝐴 𝐴𝐴• MUMPSがマルチフロンタル法による並列行列分解と求解で構築• 長所: 全体係数行列とベクトルを入力するのみで実装が容易• 短所: 台数効果が低い• ローカルSchur補元行列構築版(local-schur-ver)• 𝐴∑

𝑅𝐴𝐴𝐴0𝑅→ MUMPS → 𝑆000 𝐴

→𝑆𝐴

𝐴𝐴𝐴

• MUMPSが局所Schur補行列を構築• 長所: 台数効果が高い• 短所: 実装が少し複雑になる∑𝑅𝑆 ′𝑅依存ライブラリ等• LexADV_TryDDM依存ライブラリ• 入出力ライブラリ AdvIO• ADVENTUREプロジェクト• 反復法ライブラリ Lis• SSIプロジェクト• 直接法ライブラリ MUMPS• PARASOLプロジェクト• BLAS, BLACS, ScaLAPACKが必要 (Intel MKL等でも良い)• その他• 領域分割メッシュ作成ツール AdvMetis• ADVENTUREプロジェクト• AdvIOが必要

LexADV_TryDDMのビルド例
download AdvIO-1.2.tar.gz from adventure projecttar xvfz AdvIO-1.21.tar.gz
cd AdvIO-1.2a
AdvIO./configure
make install
git clone https://github.com/anishida/lis.git
cd lisLis./configure --enable-mpi --enable-omp
make installsudo apt install libmumps-dev

MUMPStar xvfz LexADV_TryDDM-0.1b.tar.gz
cd src/local-schur-vervi Makefile
make

TryDDM

TryDDM公開版はLis最新版への対応で少し修正が必要

Makefileの編集#-----ADVENTURE-----------------ADVIO_DIR =(AdvIOインストールディレクトリ)

ADVIO_CONFIG = $(ADVIO_DIR)/bin/advsys-config
ADVIO_INCLUDES = `$(ADVIO_CONFIG) --cflags`
ADVIO_LIBS = `$(ADVIO_CONFIG) --libs docio`#-------MUMPS---------------------------------MUMPS_INCLUDES = -I$(MUMPSインクルードディレクトリ)-L$(MUMPSリンクディレクトリ)-ldmumps -lmumps_common -lpord-ldmumps -lmumps_common -lpord(ScaLAPACKへのリンク)= (BLASへのリンク)

MUMPS_LIBS
SCALAP_LIBS==

BLAS_LIBS#-------Lis---------------------------------LIS_INCLUDE = -I
-I$(Lisインクルードディレクトリ)
LIS_LIBS = -L
-L$(Lisリンクディレクトリ) -llis-I$(MPIインクルードディレクトリ)
MPI_INCLUDE = -I

#必要に応じてMPI_LIBS = -lmpi$(MPIへのリンク) #必要に応じて
DDM_INCLUDES = -I../common-include
CC = mpicc(MPI Cコンパイラ)

CFLAGS = -O3 -DUSE_MPI -fopenmp(コンパイラフラッグ)

INCLUDES = $(DDM_INCLUDES) $(ADVIO_INCLUDES) $(LIS_INCLUDE)$(MPI_INCLUDE)$(MUMPS_INCLUDES) ¥LIBS = $(ADVIO_LIBS) $(LIS_LIBS) $(MUMPS_LIBS) $(SCALAP_LIBS) $(BLAS_LIBS) $(MPI_LIBS) -lm
LINKER = mpicc(MPIリンカ)

PROGRAMS = libtry_ddm_gl.a
...

利用方法
int try_ddm(int argc, char** argv, int *csr_ptr, int*csr_idx, double *csr_val, int A_nnz, int A_matdim, intb_dim, double *b_val, char *mesh_file, char *mesh_type);

オプション指定
係数行列右辺ベクトル
領域分割メッシュ

引数説明

int argc

引数の数

char** argv

引数リスト

int* csr_ptr

各行開始位置リスト (CSR形式)

int* csr_idx

非零要素の列番号リスト (CSR形式)

double* csr_val

非零要素の値リスト (CSR形式)

int A_nnz

非零要素数

int A_matdim

行列次元数

int b_dim

ベクトル次元数

double* b_val

値リスト (dense形式)

char* mesh_file

領域分割ファイル名

char* mesh_type

メッシュタイプ (4面体2次: tetra)サンプルプログラムとデータ
ディレクトリ構造• samples/• Makefile ... メイクファイル• sample.c ... サンプルプログラム• data/• Ab/• Ab-130k.dat ... 係数行列と右辺ベクトル• Mesh/• adv-metis/• 130k/p1d240/model/advhddm_in_0.adv ... 領域分割メッシュ• origin-mesh/• 130k.adv ... 領域分割前メッシュ(AdvIO形式)

サンプル実行例
mpiexec -n 2 ./sample ¥反復法ソルバLisへのオプション指定-i cg -p ssor -maxiter 10000 -tol 1.0e-6 -print out ¥tetra ./data/Ab/Ab-130k.dat ¥./data/Mesh/adv-metis/130k/p1d240/model/advhddm_in_0.adv####################################

TryDDMへのオプション指定

Matrix A and Vector b
matdim:128814
nnz:9726120

入力した係数行列𝐴の情報

DDM Meshfile
SubDomains:240
Global Schur matdim:36219

Schur補元行列𝑆の情報...
iteration:1617

relative residual = 9.719264E-07

linear solver status: normal end

反復法停止時の情報...
log10(||b-Ax_n|| / ||b||) = -6.040980e+00

得られた近似解の誤差[0]*** All elapsed time :82.436755 ***全体の計算時間

IsDDMとTryDDMを用いた前処理の性能評価例

BDDはILU(0)前処理などに比べても反復回数が少ないといった評価が可能に

BDD ILU ILUT
BDD-DIAG

SSOR(𝜔1)

Jacobi

領域分割数を変えたサンプル実行 (1)• 同梱サンプルは領域分割数240 (p1d240)• 領域分割数を変えた実験を行う例• AdvMetisの実行例(領域分割数120)
adventure_metis origin-mesh/130k.adv adv-metis/p1d120 120• サンプルプログラムの実行
mpiexec -n 2 ./sample ¥-i cg -p ssor -maxiter 10000 -tol 1.0e-6 -print out ¥tetra ../Ab/Ab-130k.dat ¥Mesh/adv-metis/130k/p1d120/model/advhddm_in_0.adv####################################...
DDMでは領域分割数を変えるとDDM Meshfile
SubDomains:120
Schur補元行列のサイズも変わり，そGlobal Schur matdim:25038
れが反復回数や計算時間に影響する####################################...
iteration: 1333 relative residual = 9.843737E-07
linear solver status : normal end...
[0]*** All elapsed time :60.299284 ***領域分割数を変えたサンプル実行 (2)BDDやILUはSSORに比べて領域分割数に対してロバストだな，などが分かるようになった

SSOR (TryDDM)

ILU (TryDDM)

BDD (IsDDM)

領域分割数 vs. 計算時間

LexADV_TryDDMの応用例
静磁場問題のBDD-DIAG前処理開発• 静磁場の有限要素解析における線形方程式• 磁気ベクトルポテンシャルを未知数とするA法• BDD前処理フレームワークの適用• Coarse空間の基底 𝑍

1 (BDDの理論は満たさない)

前処理なし
BDD-DIAG-like前処理
An axial model

SSOR
Jacobi前処理

まとめと今後の計画• 階層型領域分割ライブラリの概要• LexADV_IsDDMは，• 高性能・高機能だが線形ソルバ変更が困難なADVENTURE
モジュールにおいて，反復法の変更を容易に実現するもの• API設計が完了しておらず，現時点ではAdvSolidモジュールに直接組み込む形での利用となり，非公開扱いである• 名大「不老」やR-CCS「富岳」などでの最先端研究でも使われている• LexADV_TryDDMは，• DDM向けの前処理法開発を支援するもの• 自前アプリの線形方程式におけるDDM性能の予備調査にも• 今後の計画• LexADV_IsDDMの公開は，AdvSolid同梱という形で検討• LexADV_TryDDMはVer.0.2bを公開予定