CMSI配信講義B 第15回

CMSI配信講義B (2016)
第15回

1 /83

大規模MD並列化の技術2
名古屋大学 工学研究科
附属計算科学連携教育研究センター
特任講師
安藤 嘉倫

2016/7/28

CMSI配信講義B
CMSI配信講義B 第15回
第15回

目次

2 /83

・分子動力学 (MD) 法・ 分子動力学計算の並列化特性
・並列化技術 1

データ構造・並列化技術 2

MPI・並列化技術 3

OpenMP, SIMD

第一回

第二回

CMSI配信講義B 第15回

並列化技術 2

MPI

3 /83

MPI 並列化技術の要点(特に 3 次元トーラスネットワークに最適化)① 通信衝突の回避② 通信前後での配列間コピーの消去③ 通信の演算による代用
ただし簡単のため, 各項目の説明は主に二次元平面上で行いますCMSI配信講義B 第15回

3 次元トーラスネットワーク
トーラス
(円環面)

4 /83

通信ケーブル

円環状

myrank

計算ノード例) 京, FX10, Blue Gene, Anton

長所 : 空間ドメイン分割によるMPI並列化との相性が良い.
短所 : 斜め方向ノードとの通信が直にできないため通信の回数 (ホップ数) が増加
プロセス間通信の際に衝突が生じやすい.

MPI (message passing interface)

CMSI配信講義B 第15回
5 /83

nprocs=2

・分散メモリー型並列化のための規格
・MPICH, OpenMPI などのライブラリーをインストールすることで使用可
・「プロセス」と呼ばれる処理単位に仕事(タスク)が割り振られる・各種の MPI 関数をコードに挿入
・自動並列化という概念はない
表1

代表的なMPI関数

NIC

NIC

memorymemory

CPU

CPU

プロセス0
(myrank=0)

プロセス1(myrank=1)

CMSI配信講義B 第15回

MPI (message passing interface)

6 /83

サンプル hello_mpi.f:
include ‘mpif.h’

MPI 変数,関数の読み込み

CALL MPI_INIT(ierr)
MPI 起動プロセス数
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) (nprocs)取得
CALL MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr) 自プロセス番号(myrank)取得＊0 始まりWRITE(*,*) ‘Hello’, myrank
CALL MPI_FINALIZE(ierr)

MPI 終了

STOP
END

実行
> mpirun –np 4 ./a.out
コンパイル> mpif90 hello_mpi.f

Hello
Hello
Hello
Hello

0
1
2
3

CMSI配信講義B 第15回

mpi̲sendrecv 関数

7 /83

3Dトーラスネットワークでの通信に適する MPI 関数

mpi_send + mpi_recv

call MPI_SENDRECV (sendbuf,scount,stype,dest,stag,
recvbuf,rcount,rtype,source,rtag,comm,status,ierr)sendbuf : 送信データ
scount : 送信データ数
stype: 送信データの型
dest: 送信先プロセス番号(rank)stag: 送信タグ

recvbuf : 受信データ
rcount : 受信データ数
rtype : 受信データの型
source : 受信元プロセス番号
rtag: 受信タグ

dest=1source=3

dest=2source=0

dest=3source=1

dest=0source=2

rank0

rank1

rank2

rank3

s
3

s
0

s
1

s
2

d
1

d
2

d
3

d
0

送信

受信

comm : コミュニケータ
status : 状態
ierr
: 完了コード

commで指定したプロセス
グループ内で環状のシフト通信が可能.
通信相手は両隣のプロセス.

myrankデータの流れdestination : 送信先
source: 送信元

myrank からみると, call mpi_sendrecv() の結果
下流の rank のデータが受信される.

CMSI配信講義B 第15回

mpi̲sendrecv 関数
2次元

8 /83

rank12

rank13

rank14

rank15

s d
15 13

s d
12 14

s d
13 15

s d
14 12

rank8

rank9

rank10

rank11

s d
11 9

s d
8 10

s
9

s
10

rank4

rank5

rank6

rank7

s
7

s
4

s
5

s
6

+x方向

y

x
データの流れd5d6

myrank

d
11

d
7

d
8

d
4

rank0

rank1

rank2

rank3

s
3

s
0

s
1

s
2

d
1

d
2

d
3

d
0

CMSI配信講義B 第15回

mpi̲sendrecv 関数
2次元

9 /83

rank12

rank13

rank14

rank15

d0

d1

d2

d3

s8

s9

s10

s11

rank8

rank9

rank10

rank11

d12

d13

d14

d15

s4

s5

s6

s7

rank4

rank5

rank6

rank7

d8

d9

d10

d11

s0

s1

s2

s3

rank0

rank1

rank2

rank3

d4

d5

d6

d7

s12

s13

s14

s15

+y方向

y

x
データの流れMPI 並列化技術①: 通信衝突の回避例) 多極子の通信 (1スーパーセル/プロセス)
M2L 対象セル-5+4

myrank+4

M2M 対象セル

y-5x

M2M, M2L演算に必要なデータ範囲(白抜き部分は除く)

CMSI配信講義B 第15回
10/83

MPI 並列化技術①: 通信衝突の回避

CMSI配信講義B 第15回
11/83

例) 多極子の通信 (1スーパーセル/プロセス)
古い通信コード

コーディングイメージ

iy=+4, ix=-5~+4

do iy=-5,+4
do ix=-5,+4
ip_dest =送信先プロセス番号
ip_src =受信元プロセス番号
IF( (ix, iy) が白抜き部分 ) cycle
call mpi_sendrecv(..,ip_dest,..,ip_src,..)
enddo
enddo

iy=-2, ix=-5~+4

iy=-5, ix=-5~+4yx

M2M, M2L演算に必要なデータ範囲(白抜き部分は除く)

このままだと(103-53)+(23-1) =882 回のプロセス間通信 [3次元]MPI 並列化技術①: 通信衝突の回避

CMSI配信講義B 第15回
12 /83例) 多極子の通信 (1スーパーセル/プロセス)
古い通信コード

コーディングイメージ
myrank

8回

6回
2次元では
合計500回

9回

do iy=-5,+4
do ix=-5,+4
ip_dest =送信先プロセス番号
ip_src =受信元プロセス番号
IF( (ix, iy) が白抜き部分 ) cycle
call mpi_sendrecv(..,ip_dest,..,ip_src,..)
enddo
enddo

このままだと(103-53)+(23-1) =882 回のプロセス間通信 [3次元]さらに、トーラスネットワークに特有の問題:
(1) 882個の通信先ごと多数のホップ回数

用語の定義：

ホップ:
隣接ノードとの通信

MPI 並列化技術①: 通信衝突の回避

CMSI配信講義B 第15回
13 /83例) 多極子の通信 (1スーパーセル/プロセス)
古い通信コード

コーディングイメージ
other rankmyrank

8回

衝突
6回
2次元では
合計500回

用語の定義：

9回

ホップ:
隣接ノードとの通信

do iy=-5,+4
do ix=-5,+4
ip_dest =送信先プロセス番号
ip_src =受信元プロセス番号
IF( (ix, iy) が白抜き部分 ) cycle
call mpi_sendrecv(..,ip_dest,..,ip_src,..)
enddo
enddo

このままだと(103-53)+(23-1) =882 回のプロセス間通信 [3次元]さらに、トーラスネットワークに特有の問題:
(1) 882個の通信先ごと多数のホップ回数(2) 通信の衝突が至る所で発生通信時間の増加 = 並列性能の低下

MPI 並列化技術①: 通信衝突の回避例) 多極子の通信 (1スーパーセル/プロセス)

