MPIによる並列アプリケーション
開発法入門（III）
2011年5月19日・20日
2011年5月19日
20日
中島研吾
東京大学情報基盤センター
T2Kオープンスパコン（東大）並列プログラミング講習会

T2K-FVM-03

1

方 針• 〔II〕で定義した局所分散データ構造• MPIの処理をできるだけ「隠蔽」– 初期化等環境設定– 通信• hpcmw_eps_fvm_・・・という関数名hfという関数名– HPC-MW（HPC Middlewareに由来）• マルチフィジックスシミュレーション向け大規模並列計算コード開発基盤
ジ
向 大規模並 計算
ド
基盤– 並列アプリケーションにおける並列処理のソフトウェア的な隠蔽• http://hpcmw.tokyo.rist.or.jp/
htt //ht k
i t j /

T2K-FVM-03

2

マルチフィジックスシミュレーション向け
マルチフィジックスシミ
レ ション向け
大規模並列計算コード開発基盤• マルチフィジックス・マルチスケールの並列大規模連成シミュ
レ ションを円滑に実施するためのフレ ムワ ク
レーションを円滑に実施するためのフレームワーク• 連立一次方程式ソルバー，形状処理，可視化，コード間連成などの共通処理に関する機能を提供し 並列大規模連成シなどの共通処理に関する機能を提供し，並列大規模連成シ
ミュレーションコードを開発するための支援環境• PC上で開発された個別のプログラムを「plug-in」すれば，PC
C上 開発された個別 プ グ ムを「すれば Cクラスタから「地球シミュレータ」，「ペタスケール計算機」まで
様 な
様々なハードウェアに対して最適化された並列プログラムが自ドウ
対最適化された並列プ グ ムが自動的に生成される・・・というのが理想– HPC-MW（Middleware），HEC-MWT2K-FVM-03

3

大規模並列計算コード開発基盤
これまで関連して来たプロジェクト
れま 関連
来たプ ジ ク
• GeoFEM（FY.1998-FY.2002）– 固体地球シミュレーション用並列有限要素法プラットフォーム• 地球シミュレータ• HPC-MW（FY.2002-FY.2003）（プロジェクトはFY.2007迄）– 文部科学省ITプログラム「戦略的基盤ソフトウェアの開発」• RSS21「革新的シミュレーションソフトウェアの研究開発（FY.2005-）」• 観測・計算を融合した階層連結地震・津波予測システム（FY.2005-FY.2010）– 科学技術振興機構戦略的創造研究推進事業（CREST）• 「階層連結シミュレーション」，「連成」を重視• 自動チューニング機構を有するアプリケーション開発・実行環
境 ppOpen-HPC（FY.2011-FY.2015（予定））O
HPC（FY 2011 FY 2015（予定））– 科学技術振興機構戦略的創造研究推進事業（CREST）

T2K-FVM-03

GeoFEM: FY.1998
FY.1998-2002
2002
http://geofem.tokyo.rist.or.jp/• 文部科学省「科学技術振興調整費総合研究」– 「高精度の地球変動予測のための並列ソフトウェア開発に関する研究」の一部– リーダー：奥田洋司教授（東大・人工物）
ダ
奥 洋 教
東
物

• 固体地球シミュレーション用並列有限要素法プラットフォーム– 並列I/O，並列線形ソルバー，並列可視化をサポート• HPCと自然科学の緊密な協力

4

T2K-FVM-03

5

GeoFEM 「plug in」の発想
GeoFEM：「plug-in」の発想
Utilities
One-domain mesh

Partitioner

Pluggable Analysis Modules
Structure
Fluid
Wave

構造計算（Static linear）
構造計算（Dynamic
構造計算（
linear）
Contact）

Comm.
I/F
Parallel
I/O

Solver
I/F

Vis.
I/F

Platform

Equation Visualizer
solvers

P titi
Partitioned
d meshhPEs
GPPView

Visualization data

T2K-FVM-03

HPC-MW
PCクラスタから「地球シミュレータ」まで• GeoFEMにおける「Plug-in」のアイディアに基づく
GeoFEMにおける「Plug in」のアイディアに基づく• 科学技術計算（特に有限要素法）における共通プロセスの最適化チ適化チューニング，並列化等のライブラリ化による隠蔽
ング 並列化等のライブラリ化による隠蔽• 「HPC-MW」によって，PC上で開発されたコードを，PCクラス
タから「ES」まで，様々なハードウェアで最適な性能によって
「
様 な
ド
適な性能
稼動させることができる– GeoFEMで「ES」向け最適化をやったのが契機– 当時は「ベクトル」，「スカラー」が並立しそうな気配もあり• HPC-MWはGeoFEMと比較して，より広い機能をカバーし，さらに各ハードウェアへの最適化の考えを導入している

6

T2K-FVM-03

7

並列有限要素法のプロセス

Pre--Processing
Pre

Main

Post--Processing
Post

Initial Grid Data

Data Input/Output

Post Proc.
Proc.

Partitioning

Matrix Assemble

Visualization

Linear Solvers
Domain Specificp
Algorithms/Models

T2K-FVM-03

HPC-MWのサポ トする機能
HPC-MWのサポートする機能• データ入出力• 適応格子，動的負荷分散• 並列可視化• 線形ソルバー
線形ソルバ• 有限要素処理（コネクティビティ処理，係数行列生成）• カップリング関連
ティリティ（領域分割等）• 関連ユーティリティ（領域分割等）

8

T2K-FVM-03

9

HPC MWの利用イメ ジ
HPC-MWの利用イメージ
PC（単独CPU）上で開発されたFEMコード
FEM code developed on PC
I/OMatrixAssemble

LinearSolverVis.

T2K-FVM-03

10

HPC MWの利用イメ ジ
HPC-MWの利用イメージ
HPC-MWを使用する場合，共通部分は不要
FEM code developed on PC

T2K-FVM-03

11

HPC MWの利用イメ ジ
HPC-MWの利用イメージ
各H/W用に最適化されたライブラリ
FEM code developed on PC

I/OMatrixM
ti
Assemble

LinearLi
SolverVis.

HPC-MW for T2K
I/O

I/O

MatrixAssemble

LinearSolverMatrixAssemble

LinearSolverVis.

HPC-MW for Next Generation
SupercomputerVis.

HPC-MW for Earth Simulator

T2K-FVM-03

12

HPC MWの利用イメ ジ
HPC-MWの利用イメージ
各ライブラリに対して同じインタフェース
FEM code developed on PC
I/F for
I/O

I/O

MatrixM
ti
Assemble

I/F forMat.Ass.

LinearLi
SolverI/F forSolvers

I/F forVis.

I/O

MatrixAssemble

Vis.

HPC-MW for T2K

I/O

MatrixAssemble

LinearSolverLinearSolverVis.

HPC-MW for Next Generation
SupercomputerVis.

HPC-MW for Earth Simulator

T2K-FVM-03

13

HPC MWの利用イメ ジ
HPC-MWの利用イメージ「地球シミュレータ」用最適化コード
FEM code developed on PC
I/F for
I/O
I/O

I/F forMat.Ass.
MatrixAssemble

I/F forSolvers
LinearSolverI/F for
Vis.
Vis.

HPC-MW for T2K

I/O

I/O

MatrixAssemble

LinearSolverMatrixAssemble

LinearSolverVis.

HPC-MW for Next Generation
SupercomputerVis.

HPC-MW for Earth Simulator

T2K-FVM-03

14

HPC-MWの問題点• 新しいアーキテクチャが登場するごとに，最適化（チューニング）のやりなおしが必要となる場合がある。

T2K-FVM-03

15

ppOpen-HPC （FY.2011-）
http://ppopenhpc.cc.u-tokyo.ac.jp/• 各手法の限定されたプロセスに特化• メモリアクセス最適化に照準をしぼった適応的な自動チュー
ニング（A tニング（Automaticti Tuning，Auto-Tuning）
T i
A t T i ）User’s Program

ppOpen-APPL

FEM

FDMii

FVM

BEM

ppOpen-MATH
O
MATH

MG

ii
GRAPH

VIS

MP

DEM

FEM

Finite Element Method
ppOpen-AT

STATICii
DYNAMIC

ppOpen-SYS

COMM

ii
FT

FDM

Finite Difference Method

FVM

Finite Volume Method

ppOpen-HPC

Optimized Application with
Optimized ppOpen-APPL, ppOpen-MATH

BEM

Boundary Element Method

DEM

Discrete Element Method

T2K-FVM-03

16

hpcmw eps fvm
hpcmw_eps_fvm
MPI ctrl.

I/O

LinearSolverVis.

HPC-MW for FVM using MPI• 今回は，このようなミドルウェア的な機能も想定して並列プ
ログラムを開発している– ミドルウェアとして切り離し，他の並列アプリケーションにも使いまわせるような機能• SMASHのSHをカバー• hpcmw_eps_fvmで始まる関数群– MPI関連– I/O

ScienceModeling
Algorithm
SoftwareHardware

T2K-FVM-03

17

プログラム類のインスト ル
プログラム類のインストール• ソースコード– <$FVM>/src– ここで「make」すると「<$FVM>/run」 に「sol」という実行形式ができる。以下これを使用する。
きる 以
れを使用する$> cd
d <$FVM>/<$FVM>/src$> make$> ls -ll ../run/sol/run/sol
sol• チュートリアルチ
トリアル
http://nkl.cc.u-tokyo.ac.jp/tutorial/parallel_lib_tutorial/
http://nkl.cc.u-tokyo.ac.jp/tutorial/parallel_lib_tutorial.tar

T2K-FVM-03

18

ファイル準備$> cd <$FVM>/ex$> cat fvmmg.ctrl
32 32 32$> cat fvmpart.ctrl!INITIAL FILE
fvm_entire_mesh.dat!METHODRCB
X,Y,Z!REGION NUMBER
8!MESH FILE
mesh.rcb!COMMUNICATION FILE
comm.rcb!UCD
32 32 32
32-32-32-rcb-8.inp
b 8 i$> eps_fvm_mg$> eps_fvm_part$> ls -l mesh.rcb.*
mesh.rcb.0 … mesh.rcb.7$> ls -ll comm.rcb.*
comm.rcb.
comm.rcb.0 … comm.rcb.7

T2K-FVM-03

19

並列計算制御ファイル• INPUT.DAT（名称固定）
INPUT DAT（名称固定）• 実行形式「sol」と同じディレクトリになければならない（この場合は場合は<$FVM>/run）。
）• 全ての項目は省略不可。
../ex/mesh.rcb 局所分散メッシュファイルのヘッダ名../ex/comm.rcb 局所分散通信ファイルのヘッダ名../ex/result
可視化用出力ファイル名（後述）1
可視化用出力の有無（=1のとき出力）視有無きT2K-FVM-03

20

並列シミュレ ションにおけるI/O
並列シミュレーションにおけるI/O
領域分割機能割 能<$FVM>/run/
INPUT.DAT（固定）並列計算制御ファイル#GRIDout#GRIDout#GRIDout#GRIDout#D-GRID

局所分散メッシュデータ#GRIDout#GRIDout#GRIDout#GRIDout#D-COMM

局所分散通信データ

並列計算
sol
so

#M-RESULT

計算結果ファイル
（UCDファイル後半）

T2K-FVM-03

21

計算実行
$> cd <$FVM>/run$> cat INPUT
INPUT.DATDAT../ex/mesh.rcb../ex/comm.rcb../ex/result
1「go.sh」を書き換える$> qsub go.sh$> ls -l ../ex/result
resultT2K-FVM-03

22

「eps
eps_fvm」の並列化：変更点：test.f
fvm」の並列化：変更点：test.f
program eps_fvm
use hpcmw
p_eps
p _fvm_all
implicit REAL*8 (A-H,O-Z)
call hpcmw_eps_fvm_init
call hpcmw_eps_fvm_input_grid
call poi_gen
call hpcmw_eps_fvm_solver
call output_ucd
call hpcmw_eps_fvm_finalize
end program eps_fvm• 実はほとんど無い• 通信関連ファイル読み込み• 内点，外点– 内点– 内点＋外点

intNODE_tot
O
NODE_tot• MPIコールはできるだけ隠蔽– 初期化，Finalize• hpcmw_eps_fvm_util.＊– 通信用サブルーチン群• hpcmw_eps_fvm_comm.＊

T2K-FVM-03

23

program eps_fvm
use hpcmw_eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call hpcmw_eps_fvm_init
call hpcmw_eps_fvm_input_grid____
call poi_gen
call hpcmw_eps_fvm_solver___
call output_ucd
call hpcmw_eps_fvm_finalize
end program eps_fvm

T2K-FVM-03

24

変数ブロック：hpcmw eps fvm all
変数ブロック：hpcmw_eps_fvm_all!C!C***!C*** hpcmw_eps_fvm_all
hpcmw eps fvm all!C***!Cmodule hpcmw_eps_fvm_all
use hpcmw_eps_fvm_util
use hpcmw_eps_fvm_pcg
hf
use appl_cntl
end module hpcmw_eps_fvm_all
hpcmw eps fvm all

T2K-FVM-03

25

