= inside_global[0] + 2*mgn + 1
であることに注意
ソースコード（１）
サンプルコード: openacc_fdtd/
OpenACCを利用したFDTD法（電磁波解析）
158
openacc_fdtd/01_original CPU
コード。openacc_fdtd/02_openacc1 calc_ex_ey, pml_boundary_ex, pml_boundary_ey,
がOpenACC
。openacc_fdtd/03_openacc2
時間更新ループ全体がOpenACC
。openacc_fdtd/04_openacc3
初期化を含めOpenACC
。openacc_fdtd/05_openacc4
データ移動の最適化。ソースコード（２）
それぞれのファイルの内容
159
main.c
プログラムのメインコードfdtd2d.{c, h} 2次元 FDTD の 計算コード
fdtd2d_sources.{c, h}
入射光設定のための関数setup.c
計算条件の設定と変数の初期化config.{c, h}
物理定数の定義output.{cc, h}
計算結果出力のための関数bitmap* BMPファイル作成のための関数
本講習では、“main.c”、“fdtd2d.c”、” fdtd2d_sources.c”、”setup.c” のソー スコードを追記・修正していきます。
2次元波動伝搬
成分: Ex、Ey、Hz y
方向下側から平面波を入射計算条件
160
dx = lx/nx dy = ly/ny
デフォルト設定:
nx = 512 ny = 512 mgn = 8 lnx = 529 lny = 529
電磁波の境界での非物理的な反射を 防ぐための吸収境界条件
(PML)
物体を設置
入射光
プログラム中では下記の変数が使わ れているので注意
inside_global.length[0] = nx0 inside_global.length[1] = ny0 whole_global.length[0]
= nx0 + 2*mgn + 1 whole_global.length[1]
= ny0 + 2*mgn + 1
真空と異なる 媒質を設定
コード全体の流れ（main.c 内）
161
計算領域の設定
Range 構造体
inside_global, whole_global,
inside, whole
電場E の時間発展
calc_ex_ey
関数電場
E の境界条件
pml_boundary_ex
関数pml_boundary_ey
関数 電場の入射（入射光）plane_wave_incidence
関数磁場
H の時間発展
calc_hz
関数磁場
H の境界条件
pml_boundary_hz
関数 媒質と物体の初期化init_relative_permittivity関数 init_object関数
媒質と物体の設定（ユーザ定義）
set_object_er
関数(main.cで定義) CPUメモリに配列を確保
ex, ey, hz, ...
初期化と計算条件の設定
init_vars, set_initial_condition
関数init_pml_vars
関数set_pml_initial_condition
関数set_pml_rer
関数計算結果出力
MPI_Gather,write_bmp関数
F D TD
計算繰り返し時間発展計算結果出力
MPI_Gather収集, write_bmp関数
磁場
H の袖領域更新
MPI_Send/MPI_Recv
磁場
Ex の袖領域更新
MPI_Send/MPI_Recv
※
MPI
はOFF
にして あります計算領域の設定（１）
Range 構造体
計算領域の始点と大きさを保持162 // config.h
struct Range { int length[2];
int begin [2];
};
// main.c
const struct Range inside_global = { { atoi(argv[1]), atoi(argv[2]) }, { 0, 0 } };
const struct Range whole_global = { { inside_global.length[0] + 2*mgn + 1, inside_global.length[1] + 2*mgn + 1},
{ inside_global.begin[0] - mgn , inside_global.begin[1] - mgn } };
const struct Range inside = { { inside_global.length[0], inside_global.length[1]/nsubdomains }, { 0,
inside_global.length[1]/nsubdomains * rank } };
const struct Range whole = { { inside.length[0] + 2*mgn + 1, inside.length[1] + 2*mgn + 1},
{ inside.begin[0] - mgn , inside.begin[1] - mgn } };
全領域の中心領域 全領域の 全体領域
分割領域の 中心領域
分割領域の 全体領域
計算領域の設定（２）
Range 構造体
計算領域の始点と大きさを保持163 struct Range {
int length[2];
int begin [2];
};
const struct Range inside = { { inside_global.length[0], inside_global.length[1]/nsubdomains }, { 0,
inside_global.length[1]/nsubdomains * rank } };
const struct Range whole = { { inside.length[0] + 2*mgn + 1, inside.length[1] + 2*mgn + 1},
{ inside.begin[0] - mgn , inside.begin[1] - mgn } };
座標(whole.begin[0], whole.begin[1]) 座標(inside.begin[0], inside.begin[1])
inside.length[0] = nx
inside.length[1] = ny
whole.length[0] = nx + 2*mgn + 1 whole.length[1] = ny + 2*mgn + 1
プログラムでは下記の通り
分割領域の 中心領域
分割領域の 全体領域
配列の確保
物理変数配列は main.c で確保
多くの配列は whole.length[0] * whole.length[1]
ex_global, ey_global, hz_global はファイル出力に使うため、
whole_global.length[0] * whole_global.length[1]
164 // main.c
const int nelems = whole.length[0] * whole.length[1];
const int nelems_x = whole.length[0];
const int nelems_y = whole.length[1];
const size_t size = sizeof(FLOAT)*nelems;
const size_t size_x = sizeof(FLOAT)*nelems_x;
const size_t size_y = sizeof(FLOAT)*nelems_y;
const size_t size_global = sizeof(FLOAT)* whole_global.length[0] * whole_global.length[1];
FLOAT *ex = (FLOAT *)malloc(size); // 電場Ex FLOAT *ey = (FLOAT *)malloc(size); // 電場Ey FLOAT *hz = (FLOAT *)malloc(size); // 磁場Hz ...
// For output
FLOAT *ex_global = (FLOAT *)malloc(size_global);
FLOAT *ey_global = (FLOAT *)malloc(size_global);
FLOAT *hz_global = (FLOAT *)malloc(size_global);
時間発展（１）
前半
電場Eの時間発展（calc_ex_ey）、境界条件(pml_boundary_...)
入射光（plane_wave_incidence）165 while (icnt < nt) {
MPI_Status status;
const int tag = 0;
const int nhalo = whole.length[0];
const int inside_end1 = inside.begin[1] + inside.length[1];
const int src_hz = whole.length[0] * (inside_end1 - whole.begin[1] - 1);
const int dst_hz = whole.length[0] * (inside.begin[1] - whole.begin[1] - 1);
MPI_Send(&hz[src_hz], nhalo, MPI_FLOAT_T, rank_up , tag, MPI_COMM_WORLD);
MPI_Recv(&hz[dst_hz], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD, &status);
calc_ex_ey(&whole, &inside, hz, cexly, ceylx, ex, ey);
pml_boundary_ex(&whole, &inside, hz, cexy, cexyl, rer_ex, ex, exy);
pml_boundary_ey(&whole, &inside, hz, ceyx, ceyxl, rer_ey, ey, eyx);
const int j_in = 0;
plane_wave_incidence(&whole, &inside, time, j_in, wavelength, ex, ey);
time += 0.5*dt;
(後半へ）
※
MPI
はOFF
にして あります時間発展（２）
後半
磁場Hの時間発展（calc_hz）、境界条件(pml_boundary_hz)166
（前半から）
const int src_ex = whole.length[0] * (inside.begin[1] - whole.begin[1]);
const int dst_ex = whole.length[0] * (inside_end1 - whole.begin[1]);
MPI_Send(&ex[src_ex], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD);
MPI_Recv(&ex[dst_ex], nhalo, MPI_FLOAT_T, rank_up , tag, MPI_COMM_WORLD, &status);
calc_hz(&whole, &inside, ey, ex, chzlx, chzly, hz);
pml_boundary_hz(&whole, &inside, ey, ex, chzx, chzxl, chzy, chzyl, hz, hzx, hzy);
time += 0.5*dt;
icnt++;
(出力など) }
※
MPI
はOFF
にして あります167
チャレンジ課題:GPUを用いた
FDTD法による電磁波伝搬計算
の実習（C言語版のみ）
プログラムのコンパイルと実行（１）
CPUコードのコンパイルと実行
なお、pjsub ./run_no_out.sh すると出力なしで実行する。性能測定用。
168
$ module load nvidia ompi-cuda
$ cd openacc_mpi_fdtd/01_original
$ make
$ pjsub ./run.sh
$ cat run.sh.??????.out Rank 0: hostname = a090 Rank 1: hostname = a090 Rank 2: hostname = a091 Rank 3: hostname = a091 Calculation condition
nx_global = 512
（省略）
icnt = 4900, time = 2.3115e-14 [sec]
icnt = 5000, time = 2.3587e-14 [sec]
---Domain = 512 x 512
nsubdomains = 4 output_file = 1
Time = 4.103535 [sec]
---?
の数字はジョブごと に変わります。openacc_fdtd/01_original
計算領域サイズ、領域 分割数、出力の有無、
計算時間
利用したノード
MPIのmoduleが必要
プログラムのコンパイルと実行（２）
プログラムの実行時オプション
mpirun -np <nprocs> ../run <nx> <ny> <nprocs> <nt> <nout>
nprocs: 全ランク数（=分割数） ※今回は１ nx, ny: 計算領域サイズ
nt: 全時間ステップ
nout: 出力を行うタイムステップ数。50 の場合、50ステップに1回
出力する。0 を指定すると出力しない。169
$ cat run.sh
#! /bin/sh
#PBS -q h-tutorial
#PBS -l select=1:mpiprocs=1:ompthreads=0
（省略）
mkdir -p sim_run cd sim_run
nprocs=1
mpirun -np $nprocs ../run 512 512 $nprocs 5000 50
openacc_fdtd/01_original
計算結果の表示
計算結果は sim_run に BMP として出力される
$ cd sim_run/
計算結果の表示
1枚のBMPを見る
$ display e05000.bmp
複数のBMPファイルをアニメーションで表示$ animate *.bmp
なお
ssh -Y [email protected]
と
–Y をつけていないと表示されない。うまく表示できない場合は画像を
手元にコピーして表示してください。
170
openacc_fdtd/01_original
計算結果の例
出力されたBMPファイルの一例
Ex （電場の x 成分）の出力
171
x y
z
入射光
実習1
calc_ex_ey, pml_boundary_ex, pml_boundary_ey を OpenACC化しましょう。
Makefile
コンパイルオプションの修正 main.c
OpenACCヘッダーの追加
data 指示文の追加
fdtd2d.c
kernels 指示文、loop 指示文の追加
172
openacc_fdtd/02_openacc1
解答例は、
実行速度が遅くても、動 くプログラムである状態 を保ちながらOpenACC 化します。
末端の関数から
OpenACC化するのがよ
いでしょう。kernels, loop指示文
fdtd2d.c 内の関数
173 void calc_ex_ey(const struct Range *whole, const struct Range *inside,
const FLOAT *hz, const FLOAT *cexly, const FLOAT *ceylx, FLOAT *ex, FLOAT *ey) {
const int nx = inside->length[0];
const int ny = inside->length[1];
const int mgn[] = { inside->begin[0] - whole->begin[0], inside->begin[1] - whole->begin[1] };
const int lnx = whole->length[0];
#pragma acc kernels present(hz, cexly, ex)
#pragma acc loop independent for (int j=0; j<ny+1; j++) {
#pragma acc loop independent for (int i=0; i<nx; i++) {
const int ix = (j+mgn[1])*lnx + i+mgn[0];
const int jm = ix - lnx;
//ex[ix] += cexly[ix]*(hz[ix]-hz[jm]) - cexlz[ix]*(hy[ix]-hy[km]);
ex[ix] += cexly[ix]*(hz[ix]-hz[jm]);
} } (省略）
}
実習2
main 関数内の while 内をすべて OpenACCにしましょう。
main.c
data 指示文の移動と copyin
などの最適化 fdtd2d.c
残りの関数にkernels 指示文、loop 指示文の追加 fdtd2d_sources.c
kernels 指示文、loop 指示文の追加
174
openacc_fdtd/03_openacc2
解答例は、
data 指示文
main関数のwhile 外に data を移動
175
#pragma acc data ¥
copyin(ex[0:nelems], ey[0:nelems], hz[0:nelems]) ¥
copyin(cexly[0:nelems], ceylx[0:nelems], chzlx[0:nelems], chzly[0:nelems]) ¥ copyin(exy[0:nelems], eyx[0:nelems], hzx[0:nelems], hzy[0:nelems]) ¥
copyin(cexy[0:nelems_y], ceyx[0:nelems_x], chzx[0:nelems_x], chzy[0:nelems_y]) ¥ copyin(cexyl[0:nelems_y], ceyxl[0:nelems_x], chzxl[0:nelems_x], chzyl[0:nelems_y]) ¥ copyin(obj[0:nelems], er[0:nelems]) ¥
copyin(rer_ex[0:nelems], rer_ey[0:nelems]) {
while (icnt < nt) { MPI_Status status;
const int tag = 0;
const int nhalo = whole.length[0];
const int inside_end1 = inside.begin[1] + inside.length[1];
const int src_hz = whole.length[0] * (inside_end1 - whole.begin[1] - 1);
const int dst_hz = whole.length[0] * (inside.begin[1] - whole.begin[1] - 1);
#pragma acc host_data use_device(hz) {
MPI_Send(&hz[src_hz], nhalo, MPI_FLOAT_T, rank_up , tag, MPI_COMM_WORLD);
MPI_Recv(&hz[dst_hz], nhalo, MPI_FLOAT_T, rank_down, tag, MPI_COMM_WORLD, &status);
}
実習3
初期化を含めて全てOpenACCにします。ただし、set_object_er がCPU上のユーザ 定義関数のため、これ以降の初期化関数をOpenACCにします。
main.c
data 指示文の移動と最適化（多くが create になるはずです）
setup.c
kernels 指示文、loop 指示文の追加
176
openacc_fdtd/04_openacc3
解答例は、
実習4
計算領域のサイズなどを変更して性能測定してみましょう。
OpenACCコードをさらに最適化しましょう。
NVCOMPILER_ACC_TIMEも活用しましょう。
実は単純にfdtd2d.c に kernels と loop を入れても、いくつかの関数で暗
黙の
copyin
が発生します。これも修正していきましょう。openacc_fdtd/05_openacc4
177解答例は、
$ make calc_ex_ey:
25, Generating present(ex[:],cexly[:]) Generating implicit copyin(mgn[:]) Generating present(hz[:])
27, Loop is parallelizable 29, Loop is parallelizable
Accelerator kernel generated Generating Tesla code
27, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */
29, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */
37, Generating present(ey[:],ceylx[:]) Generating implicit copyin(mgn[:])
実習5
実はUnified Memoryを利用するとずっと簡単に実装できます。
実習2・3で行ったdata 指示文に関する実装する必要がないため -ta=tesla,cc80,managed として性能を比較してみましょう
性能がだいぶ違うと思います。改善するならどうすべきでしょうか。 managed memory はページ単位（結構大きい）でデータ転送を行います
05_openacc4
で使っているupdate 指示文は、data
指示文で確保済みのデータを
CPU-GPU間でコピーするための指示文です。
openacc_fdtd/06_openacc5
178性能改善済みの解答例は、
#pragma acc update host(ex[src:sendnelems],ey[src:sendnelems],hz[src:sendnelems]) for(i = 0;i < sendnelems;i++){
ex_global[dst+i] = ex[src+i];
ey_global[dst+i] = ey[src+i];
hz_global[dst+i] = hz[src+i];
}
src
から始まりsendnelems 個の 要素をhost(CPU)にコピーするQ & A
アカウントは1ヶ月有効です。
資料のPDF版はWEBページに掲載します。
https://www.cc.u-tokyo.ac.jp/events/lectures/157/
アンケートへの協力をお願いします。
179
補足スライド
180
性能を出すためにはスレッド数>>コア数
推奨スレッド数
CPU： スレッド数=コア数 (高々数十スレッド)
GPU： スレッド数>=コア数*4~ (数万~数百万スレッド)
最適値は他のリソースとの兼ね合いによる 理由：高速コンテキストスイッチによるメモリレイテンシ隠し
CPU : レジスタ・スタックの退避はOSがソフトウェアで行う(遅い)
GPU : ハードウェアサポートでコストほぼゼロ
メモリアクセスによる暇な時間(ストール)に他のスレッドを実行181
1core=1
スレッドのときメモリ
read
開始 メモリread
終了1core=N
スレッドのとき階層的スレッド管理とコミュニケーション
階層的なコア/スレッド管理
P100は56 SMを持ち、1 SMは64 CUDA coreを持つ。トータル3584 CUDA core
1 SMが複数のスレッドブロックを担当し、1
CUDA core が複数スレッドを担当
スレッド間のコミュニケーション
同一スレッドブロック内のスレッドは高速コ ミュニケーション可能
異なるスレッドブロックに属するスレッド間 はコミュニケーションが低速
いったんメモリに書き出したり、CPUに処理を戻 さなくてはならない182
cited from :
http://cuda- programming.blogspot.jp/2012/12/thread-hierarchy-in-cuda-programming.html
スレッドブロック
Warp 単位の実行
連続した32スレッドを1単位 = Warp と呼ぶ
このWarpは足並み揃えて動く
実行する命令は32スレッド全て同じ
データは違ってもいい
183
4 3 5 … 8 0
スレッド
1 2 3 … 31 32
配列
A
配列
B 2 3 1 … 1 9
× × × … × ×
4 3 5 … 8 0
スレッド
1 2 3 … 31 32
配列
A
配列
B 2 3 1 … 1 9
÷ ×
＋… − ×
OK ！ NG ！
Warp内分岐
Divergent Branch
Warp 内で分岐すること。Warp単位の分岐ならOK。
184
: :
if ( TRUE ) { : :
} else { : : } :
:
: :
if (
奇数スレッド) { : :
} else { : : } :
:
一部スレッドを眠らせて全分岐を実行
最悪ケースでは 32 倍のコスト
else
部分は実行せずジャンプCUDA 8
以前のバージョンCUDA 9
以上では多少マシになるが、ペナルティがあることに変わりはない
Divergent branch
なしDivergent branch
ありコアレスドアクセス
同じWarp内のスレッド（連続するスレッド）は近いメモリアドレスへ アクセスすると効率的
コアレスドアクセス（coalesced access）と呼ぶ
メモリアクセスは128 Byte 単位で行われる。128 Byte に収まれば1回の アクセス、超えれば128 Byte アクセスをその分繰り返す。185