other rankmyrank

8回

衝突
6回
2次元では
合計500回

9回

トーラスネットワークに特有の問題:
(1) 882個の通信先ごと多数のホップ回数(2) 通信の衝突が至る所で発生

CMSI配信講義B 第15回
14 /83

新しい通信コード

1) ± x 軸方向通信

9回

全てのプロセスが・同じタイミング 衝突回避
・同じ方向

MPI 並列化技術①: 通信衝突の回避例) 多極子の通信 (1スーパーセル/プロセス)

myrank8回

15 /83

新しい通信コード

2) ± y 軸方向通信
other rank

CMSI配信講義B 第15回

9回

全てのプロセスが・同じタイミング
・同じ方向

衝突回避

衝突
6回
2次元では
合計500回

9回

トーラスネットワークに特有の問題:
(1) 882個の通信先ごと多数のホップ回数(2) 通信の衝突が至る所で発生

MPI 並列化技術①: 通信衝突の回避例) 多極子の通信 (1スーパーセル/プロセス)

CMSI配信講義B 第15回
16 /83

新しい通信コード

冗長にデータが送られてくる領域

other rankmyrank

8回

衝突

衝突回避

合計9+9=18 回
6回

2次元では
合計500回

9回

トーラスネットワークに特有の問題:
(1) 882個の通信先ごと多数のホップ回数(2) 通信の衝突が至る所で発生

３次元では
合計27回

(1) ホップ回数を最小限に抑える(2) 通信の衝突を回避

MPI 並列化技術①: 通信衝突の回避

CMSI配信講義B 第15回
17 /83x方向通信コーディングイメージ
ipx_dest =送信先プロセス番号(-x)
ipx_src =受信元プロセス番号(-x)
do ix= -1, -4, -1
call mpi_sendrecv(..,ipx_dest,..,ipx_src,..)enddo
袖部つき局所化されたデータ構造の所定位置に受信データを格納. 次回通信では前回の受信データを送受信.
+x
myrank
ix=-1

ix=-2

ix=-3

ix=-4

:送信データ(-x)

注) call mpi_sendrecv() の結果,
myrank からみた -x 方向通信により +x 方向のデータを受信

MPI 並列化技術①: 通信衝突の回避
x方向通信コーディングイメージ
ipx_dest =送信先プロセス番号(-x)
ipx_src =受信元プロセス番号(-x)
do ix= -1, -4, -1
call mpi_sendrecv(..,ipx_dest,..,ipx_src,..)enddo
袖部つき局所化されたデータ構造の所定位置に受信データを格納. 次回通信では前回の受信データを送受信.
+x
myrank
ix=-1

ix=-2

:送信データ(-x)

CMSI配信講義B 第15回
18 /83

ipx_dest =送信先プロセス番号(+x)
ipx_src =受信元プロセス番号(+x)
do ix= +1, +5, +1
call mpi_sendrecv(..,ipx_dest,..,ipx_src,..)
enddo
ix=+1

ix=+2

ix=+3

ix=-3

ix=+4

ix=-4

ix=+5

+x

:送信データ(+x)

MPI 並列化技術①: 通信衝突の回避

CMSI配信講義B 第15回
19 /83y方向通信コーディングイメージ
ipy_dest =送信先プロセス番号(-y)
ipy_src =受信元プロセス番号(-y)
do iy= -1, -4, -1
call mpi_sendrecv(..,ipy_dest,..,ipy_src,..)
enddo

iy=-1

+y

myrank:送信データ(-y)+x

袖部つき局所化されたデータ構造の所定位置に棒状の受信データを格納. 次回通信では前回の棒状の受信データを送受信.

iy=-2

MPI 並列化技術①: 通信衝突の回避

CMSI配信講義B 第15回
20 /83y方向通信コーディングイメージ
ipy_dest =送信先プロセス番号(-y)
ipy_src =受信元プロセス番号(-y)
do iy= -1, -4, -1
call mpi_sendrecv(..,ipy_dest,..,ipy_src,..)
enddo

+y

iy=-3

+x

袖部つき局所化されたデータ構造の所定位置に棒状の受信データを格納. 次回通信では前回の棒状の受信データを送受信.

iy=-4

MPI 並列化技術①: 通信衝突の回避

CMSI配信講義B 第15回
21 /83y方向通信コーディングイメージ
ipy_dest =送信先プロセス番号(-y)
ipy_src =受信元プロセス番号(-y)
do iy= -1, -4, -1
call mpi_sendrecv(..,ipy_dest,..,ipy_src,..)enddo

ipy_dest =送信先プロセス番号(+y)
ipy_src =受信元プロセス番号(+y)
do iy= +1, +5, +1
call mpi_sendrecv(..,ipy_dest,..,ipy_src,..)
enddo

+y

iy=+1

+x

:送信データ(+y)

MPI 並列化技術①: 通信衝突の回避

CMSI配信講義B 第15回
22 /83y方向通信コーディングイメージ
ipy_dest =送信先プロセス番号(+y)
ipy_src =受信元プロセス番号(+y)
do iy= +1, +5, +1
call mpi_sendrecv(..,ipy_dest,..,ipy_src,..)
enddo
+y

iy=+2

iy=+3

+x

:送信データ(+y)

MPI 並列化技術①: 通信衝突の回避+y y方向通信コーディングイメージ

CMSI配信講義B 第15回
23 /83
iy=+5

iy=+4

:送信データ(+y)+xMPI 並列化技術①: 通信衝突の回避

ipx_dest =送信先プロセス番号(-x)
固定
ipx_src =受信元プロセス番号(-x)
do ix= -1, -4, -1
call mpi_sendrecv(..,ipx_dest,..,ipx_src,..)enddo
ipx_dest =送信先プロセス番号(+x)
固定
ipx_src =受信元プロセス番号(+x)
do ix= +1, +5, +1
call mpi_sendrecv(..,ipx_dest,..,ipx_src,..)enddo
ipy_dest =送信先プロセス番号(-y)
固定
ipy_src =受信元プロセス番号(-y)
do iy=-1, -4, -1
call mpi_sendrecv(..,ipy_dest,..,ipy_src,..)enddo
ipy_dest =送信先プロセス番号(+y)
固定
ipy_src =受信元プロセス番号(+y)
do iy=+1, +5, +1
call mpi_sendrecv(..,ipy_dest,..,ipy_src,..)
enddo
-z

24 /83

新しい通信コード

コーディングイメージ

残るz方向も同様

CMSI配信講義B 第15回

+z

冗長にデータが送られてくる領域-x
myrank

衝突回避+x

合計9+9=18 回

-y

+y

３次元では
合計27回

３軸逐次隣接通信による局所的
収集処理

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去

演算部

旧データ構造新データ構造（メタデータ•袖部付き局所化）

wk_x

meta_x

通信対象を抽出/頭詰め
通信バッファー配列へコピー

通信部

trans_x

受信バッファー配列からのコピー
元の配列末尾へ格納

演算部

25 /83

wk_x

最初に通信する1軸方向のみ最小限のコピー

meta_x
最初に通信する1軸方向のみ最小限のコピー

meta_xCMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去-x 軸方向通信

ipx_mdest =送信先プロセス番号(-x)
ipx_msrc =受信元プロセス番号(-x)
do ix=-1,-2,-1c 原子数情報の送受信
icbufm( ) = na_per_cell( 3, 3 )call mpi_sendrecv(icbufm, ...,ipx_mdest, ...,
ircbufm,..., ipx_msrc, ...)
na_per_cell( 4, 3 )=ircbufm( )受信データのtag(4,3)を作成
送信バッファーへc 座標情報の送受信
コピー
buffm( ) = meta_x( )call mpi_sendrecv(buffm,..,ipx_mdest,..,rbuffm, ...,ipx_msrc,..)meta_x( ) = rbuffm( )enddo
データを左づめにしつつ
元配列にコピー注) myrank からみると -x 方向通信により+x 方向の座標情報を受信[サイズ不定のセル単位データを送受信]