module hpcmw_eps_fvm_util
hpcmw eps fvm util
• 変数ブロック
変数ブロ ク– メッシュ
– 通信• MPI初期化・終了等に関連したサブルーチン群– hpcmw_eps_fvm_init• MPI_Init– hpcmw_eps_fvm_finalize• MPI_Finalize– hpcmw_eps_fvm_abortf• MPI_Abort– hpcmw_eps_fvm_define_file_name
h
f
d fi
fil
• 分散ファイル名定義

T2K-FVM-03

26

hpcmw eps fvm util（1/3）メッシュ関連
hpcmw_eps_fvm_util（1/3）メッシュ関連
変数名

型

配列サイズ

内 容

NODE_tot

I-

内点数＋外点数

intNODE_tot

I

-

内点数

NODE_GLOBAL(:)

INODE_tot

グローバル要素番号

NODE_VOL(:)

R

NODE_tot

要素体積

NODE_COND(:)

R

NODE_tot

要素熱伝導率NODE_XYZ(:)

R

3*NODE_tot

要素重心座標（3次元）

CONN_totI-

コネクティビティ総数

CONN_node(:)

I

2*CONN_tot

コネクティビティ構成要素ビ
構成 素

CONN_COEF(:)

R

CONN_tot

コネクティビティ係数

FIX_NODE_tot

I-

ディリクレ境界条件適用要素数

FIX NODE ID( )FIX_NODE_ID(:)

I

FIX NODE t t
FIX_NODE_tot

ディリクレ境界条件適用要素番号

FIX_NODE_COEF(:)

R

FIX_NODE_tot

ディリクレ境界条件係数

FIX_NODE_VAL(:)

R

FIX_NODE_tot

ディリクレ境界条件値

SURF NODE tot
SURF_NODE_tot

I-

ノイマン境界条件適用要素数

SURF_NODE_ID(:)

I

SURF_NODE_tot

ノイマン境界条件適用要素番号

SURF_NODE_FLUX(:)

R

SURF_NODE_tot

ノイマン境界条件フラックス

BODY NODE tot
BODY_NODE_tot

I-

体積発熱境界条件適用要素数

BODY_NODE_ID(:)

I

BODY_NODE_tot

体積発熱境界条件適用要素番号

BODY_NODE_FLUX(:)

R

BODY_NODE_tot

体積発熱境界条件フラックス

T2K-FVM-03

27

hpcmw_eps_fvm_util（2/3）通信関連
p
_ p _
_
変数名

型

配列サイズ

内 容

PETOT

I

-

プロセッサ数

errno

I

-

エラーもどり値

my_rank

I

-

ランク番号

g
_p
pe
n_neighbor

I

-

隣接領域数

neighbor_pe(:)

In_neighbor_pe

隣接領域ID

import_index(:)

I

0:n_neighbor_pe

受信テーブル用インデックス

import_item(:)

I

import_index(n_neighbor_pe)

export_index(:)

I

0:n_neighbor_pe

export_item(:)

I

export_index(n_neighbor_pe)

HPCMW_NAME_LEN

I-

NAME lengthパラメータ（＝63）

HPCMW_HEADER_LEN I

I-

ヘッダー長さパラメータ（＝127）

HPCMW_MSG_LEN

I-

メッセージ長さパラメータ（＝255）

HPCMW_FILNAME_LEN

I-

ファイル名長さパラメータ（＝1023）

hpcmw_eps_fvm_files(:)

C4

分散ファイル名，（1）メッシュファイル，（2）結果ファイル，（4）通信ファイル

受信テーブル
送信テーブル用インデックス
送信テーブルT2K-FVM-03

28

hpcmw_eps_fvm_util（3/3）
p
_ p _
_
MPI用パラメータ
変数名型

パラメータ値

内 容

hpcmw_sum

I

46801

MPI_SUM

hpcmw_prod

I

46802

MPI_PROD

hpcmw_max

I

46803

MPI_MAX

hpcmw_min

I

46804

MPI_MIN

hpcmw_integer

I

53951

MPI_INTEGER

hpcmw_single_precision

I

53952

MPI_SINGLE_PRECISION

hpcmw_double_presision

I

53953

MPI_DOUBLE_PRECISION

hpcmw_character

I

53954

MPI_CHARACTERT2K-FVM-03

29

マトリクス関連（hpcmw eps fvm pcg）
マトリクス関連（hpcmw_eps_fvm_pcg）
変数名

型

サイズ

内容
連立一次方程式係数マトリクス非対角成分総数

NPLU

I

-

D(:)

R

NODE_tot

連立一次方程式係数マトリクス対角成分

PHI(:)

R

NODE_tot

連立一次方程式未知数ベクトル

BFORCE(:)

R

NODE_tot

連立一次方程式右辺ベクトル

index(:)

I

0:NODE tot
0:NODE_tot

係数マトリクス非対角成分要素番号用一次元圧縮配列（非対角成分数）
係数マトリクス非対角成分要素番号用
次元圧縮配列（非対角成分数）

item(:)I

NPLU

係数マトリクス非対角成分要素番号用一次元圧縮配列（非対角成分要素番号）

AMAT(:)R

NPLU

係数マトリクス非対角成分要素番号用一次元圧縮配列（非対角成分）
係数マトリクス非対角成分要素番号用
次元圧縮配列（非対角成分）

do i= 1, N
q(i)=(i) D(i)*p(i)(i)* (i)do k= index(i-1)+1, index(i)q(i)= q(i) + AMAT(k)*p(item(k))
enddo
dd
enddo

T2K-FVM-03

30

変数ブロック：appl cntl
変数ブロック：appl_cntl!C!C***!C*** appl_cntl!C***!C
module appl_cntl
use hpcmw_eps_fvm_util!C!C-- FILE NAME
character(len=HPCMW_HEADER_LEN):: HEADERgrid
character(len=HPCMW_HEADER_LEN):: HEADERresult
character(len=HPCMW_HEADER_LEN):: HEADERcomm
character(len=HPCMW_HEADER_LEN):: AVSfile!C!C-- MESH info.
integer(kind=kint) :: NX, NY, NZ, NXP1, NYP1, NZP1
integer(kind=kint) :: PVISFLAG（略）
end module appl_cntl

T2K-FVM-03

31

program eps_fvm
use hpcmw_eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call hpcmw_eps_fvm_init
call hpcmw_eps_fvm_input_grid____
call poi_gen
call hpcmw_eps_fvm_solver___
call output_ucd
call hpcmw_eps_fvm_finalize
end program eps_fvm

T2K-FVM-03

32

初期化（hpcmw eps fvm util f）
初期化（hpcmw_eps_fvm_util.f）!C!C***!C*** HPCMW_EPS_FVM_INIT
!C***
!C
!C
INIT. HPCMW-FEM process's!Csubroutine HPCMW_EPS_FVM_INIT
integerg:: ierrcall MPI_INIT(ierr)
call MPI_COMM_SIZE (MPI_COMM_WORLD, PETOT, ierr)call MPI_COMM_RANK (MPI_COMM_WORLD,, myy_rank , ierr)
if (my_rank.eq.0) thenzero= 1else
zero= 0endif
end subroutine hpcmw
p_eps
p _fvm_init

T2K-FVM-03

33

Finalize （hpcmw_eps_fvm_util.f）（hpcmw eps fvm util f）!C!C***!C*** HPCMW_EPS_FVM_FINALIZE!C***!Csubroutine HPCMW_EPS_FVM_FINALIZE
integer :: ierrcall MPI_FINALIZE ((ierr))
if (my_rank.eq.0) stop ' * normal termination'end subroutine hpcmw_eps_fvm_finalize

T2K-FVM-03

34

Abort （hpcmw_eps_fvm_util.f）（hpcmw eps fvm util f）!C!C***!C*** HPCMW_EPS_FVM_ABORT!C***!C
subroutine HPCMW_EPS_FVM_ABORT
integer :: ierrcall MPI_BARRIER ((MPI_COMM_WORLD,, ierr))call MPI_ABORT(MPI_COMM_WORLD, ierr)
end subroutine hpcmw_eps_fvm_abort

T2K-FVM-03

35

program eps_fvm
use hpcmw_eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call hpcmw_eps_fvm_init
call hpcmw_eps_fvm_input_grid____
call poi_gen
call hpcmw_eps_fvm_solver___
call output_ucd
call hpcmw_eps_fvm_finalize
end program eps_fvm

T2K-FVM-03

36

デ タ入力部分（1/5）
データ入力部分（1/5）
subroutine hpcmw_eps_fvm_input_grid
use hpcmw_eps_fvm_all
implicit REAL*8 (A(A-H,O-Z)
H O Z)
character(len=HPCMW_NAME_LEN) :: member
character(len=80) :: LINE!C!C +-------+!C | FILES |!C +-------+
制御ファイル 「INPUT.DAT」!C===
open (11,(11 file='INPUT.DAT',
file='INPUT DAT' status='unknown')
read (11,'(a127)') HEADERgrid
read (11,'(a127)') HEADERcomm
read (11,‘(a127)’) AVSfile
read (11,*)
PVISFLAG
close (11)allocate (hpcmw_eps_fvm_files(4))
member= 'gridfile'call hpcmw_eps_fvm_define_file_name
hpcmw eps fvm define file name (member,(member HEADERgrid)#D-GRID

member= 'commfile'call hpcmw_eps_fvm_define_file_name (member, HEADERcomm)#D-COMM!C===

ヘッダーを与えて分散ファイル名を生成する。

T2K-FVM-03

37

並列シミュレ ションにおけるI/O
並列シミュレーションにおけるI/O
領域分割機能割 能<$FVM>/run/
INPUT.DAT（固定）並列計算制御ファイル#GRIDout#GRIDout#GRIDout#GRIDout#D-GRID

局所分散メッシュデータ#GRIDout#GRIDout#GRIDout#GRIDout#D-COMM

局所分散通信データ

並列計算
sol
so

#M-RESULT

計算結果ファイル
（UCDファイル後半）

T2K-FVM-03

局所分散ファイル名生成（1/2）（h（hpcmw_eps_fvm_util.f）
f
til f）!C!C***!C*** HPCMW_EPS_FVM_DEFINE_FILE_NAME!C***!C, HEADERo)
subroutine HPCMW_EPS_FVM_DEFINE_FILE_NAME (member,character (len=HPCMW_HEADER_LEN) :: HEADERo, FILENAME
character (len=HPCMW_NAME_LEN):: membercharacter (len=HPCMW_NAME_LEN):: HEADER
character (len= 1) :: SUBindex1
character (len= 2) :: SUBindex2
character (len= 3) :: SUBindex3
character (len= 4) :: SUBindex4
character (len= 5) :: SUBindex5
character (len= 6) :: SUBindex6
integer:: LENGTH, ID
HEADER= adjustL (HEADERo)LENGTH= len_trim(HEADER)38

T2K-FVM-03

局所分散ファイル名生成（2/2）（h（hpcmw_eps_fvm_util.f）
f
til f）
if (my_rank.le.9) thenID= 1write(SUBindex1 ,'(i1.1)') my_rank
else if (my_rank.le.99) thenID= 2write(SUBindex2 ,'(i2.2)') my_rank
else if (my(my_rank.le.999)
rank le 999) thenID= 3write(SUBindex3 ,'(i3.3)') my_rank
else if (my_rank.le.9999) thenID= 4write(SUBindex5 ,'(i4(i4.4)
4)')) mymy_rankrankelse if (my_rank.le.99999) thenID= 5write(SUBindex6 ,'(i5.5)') my_rank
else if (my_rank.le.999999) thenID= 6IDwrite(SUBindex4 ,'(i6.6)') my_rank
endif

領域数1,000,000まで対応可能

if (ID.eq.1) filename= HEADER(1:LENGTH)//'.'//SUBindex1
if (ID.eq.2)(ID eq 2) filename
filename= HEADER(1:LENGTH)//HEADER(1:LENGTH)//'.'//SUBindex2//SUBindex2
if (ID.eq.3) filename= HEADER(1:LENGTH)//'.'//SUBindex3
if (ID.eq.4) filename= HEADER(1:LENGTH)//'.'//SUBindex4
if (ID.eq.5) filename= HEADER(1:LENGTH)//'.'//SUBindex5
if (ID.eq.6) filename= HEADER(1:LENGTH)//'.'//SUBindex6
if (member.eq. 'gridfile') hpcmw_eps_fvm_files(1)= filename
if (member.eq. 'commfile') hpcmw_eps_fvm_files(4)= filename
end subroutine hpcmw_eps_fvm_define_file_name39T2K-FVM-03

40

デ タ入力部分（2/5）
データ入力部分（2/5）!C!C +------------+!C | MESH INPUT |!C ++------------++!C===
IUNIT= 11
open (IUNIT,file= hpcmw_eps_fvm_files(1), status='unknown')#D-GRID, ここらへんはserialと同じ各PEでフ イル名を書き出して見よ
各PEでファイル名を書き出して見よ!C!C-- NODE
read (IUNIT, '(10i10)') NODE_tot
allocate (NODE_VOL(NODE_tot), NODE_COND(NODE_tot),&
NODE
NODE_XYZ(3*NODE_tot))
XYZ(3*NODE tot))&

do i= 1, NODE_tot
read (IUNIT,'(i10,5e16.6)') ii, NODE_VOL(i), NODE_COND(i),&(NODE_XYZ(3*i-3+k), k=1, 3)
enddo

&

!C
!C-- CONNECTION
read (IUNIT,'(10i10)') CONN_tot
allocate (CONN_NODE(2(CONN NODE(2*CONNCONN_tot),tot) CONNCONN_COEF(CONN_tot))
COEF(CONN tot))
do i= 1, CONN_tot
read (IUNIT,'( 2i10, 3e16.6)') (CONN_NODE(2*i-2+k), k= 1, 2),&
AREA, D1, D2
in1= CONN_NODE(2*i-1)
in2= CONN_NODE(2CONN NODE(2*ii )C1 = NODE_COND(in1)
C2 = NODE_COND(in2)
CONN_COEF(i)= AREA / ( D1/C1 + D2/C2 )
enddo

&

T2K-FVM-03

41

有限体積法：隣接メッシュの情報が必要
熱流束に関するつりあい式
隣接要素との熱伝導

∑d
k

Sik
+

ik

λi

d ki

温度固定境界(Tk − Ti ) + ∑
e

λk

∑

Sie
(
TiBe − Ti ) +Sid q&id + Vi Q& i = 0
d ie
d

λi

要素境界面
通過熱流束

体積発熱

a
b
dai
dbi

Sib
Sia

dia

die

TiBe

Sie
Sid

dib

i

dic
dci

Sic
qid

c

λ ：熱伝導率
Vi ：要素体積
S ：表面面積
表面面積
dij ：要素中心から表面までの距離
q ：表面フラックス
Q ：体積発熱
体積発熱
TiB ：境界温度

T2K-FVM-03

42

デ タ入力部分（3/5）
データ入力部分（3/5）!C!C-- DIRICHLET
read (IUNIT,'(10i10)') FIX_NODE_tot
allocate (FIX
(FIX_NODE_ID(FIX_NODE_tot),
NODE ID(FIX NODE tot) FIX
FIX_NODE_COEF(FIX_NODE_tot))
NODE COEF(FIX NODE tot))
allocate (FIX_NODE_VAL(FIX_NODE_tot))
do i= 1, FIX_NODE_tot
read (IUNIT, '(i10, 3e16.6)')&
FIX
FIX_NODE_ID(i),NODE ID(i) AREA,
AREA DIST
DIST, FIXFIX_NODE_VAL(i)NODE VAL(i)icel= FIX_NODE_ID(i)COND= NODE_COND(icel)
FIX_NODE_COEF(i)= AREA / (DIST/COND)
enddo!C!C-- NEUMANN
read (IUNIT,'(10i10)') SURF_NODE_tot
allocate& (SURF_NODE_ID (SURF_NODE_tot), SURF_NODE_FLUX(SURF_NODE_tot))
do i= 1, SURF_NODE_tot
read (IUNIT, '(i10, 3e16.6)') SURF_NODE_ID(i), AREA, FLUXSURF_NODE_FLUX(i)= AREA*FLUX
enddo!C!C-- BODY FLUX
read (IUNIT,'(10i10)') BODY_NODE_tot
allocate (BODY_NODE_FLUX(NODE_tot))
do i= 1,
1 BODY_NODE_tot
BODY NODE tot
read (IUNIT, '(i10, 3e16.6)') icel, FLUXBODY_NODE_FLUX(icel)= FLUX * NODE_VOL(icel)
enddo
close (IUNIT)&&T2K-FVM-03

43

局所分散メッシュファイル（mesh.0）
要
素

コ
ネ
ク
テ
ィ
ビ
テ
ィ
境
界
条
件

8
1
2
3
4
5
6
7
8
8
5
6
7
1
1
2
8
3
2
2
4
0
1
1

1.000000E+00
1.000000E+00
11.000000E+00
000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
11.000000E+00
000000 00
1
2
1
2
3
4
3
4

1.000000E+00
1.000000E+00
1.000000E+00
1 000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1 000000 00

1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
11.000000E+00
000000E+00

1.000000E+00
1.000000E+00

1.000000E+00

2.500000E+00
3.500000E+00
2.500000E+00
2 500000E+00
3.500000E+00
2.500000E+00
3.500000E+00
1.500000E+00
1.500000E+00
1 500000 00

5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
55.000000E
000000E-01
01

5.000000E-01
5.000000E-01

2.500000E+00
2.500000E+00
33.500000E+00
500000E+00
3.500000E+00
1.500000E+00
1.500000E+00
2.500000E+00
33.500000E+00
500000 00

5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E
5 000000E-01
01

0.000000E+00
0.000000E+00

5.000000E-01
5.000000E-01
55.000000E
000000E-01
01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
55.000000E-01
000000 01

PE#0
14

15

16

8

3

4

10

11

12

7

1

2

PE#1

7
8
基本的に初期全体メッシュファイルと同じ
5
6
局所要素番号による記述
境界条件（ディリクレ，ノイマン，体積発熱）：「内点」のみの情報 PE#2

T2K-FVM-03

44

局所分散メッシュファイル（mesh.0）
要
素

コ
ネ
ク
テ
ィ
ビ
テ
ィ
境
界
条
件

8
1
2
3
4
5
6
7
8
8
5
6
7
1
1
2
8
3
2
2
4
0
1
1

1.000000E+00
1.000000E+00
11.000000E+00
000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
11.000000E+00
000000 00
1
2
1
2
3
4
3
4

1.000000E+00
1.000000E+00
1.000000E+00
1 000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1 000000 00

1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
11.000000E+00
000000E+00

1.000000E+00
1.000000E+00

1.000000E+00

2.500000E+00
3.500000E+00
2.500000E+00
2 500000E+00
3.500000E+00
2.500000E+00
3.500000E+00
1.500000E+00
1.500000E+00
1 500000 00

5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
55.000000E
000000E-01
01

5.000000E-01
5.000000E-01

2.500000E+00
2.500000E+00
33.500000E+00
500000E+00
3.500000E+00
1.500000E+00
1.500000E+00
2.500000E+00
33.500000E+00
500000 00

5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
55.000000E
000000E-01
01

0.000000E+00
0.000000E+00

①
②
③
④
⑤
⑥
⑦
⑧

5.000000E-01
5.000000E-01
55.000000E
000000E-01
01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
55.000000E-01
000000 01

PE#0
14

⑦

8

16

3

4

⑤

10

7

15 ⑧

11

12

1 ④

2

③

PE#1

基本的に初期全体メッシュファイルと同じ
局所要素番号による記述
コネクティビティ：「内点～内点」，「内点～外点」のみの情報

⑥

①

②
7

8

5

6
PE#2

T2K-FVM-03

45

デ タ入力部分（4/5）
データ入力部分（4/5）!C!C +------------+!C | COMM INPUT |!C ++------------++!C===
IUNIT= 12
open (IUNIT,file= hpcmw_eps_fvm_files(4), status='unknown')

分散通信データ#D COMM#D-COMM

read (IUNIT,(IUNIT '(a)')(a) ) LINE
read (IUNIT,'(6i12)') n_neighbor_pe
allocate (neighbor_pe(n_neighbor_pe))
allocate (import_index(0:n_neighbor_pe))
allocate (export(export_index(0:n_neighbor_pe))
index(0:n neighbor pe))import_index= 0
export_index= 0
read (IUNIT,(IUNIT '(a)')(a) ) LINE
read (IUNIT,'(6i12)') (neighbor_pe(k), k= 1, n_neighbor_pe)
read (IUNIT,'(a)') LINE
read (IUNIT,(IUNIT '(6i12)')(6i12) ) (import(import_index(k),index(k) k= 1,
1 n_neighbor_pe)
n neighbor pe)
nn= import_index(n_neighbor_pe)
allocate (import_item(nn))
read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') (import_item(k), k= 1, nn)

T2K-FVM-03

46

局所分散通信ファイル（comm.0）
隣接領域#NEIBPEtot
2#NEIBPE
1
2#IMPORT index
2
4#IMPORT items
7856#EXPORT index
24#EXPORT items
13
12#INTERNAL NODE
4#TOTAL NODE8#GLOBAL NODE ID
11 12 15 16
10 14

PE#0

7

14

15

16

8

3

4

10

11

12

7

1

2

7

8

5

6

8
PE#1

PE#2

T2K-FVM-03

47

局所分散通信ファイル（comm.0）
受信テーブル，外点情報#NEIBPEtot
2#NEIBPE
1
2#IMPORT index
2
4#IMPORT items
7856#EXPORT index
24#EXPORT items
13
12#INTERNAL NODE
4#TOTAL NODE8#GLOBAL NODE ID
11 12 15 16
10 14

PE#0

7

14

15

16

8

3

4

10

11

12

7

1

2

7

8

5

6

8
PE#1

PE#2

T2K-FVM-03

48

デ タ入力部分（5/5）
データ入力部分（5/5）
read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') (export_index(k), k= 1, n_neighbor_pe)
nn= export_index(n_neighbor_pe)
allocate (export(export_item(nn))
item(nn))read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') (export_item(k), k= 1, nn)
read (IUNIT,'(a)') LINE
read (IUNIT,(IUNIT '(6i12)')(6i12) ) intNODE
intNODE_tot
tot
read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') nn
allocate (NODE_GLOBAL(nn))(NODE GLOBAL(nn))
read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') (NODE_GLOBAL(k), k= 1, nn)close (IUNIT)!C===
end subroutine hpcmw_eps_fvm_input_gridT2K-FVM-03

49

局所分散通信ファイル（comm.0）
送信テーブル，境界点情報#NEIBPEtot
2#NEIBPE
1
2#IMPORT index
2
4#IMPORT items
7856#EXPORT index
24#EXPORT items
13
12#INTERNAL NODE
4#TOTAL NODE8#GLOBAL NODE ID
11 12 15 16
10 14

PE#0

7

14

15

16

8

3

4

10

11

12

7

1

2

7

8

5

6

8
PE#1

PE#2

T2K-FVM-03

50

デ タ入力部分（5/5）
データ入力部分（5/5）
read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') (export_index(k), k= 1, n_neighbor_pe)
nn= export_index(n_neighbor_pe)
allocate (export(export_item(nn))
item(nn))read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') (export_item(k), k= 1, nn)
read (IUNIT,'(a)') LINE
read (IUNIT,(IUNIT '(6i12)')(6i12) ) intNODE
intNODE_tot
tot

内点数

read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') nn

内点＋外点数（NODE_tot）

allocate (NODE_GLOBAL(nn))(NODE GLOBAL(nn))
read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') (NODE_GLOBAL(k), k= 1, nn)close (IUNIT)!C===
end subroutine hpcmw_eps_fvm_input_gridT2K-FVM-03

51

局所分散通信ファイル（comm.0）
内点数，総要素数（内点＋外点），全体要素番号#NEIBPEtot
2#NEIBPE
1
2#IMPORT index
2
4#IMPORT items
7856#EXPORT index
24#EXPORT items
13
12#INTERNAL NODE
4#TOTAL NODE8#GLOBAL NODE ID
11 12 15 16
10 14

PE#0

全体要素番号（局所番号順）
7

14

15

16

8

3

4

10

11

12

7

1

2

7

8

5

6

8
PE#1

PE#2

T2K-FVM-03

52

program eps_fvm
use hpcmw_eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call hpcmw_eps_fvm_init
call hpcmw_eps_fvm_input_grid____
call poi_gen
call hpcmw_eps_fvm_solver___
call output_ucd
call hpcmw_eps_fvm_finalize
end program eps_fvm

T2K-FVM-03

53

poi_gen（1/2）• 係数行列生成部は1PEの場合とほとんど同じ• 理由– 係数行列生成に必要な情報が全て局所分散メッシュファイルに含まれている。
含まれている局所分散メ シ
局所分散メッシュファイル
イル
8
1
2
3
4
5
6
7
8
8
5
6
7
1
1
2
8
3
2
2
4
0
1
1

1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00

1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1
2
1
2
3
4
3
4

2.50E+00
3.50E+00
2.50E+00
3.50E+00
2.50E+00
3.50E+00
1.50E+00
1.50E+00

1.000000E+00
1.000000E+00
11.000000E+00
000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00

1.000000E+00
1.000000E+00
1.000000E+00

2.50E+00
2.50E+00
3.50E+00
3.50E+00
1.50E+00
1.50E+00
2.50E+00
3.50E+00
5.000000E-01
5.000000E-01
55.000000E-01
000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01

5.000000E-01
5.000000E-01

5.00E-01
5.00E-01
5.00E-01
5.00E-01
5.00E-01
5.00E-01
5.00E-01
5.00E-01
5.000000E-01
5.000000E-01
55.000000E-01
000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01

0.000000E+00
0.000000E+00

T2K-FVM-03

54

poi_gen（2/2）• 例えば1番の要素（全体番号11番）における係数行列生成に必要な隣接要素の情報（2 3 5 7番（全体番号：
成に必要な隣接要素の情報（2,3,5,7番（全体番号：
12,15,7,10））は全て局所分散メッシュデータに全て含ま
れている。
れている
– うち，5番，7番は外点
局所分散メ シ
局所分散メッシュファイル
イル
8
1
2
3
4
5
6
7
8
8
5
6
7
1
1
2
8
3
2
2
4
0
1
1

1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00

1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1
2
1
2
3
4
3
4

2.50E+00
3.50E+00
2.50E+00
3.50E+00
2.50E+00
3.50E+00
1.50E+00
1.50E+00

1.000000E+00
1.000000E+00
11.000000E+00
000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00

1.000000E+00
1.000000E+00
1.000000E+00

2.50E+00
2.50E+00
3.50E+00
3.50E+00
1.50E+00
1.50E+00
2.50E+00
3.50E+00
5.000000E-01
5.000000E-01
55.000000E-01
000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01

5.000000E-01
5.000000E-01

5.00E-01
5.00E-01
5.00E-01
5.00E-01
5.00E-01
5.00E-01
5.00E-01
5.00E-01
5.000000E-01
5.000000E-01
55.000000E-01
000000E 01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01
5.000000E-01

0.000000E+00
0.000000E+00

T2K-FVM-03

55

program eps_fvm
use hpcmw_eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call hpcmw_eps_fvm_init
call hpcmw_eps_fvm_input_grid____
call poi_gen
call hpcmw_eps_fvm_solver___
call output_ucd
call hpcmw_eps_fvm_finalize
end program eps_fvm

T2K-FVM-03

56

前処理付き共役勾配法の並列化
Preconditioned Conjugate Gradient Method （CG）
0 = b-[A]x(0)0Compute r(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)(i 1) z(i-1)(i 1)
ρi-1= r(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

並列計算，領域間通信が必要な部分• 行列ベクトル積• 内積

T2K-FVM-03

57

hpcmw eps fvm comm.＊
hpcmw_eps_fvm_comm.
• MPIに関連した共通サブルーチン群– 通信• MPIサブルーチンの代替
MPIサブル チンの代替• 「eps_fvm」における1対1通信– 一般化された通信テーブル• 利点– MPIサブルーチンを直接呼ぶよりも引数を減らすことができる。
MPIサブル チンを直接呼ぶよりも引数を減らす とができる。
– 1対1通信においては「送信バッファへの代入，送信，受信，受信バッファからの読み出し」という操作をそのたびに記述する必要が無くなる。

T2K-FVM-03

58

hpcmw_eps_fvm_commの内容h
http://nkl.cc.u-tokyo.ac.jp/tutorial/parallel_lib_tutorial/
// kl
k
j /
i l/
ll l lib
i l/
MPI_BARRIER
subroutine hpcmw
hpcmw_eps_fvm_barrier
eps fvm barrier
MPI_ALLREDUCE(スカラー）
subroutine hpcmw_eps_fvm_allreduce_R ( VAL, ntag)
subroutine hpcmw
hpcmw_eps_fvm_allreduce_I
eps fvm allreduce I ( VAL
VAL, ntag)
MPI_BCAST（スカラー）
subroutine hpcmw_eps_fvm_bcast_R ( VAL, nbase)subroutinebtih
hpcmw_eps_fvm_bcast_Ifbt I ( VAL
VAL, nbase)b)subroutine hpcmw_eps_fvm_bcast_C ( VAL, n, nbase)
MPI_ALLREDUCE(ベクトル）
subroutine hpcmw_eps_fvm_allreduce_RV ( VAL, n, ntag)
subroutine hpcmw_eps_fvm_allreduce_IV ( VAL, n, ntag)
MPI_BCAST（ベクトル）
subroutine hpcmw_eps_fvm_bcast_RV ( VAL, n, nbase)subroutine hpcmw_eps_fvm_bcast_IV ( VAL, n, nbase)subroutine hpcmw
hpcmw_eps_fvm_bcast_CV
eps fvm bcast CV ( VAL, n, nn, nbase)
一対一通信用（ベクトル）
subroutine hpcmw_eps_fvm_update_1_R (X, n)T2K-FVM-03

59

hpcmw
p
_eps
p _fvm_allreduce_R!C!C***!C*** hpcmw
p_eps
p _fvm_allREDUCE_R!C***!C
subroutine hpcmw_eps_fvm_allreduce_R ( VAL, ntag)
use hpcmw_eps_fvm_util
implicit REAL*8 (A-H,O-Z)
i tinteger:: ntag,tiierr
real(kind=kreal) :: VAL, VALM
if (ntag .eq. hpcmw_sum) thencall MPI_allREDUCE
&
(VAL,
(
, VALM,, 1,, MPI_DOUBLE_PRECISION,, MPI_SUM,,&
MPI_COMM_WORLD, ierr)
endif
if (ntag .eq. hpcmw_max) thencall MPI_allREDUCE&(VAL(VAL, VALM
VALM, 1
1, MPIMPI_DOUBLE_PRECISION,DOUBLE PRECISION MPIMPI_MAX,MAX&MPI_COMM_WORLD, ierr)
endif
if (ntag .eq. hpcmw_min) thencall MPI_allREDUCE&(VAL, VALM, 1, MPI_DOUBLE_PRECISION, MPI_MIN,&
MPI_COMM_WORLD, ierr)
endif
VAL= VALM
end subroutine hpcmw_eps_fvm_allreduce_R

&
&

&
&

&
&

T2K-FVM-03

60

hpcmw
p
_eps
p _fvm_allreduce_RV!C!C***!C*** hpcmw_eps_fvm_allREDUCE_RV
!C
!C***
!C
subroutine hpcmw_eps_fvm_allreduce_RV ( VAL, n, ntag)
use hpcmw_eps_fvm_util
implicit REAL*8 (A-H,O-Z)
integer :: n, ntag, ierrreal(kind=kreal), dimension(n) :: VAL
real(kind=kreal), dimension(:), allocatable :: VALM
allocate (VALM(n))
if (ntag .eq. hpcmw_sum) thencall MPI_allREDUCE
MPI allREDUCE&(VAL, VALM, n, MPI_DOUBLE_PRECISION, MPI_SUM,&
MPI_COMM_WORLD, ierr)
endif
if (ntag .eq. hpcmw_max) thencall MPI_allREDUCE&(VAL, VALM, n, MPI_DOUBLE_PRECISION, MPI_MAX,&MPI_COMM_WORLD, ierr)
endif
if (ntag .eq.
eq hpcmw_min)
hpcmw min) thencall MPI_allREDUCE&(VAL, VALM, n, MPI_DOUBLE_PRECISION, MPI_MIN,&
MPI_COMM_WORLD, ierr)
endif
VAL= VALM
deallocate (VALM)
end subroutine hpcmw_eps_fvm_allreduce_RV

&
&

&
&

&
&

T2K-FVM-03

61

hpcmw
pc
_eps_
eps fvm_update_
update 1_R（1/2）（ / ）subroutine hpcmw_eps_fvm_update_1_R (X, n)
use hpcmw_eps_fvm_util
implicit REALREAL*88 (A-H(A H,O
O-Z)
Z)
integer :: n, nn, ierrreal(kind=kreal), dimension(n) :: Xreal(kind=kreal), dimension(:), allocatable :: WS, WR
integer(kind=kint ), dimension(:,:), allocatable :: sta1
integer(kind=kint ), dimension(:,:), allocatable :: sta2
integer(kind=kint ), dimension(: ), allocatable :: req1
integer(kind=kint ), dimension(: ), allocatable :: req2
nn= max (n, import_index(n_neighbor_pe),&
exportexport_index(n_neighbor_pe))
index(n neighbor pe))allocate (WS(nn), WR(nn))

送信，受信バッファの定義!C!C-- INIT.
allocate (sta1(MPI_STATUS_SIZE,n_neighbor_pe))
allocate (sta2(MPI_STATUS_SIZE,n_neighbor_pe))
allocate (req1(n_neighbor_pe))
allocate (req2(n_neighbor_pe))&

T2K-FVM-03

62

hpcmw
pc
_eps_
eps fvm_update_
update 1_R（2/2）（ / ）!C!C-- SEND
do neib= 1, n_neighbor_pe
istart= export_index(neib-1)
iinum= export_index(neibt i d ( ib ) - istart
i t t

送信バッファへの代入

do k= istart+1, istart+inum
WS(k)= X(export_item(k))
enddocall MPI_ISENDMPI ISEND (WS(istart+1), inum, MPI
MPI_DOUBLE_PRECISION,DOUBLE PRECISION,&
neighbor_pe(neib), 0, MPI_COMM_WORLD,&
req1(neib), ierr)
enddo

送信
!C
!C-- RECEIVE
d neib=do
ib 1
1, n_neighbor_pe
i hb
istart= import_index(neib-1)
inum = import_index(neib ) - istartcall MPI_IRECV (WR(istart+1), inum, MPI_DOUBLE_PRECISION,&
neighbor_pe(neib), 0, MPI_COMM_WORLD,&
req2(neib),eq ( e b), ierr)
e )
enddo
call MPI_WAITALL (n_neighbor_pe, req2, sta2, ierr)do neib= 1, n_neighbor_pe
i t t i
istart=import_index(neib-1)
t i d ( ib 1)
inum = import_index(neib ) - istart
do k= istart+1, istart+inum
X(import_item(k))= WR(k)
enddo
enddo
call MPI_WAITALL (n_neighbor_pe, req1, sta1, ierr)deallocate (sta1, sta2, req1, req2, WS, WR)]end subroutine hpcmw_eps_fvm_update_1_R

&
&

&
&

T2K-FVM-03

63

送信（MPI_Isend/Irecv/Waitall）_
SENDbuf

neib#1neib#2neib#3neib#4BUFlength_e

BUFlength_e

BUFlength_e

BUFlength_e

export_index(0)+1

export_index(1)+1 export_index(2)+1

export_index(3)+1

export_index(4)do neib= 1, NEIBPETOT
do k= export
export_index(neib-1)+1,
index(neib 1)+1 export
export_index(neib)
index(neib)
kk= export_item(k)SENDbuf(k)= VAL(kk)送信バッファへの代入
enddo
enddo
温度などの変数を直接送信，受信に使
do neib= 1, NEIBPETOTiS_e= export_index(neib-1) + 1
iE_e= export_index(neib )BUFlengthg _e= iE_e + 1 - iS_e&&

うのではなく，このようなバッファへ一回代入して計算することを勧める。

call MPI_Isend&(SENDbuf(iS_e), BUFlength_e, MPI_INTEGER, NEIBPE(neib), 0,&
MPI_COMM_WORLD, request_send(neib), ierr)enddo
call MPI_Waitall (NEIBPETOT, request_send, stat_recv, ierr)

T2K-FVM-03

64

hpcmw
pc
_eps_
eps fvm_update_
update 1_R（2/2）（ / ）!C!C-- SEND
do neib= 1, n_neighbor_pe
istart= export_index(neib-1)
iinum= export_index(neibt i d ( ib ) - istarti t t
do k= istart+1, istart+inum
WS(k)= X(export_item(k))
enddocall MPI_ISENDMPI ISEND (WS(istart(WS(istart+1),1), inum, MPI_DOUBLE_PRECISION,MPI DOUBLE PRECISION,&
neighbor_pe(neib), 0, MPI_COMM_WORLD,&
req1(neib), ierr)enddo!C!C-- RECEIVE
d neib=do
ib 1
1, n_neighbor_pe
i hb
受信istart= import_index(neib-1)
inum = import_index(neib ) - istartcall MPI_IRECV (WR(istart+1), inum, MPI_DOUBLE_PRECISION,&
neighbor_pe(neib), 0, MPI_COMM_WORLD,&
req2(neib),eq ( e b), ierr)
e )
enddo
call MPI_WAITALL (n_neighbor_pe, req2, sta2, ierr)do neib= 1, n_neighbor_pe
i t t i
istart=import_index(neib-1)
t i d ( ib 1)
inum = import_index(neib ) - istart
do k= istart+1, istart+inum
X(import_item(k))= WR(k)
enddo
enddo
call MPI_WAITALL (n_neighbor_pe, req1, sta1, ierr)deallocate (sta1, sta2, req1, req2, WS, WR)]end subroutine hpcmw_eps_fvm_update_1_R

&
&

&
&

T2K-FVM-03

65

hpcmw
pc
_eps_
eps fvm_update_
update 1_R（2/2）（ / ）!C!C-- SEND
do neib= 1, n_neighbor_pe
istart= export_index(neib-1)
iinum= export_index(neibt i d ( ib ) - istarti t t
do k= istart+1, istart+inum
WS(k)= X(export_item(k))
enddocall MPI_ISENDMPI ISEND (WS(istart(WS(istart+1),1), inum, MPI_DOUBLE_PRECISION,MPI DOUBLE PRECISION,&
neighbor_pe(neib), 0, MPI_COMM_WORLD,&
req1(neib), ierr)enddo!C!C-- RECEIVE
d neib=do
ib 1
1, n_neighbor_pe
i hb
istart= import_index(neib-1)
inum = import_index(neib ) - istartcall MPI_IRECV (WR(istart+1), inum, MPI_DOUBLE_PRECISION,&
neighbor_pe(neib), 0, MPI_COMM_WORLD,&
req2(neib),eq ( e b), ierr)
e )
enddo
call MPI_WAITALL (n_neighbor_pe, req2, sta2, ierr)do neib= 1, n_neighbor_pe
i t t i
istart=import_index(neib-1)
t i d ( ib 1)
inum = import_index(neib ) - istart
do k= istart+1, istart+inum
X(import_item(k))= WR(k)
enddo
enddo

受信バッファからの受信バ
フ からのデータ取り出しcall MPI_WAITALL (n_neighbor_pe, req1, sta1, ierr)deallocate (sta1, sta2, req1, req2, WS, WR)]end subroutine hpcmw_eps_fvm_update_1_R

&
&

&
&

T2K-FVM-03

66

受信（MPI_Isend/Irecv/Waitall）_
do neib= 1, NEIBPETOTiS_i= import_index(neib-1) + 1
iE_i= import_index(neib )BUFlength i= iE
BUFlength_i=iE_ii + 1 - iS_i
iS icall MPI_Irecv&&(RECVbuf(iS_i), BUFlength_i, MPI_INTEGER, NEIBPE(neib), 0,&&
MPI_COMM_WORLD, request_recv(neib), ierr)enddo
ddcall MPI_Waitall (NEIBPETOT, request_recv, stat_recv, ierr)do neib= 1,, NEIBPETOT
do k= import_index(neib-1)+1, import_index(neib)
kk= import_item(k)VAL(kk)= RECVbuf(k)
enddo
enddo

受信バッファから代入

RECVbuf

neib#1neib#2neib#3neib#4BUFlength_i

BUFlength_i

BUFlength_i

BUFlength_i

import_index(0)+1

import_index(1)+1 import_index(2)+1

import_index(3)+1

import_index(4)T2K-FVM-03

67

hpcmw solver.f
hpcmw_solver.f
subroutine hpcmw
p_eps
p _fvm_solver
use hpcmw_eps_fvm_all
implicit REAL*8 (A-H,O-Z)
EPS = 1.d-8
ITR = NODE_tot
call hpcmw_eps_fvm_solver_CG&( i
intNODE_tot,
tNODE t t NODENODE_tot,t t NPLU
NPLU, DD, BFORCEBFORCE, DELPHI
DELPHI, EPS
EPS,&
ITR, IER, index, item, AMAT, COMMtime)
open (11, file='fvmmg.ctrl', status='unknown')
read
d (11
(11,*)*) NX
NX, NY
NY, NZ
close (11)iS= NX*NY*NZ/2 + NX*NY/2
do i= iS+1,
iS+1 iS+NX
write (*,'(i8,3(1pe16.6))') i, DELPHI(i)enddo
end subroutine hpcmw
hpcmw_eps_fvm_solver
eps fvm solver&&

T2K-FVM-03

68

共役勾配法
!C
!C***
!C*** CG!C***!C
subroutine hpcmw_eps_fvm_solver_CG&( N, NP, NPLU, D, B, X, EPS, ITR, IER,&
index
index, item,
item COEF,
COEF Tcomm)use hpcmw_eps_fvm_util
implicit REAL*8 (A-H,O-Z)real(kind=kreal),l(ki d kl) dimension(Ndii (N ) :: Dreal(kind=kreal), dimension(NP) :: Breal(kind=kreal), dimension(NP) :: X
integer, dimension(0:N) :: index
integer, dimension(NPLU):: item
real(kind=kreal), dimension(NPLU):: COEF
real(kind=kreal) :: EPS, Tcomm
integer :: ITR, IER
integer :: P, Q, R, Z, DD
real(kind=kreal), dimension(:,:), allocatable, save :: W

N :
NP:

intNODE_tot
NODE_tot

COEF: AMAT
B
: RHS

&
&

T2K-FVM-0369

共役勾配法の 並列化」（1/4）
共役勾配法の「並列化」（1/4）• 行列ベクトル積!C!C +-------------+!C | {q}= [A]{p} |!C +-------------+!C===

exchange W(i,P)

do i= 1, N
W(i,Q) = D(i) * W(i,P)
do j= index(i-1)+1, index(i)W(i,Q) = W(i,Q) + COEF(j) * W(item(j),P)
enddo
enddo
!C===

T2K-FVM-0370

共役勾配法の 並列化」（１/4）
共役勾配法の「並列化」（１/4）• 行列ベクトル積!C!C ++-------------++!C | {q}= [A]{p} |!C +-------------+!C===
call hpcmw_eps_fvm_update_1_R (W(1,P), NP)
do i= 1, N
W(i,Q) = D(i) * W(i,P)
do j= index(i-1)+1, index(i)W(i,Q) = W(i,Q) + COEF(j) * W(item(j),P)
enddo
enddo
!C===

T2K-FVM-0371

共役勾配法の「並列化」（2/4）• 内積：MPI_ALLREDUCE!C!C +-------------+!C | RHO= {r}{z} |!C +-------------+
+
+
!C===
RHO= 0.d0

!C===

!C
!C +-------------+!C | RHO= {r}{z} |!C +-------------+
+
+
!C===
RHO= 0.d0

do i= 1
1, N
RHO= RHO + W(i,R)*W(i,Z)
enddo

do i= 1
1, NP
RHO= RHO + W(i,R)*W(i,Z)
enddo

allreduce RHO

allreduce RHO!C===
こうしてはいけない・・・何故か?
うな
何故か（N：内点数，NP：内点＋外点）
N : i
intNODE_tot
tNODE t t
NP: NODE_tot

T2K-FVM-0372

共役勾配法の 並列化」（2/4）
共役勾配法の「並列化」（2/4）• 内積：MPI_ALLREDUCE
内積!C!
!C +-------------+!C | RHO= {r}{z} |!C +-------------+!C===
RHO= 0.d0
do i= 1
1, N
RHO= RHO + W(i,R)*W(i,Z)
enddo
call
ll hpcmw_eps_fvm_allreduce_R
h
f
ll d
(
(RHO,
O h
hpcmw_sum))!C===T2K-FVM-0373

共役勾配法の 並列化」（3/4）
共役勾配法の「並列化」（3/4）• N（intNODE_tot）とNP（NODE_tot）の違いに注意。
N（）とNP（）の違いに注意– 基本的に計算はN個分やればよい（ループはdo i=1,N）。
– 外点の値を変えるような（左辺にくる）計算はしない。
• 外点の値が必要なときは行列ベクトル積のみ。
• そのときに通信して外点の値をもらってくれば良い。
そ とき 通信外点 値をもらくれば良T2K-FVM-0374

共役勾配法の「並列化」（4/4）• 最後に，従属変数（X）の外点における最新値をもらっておくことを忘れないように。
を忘れなう 。
DNRM20= 0.d0
do i= 1, N
X(i) = X(i)+ ALPHA * W(i,P)
W(i,R)=(i ) W(i,R)(i ) - ALPHA * W(i,Q)
(i )
enddo
DNRM2 = 0.0
do i= 1, N
DNRM2= DNRM2 + W(i,R)**2
enddo
call hpcmw
p_eps
p _fvm_allreduce_R (DNRM2, hpcmw
p_sum)
RESID= dsqrt(DNRM2/BNRM2)
if ( RESID.le.EPS) goto 900
RHO1 = RHO
enddo
900 continue
call hpcmw
hpcmw_eps_fvm_update_1_R
eps fvm update 1 R (X, NP)
returnend subroutine hpcmw_eps_fvm_solver_CG

T2K-FVM-03

75

program eps_fvm
use hpcmw_eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call hpcmw_eps_fvm_init
call hpcmw_eps_fvm_input_grid____
call poi_gen
call hpcmw_eps_fvm_solver___
call output_ucd
call hpcmw_eps_fvm_finalize
end program eps_fvm

T2K-FVM-03

76

可視化• いわゆる，「並列可視化」をやってもらうつもりだったのであいわゆる 「並列可視化」をやってもらうつもりだったのであるが，プログラムが完全構造格子用に作られたものであるた
め 「eps fvm」のような非構造格子的なものでは対応が難
め，「eps_fvm」のような非構造格子的なものでは対応が難しい。
• 幸い，MicroAVSが各自で使えるので，少し違ったアプロー
チを試みることにした。
チを試みることにしたT2K-FVM-03

77

やる きことやるべきこと初期全体メッシュ

領域分割

結果をもう一回
一つのファイルにする

分散処理

T2K-FVM-03

78

可視化• いわゆる，「並列可視化」をやってもらうつもりだったのであいわゆる 「並列可視化」をやってもらうつもりだったのであるが，プログラムが完全構造格子用に作られたものであるた
め 「eps fvm」のような非構造格子的なものでは対応が難
め，「eps_fvm」のような非構造格子的なものでは対応が難しい。
• 幸い，MicroAVSが各自で使えるので，少し違ったアプロー
チを試みることにした。
チを試みることにした• 本当の「並列可視化」とは言えないが，「 MPI_Gatherv」等
関数
す 機
，
う
の関数を利用する機会にもなるので，このようなやりかたを試みる。

T2K-FVM-03

79

AVS UCDファイル• 最終出力はUCD（Unstructured Cell Data）フォーマットとしData）フォ マットとして，これをMicroAVSで読み込む。
ファイル名は「* inp」でなければならない• ファイル名は「*.inp」でなければならない。
• UCDファイルは以下の2部分から構成されているのであるが，今回は「形状」については 初期全体メ シ ジ ネレ タで
今回は「形状」については，初期全体メッシュジェネレータで自動的に生成されている（fvm_entire_mesh.inpgeo）。
）– 形状（節点，要素）– 結果• 結果のみについてアペンドすればよい。

T2K-FVM-03

80#MGCTRL<$FVM>/ex/
fvmmg.ctrl

メッシュジェネレータ
eps fvm mg
eps_fvm_mg

各座標軸方向の要素数（名称固定）

初期全体
メッシュ生成生成#S-GRID<$FVM>/ex/
fvm_entire_mesh.dat#S-GRID-UCD<$FVM>/ex/
fvm_entire_mesh.inp#S-GRID-GEO<$FVM>/ex/
fvm_entire_mesh.inp-geo

初期全体メッシュデータ（名称固定）

初期全体メッシュデータ
UCDファイル（AVS用）（名称固定）

初期全体メッシュデータ
UCDファイル形状部（AVS用）（名称固定）

節点数，要素数等

節点数，要素数等

節点座標

節点座標

要素コネクティビティ

要素コネクティビティ

要素の値，計算結果

T2K-FVM-03

81

AVS UCDファイルの例（1/2）
45節点 16要素の六面体
45節点，16要素の六面体
45
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

16
0.000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
0 000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
1.000000E
1 000000E+00
00
2.000000E+00
0.000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
1.000000E+00
1 000000E+00
2.000000E+00
0.000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
1.000000E+00
2.000000E
2.000000E+00
00
0.000000E+00
1.000000E+00
2.000000E+00
0.000000E+00
1.000000E+00
2.000000E+00

00.000000E+00
0.000000E+00
0.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00
0.000000E+00
0 000000E+00
0.000000E+00
0.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
2.000000E+00
2 000000E+00
2.000000E
00
2.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.000000E+00
1.000000E+00
1 000000E+00
1.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E
00
1.000000E+00
1.000000E+00
1.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00

1

00.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1 000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1.000000E+00
1 000000E+00
1.000000E
00
1.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00
2.000000E+00
3.000000E+00
3.000000E+00
3.000000E+00
3.000000E+00
3 000000E+00
3.000000E+00
3.000000E+00
3.000000E+00
3.000000E+00
3.000000E+00
4.000000E+00
4.000000E+00
4.000000E+00
4.000000E
00
4.000000E+00
4.000000E+00
4.000000E+00
4.000000E+00
4.000000E+00
4.000000E+00

T2K-FVM-03

82

AVS UCDファイルの例（2/2）
45節点 16要素の六面体
45節点，16要素の六面体
1 1 hex
2 1 hex
3 1 hexh
4 1 hex
5 1 hex
6 1 hex
7 1 hex
8 1 hex
9 1 hex
10 1 hex
11 1 hex
12 1 hex
13 1 hex
14 1 hex
15 1 hex
16 1 hex
1 1
COLOR, color
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

1
2
4
5
10
11
13
14
19
20
22
23
28
29
31
32
2.00e+00
2.00e+00
2.00e+00
2.00e+00
2.00e+00
2.00e+00
2.00e+00
2 00 00
2.00e+00
1.00e+00
1.00e+00
1.00e+00
1.00e+00
1.00e+00
1.00e+00
1 00e+00
1.00e+00
1.00e+00

2
3
5
6
11
12
14
15
20
21
23
24
29
30
32
33

5
6
8
9
14
15
17
18
23
24
26
27
32
33
35
36

4
5
7
8
13
14
16
17
22
23
25
26
31
32
34
35

10
11
13
14
19
20
22
23
28
29
31
32
37
38
40
41

この部分のみ自分で作成す
れば良い

11
12
14
15
20
21
23
24
29
30
32
33
38
39
41
42

14
15
17
18
23
24
26
27
32
33
35
36
41
42
44
45

13
14
16
17
22
23
25
26
31
32
34
35
40
41
43
44

T2K-FVM-03

83

手順• 各プロセッサにおける計算結果（PHI(:)に格納される），グ
ロ バル要素番号を あるプロセッサに集める
ローバル要素番号を，あるプロセッサに集める：
MPI_Gatherv使用– 「内点」の結果のみで良い。
「内点」の結果のみで良い• そのプロセッサから，結果を以下の書式で，グローバル番号順に書き出す（「 は入れない）
順に書き出す（「,」は入れない）：<グローバル要素番号><計算結果>• ファイルを「fvm_entire_mesh.inp-geo」にアペンドしてUCDファイルを生成する。
UCDファイルを生成する• ヒント– MPI_Gathervを使う– グローバルIDは局所分散通信ファイルより得る

T2K-FVM-03

MPI_Gatherv

MPI_Gatherv：
局所ベクト から全体ベクト 生成
局所ベクトルから全体ベクトル生成• cd <$FVM>/gatherv• 「a2.0~a2.3」から局所ベクトル情報を読み込み，「全体ベクトクトル」情報を各プロセッサに生成するプログラムを作成
情報を各プ セ サ 生成するプ グ ムを作成する。MPI_Gathervを使用する。
PE#0

PE#1

PE#2

PE#3

8
101.0
103.0
105 0
105.0
106.0
109.0
111.0
121.0
151.0

5
201.0
203.0
205 0
205.0
206.0
209.0

7
301.0
303.0
305 0
305.0
306.0
311.0
321.0
351.0

3
401.0
403.0
405 0
405.0

84

T2K-FVM-03

MPI_Gatherv

局所ベクトルから全体ベクトル生成
MPI_Gathervを使う場合（1/5）

PE#0

PE#0

PE#1

PE#1

PE#2

PE#2
MPI Gatherv
MPI_Gatherv

PE#3

PE#3

85

T2K-FVM-03

MPI_Gatherv

MPI Gatherv
MPI_Gatherv•

MPI_Gather の可変長さベクトル版– 「局所データ」から「全体データ」を生成する• call MPI_Gatherv (sendbuf, scount, sendtype, recvbuf, rcounts,
displs, recvtype, root, comm, ierr)
–
–
–
–
–
–
–
–
–
–

sendbufscountsendtype
recvbufrcounts
displsrecvtype
root
comm
ierr

任意
整数
整数
任意
整数
整数
整数
整数
整数
整数

I
I
I
O
I
I
I
I
I
O

送信バッファの先頭アドレス，
送信バッファの先頭アドレス
送信メッセージのサイズ
送信メッセージのデータタイプ
受信バッファの先頭アドレス，
受信メ セ ジのサイズ（配列 サイズ
受信メッセージのサイズ（配列：サイズ＝PETOT））
受信メッセージのインデックス（配列：サイズ＝PETOT+1）
受信メッセージのデータタイプ
受信メッセージ受信元（ランク）
コミュニケータを指定する
完了コード

86

T2K-FVM-03

MPI_Gatherv

MPI Gatherv（続き）
MPI_Gatherv（続き）• call MPI_Gatherv (sendbuf, scount, sendtype, recvbuf, rcounts,
di ldispls,recvtype,troot,
t comm, iierr))
整数 I 受信メッセージのサイズ（配列：サイズ＝PETOT）
整数 I 受信メッセージのインデックス（配列：サイズ＝PETOT+1）––

rcounts
displs–

この2つの配列は，最終的に生成される「全体データ」のサイズに関する配列であるため，各プロセスで配列の全ての値が必要になる：• もちろん各プロセスで共通の値を持つ必要がある。

–

通常はstride(i)=rcounts(i)

PE#0

PE#1

stride(1)rcounts(1)displs(1)=0
displs(1)0stride(2)rcounts(2)PE#2

PE#(m-2) PE#(m-1)

stride(3)stride(m-1)

stride(m)rcounts(3)rcounts(m-1)

rcounts(m)displs(2)displs(2)=displs(1) + stride(1)size(recvbuf)= displs(PETOT+1)= sum(stride)

displs(m 1)displs(m+1)=
displs(m) + stride(m)87

T2K-FVM-03

PE#0

PE#2

N

N

rcounts(3
3)

PE#1

N

stride(1)displs(2)p ( )rcounts((2)

局所データから全体データを局所デタから全体デ タを生成する

rcounts(1)MPI_Gathervでやっていること

displs(1) MPI_
Gatherv

stride(2)displs(3)p ( )stride(3)displs(4)Nrccounts(4)

PE#3

stride(4)displs(5)局所データ：sendbuf

全体データ：recvbuf

88

T2K-FVM-03

MPI_Gatherv

MPI_Gathervでやって
いること
PE#0

PE#1

N

N

displs(2)stride(2)= rcounts(2)displs(3)stride(3)= rcounts(3)rcountss
(4)

PE#3

N

stride(1)= rcounts(1)rccounts(3))

PE#2

N

rccounts(1)) rcountts(2)局所デ タから全体デ タを生成する
局所データから全体データを生成する

displs(1)p ( )stride(4)= rcounts(4)di l (4)displs(4)displs(5)局所データ：sendbuf

全体データ：recvbuf

89

T2K-FVM-03

MPI_Gatherv

MPI_Gatherv詳細（1/2）_• call MPI_ALLGATHERV (sendbuf, scount, sendtype, recvbuf, rcounts,
displs, recvtype, root, comm, ierr)––rcounts
displs

整数
整数

I
I

受信メッセージのサイズ（配列：サイズ＝PETOT）
受信メッセージのインデックス（配列：サイズ＝PETOT+1）• rcounts– 各PEにおけるメッセージサイズ：局所データのサイズ• displs– 各局所データの全体データにおけるインデックス– displs(PETOT+1)が全体データのサイズ

PE#0

PE#1

stride(1)rcounts(1)displs(1)=0( )stride(2)rcounts(2)PE#2

PE#(m-2) PE#(m-1)

stride(3)stride(m-1)

stride(m)rcounts(3)rcounts(m-1)

rcounts(m)displs(2)=p ( )displs(1) + stride(1)size(recvbuf)= displs(PETOT+1)= sum(stride)

displs(m+1)=
p ()
displs(m) + stride(m)90

T2K-FVM-03

MPI_Gatherv

MPI_Gatherv詳細（2/2）_• rcountsとdisplsは各プロセスで共通の値が必要– 各プロセスのベクトルの大きさ N をallgatherして，rcounts
をしてに相当するベクトルを作る。
– rcountsから各プロセスにおいてdisplsを作る（同じものができる）。
• stride(i)=( ) rcounts(i)( ) とする– rcountsの和にしたがってrecvbufの記憶領域を確保する。
PE#0

PE#1

stride(1)rcounts(1)displs(1)=0( )stride(2)rcounts(2)PE#2

PE#(m-2) PE#(m-1)

stride(3)stride(m-1)

stride(m)rcounts(3)rcounts(m-1)

rcounts(m)displs(2)=p ( )displs(1) + stride(1)size(recvbuf)= displs(PETOT+1)= sum(stride)

displs(m+1)=
p ()
displs(m) + stride(m)91

T2K-FVM-03

MPI_Gatherv

MPI G thMPI_Gatherv使用準備
使用準備• “a2.0”~”a2.3”から，全体ベクトルを生成する。
“ 2 0” ” 2 3”から 全体ベクトルを生成する• 各ファイルのベクトルのサイズが，8,5,7,3であるから，長さ23（=8+5+7+3）のベクトルができることになる。

S1-2

PE#0

PE#1

PE#2

PE#3

8
101.0
103.0
105 0
105.0
106.0
109.0
111.0
121.0
151.0

5
201.0
203.0
205 0
205.0
206.0
209.0

7
301.0
303.0
305 0
305.0
306.0
311.0
321.0
351.0

3
401.0
403.0
405 0
405.0

92

T2K-FVM-03

MPI_Gatherv

局所 全体ベクトル生成 手順
局所⇒全体ベクトル生成：手順
PE#0

PE#1

stride(1)rcounts(1)displs(1)=0

stride(2)rcounts(2)PE#2

PE#(m-2) PE#(m-1)

stride(3)stride(m-1)
stride(m1)

stride(m)rcounts(3)rcounts(m-1)

rcounts(m)displs(2)=displs(1) + stride(1)displs(m+1)=
displs(m) + stride(m)size(recvbuf)= displs(PETOT+1)= sum(stride)

•
•
•
•
S1-2

局所ベクトル情報を読み込む「rcounts」， 「displs」を作成する「recvbuf」を準備する
Gatherv
93

T2K-FVM-03

MPI_Gatherv

局所⇒全体ベクトル生成（1/2）<$FVM>/gatherv/mpigathertest.*
implicit REAL*8 (A-H,O-Z)
include 'mpif.h'
integer :: PETOT, my_rank, SOLVER_COMM, ierrreal(kind=8), dimension(:), allocatable :: VEC, VEC2, VECg
integer (kind=4), dimension(:), allocatable :: COUNT, COUNTindex
character(len=80)ca acte ( e 80):: filenamee a ecall MPI_INIT(ierr)
call MPI_COMM_SIZE (MPI_COMM_WORLD, PETOT, ierr )call MPI_COMM_RANK (MPI_COMM_WORLD, my_rank, ierr )call MPI_COMM_DUP (MPI_COMM_WORLD, SOLVER_COMM, ierr)
if (my_rank.eq.0) filename= 'a2.0'
if (my_rank.eq.1) filename= 'a2.1'
if (my_rank.eq.2) filename= 'a2.2'
if (my_rank.eq.3) filename= 'a2.3'
open (21,(21 file= filename,filename status= 'unknown')
unknown )read (21,*) N
allocate (VEC(N))
do i= 1, N
read (21,*) VEC(i)enddo

中身を書き出して見よう

allocate (COUNT(PETOT), COUNTindex(PETOT+1))
call MPI_allGATHER ( N, 1, MPI_INTEGER,&
COUNT, 1, MPI_INTEGER,&MPI_COMM_WORLD, ierr)
COUNTindex(1)= 0&
各PEにおけるベクトル長さの情報が&「COUNT」に入る（「rcounts」）

do ip= 1, PETOTCOUNTindex(ip+1)= COUNTindex(ip) + COUNT(ip)enddo

中身を書き出して見よう
94

T2K-FVM-03

MPI_Gatherv

MPI Allgather
MPI_Allgather•

P#0 A0P#1 B0P#0 A0 B0 C0 D0
All g
gatherP#1 A0 B0 C0 D0P#2 C0P#2 A0 B0 C0 D0P#3 D0P#3 A0 B0 C0 D0

MPI_Gather＋MPI_Bcast• call MPIMPI_Allgather
Allgather (sendbuf(sendbuf, scountscount, sendtype
sendtype, recvbuf
recvbuf, rcountrcount,recvtype, comm, ierr)
–
–
–
–
–
–
–
–

sendbufscountsendtype
recvbufrcountrecvtype
t
comm
ierr

任意
整数
整数
任意
整数
整数
整数
整数

I
I
I
O
I
I
I
O

送信バッファの先頭アドレス，
送信メッセージのサイズ
送信メッセ
ジのサイズ
送信メッセージのデータタイプ
受信バッファの先頭アドレス，
受信メッセージのサイズ
受信メッセ ジのデ タタイプ
受信メッセージのデータタイプ
コミュニケータを指定する
完了コード

95

T2K-FVM-03

MPI_Gatherv

局所⇒全体ベクトル生成（2/2）<$FVM>/gatherv/mpigathertest.*
do ip= 1, PETOTCOUNTindex(ip+1)= COUNTindex(ip) + COUNT(ip)
enddo「displs」に相当するものを生成。

allocate (VECg(COUNTindex(PETOT+1)))
VECg= 0.d0&&&

call MPI_Gatherv
MPI Gatherv( VEC , N, MPI_DOUBLE_PRECISION,
VECg, COUNT, COUNTindex, MPI_DOUBLE_PRECISION,
0, MPI_COMM_WORLD, ierr)&&&

if (my_rank.eq.0) then
d i
do
i= 1
1, COUNTi
COUNTindex(PETOT+1)d (PETOT+1)
write (*,'(2i8,f10.0)') my_rank, i, VECg(i)
enddo
endif
call MPI_FINALIZE (ierr)
stop
end

PE#0

PE#1

stride(1)t id (1)rcounts(1)displs(1)=0

stride(2)t id (2)rcounts(2)PE#2

PE#(m-2) PE#(m-1)

stride(3)t id (3)stride(m-1)
t id ( 1)stride(m)t id ( )rcounts(3)rcounts(m-1)

rcounts(m)displs(2)=displs(1) + stride(1)size(recvbuf)= displs(PETOT+1)= sum(stride)

displs(m+1)=
displs(m) + stride(m)96

T2K-FVM-03

MPI_Gatherv

局所⇒全体ベクトル生成（2/2）<$FVM>/gatherv/mpigathertest.*
do ip= 1, PETOTCOUNTindex(ip+1)= COUNTindex(ip) + COUNT(ip)enddoallocate (VECg(COUNTindex(PETOT+1)))
VECg= 0.d0

「recbuf」のサイズ

call MPI_Gatherv
MPI Gatherv&( VEC , N, MPI_DOUBLE_PRECISION,&
VECg, COUNT, COUNTindex, MPI_DOUBLE_PRECISION,&
0, MPI_COMM_WORLD, ierr)&&&

if (my_rank.eq.0) then
d i
do
i= 1
1, COUNTi
COUNTindex(PETOT+1)d (PETOT+1)
write (*,'(2i8,f10.0)') my_rank, i, VECg(i)
enddo
endif
call MPI_FINALIZE (ierr)
stop
end

PE#0

PE#1

stride(1)t id (1)rcounts(1)displs(1)=0

stride(2)t id (2)rcounts(2)PE#2

PE#(m-2) PE#(m-1)

stride(3)t id (3)stride(m-1)
t id ( 1)stride(m)t id ( )rcounts(3)rcounts(m-1)

rcounts(m)displs(2)=displs(1) + stride(1)size(recvbuf)= displs(PETOT+1)= sum(stride)

displs(m+1)=
displs(m) + stride(m)97

T2K-FVM-03

MPI_Gatherv

局所⇒全体ベクトル生成（2/2）<$FVM>/gatherv/mpigathertest.*
do ip= 1, PETOTCOUNTindex(ip+1)= COUNTindex(ip) + COUNT(ip)enddoallocate (VECg(COUNTindex(PETOT+1)))
VECg= 0.d0&&&

call MPI_Gatherv
MPI Gatherv( VEC , N, MPI_DOUBLE_PRECISION,
VECg, COUNT, COUNTindex, MPI_DOUBLE_PRECISION,
0, MPI_COMM_WORLD, ierr)&&&call MPI_Gatherv
if (my_rank.eq.0) then
d i
do
i= 1
1, COUNTi
COUNTindex(PETOT+1)d (PETOT+1)((sendbuf,
db f scount,t sendtype,dt
recvbuf,b f rcounts,
t displs,
di lwrite (*,'(2i8,f10.0)') my_rank, i, VECg(i) recvtype, root, comm, ierr)
enddo
endif
call MPI_FINALIZE (ierr)
stop
end

PE#0

PE#1

stride(1)t id (1)rcounts(1)displs(1)=0

stride(2)t id (2)rcounts(2)PE#2

PE#(m-2) PE#(m-1)

stride(3)t id (3)stride(m-1)
t id ( 1)stride(m)t id ( )rcounts(3)rcounts(m-1)

rcounts(m)displs(2)=displs(1) + stride(1)size(recvbuf)= displs(PETOT+1)= sum(stride)

displs(m+1)=
displs(m) + stride(m)98

T2K-FVM-03

実行
FORTRAN$ cd
d <$FVM>/<$FVM>/gathervth$ mpif90 –Oss –noparallel mpigathertest.f（go.sh修正）
修$ qsub go.shC$ cd <$FVM>/gatherv$ mpiccp–Os –noparallelp
mpigathertest.cp g（go.sh修正）$ qsub go
go.sh
sh
99

T2K-FVM-03

100

手順• 各プロセッサにおける計算結果（PHI(:)に格納される）を，各プ セッサにおける計算結果（(:)に格納される）を，あるプロセッサに集める：MPI_Gatherv使用– 「内点」の結果のみで良い。

• そのプロセッサから，結果を以下の書式で，グローバル番号順に書き出す（「 」は入れない）：
順に書き出す（「,」は入れない）：<グローバル要素番号><計算結果>• ファイルを「fvm_entire_mesh.inp-geo」にアペンドしてファイルを「fvm entire mesh inp geo」にアペンドしてUCDファイルを生成する。
• ヒント– MPI_Gathervを使う– グローバルIDは局所分散通信ファイルより得る

T2K-FVM-03

101

局所分散通信ファイル（comm.0）
内点数，総要素数（内点＋外点），全体要素番号#NEIBPEtot
2#NEIBPE
1
2#IMPORT index
2
4#IMPORT items
7856#EXPORT index
24#EXPORT items
13
12#INTERNAL NODE
4#TOTAL NODE8#GLOBAL NODE ID
11 12 15 16
10 14

PE#0

全体要素番号（局所番号順）
7

14

15

16

8

3

4

10

11

12

7

1

2

7

8

5

6

8
PE#1

PE#2

T2K-FVM-03

102

要素のグローバルID
要素のグロ
バルID
hpcmw_eps_fvm_input_grid
read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') (export_index(k), k= 1, n_neighbor_pe)p_index(n( _neighbor
g
_p
pe)
)
nn= export
allocate (export_item(nn))
read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') (export_item(k), k= 1, nn),LINEread (IUNIT,'(a)')
read (IUNIT,'(6i12)') intNODE_tot
read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') nn
allocate (NODE_GLOBAL(nn))
read (IUNIT,'(a)') LINE
read (IUNIT,'(6i12)') (NODE_GLOBAL(k), k= 1, nn)close (IUNIT)!C===
end subroutine hpcmw_eps_fvm_input_gridT2K-FVM-03

103

例題
PE#2
13

14

15

PE#3

16

13
3

14
4

15
3

16
4

9
1

10
2

11
1

12
2

9

10

11

12

5

6

7

8

5
3

6
4

7
3

8
4

1

2

3

4

1
1

2
2

3
1

4
2

PE#0

PE#1

内点

T2K-FVM-03

104

例題：計算結果例：PHIの中身
PE#
PE#2

PE#3

13
3

14
4

15
3

16
4

9
1

10
2

11
1

12
2

5
3

6
4

7
3

8
4

1
1

2
2

3
1

4
2

PE#0

PE#1

要素番号
（ロ カル）（グロ バル）（ローカル）（グローバル）

PHI(i)

0
0
0
0

1
2
3
4

1
2
5
6

1.150000E+01
9.500000E+00
1.150000E+01
9.500000E+00

2
2
2
2

1
2
3
4

3
4
7
8

6.500000E+00
6500000E+00
2.500000E+00
6.500000E+00
2.500000E+00

1
1
1
1

1
2
3
4

9
10
13
14

1.150000E+01
9.500000E+00
1.150000E+01
9.500000E+00

3
3
3
3

1
2
3
4

11
12
15
16

6.500000E+00
2 500000E+00
2.500000E+00
6.500000E+00
2.500000E+00

T2K-FVM-03

105

hpcmw eps fvm global output（1/2）
hpcmw_eps_fvm_global_output（1/2）
subroutine hpcmw_eps_fvm_global_output
use hpcmw_eps_fvm_all
implicit REAL*8 (A-H,O-Z)
integer, dimension(:), allocatable :: rcounts, displs
integer, dimension(:), allocatable :: NODE_ID_G
integer, dimension(:), allocatable :: NEWtoOLD
real(kind=kreal), dimension(:), allocatable :: VAL!C!C-- INIT.
allocate (rcounts(PETOT), displs(0:PETOT))
rcounts= 0
displs = 0

MPI_Gathervの準備

call MPI_Allgather&(intNODE_tot, 1, MPI_INTEGER, rcounts, 1, MPI_INTEGER,&MPI COMM WORLD i
MPI_COMM_WORLD,
ierr))
do ip= 1, PETOT
displs(ip)= displs(ip-1) + rcounts(ip)enddo
NODE_tot_G= displs(PETOT)

領域全体の要素数（内点数の和）allocate (NODE_ID_G(NODE_tot_G), NEWtoOLD(NODE_tot_G))
allocate (VAL(NODE(VAL(NODE_tot_G))
tot G))&&

T2K-FVM-03

106

hpcmw eps fvm global output（2/2）
hpcmw_eps_fvm_global_output（2/2）!C!C-- GLOBAL ARRAY
call MPI_GATHERv&(NODE_GLOBAL, intNODE_tot,
MPI_INTEGER,&NODE_ID_G , rcounts, displs(0), MPI_INTEGER,&
0, MPI_COMM_WORLD, ierr)&&&

グローバル要素番号
グロ
バル要素番号 NODE_ID_G
NODE ID G

call MPI_GATHERv&(PHI, intNODE_tot,
MPI_DOUBLE_PRECISION,&
VAL, rcounts, displs(0), MPI_DOUBLE_PRECISION,&
0, MPI_COMM_WORLD, ierr)&&&

if (my_rank.eq.0) thendo i= 1, NODE_tot_G
j= NODE_ID_G(i)NEWtoOLD(j)= i
enddo
dd
IUNIT= 12
open (IUNIT, file= AVSfile, status='unknown', position='append')
do j
j= 1, NODE_tot_G
ii= NEWtoOLD(j)write (IUNIT,'(i8,1pe16.6)') j, VAL(ii)enddo
close (IUNIT)
endif
end subroutine hpcmw_eps_fvm_global_output

T2K-FVM-03

107

hpcmw eps fvm global output（2/2）
hpcmw_eps_fvm_global_output（2/2）!C!C-- GLOBAL ARRAY
call MPI_GATHERv&(NODE_GLOBAL, intNODE_tot,
MPI_INTEGER,&NODE_ID_G , rcounts, displs(0), MPI_INTEGER,&
0, MPI_COMM_WORLD, ierr)&&&call MPI_GATHERv&(PHI, intNODE_tot,
MPI_DOUBLE_PRECISION,&
VAL, rcounts, displs(0), MPI_DOUBLE_PRECISION,&
0, MPI_COMM_WORLD, ierr)
if (my_rank.eq.0) thendo i= 1, NODE_tot_G
j= NODE_ID_G(i)NEWtoOLD(j)= i
enddo
dd

&
&
&

計算結果ベクトル VAL

IUNIT= 12
open (IUNIT, file= AVSfile, status='unknown', position='append')
do j
j= 1, NODE_tot_G
ii= NEWtoOLD(j)write (IUNIT,'(i8,1pe16.6)') j, VAL(ii)enddo
close (IUNIT)
endif
end subroutine hpcmw_eps_fvm_global_output

T2K-FVM-03

108

この状態のVALの中身
PE番号順に格納されており，
グローバル要素番号順になっていない
PE#2

i

0
0
0
0
1
1
1
1
2
2
2
2
3
3
3
3

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

PE#3

13
3

14
4

15
3

16
4

9
1

10
2

11
1

12
2

5
3

6
4

7
3

8
4

1
1

2
2

3
1

4
2

PE#0

PE#

PE#1

要素番号
VAL(i)（グローバル）NODE_ID_G(i)( )11.150000E+01
29.500000E+00
51.150000E+01
699.500000E+00
500000E+00
36.500000E+00
42.500000E+00
76.500000E+00
82.500000E+00
91.150000E+01
109.500000E+00
13
1 150000E+01
1.150000E+01
149.500000E+00
116.500000E+00
122.500000E
2.500000E+0000
156.500000E+00
162.500000E+00

T2K-FVM-03

109

hpcmw eps fvm global output（2/2）
hpcmw_eps_fvm_global_output（2/2）!C!C-- GLOBAL ARRAY
call MPI_GATHERv&(NODE_GLOBAL, intNODE_tot,
MPI_INTEGER,&NODE_ID_G , rcounts, displs(0), MPI_INTEGER,&
0, MPI_COMM_WORLD, ierr)call MPI_GATHERv&(DELPHI, intNODE_tot,
MPI_DOUBLE_PRECISION,&
VAL, rcounts, displs(0), MPI_DOUBLE_PRECISION,&
0, MPI_COMM_WORLD, ierr)
if (my_rank.eq.0) thendo i= 1, NODE_tot_G
j= NODE_ID_G(i)NEWtoOLD(j)= i
enddo
dd

&
&
&
&
&
&

Global番号順に並べ替え（my_rank=0のみ）

IUNIT= 12
open (IUNIT, file= AVSfile, status='unknown', position='append')
do j
j= 1, NODE_tot_G
ii= NEWtoOLD(j)write (IUNIT,'(i8,1pe16.6)') j, VAL(ii)enddo
close (IUNIT)
endif
end subroutine hpcmw_eps_fvm_global_output

T2K-FVM-03

110

hpcmw eps fvm global output（2/2）
hpcmw_eps_fvm_global_output（2/2）!C!C-- GLOBAL ARRAY
call MPI_GATHERv&(NODE_GLOBAL, intNODE_tot,
MPI_INTEGER,&NODE_ID_G , rcounts, displs(0), MPI_INTEGER,&
0, MPI_COMM_WORLD, ierr)call MPI_GATHERv&(DELPHI, intNODE_tot,
MPI_DOUBLE_PRECISION,&
VAL, rcounts, displs(0), MPI_DOUBLE_PRECISION,&
0, MPI_COMM_WORLD, ierr)
if (my_rank.eq.0) thendo i= 1, NODE_tot_G
j= NODE_ID_G(i)NEWtoOLD(j)= i
enddo
dd

&
&
&

書き出し（my_rank=0のみ）

IUNIT= 12
open (IUNIT, file= AVSfile, status='unknown')
do jj= 1, NODE_tot_G
ii= NEWtoOLD(j)write (IUNIT,'(i8,1pe16.6)') j, VAL(ii)enddo
close (IUNIT)
endif
end subroutine hpcmw_eps_fvm_global_output&&&

T2K-FVM-03

111

並び替えた後のVALの中身
PE#2

PE#3

13
3

14
4

15
3

16
4

9
1

10
2

11
1

12
2

5
3

6
4

7
3

8
4

1
1

2
2

3
1

4
2

PE#0

PE#1

PE#

i

0
0
1
1
0
0
1
1
2
2
3
3
2
2
3
3

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

要素番号
VAL(i)（グローバル）
1
11.150000E+01
150000E+01
29.500000E+00
36.500000E+00
42.500000E+00
51.150000E+01
69.500000E+00
76.500000E+00
82 500000E+00
2.500000E+00
91.150000E+01
109.500000E+00
116.500000E+00
122.500000E+00
131.150000E+01
149.500000E+00
15
6 500000E+00
6.500000E+00
162.500000E+00

T2K-FVM-03

112

並列シミュレ ションにおけるI/O
並列シミュレーションにおけるI/O
領域分割機能割 能<$FVM>/run/
INPUT.DAT（固定）並列計算制御ファイル#GRIDout#GRIDout#GRIDout#GRIDout#D-GRID

局所分散メッシュデータ#GRIDout#GRIDout#GRIDout#GRIDout#D-COMM

局所分散通信データ

並列計算
sol
so

#M-RESULT

計算結果ファイル
（UCDファイル後半）

T2K-FVM-03

113

並列計算制御ファイル• INPUT.DAT（名称固定）
INPUT DAT（名称固定）• 実行形式「sol」と同じディレクトリになければならない（この場合は場合は<$FVM>/run）。
）• 全ての項目は省略不可。
../ex/mesh.rcb 局所分散メッシュファイルのヘッダ名../ex/comm.rcb 局所分散通信ファイルのヘッダ名../ex/result
可視化用出力ファイル名（後述）1
可視化用出力の有無（=1のとき出力）視有無きT2K-FVM-03

114

計算実行，ファイル生成$> cd <$FVM>/run$> cat INPUT
INPUT.DATDAT../ex/mesh.rcb../ex/comm.rcb../ex/result/ /1$> qsub go.sh$> ls -l ../ex/result これが書きだされた
result
結果ファイル$> cd ../ex$> cat fvm
fvm_entire_mesh.inp-geo
entire mesh inp geo result > test
test.inp
inp
形状部分に結果ファイルをアペンドして AVS出力用のアペンドして，AVS出力用のファイルを生成する。

T2K-FVM-03

115「eps fvm」の並列化「eps_fvm」の並列化
program eps_fvm
fuse hpcmw_eps_fvm_all
implicit REAL*8 (A-H,O-Z)
call hpcmw_eps_fvm_init
call hpcmw_eps_fvm_input_grid
call poi_gen
call
ll h
hpcmw_eps_fvm_solverflcall hpcmw_eps_fvm_global_output
call hpcmw_eps_fvm_finalize

可視化処理の部分を除くと，並列計
可視化処理部分を除くと 並列計
算特有の処理は，共役勾配法の内積 行列ベクト 積程度 それも積，行列ベクトル積程度。それも，hpcmw_eps_fvm_comm.＊のようなサブなサブルーチン群を使用すれば，一チ 群を使用すれば行の挿入で済んでしまう。

end program eps_fvm

並列計算の「隠蔽」
並列分散データ構造の適切な設計が全てが全て：一般化された通信テーブル
般化された通信テ ブル

T2K-FVM-03

116

並列化した eps fvm」
並列化した「eps-fvm」• 並列分散データ入力
並列分散デ タ入力
• 共役勾配法• 可視化• 課題

T2K-FVM-03

117

課題1
• Strong
g Scalingg– 問題サイズを固定，PE数を変化させて時間（全体，各部分）を測定。

• Weak Scaling– PEあたりの問題サイズを固定，1反復あたりの計算時間を求める。

• 考慮すべき項目– 問題サイズ– 領域分割手法（RCB，METIS，1D～3D）の影響。
領域分割手法（RCB METIS 1D 3D）の影響• 注意– 通信時間を計測する工夫をして見よ。

T2K-FVM-03

118

反復回数を固定したい場合• hpcmw_eps_fvm_solver!C!C***!C*** hpcmw_solver
!C
!C***
!C
subroutine hpcmw_eps_fvm_solver
use hpcmw_eps_fvm_all
hpcmw eps fvm all
implicit REAL*8 (A-H,O-Z)
ISET= 0
EPS = 1
1.d-8
d-8
この部分を例えば「ITR=50」のようにする。
ITR = intNODE_tot
call hpcmw_eps_fvm_allreduce_R (ITR, hpcmw_max)
ITR = PETOT * ITR
S_TIME= MPI_WTIME()
call hpcmw_eps_fvm_solver_CG&&( intNODE_tot, NODE_tot, N2, D, BFORCE, DELPHI, EPS, &&ITR
ITR, IER,
IER WAindex,
WAi dWA
WAcoef,f COMMti
COMMtime))
E_TIME= MPI_WTIME()ISET= ISET + 1