ix=-1

meta_x(n)Y

コーディングイメージ

26 /83

myrank(3,3) (4,3)

na1cell

meta_x(1):送信データ(-x)

X

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去-x 軸方向通信

ipx_mdest =送信先プロセス番号(-x)
ipx_msrc =受信元プロセス番号(-x)
do ix=-1,-2,-1c 原子数情報の送受信
icbufm( ) = na_per_cell( 4, 3 )call mpi_sendrecv(icbufm, ...,ipx_mdest, ...,
ircbufm,..., ipx_msrc, ...)
na_per_cell( 5, 3 )=ircbufm( )受信データのtag(5,3)を作成
送信バッファーへc 座標情報の送受信
コピー
buffm( ) = meta_x( )call mpi_sendrecv(buffm,..,ipx_mdest,..,rbuffm, ...,ipx_msrc,..)meta_x( ) = rbuffm( )enddo
データを左づめにしつつ
元配列にコピー注) myrank からみると -x 方向通信により+x 方向の座標情報を受信[サイズ不定のセル単位データを送受信]

ix=-2

meta_x(n)Y

コーディングイメージ

27 /83

myrank(4,3) (5,3)

na1cell

meta_x(1):送信データ(-x)

X

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去+x 軸方向通信

ipx_pdest =送信先プロセス番号(+x)
ipx_psrc =受信元プロセス番号(+x)
do ix=+1,+2,+1c 原子数情報の送受信
icbufp( ) = na_per_cell( 3, 3 )call mpi_sendrecv(icbufp, ...,ipx_pdest, ...,
ircbufp,..., ipx_psrc, ...)
na_per_cell( 2, 3 )=ircbufp( )受信データのtag(2,3)を作成
送信バッファーへc 座標情報の送受信
コピー
buffp( ) = meta_x( )call mpi_sendrecv(buffp,..,ipx_pdest,..,rbuffp, ...,ipx_psrc,..)meta_x( ) = rbuffp( )enddo
データを右づめにしつつ
元配列にコピー注) myrank からみると +x 方向通信により–x 方向の座標情報を受信[サイズ不定のセル単位データを送受信]

ix=+1

meta_x(n)Y

コーディングイメージ

28 /83

myrank(2,3) (3,3)

na1cell

meta_x(1):送信データ(+x)

X

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去+x 軸方向通信

ipx_pdest =送信先プロセス番号(+x)
ipx_psrc =受信元プロセス番号(+x)
do ix=+1,+2,+1c 原子数情報の送受信
icbufp( ) = na_per_cell( 2, 3 )call mpi_sendrecv(icbufp, ...,ipx_pdest, ...,
ircbufp,..., ipx_psrc, ...)
na_per_cell( 1, 3 )=ircbufp( )受信データのtag(1,3)を作成
送信バッファーへc 座標情報の送受信
コピー
buffp( ) = meta_x( )call mpi_sendrecv(buffp,..,ipx_pdest,..,rbuffp, ...,ipx_psrc,..)meta_x( ) = rbuffp( )enddo
データを右づめにしつつ
元配列にコピー注) myrank からみると +x 方向通信により–x 方向の座標情報を受信[サイズ不定のセル単位データを送受信]

ix=+2

meta_x(n)Y

コーディングイメージ

29 /83

myrank(1,3) (2,3)

na1cell

meta_x(1):送信データ(+x)

X

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去-y 軸方向通信

ipy_mdest =送信先プロセス番号(-y)
ipy_msrc =受信元プロセス番号(-y)
棒状データを構成する
do iy=-1,-2,-1
サブセル数を指定
c 原子数情報の送受信[サイズ一定の棒状データを送受信]

iy=-1

call mpi_sendrecv(na_per_cell(1,3),5,...,ipy_mdest, ...,
na_per_cell(1,4),5,..., ipy_msrc, ...)
受信データのtag(1:5,4)を作成(1,3)セルの先頭
c 座標情報の送受信
アドレスを指定(1,4)

call mpi_sendrecv(meta_x( ),naline,..,ipy_mdest,..,meta_x( ),naline,...,ipy_msrc,..)

(1,3)

enddo

meta_x(n)Y

コーディングイメージ

30 /83(3,3)(1,4)セルの先頭アドレスを指定

na1cell
注) myrank からみると -y 方向通信により+y 方向の座標情報を受信

naline=5*na1cell

meta_x(1):送信データ(-y)

X

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去-y 軸方向通信

ipy_mdest =送信先プロセス番号(-y)
ipy_msrc =受信元プロセス番号(-y)
棒状データを構成する
do iy=-1,-2,-1
サブセル数を指定
c 原子数情報の送受信
call mpi_sendrecv(na_per_cell(1,4),5,...,ipy_mdest, ...,
na_per_cell(1,5),5,..., ipy_msrc, ...)
受信データのtag(1:5,5)を作成(1,4)セルの先頭
c 座標情報の送受信
アドレスを指定[サイズ一定の棒状データを送受信]

iy=-2

(1,5)
(1,4)

call mpi_sendrecv(meta_x( ),naline,..,ipy_mdest,..,meta_x( ),naline,...,ipy_msrc,..)enddo

meta_x(n)Y

コーディングイメージ

31 /83(3,3)(1,5)セルの先頭アドレスを指定注) myrank からみると -y 方向通信により+y 方向の座標情報を受信

naline=5*na1cell

meta_x(1):送信データ(-y)

X

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去+y 軸方向通信

ipy_pdest =送信先プロセス番号(+y)
ipy_psrc =受信元プロセス番号(+y)
棒状データを構成する
do iy=+1,+2,+1
サブセル数を指定
c 原子数情報の送受信[サイズ一定の棒状データを送受信]

iy=+1

meta_x(n)Y

コーディングイメージ

32 /83

call mpi_sendrecv(na_per_cell(1,3),5,...,ipy_pdest, ...,
na_per_cell(1,2),5,..., ipy_psrc, ...)
受信データのtag(1:5,2)を作成(1,3)セルの先頭
c 座標情報の送受信
アドレスを指定
call mpi_sendrecv(meta_x( ),naline,..,ipy_pdest,..,meta_x( ),naline,...,ipy_psrc,..)enddo(1,2)セルの先頭アドレスを指定注) myrank からみると +y 方向通信により-y 方向の座標情報を受信

(1,3)

(3,3)

(1,2)
naline=5*na1cell

meta_x(1):送信データ(+y)

X

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去+y 軸方向通信

ipy_pdest =送信先プロセス番号(+y)
ipy_psrc =受信元プロセス番号(+y)
棒状データを構成する
do iy=+1,+2,+1
サブセル数を指定
c 原子数情報の送受信[サイズ一定の棒状データを送受信]

iy=+2

meta_x(n)Y

コーディングイメージ

33 /83

call mpi_sendrecv(na_per_cell(1,2),5,...,ipy_pdest, ...,
na_per_cell(1,1),5,..., ipy_psrc, ...)
受信データのtag(1:5,1)を作成(1,2)セルの先頭
c 座標情報の送受信
アドレスを指定
call mpi_sendrecv(meta_x( ),naline,..,ipy_pdest,..,meta_x( ),naline,...,ipy_psrc,..)enddo(1,1)セルの先頭アドレスを指定注) myrank からみると +y 方向通信により-y 方向の座標情報を受信