T2K-FVM-03

119

１D 3D分割
１D～3D分割
どの方法が良いか考えて見よ

1D型

2D型

3D型

T2K-FVM-03

120

１D 3D分割
１D～3D分割
通信量の総和（各辺4N，8領域とする）

1D型
16 N2 x 7 = 112 N2

2D型
16 N2 x 4 = 64 N2

3D型
16 N2 x 3 = 48 N2

T2K-FVM-03

121

１D 3D分割
１D～3D分割
80x80x80要素，8領域の計算例

1D型
16 N2 x 7 = 112 N2### solver time### comm. time### work ratio

8.867188E+00
3.488281E+00
6.066079E+01

2D型

3D型

16 N2 x 4 = 64 N2### solver time### comm. time### work ratio

6.574219E+00
1.285156E+00
8.045157E+01

16 N2 x 3 = 48 N2### solver time### comm. time### work ratio

6.230469E+00
1.031250E+00
8.344828E+01

T2K-FVM-03

122

１D 3D分割
１D～3D分割
fvmpart.ctrl

1D型

2D型

3D型

!METHOD
RCB
X,X,X!METHOD
!
O
RCB
X,Y,X

!METHOD
!
RCB
X,Y,Z

!REGION NUMBER
8!REGION NUMBER
8!REGION NUMBER
8

T2K-FVM-03

123

ジョブスクリプト：go.sh
#@$-r
$
test
#@$-q tutorial#@$-N 1#@$-J T8#@$-e err#@$-o aaa
aaa.lstlst#@$-lM 28GB#@$-lT 00:05:00#@$-s /bin/sh#@$
cd $PBS_O_WORKDIR
mpirun numactl --localalloc ./sol/sol

T2K-FVM-03

124

計算例：Strong Scaling：T2K東大• 1003要素– 問題規模固定• 4コアの場合を基準

T2K-FVM-03

125

計算例（三次元弾性問題）（T2K）• 8,957,952自由度（3×1443），Strong Scaling• 32
32～1024コア（32コア
1024コア（32コア Flat MPIの場合を基準）

T2K-FVM-03

126

Flat MPI vs. Hybrid
Flat-MPI：Each PE -> Independent

core
core

core
core

core
memorry
m

core

core
memorrymmemorry
m

core

core

core
core
core

core
core
core
core

mem
mory

core
core
core
core

mem
mory

mem
mory

Hybrid：Hierarchal Structure (OpenMP+MPI)
core
core
core
core

T2K-FVM-03

127

Flat MPI, Hybrid (4x4, 8x2, 16x1)
Higher Performance of HB16x1 is important0
1
2
3
Flat MPI