(3,3)
(1,2)
(1,1)
meta_x(1)naline=5*na1cell

X:送信データ(+y)

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去
コーディングイメージ
ipz_mdest =送信先プロセス番号(-z)
ipz_msrc =受信元プロセス番号(-z)
do iz=-1,-2,-1c 原子数情報の送受信-z 軸方向通信[サイズ一定の面状データを送受信]

z
面状データを構成する
サブセル数を指定

call mpi_sendrecv(na_per_cell( ),25,...,ipz_mdest, ...,
na_per_cell( ),25,..., ipz_msrc, ...)
受信データのtagを作成
c 座標情報の送受信(1,5,5)(1,1,5)(5,5,5)(5,1,5)(1,5,4)(1,1,4)(5,5,4)(5,1,4)(1,5,3)

call mpi_sendrecv(meta_x( ),narea,..,ipz_mdest,..,meta_x( ),narea,...,ipz_msrc,..)enddo(5,5,3)

myrank(1,1,3)(5,1,3)(1,5,2)注) myrank からみると -z 方向通信により+z 方向の座標情報を受信

34 /83(1,1,2)

narea=25*na1cell(1,1,1)(5,5,2)(5,1,2)(1,5,1)(5,5,1)(5,1,1)x

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去
コーディングイメージ
ipz_mdest =送信先プロセス番号(-z)
ipz_msrc =受信元プロセス番号(-z)
do iz=-1,-2,-1c 原子数情報の送受信-z 軸方向通信[サイズ一定の面状データを送受信]

z
面状データを構成する
サブセル数を指定

call mpi_sendrecv(na_per_cell(1,1,3),25,...,ipz_mdest, ...,
na_per_cell(1,1,4),25,..., ipz_msrc, ...)
受信データのtag(1:5,1:5,4)を作成(1,1,3)セルの先頭
c 座標情報の送受信
アドレスを指定35 /83(1,5,5)(1,1,5)(5,5,5)iz=-1(5,1,5)(1,5,4)(1,1,4)(5,5,4)(5,1,4)(1,5,3)(5,5,3)

call mpi_sendrecv(meta_x( ),narea,..,ipz_mdest,..,meta_x( ),narea,...,ipz_msrc,..)enddo(1,1,4)セルの先頭アドレスを指定注) myrank からみると -z 方向通信により+z 方向の座標情報を受信(1,1,3)(5,1,3)(1,5,2)(1,1,2)

narea=25*na1cell(1,1,1):送信データ(-z)(5,5,2)(5,1,2)(1,5,1)(5,5,1)(5,1,1)x

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去
コーディングイメージ
ipz_mdest =送信先プロセス番号(-z)
ipz_msrc =受信元プロセス番号(-z)
do iz=-1,-2,-1c 原子数情報の送受信-z 軸方向通信[サイズ一定の面状データを送受信]

z
面状データを構成する
サブセル数を指定

call mpi_sendrecv(na_per_cell(1,1,4),25,...,ipz_mdest, ...,
na_per_cell(1,1,5),25,..., ipz_msrc, ...)
受信データのtag(1:5,1:5,5)を作成(1,1,4)セルの先頭
c 座標情報の送受信
アドレスを指定36 /83(1,5,5)(1,1,5)(5,5,5)(5,1,5)(1,5,4)(1,1,4)(5,5,4)(5,1,4)(1,5,3):送信データ(-z)(5,5,3)

call mpi_sendrecv(meta_x( ),narea,..,ipz_mdest,..,meta_x( ),narea,...,ipz_msrc,..)enddo(1,1,5)セルの先頭アドレスを指定注) myrank からみると -z 方向通信により+z 方向の座標情報を受信(1,1,3)(5,1,3)(1,5,2)(1,1,2)

narea=25*na1cell(1,1,1)iz=-2(5,5,2)(5,1,2)(1,5,1)(5,5,1)(5,1,1)x

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去
コーディングイメージ
ipz_mdest =送信先プロセス番号(-z)
ipz_msrc =受信元プロセス番号(-z)
do iz=-1,-2,-1c 原子数情報の送受信-z 軸方向通信[サイズ一定の面状データを送受信]

z
面状データを構成する
サブセル数を指定

call mpi_sendrecv(na_per_cell(1,1,3),25,...,ipz_mdest, ...,
na_per_cell(1,1,2),25,..., ipz_msrc, ...)
受信データのtag(1:5,1:5,2)を作成(1,1,3)セルの先頭
c 座標情報の送受信
アドレスを指定37 /83(1,5,5)(1,1,5)(5,5,5)iz=+1(5,1,5)(1,5,4)(1,1,4)(5,5,4)(5,1,4)(1,5,3)(5,5,3)

call mpi_sendrecv(meta_x( ),narea,..,ipz_mdest,..,meta_x( ),narea,...,ipz_msrc,..)enddo(1,1,2)セルの先頭アドレスを指定注) myrank からみると +z 方向通信により-z 方向の座標情報を受信(1,1,3)(5,1,3)(1,5,2)(1,1,2)

narea=25*na1cell(1,1,1):送信データ(+z)(5,5,2)(5,1,2)(1,5,1)(5,5,1)(5,1,1)x

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去
コーディングイメージ
ipz_mdest =送信先プロセス番号(-z)
ipz_msrc =受信元プロセス番号(-z)
do iz=-1,-2,-1c 原子数情報の送受信-z 軸方向通信[サイズ一定の面状データを送受信]

z
面状データを構成する
サブセル数を指定

call mpi_sendrecv(na_per_cell(1,1,2),25,...,ipz_mdest, ...,
na_per_cell(1,1,1),25,..., ipz_msrc, ...)
受信データのtag(1:5,1:5,2)を作成(1,1,2)セルの先頭
c 座標情報の送受信
アドレスを指定38 /83(1,5,5)(1,1,5)(5,5,5)iz=+2(5,1,5)(1,5,4)(1,1,4)(5,5,4)(5,1,4)(1,5,3)(5,5,3)

call mpi_sendrecv(meta_x( ),narea,..,ipz_mdest,..,meta_x( ),narea,...,ipz_msrc,..)enddo(1,1,3)(1,1,1)セルの先頭アドレスを指定(5,1,3)(1,5,2)(1,1,2)注) myrank からみると +z 方向通信により-z 方向の座標情報を受信

narea=25*na1cell→ すべての袖部のデータが受信された(1,1,1)(5,5,2)(5,1,2)(1,5,1):送信データ(+z)(5,5,1)(5,1,1)x

CMSI配信講義B 第15回

MPI 並列化技術②:通信前後での配列間コピーの消去

39 /83汎用性のため実際のコードでは以下の点にも対応.
・任意のmyrank所持サブセルブロック厚さ2n (n≥0)
・直方体形状のmyrank所持サブセルブロック
narea=96*na1cell
narea=36*na1cell

21
myrank23
22

21
myranknaline=6*na1cell

allocate( tag(1:6,1:6) )allocate( na_per_cell(1:6,1:6) )naline=12*na1cell

allocate( tag(1:12,1:8) )allocate( na_per_cell(1:12,1:8) )注) 開発経緯の都合上, 実際のコードでは z → y → x の順に通信しています.
comm_3.f, comm_fmm.f

CMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

40 /83

代用の目的 (1): プログラム全体の計算時間の短縮
冗長な多重演算時間 < 単一演算結果の通信時間

代用
現在のスパコン構成では, しばしば生じる関係.