Hybrid
4x4

0

1

2

3

Hybrid
H
b id
8x2

0

1

2

3

Hybrid
y
16x1

0

1

2

3

T2K-FVM-03

128

D
Domaini D
Decomposition
itiexample: 6 nodes, 24 sockets, 96 cores

Flat MPI

HB 4x4

HB 16x1

T2K-FVM-03

理想値からのずれ• MPI通信そのものに要する時間– データを送付している時間– ノード間においては通信バンド幅によって決まる– 通信時間は送受信バッファのサイズに比例• MPIの立ち上がり時間– latency– 送受信バッファのサイズによらない• 呼び出し回数依存，プロセス数が増加すると増加する傾向– 通常，数～数十μsecのオーダー
通常 数 数十のオ ダ• MPIの同期のための時間– プロセス数が増加すると増加する傾向
プ
数が増加すると増加する傾向• 計算時間が小さい場合（問題規模が小さい場合）はこれらの効果を無視 きな
効果を無視できない。
– 特に，送信メッセージ数が小さい場合は，「Latency」が効く。

129

T2K-FVM-03

研究例• Weak Scaling– プロセッサあたりの問題規模を固定– プロセッサ数を増加させる– 「計算性能」としては不変のはずであるが，普通はプロセッサ数を増やすと性能は悪化• Nakajima,j, K. ((2007),), The Impact
pof Parallel
Programming Models on the Linear Algebra
Performance for Finite Element Simulations,, Lecture
Notes in Computer Science 4395, 334-348.
– 並列有限要素法（特並列有限要素法（特に latencyy の影響大）
影響大）
130