代用の目的 (2): デバッグのしやすいコードの作成
・多重演算だが, 分かりやすいコード

代用・単一演算だが, 複雑なコード
上のコードが完成した上で, 下のコードを派生させる.
最初から下のコードを作ることは難易度高. (腕に覚えある人はどうぞ)CMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

41 /83

冗長な多重演算の例
・分子内相互作用計算[目的(2)]

・分子間近距離相互作用計算[目的(1),(2)]・FMM 上位階層の M2M/L2L[目的(1)]CMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

42 /83

冗長な多重演算: 分子内相互作用計算
dihedral	

j	

i	

n=2

j	

k	

K! "#1+ cos(n! (ri , rj , rk , rl ) ! " )$%

l	
n=1

n=3

k	
l	

do n=1,ndihedrals
オリジナルコード
phi=phi(ri, rj, rk, rl)
ポテンシャルの計算
Fi, Fj, Fk, Fl の計算
f(i)=f(i)+Fi
f(j)=f(j)+Fj
f(k)=f(k)+Fk
f(l)=f(l)+Fl
enddo 	

k	

j	

l	

i	

i	

n=4
n=6

n=5

CMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

43 /83

冗長な多重演算: 分子内相互作用計算
myrank (rank=0)
dihedral	

i	
j	

n=2

i	

j	

k	
l	

K! "#1+ cos(n! (ri , rj , rk , rl ) ! " )$%＊ELSEの場合, 実際には合力をまとめて
該当プロセスへ通信

rank=1 へ Fk, Fl の通信が必要な原子

k	
l	

n=1

n=3
rank=1 へ Fl の通信が必要な原子MPI並列化コード
do n=1,ndihedrals(myrank)
phi=phi(ri, rj, rk, rl)
ポテンシャルの計算
Fi, Fj, Fk, Fl の計算
IF(ri in myrank) f(i)=f(i)+Fi, ELSE Fiを通信
IF(rj in myrank) f(j)=f(j)+Fj, ELSE Fjを通信
IF(rk in myrank) f(k)=f(k)+Fk, ELSE Fkを通信
IF(rl in myrank) f(l)=f(l)+Fl, ELSE Flを通信
enddo 	

rank=1

rank=3 へ Fj の通信が必要な原子rank=3

rank=4 へ Fi の通信が必要な原子rank=4CMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

44 /83

冗長な多重演算: 分子内相互作用計算

i0 : myrank内の原子通し番号

myrankdihedral	

i0=1

i	
j	

i	

rank=1
i0=5

j	

k	
l	

i0=6

i0=2

演算冗長MPI並列化コード

md_charmm_f90.f

l	

i0=3

K! "#1+ cos(n! (ri , rj , rk , rl ) ! " )$%do i0=1,natom(myrank)dihedral番号を i0 より逆引き
phi=phi(ri, rj, rk, rl)
ポテンシャルの計算
Fi の計算f(i)=f(i)+Fi
enddo 	

k	

i0=7
i0=4

4重の多重計算. しかし・分かりやすいため多重足し込み等の間違いをしない
・i0ごとの演算独立性からOpenMP, SIMD並列が容易
・力の通信なし
rank=3

rank=4CMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

45 /83

冗長な多重演算: 分子間近距離相互作用計算
二体力 Fij = !Fji

例) Lennard-Jones, Coulombの粒子対計算

do icell(myrank)
do jcell_list(myrank or otherranks)
do i=na_per_cell(icell)
do j=na_per_cell(jcell)
rij=rij(ri, rj)
カットオフ判定
ポテンシャルの計算
Fij の計算＊実際にはここでIF文は使
わず, DOループ外で判定
f(i)=f(i)+Fij
IF(rj in myrank) f(j)=f(j)-Fij
ELSE storeF(j)=storeF(j)-Fij
enddo ! jenddo ! ienddo ! jcell
Fij の通信
enddo ! icellが必要な原子storeF(j)を対象プロセスへ通信

jcell
myrankicellother ranks

CMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

46 /83

冗長な多重演算: 分子間近距離相互作用計算
二体力 Fij = !Fji

例) Lennard-Jones, Coulombの粒子対計算

do icell(myrank)
do jcell_list(myrank or otherranks)
do i=na_per_cell(icell)
do j=na_per_cell(jcell)
rij=rij(ri, rj)
カットオフ判定
ポテンシャルの計算
Fij の計算f(i)=f(i)+Fij！i 原子のみ足し込み
enddo ; enddo ; enddo ; enddo

jcell
myrankicell
2重の多重計算. しかし・力の通信なし
・OpenMP, SIMD並列が容易 = 演算効率良other ranks

実演算時間 < 実通信時間の範囲で有効 (例えばノード数が多く演算時間が短いとき)

CMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

47 /83

冗長な多重演算: FMM 上位階層の M2M

レベル 2

myrankM2M (L1→L2)

M2M 対象セル

レベル 1

M2M (L0→L1)
レベル０
myrankM2M 対象セルCMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

48 /83

冗長な多重演算: FMM 上位階層の M2M
other ranks

レベル 2

myrankM2M (L1→L2)

M2M 対象セル

レベル 1

M2M (L0→L1)
レベル０
myrankM2M 対象セルCMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

49 /83

冗長な多重演算: FMM 上位階層の M2M
mpi_bcast

単一演算の場合:

レベル 2

通信6+演算1M2M (L1→L2)

M2M 対象セル

レベル 1
mpi_gather

M2M (L0→L1)
レベル０
myrankM2M 対象セルCMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

50 /83

冗長な多重演算: FMM 上位階層の M2M
冗長演算の場合:

レベル 2

通信3+演算4M2M (L1→L2)
レベル 1
mpi_sendrecv

M2M (L0→L1)
レベル０
myrankM2M 対象セルCMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

51 /83

冗長な多重演算: FMM 上位階層の L2L
FMM-ewald

myrank

L2L (L2→L1)

レベル 2

L2L 対象セル

レベル 1

L2L (L1→L0)
L2L 対象セル

レベル０

CMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

52 /83

冗長な多重演算: FMM 上位階層の L2L

other ranks

mpi_bcast

単一演算の場合:

レベル 2

通信3+演算9
FMM-ewald

L2L (L2→L1)
レベル 1

L2L (L1→L0)
レベル０

CMSI配信講義B 第15回

MPI 並列化技術③: 通信の演算による代用

53 /83

冗長な多重演算: FMM 上位階層の L2L
冗長演算の場合:

FMM-ewald

FMM-ewald

レベル 2

通信0+演算12
FMM-ewald

注) 最上位層の wm は
M2Mの過程でプロセス間で冗長に保持FMM-ewald

L2L (L2→L1)
レベル 1

L2L (L1→L0)
レベル０

CMSI配信講義B 第15回

並列化技術 3
演算効率化の前提
・データの連続化
・ブロック化によるキャッシュの有効利用
OpenMP 並列化技術
・スレッド間のロードバランス調整・スレッド並列前後処理の削減
SIMD 並列化技術
・IF文の削除・ベクトル長の確保

54 /83CMSI配信講義B 第15回

演算効率化のために京) 16GB

Main memory

55 /83

CPUを効率よく動作させるためには
L1 キャッシュの有効利用が必要.

100京) 6MBL2 cache

MD計算では, 例えば・L1上相手 j 原子座標の再利用促進
・L1上 wm, M2L 変換行列の再利用促進

10京) 32KB

L1 cache1京) 128GFLOPS

前提(1) メモリ上の演算対象データの連続化(2) ブロック化CPU

相対アスセス速度CMSI配信講義B 第15回

データの連続化 [1] 座標