T2K-FVM-03● Flat-MPI DJDS
○ Hybrid DJDS
Flat-MPI(ideal)Hybrid (ideal)

Weak Scaling: LARGE
Earth Simulator

IBM SPSP-3 (Seaborg at LBNL)

1,572,864 DOF/PE(=3x128x64x64)

375,000 DOF/PE(=3x503)

4000

150

2000

2.2G
G DOF
O
3.80 TFLOPS
33.7% of peak(176 nodes,
1 408 PE’s)
1,408
PE’ )

1000

0

GFLOPS

GF
FLOPS

3000

100

0 38G DOF
0.38G
O
110 GFLOPS
7.16% of peak(128 nodes,
1 024 PE’s)
1,024
PE’ )

50

0
0

250

500

750

PE#

1000

1250

1500

0

250

500

750

1000

1250

1500

PE#

131

T2K-FVM-03● Flat-MPI DJDS
○ Hybrid DJDS
Flat-MPI(ideal)Hybrid (ideal)

Weak Scaling: SMALL
Earth Simulator

IBM SPSP-3 (Seaborg at LBNL)

98,304 DOF/PE(=3x323)

98,304 DOF/PE(=3x323)

4000

150

Hybrid is much faster
if PE# is large !!
GFLOPSGFLOPS