56 /83

従来の配列
wk_x(n)キャッシュの有効利用不可j 原子座標へのランダムアクセス
rcutjother ranksからの転送データ

原子 j	

原子 i	
myrankの所持データ

wk_x(1)iデータの連続化 [1] 座標

CMSI配信講義B 第15回
57 /83

メタデータ配列

前回説明

meta_x: X-Y 平面上での原子の相対位置関係 (= メタデータ ) を保持

Y

meta_x(n)
(1,5)

(2,5)

(3,5)

(4,5)

(5,5)

(1,4)

(2,4)

(3,4)

(4,4)

(5,4)

(1,3)

(2,3)

(3,3)

(4,3)

(5,3)
原子 i	

(1,2)

(2,2)

(3,2)

(4,2)

(5,2)

(1,1)

(2,1)

(3,1)

(4,1)

(5,1)

meta_x(1)X1) 自プロセスを中心にデータを局所化
2) 袖部に直にデータ装填

データの連続化 [1] 座標

CMSI配信講義B 第15回
58 /83

メタデータ配列

前回説明

meta_x: X-Y 平面上での原子の相対位置関係 (= メタデータ ) を保持

Y

meta_x(n)
(1,5)

(2,5)

(3,5)

(4,5)

(5,5)

(1,4)

(2,4)

(3,4)

(4,4)

(5,4)

(1,3)

(2,3)

(3,3)

(4,3)

(5,3)

帯単位での j 原子座標への連続アクセス

原子 j	

原子 i	
(1,2)

(2,2)

(3,2)

(4,2)

(5,2)

(1,1)

(2,1)

(3,1)

(4,1)

(5,1)

meta_x(1)X1) 自プロセスを中心にデータを局所化
2) 袖部に直にデータ装填

データの連続化 [2] 多極子

CMSI配信講義B 第15回
59 /83

前回説明

従来の配列

袖部付き局所化配列

M2L 対象セル
myrankM2M 対象セル

データへのアクセス
に不連続性

1) 自プロセスを中心にデータを局所化
2) 袖部に直にデータ装填

データの連続化 [2] 多極子

CMSI配信講義B 第15回
60 /83

前回説明

従来の配列

袖部付き局所化配列

M2L 対象セル
myrankM2M 対象セル

データへのアクセス
に不連続性

データへの
連続アクセス

1) 自プロセスを中心にデータを局所化
2) 袖部に直にデータ装填CMSI配信講義B 第15回

ブロック化によるキャッシュの有効利用 (1)61 /83

ブロック化前のループ構造:
do icell(myrank)
do jcelldo iatm=tag(icell),tag(icell)+na_per_cell(icell)-1
do jatm=tag(jcell),tag(jcell)+na_per_cell(jcell)-1
ポテンシャルの計算
力の計算
enddo
enddo
enddo
enddo

jcell 範囲
myrankicell

座標データ転送済範囲

CMSI配信講義B 第15回

ブロック化によるキャッシュの有効利用 (1)62 /83

ブロック化前のループ構造:
do icell(myrank)
do jcelldo iatm=tag(icell),tag(icell)+na_per_cell(icell)-1
do jatm=tag(jcell),tag(jcell)+na_per_cell(jcell)-1
ポテンシャルの計算
力の計算
Main memory
enddo
100
enddo
L2 cache10
enddo
L1 cacheenddo
1

jcell 範囲
icell

myrankicell’
jcell’ 範囲

CPU

問題点:

座標データ転送済範囲(1) 右図灰色部分のデータは icell → icell’ と変わった結果メモリから再ロードされる
よって、いったんキャッシュに乗ったデータを使い切っていない(2) jcell 内原子数が少なく最内ベクトル長が稼げない. よって, SIMDの性能が出ない

CMSI配信講義B 第15回

ブロック化によるキャッシュの有効利用 (1)63 /83

ブロック化後のループ構造:
jcellを外側へ移動

do jcell_line
do icell(myrank)[along icell_line]
do iatom=tag(icell),tag(icell)+na_per_cell(icell)-1
do jatom=tag(jcell),tag(jcell+4)+na_per_cell(jcell+4)-1
ポテンシャルの計算
力の計算
Main memory
enddo
100
L2 cache
enddo
10
enddo
L1 cacheenddo
1tag(jcell)tag(jcell+4)+na_per_cell(jcell+4)

jcell_line

icell_line myrankCPU

座標データ転送済範囲(1) jcell_lineの原子座標データは 1 回のみロードされる.
すなわち、いったんキャッシュに乗ったデータを使い切る.
(2) jcell_line 内原子数は５倍のベクトル長, かつ連続. よってSIMDの性能向上

CMSI配信講義B 第15回

ブロック化によるキャッシュの有効利用 (1)64 /83

最外DOループによる jcell_line シフトの概念図

icellのシフト

icell=jcellでは
IF分岐

以下同様

実際のコードは3次元のシフトに対応

md_direct_f90.f

CMSI配信講義B 第15回

ブロック化によるキャッシュの有効利用 (2)65 /83

ブロック化前のループ構造 (M2L):
do icell(myrank)
do icy=1,10
do icx=1,10
if(icx,icyが近接2セル以遠の領域にあれば)
do m1=1,(nmax+1)2
do m2=1,(nmax+1)2
wl_local=wl_local+m2l*wm_local
enddo
enddo
endif
enddo
enddo
enddo

wm_local((nmax+1)2,10,10)
M2L 対象セル

myrankicell

多極子データ転送済範囲

CMSI配信講義B 第15回

ブロック化によるキャッシュの有効利用 (2)66 /83

ブロック化前のループ構造 (M2L):
do icell(myrank)
do icy=1,10
do icx=1,10
if(icx,icyが近接2セル以遠の領域にあれば)
do m1=1,(nmax+1)2
do m2=1,(nmax+1)2
wl_local=wl_local+m2l*wm_local
enddo
enddo
endif
enddo
enddo
enddo

wm_local((nmax+1)2,10,10)
M2L 対象セル

myrankicell

多極子データ転送済範囲

CMSI配信講義B 第15回

ブロック化によるキャッシュの有効利用 (2)67 /83

ブロック化前のループ構造 (M2L):
do icell(myrank)
do icy=1,10
do icx=1,10
if(icx,icyが近接2セル以遠の領域にあれば)
do m1=1,(nmax+1)2
do m2=1,(nmax+1)2
wl_local=wl_local+m2l*wm_local
enddo
enddo
endif
enddo
enddo
enddo

wm_local((nmax+1)2,10,10)

myrankicell

問題点・ icell が変わるごと右図灰色の領域のwm_local, 変換行列m2lはメモリから再ロード.
よって、いったんキャッシュに乗ったデータを使い切っていない.

CMSI配信講義B 第15回

ブロック化によるキャッシュの有効利用 (2)68 /83

ブロック化後のループ構造 (M2L):
iblk=1
do iblk=1,nblockdo icy=icyblkst(iblk),icyblkend(iblk)
do icx=icxblkst(iblk),icxblkend(iblk)
do icell(myrank)
if(icx,icyが近接2セル以遠の領域にあれば)
do m1=1,(nmax+1)2
do m2=1,(nmax+1)2
wl_local=wl_local+m2l*wm_local
enddo
enddo
endif
enddo
enddo
enddo
enddo

icxblkst(1)icxblkend(1)icyblkend(1)M2L

iblk=1

iblk=2
icyblkst(1)myrankiblk=3

iblk=4・iblk番目のブロック内wm_local, 変換行列m2lは icell に対し 1 回のみロードされる.
すなわち、いったんキャッシュに乗ったデータを使い切る.

CMSI配信講義B 第15回

ブロック化によるキャッシュの有効利用 (2)69 /83

do iblk=1,nblockdo icy=icyblkst(iblk),icyblkend(iblk)
do icx=icxblkst(iblk),icxblkend(iblk)
do icell(myrank)
if(icx,icyが近接2セル以遠の領域にあれば)
do m1=1,(nmax+1)2
do m2=1,(nmax+1)2
wl_local=wl_local+m2l*wm_local
enddo
enddo
endif
enddo
enddo
enddo
enddo

icxblkst(2)icxblkend(2)M2L

iblk=1

iblk=2

icyblkst(2)iblk=2

icyblkend(2)ブロック化後のループ構造 (M2L):

myrankiblk=3

iblk=4・iblk番目のブロック内wm_local, 変換行列m2lは icell に対し 1 回のみロードされる.
すなわち、いったんキャッシュに乗ったデータを使い切る.

OpenMP 並列化技術

CMSI配信講義B 第15回

OpenMP 言語拡張機能

70 /83NIC

memory・共有メモリー型並列化のための規格
3
2
core
core
・商用コンパイラ (ifort, pgf90, frtpx) に標準的に備わっている
スレッド
core
core
CPU 1
0
→ 表2のコンパイルオプションを入れることで有効化
・「スレッド」と呼ばれる単位にタスクが割り振られる・!$omp (Fortran) ないし#pragma omp (C言語) ではじまる各種の指示文(ディレクテブ) をプログラムに追記する [do/for ループごと]・コンパイラに自動 OpenMP 並列化のオプションがある→ ただし実用上は配列への初期値代入程度にしか使えない
表2

OpenMP, SIMD並列化に関わるコンパイルオプションCMSI配信講義B 第15回

OpenMP 並列化技術

71 /83

サンプル hello_omp.f:
include ‘omp_lib.h’

OpenMP変数,関数の読み込み

nomp = omp_get_max_threads() スレッド数(nomp)取得 ← 環境変数 OMP_NUM_THREADS!$omp paralleliam = omp_get_thread_num()!$omp doDo i=1,nomp
WRITE(*,*) ‘Hello’, iam
Enddo!$omp end do!$omp end parallel

自スレッド番号(iam)取得 ＊0 始まり

STOP
END
コンパイル
> ifort -openmp hello_omp.f

実行> export OMP_NUM_THREADS=4> ./a.out
Hello
Hello
Hello
Hello

0
1
2
3

CMSI配信講義B 第15回

OpenMP 並列化技術

72 /83・ロードインバランス調整 (M2L)!$omp parallel!$omp dodo iblk=1,nblockdo icy=icyblkst(iblk),icyblkend(iblk)
do icx=icxblkst(iblk),icxblkend(iblk)
do icell(myrank)
if(icx,icyが近接2セル以遠の領域にあれば)
do m1=1,(nmax+1)2
do m2=1,(nmax+1)2
wl_local=wl_local+m2l*wm_local
enddo
enddo
endif
enddo
enddo
enddo
enddo
!$omp end do!$omp end parallel

あらかじめM2L対象セルをスレッド数にあわせ均等にブロック化しておく.

任意のスレッド数に対応できないという問題

CMSI配信講義B 第15回

OpenMP 並列化技術

73 /83・ロードインバランス調整 (M2L)!$omp parallel!$omp do schedule(static,nchunk)
do load=1,nload
あらかじめM2L対象セル番地をlddirに
icx=lddir(1,load)
登録しておく. nchunk個ごとスレッドに割り当てる.
icy=lddir(2,load)
do icell(myrank)
if(icx,icyが近接2セル以遠の領域にあれば)
do m1=1,(nmax+1)2・任意のスレッド数に対応
do m2=1,(nmax+1)2
・異方的セル分割(均等2ベキ以外)にも対応可能
wl_local=wl_local+m2l*wm_local
enddo
enddo
endif
enddo
enddo
!$omp end do!$omp end parallel
md_fmm_f90.f

CMSI配信講義B 第15回

OpenMP 並列化技術
・スレッド並列前後処理の削減

do jcell_line
do icell[along icell_line]!$omp parallel!$omp dodo iatom=tag(icell),tag(icell)+na_per_cell(icell)-1
do jatom=tag(jcell),tag(jcell+4)+na_per_cell(jcell+4)-1
ポテンシャルの計算
力の計算
enddo
enddo
!$omp enddo!$omp end parallel
enddo
enddo

74 /83CMSI配信講義B 第15回

OpenMP 並列化技術

75 /83・スレッド並列前後処理の削減$!omp parallelを大外に出す

do jcell_line
do icell[along icell_line]!$omp parallel!$omp dodo iatom=tag(icell),tag(icell)+na_per_cell(icell)-1
do jatom=tag(jcell),tag(jcell+4)+na_per_cell(jcell+4)-1
ポテンシャルの計算
力の計算
enddo
enddo
!$omp enddo!$omp end parallel
enddo
enddo

jcell_line数*icell数回のparallel領域open/closeのオーバヘッドが削減!$omp parallel
do jcell_line
do icell[along icell_line]!$omp dodo iatom=tag(icell),tag(icell)+na_per_cell(icell)-1
do jatom=tag(jcell),tag(jcell+4)+na_per_cell(jcell+4)-1
ポテンシャルの計算
力の計算
enddo
enddo
!$omp enddo
enddo
enddo
!$omp end parallel

SIMD 並列化技術

CMSI配信講義B 第15回
76 /83

data x 4

SIMD (single instruction multiple data)
core

・複数のデータに対する単一命令実行
・ハードウェアに固有の拡張命令セットを利用 (/proc/cpuinfo参照)例) intelSSE pentium III~
128bit (32bit X 4, 64 bit X 2)
AVX
Sandy Bridge~ 256bit (32bit X 8, 64 bit X 4)・SIMD 化を行なう方法:
・コンパイラの SIMD 自動並列化機能にまかせる
コンパイルメッセージ, プロファイル結果を読みながら,
Fortran/C コードを調整し最適化→ 一般人
・intrinsic 関数でコーディング→ プロフェッショナル＊SIMDの高効率化はデータの連続化、キャッシュの有効利用がされている前提
表2

OpenMP, SIMD並列化に関わるコンパイルオプション+,×CMSI配信講義B 第15回

SIMD 並列化技術

77 /83

SIMD自動並列化の例
SIMD自動並列化無しコンパイル>ifort -O0 simd.f

サンプル simd.f:>time ./a.out
real
0m1.782s
user
0m1.274s
sys
0m0.508s

real(8)::a(10000,10000),b(10000)
real(8)::c(10000)
a=1d0
b=1d0
SIMD自動並列化コンパイル>ifort -xHOST -vec-report simd.f
do i=1,10000
simd.f(3): (col. 2) remark: LOOP WAS VECTORIZED.
行列ベクトル積c(i)=0d0
simd.f(4): (col. 2) remark: LOOP WAS VECTORIZED.
C(i)=ΣA(i,j)*B(j)do j=1,10000
simd.f(5): (col. 2) remark: PERMUTED LOOP WAS
VECTORIZED.
c(i)=c(i)+a(i,j)*b(j)このコンパイルメッセージを読み
enddo
ながらコードを調整し最適化して>time ./a.out
いく.
enddo
real
0m0.718s
stop
end

user
sys

0m0.397s
0m0.321s

・計算の多い部分がVECTORIZEDされているか？

CMSI配信講義B 第15回

SIMD 並列化技術
・IF文の削除, ベクトル長の確保