3000

2000

100

0 0G DOF
0.10G
O
115 GFLOPS
7.51% of peak(128 nodes,
1 024 PE’s)
1,024
PE’ )

50

1000

0

0
0

250

500

750

PE#

1000

1250

1500

0

250

500

750

1000

1250

1500

PE#

132

T2K-FVM-03

通信：メモリーコピー
実は意外にメモリの負担もかかる
do neib=d
ib 1
1, NEIBPETOT
do k= export_index(neib-1)+1, export_index(neib)
kk= export_item(k)SENDbuf(k)= VAL(kk)
enddo
enddo
do neib= 1, NEIBPETOTiS_e= export_index(neib-1) + 1
iE e= export_index(neibiE_eexport index(neib )BUFlength_i= iE_e + 1 - iS_e
iS_i= import_index(neib-1) + 1
iE_i= import_index(neib )BUFlength_i= iE_i + 1 - iS_i&&&

call MPI_SENDRECV&(SENDbuf(iS_e), BUFlength_e, MPI_INTEGER, NEIBPE(neib), 0,&
RECVbuf(iS_i), BUFlength_i, MPI_INTEGER, NEIBPE(neib), 0,&
MPI_COMM_WORLD, stat_sr, ierr)enddo
do neib= 1, NEIBPETOT
do k= import_index(neib-1)+1, import_index(neib)
kk= import
po t_item(k)te ( )
VAL(kk)= RECVbuf(k)
enddo
enddo
133

T2K-FVM-03

134

課題2
• 更にnumactlの影響について検討せよ– 問題サイズの影響
ズ• 全体の問題サイズ，1コアあたりの問題サイズ– 4ノードまで使用して見よ
4ノ ドまで使用して見よ– 「--physcpubind」の影響についても検討せよ：コアを指定できる– 1ノードのうち8コアを選択して使用する場合について検討してみよ
ド うちを選択し 使用する場合に検討しよ– プロセス番号～コア番号の関係についても検討してみよ
STREAM

T2K-FVM-03

135

計算例：Strong Scaling：T2K東大• 1003要素• 4コアの場合を基準• Policy-3適用

Policy:3# /#!/bin/bash/
MYRANK=$MXMPI_ID
MYVAL=$(expr $MYRANK / 4)
SOCKET=$(expr $MYVAL % 4)numactl --cpunodebind=$SOCKET --membind=$SOCKET $@T2K-FVM-03

136

終わりに
• 「並列化」– 分散データ構造の重要性
分散デ タ構造の重要性– 元のアプリケーションに対する深い理解が必要– アプリケーションによって並列化の戦略も異なる
アプリケ ションによ て並列化の戦略も異なる• 5月27日まで利用できます。
5月27日まで利用できます• 質問，相談等は随時（5月27日以降もOK，ただ最近レスポン
スが悪くて不評）協• アンケートにご協力ください• 領域分割ツール等– /home/t00000/pre.tar

T2K-FVM-03

137

今後の技術的動向• マルチコア• GPU• ハイブリッド並列プログラミングモデル– 例：OpenMP+MPI
例：OpenMP MPI– MPIの部分は本講習会のものと基本的に変わらない