78 /83分子のループではOpenMP並列の粒度,
および最内ベクトル長が確保できない
分子間話をだいぶ遡ぼれば・・・

分子内

1-4 scale	

分子間1-2 void	

do imol=1,nmol-1
do jmol=imol+1,nmol
do i=1,natom(imol)
do j=1,natom(jmol)
rij=rij(ri, rj)
LJカットオフ判定
φnonbond=φnonbond+φij
f(i)=f(i)+Fi
f(j)=f(j)+Fj
enddo
enddo
enddo
enddo	

分子内 nonbond 項計算の注意点:

分子内do imol=1,nmol
do i=1,natom(imol)-1
do j=i+1,natom(imol)
rij=rij(ri, rj)
LJカットオフ判定
0 if 1-2,-3 voidx = s if 1-4 scale
1 elseφnonbond=φnonbond+x*φij
f(i)=f(i)+x*Fi
f(j)=f(j)+x*Fj
enddo
enddo
enddo	

・1番目および2番目の隣接原子とは相互作用しない (1-2, 1-3 void)・3番目の隣接原子との相互作用は因子 s でスケールする (1-4 scale) [s は LJ, Coulomb べつ]・4番目以降の隣接原子とは通常の相互作用

CMSI配信講義B 第15回

SIMD 並列化技術

79 /83分子間と分子内の統合

・IF文の削除, ベクトル長の確保

話をだいぶ遡ぼれば・・・

分子内

1-4 scale	

分子間1-2 void	

分子内 nonbond 項計算の注意点:

原子ループでOpenMP
並列の粒度を確保

do jcell_line
do icell [along icell_line]
do iatom=tag(icell),tag(icell)+na_per_cell(icell)-1
do jatom=tag(jcell),tag(jcell+4)+na_per_cell(jcell+4)-1
rij=rij(ri,rj)
LJカットオフ判定
jcell_line原子ループ
でベクトル長を確保
0 if 1-2,-3 voidx=s if 1-4 scale
1 elseφnonbond=φnonbond+x*φij
f(i)=f(i)+x*Fi
代償として if 文がf(j)=f(j)+x*Fj
ループ最内側に移動
enddo
enddo
enddo
enddo

・1番目および2番目の隣接原子とは相互作用しない (1-2, 1-3 void)・3番目の隣接原子との相互作用は因子 s でスケールする (1-4 scale), s は LJ, Coulomb べつ・4番目以降の隣接原子とは通常の相互作用

CMSI配信講義B 第15回

SIMD 並列化技術

80 /83

・IF文の削除, ベクトル長の確保
jcell_line原子ループ do iatom=tag(icell),do jcell_line
でベクトル長を確保
do icell [along icell_line]
tag(icell)+na_per_cell(icell)-1
do iatom=tag(icell),do jatom=1,voidpair123(iatom)tag(icell)+na_per_cell(icell)-1
rij=rij(ri,rj)
do jatom=tag(jcell),
φnonbond=φnonbond–φij
tag(jcell+4)+na_per_cell(jcell+4)-1
f(i)=f(i) – Fi
voidを引き算
rij=rij(ri,rj)
f(j)=f(j) – Fjif(rij≥rcut) LJ_epsilon=0d0
マスク処理
enddo
φnonbond=φnonbond+φij
を利用
do jatom=1,scalepair14(iatom)f(i)=f(i)+Fi
rij=rij(ri,rj)
f(j)=f(j)+Fj
x=1-s
ひとまずvoid,scaleの
enddo
φnonbond=φnonbond–x*φij
区別無くすべて計算
scale との差分を
enddo
f(i)=f(i) – x*Fi
引き算
enddo
f(j)=f(j) – x*Fj
enddo
enddo
enddo

赤字：専用計算機 (MDGRAPE) の分野で使われてきたテクニック

問題点
桁落ちによる精度低下 (倍精度なら実用上問題なし).

SIMD並列化技術
コンパイラメッセージ

CMSI配信講義B 第15回
81 /83frtpx -Qt md_direct_f90.f<<< Loop-information Start >>><<< [OPTIMIZATION]
<<<
SIMD
<<<
SOFTWARE PIPELINING<<< Loop-information End >>>
140
9
p
v
do j0=tag(jzb-2,jyb,jxb),
141
9&tag(jzb+2,jyb,jxb)
142
9&+ na_per_cell(jzb+2,jyb,jxb)-1
143
9
p
v
rx=xi-wkxyz(1,j0)
144
9
p
v
ry=yi-wkxyz(2,j0)
145
9
p
v
rz=zi-wkxyz(3,j0)
146
9
p
v
r2=rx*rx+ry*ry+rz*rz
149
9! ^^^ spherical cut-off ^^^
150
10
p
v
if(r2<=cutrad2) then
151
10
p
v
eps=epsilon_sqrt_i0
15210&*epsilon_sqrt_table(ic,iam)
153
10
p
v
else
154
10
p
v
eps=0d0155
10pv
endif !cut-off
167
9
p
v
sUlj12=sUlj12+Ulj12
168
9
p
v
sUlj6 =sUlj6 +Ulj6
184
9
p
v
sUcoulomb=sUcoulomb+Ucoulomb
170
9
p
v
stlcx=stlcx+tlx
171
9
p
v
stlcy=stlcy+tly
172
9
p
v
stlcz=stlcz+tlz
185
9
p
v
stlcx=stlcx+tcx
186
9
p
v
stlcy=stlcy+tcy
187
9
p
v
stlcz=stlcz+tcz
188
9
p
v
ic=ic+1
189
9pv
enddo !j0

SIMD並列化技術
コンパイラメッセージ
1050
1051
1052
1053
1091

1
1
1
1
4

p
p
p
p

1092

5

p

1093
1094
1095
1096
1097
1098
1105
1106

6
6
6
6
6
5
1
1

p
p

6v
6v

p
p

6v

p

CMSI配信講義B 第15回
82 /83frtpx -Qt md_fmm_f90.f
DO load = 1, nload
ic = lddir(1,load)
jc = lddir(2,load)
kc = lddir(3,load)**** multipole to local translation<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
PREFETCH: 6<<<
wwl_localx: 6<<< Loop-information End >>>do m1=1,(nmax+1)*(nmax+1)<<< Loop-information Start >>><<< [OPTIMIZATION]
<<<
SIMD
<<<
SOFTWARE PIPELINING<<< Loop-information End >>>do m2=1,(nmax+1)*(nmax+1)
wwl_localx(m1,icz0,icy0,icx0,iam)$= wwl_localx(m1,icz0,icy0,icx0,iam)$+ wm_localx(m2,icz1,icy1,icx1)*shml(m2,m1,kc,jc,ic,nl)
enddo
enddo
ENDDO ! load

まとめ

CMSI配信講義B 第15回
83 /83・3次元トーラスネットワークに最適化したMPI並列化手法について, 通信間衝突の回避,通信前後での配列間コピーの消去,および通信の演算による代用をキーワードに解説した.
・演算効率化の前提となる,データの連続化およびブロック化によるキャッシュの有効利用について解説した.
・OpenMP並列化技術 (スレッド間のロードバランス調整,スレッド並列前後処理の削減) およびSIMD並列化技術 (IF文の削除,ベクトル長の確保) について解説した.

CMSI配信講義B 第15回

付録

/83

www.modylas.org
・マニュアル(PDF)
・チュートリアル(PDF)のダウンロード

登録後このようなメールが送られてくる.
リンク先からソースコードをダウンロード	

ユーザー登録

使い方については
第21回CMSI神戸ハンズオン:
MODYLASチュートリアル 資料なども参照。