OpenMPによる並列化実装
八木 学（理化学研究所 計算科学研究センター）

KOBE HPC Spring School 20192019年3月14日

スレッド並列とプロセス並列
スレッド並列OpenMP、自動並列化
プロセス
スレッド

スレッド

スレッド

プロセス並列
MPI
プロセス

プロセス

メモリ

メモリ

スレッド

プロセス間通信
PrivatePrivatePrivateGlobal
メモリ空間Private

プロセス

プロセス

メモリ

メモリ

並列計算について・CPUをN個使って並列計算した時、計算速度がN倍になるのが理想だが・・・- 並列化率の問題（アムダールの法則）- 通信時間ボトルネック

・京のような大型計算機を有効利用するためには、如何に計算速度をN倍に近づけられるかが重要アムダールの法則
プログラムの並列化できる割合を P とし、プロセッサ数をn とすると、並列計算した時の性能向上率は1𝑃
1−𝑃 +𝑛で与えられる。
９割並列化できるが、１割逐次処理が残ってしまうような場合、どれだけプロセッサを投入しても計算速度は１０倍以上にはならない。

OpenMPとは・共有メモリ型計算機用の並列計算API→ノード内の並列（ノード間は不可）
・ユーザーが明示的に並列のための指示を与える→コンパイラの自動並列とは異なる
・標準化された規格であり、広く使われている・指示行の挿入で並列化できるため、比較的手軽

OpenMPによるスレッド並列
Fork-Joinモデル
スレッド０
Fork

スレッド０

スレッド1

スレッド2

スレッド3

並列リージョン

Join

スレッド０

...
非並列
Fork
!$omp parallel...
並列リージョン
!$omp end parallel
...
非並列

Join

...
非並列
#pragma omp parallel{...
並列リージョン}
Join
...
非並列処理

Fork

スレッド数の指定
・シェルの環境変数で与える（推奨）export OMP_NUM_THREADS=16 （bashの場合)
setenv OMP_NUM_THREADS 16 （tcshの場合）・プログラム内部で設定することも可能#include <omp.h>
omp_set_num_threads(16);スレッド数の指定
・シェルの環境変数で与える（推奨）export OMP_NUM_THREADS=16 （bashの場合)
setenv OMP_NUM_THREADS 16 （tcshの場合）・プログラム内部で設定することも可能!$use omp_lib
call omp_set_num_threads(16)コンパイル・コンパイルオプションでOpenMPを有効にする
gcc -fopenmp test.f90

icc -qopenmp test.f90
・オプションを指定しない場合はOpenMPの指示行は無視される。
#pragma omp parallel forfor (i=0; i<100; i++) {a[i] = b[i] + c;}

指示行はFortranの場合 !$OMP から始まる。
行頭の ! は通常はコメントを意味する。
Cの場合は #pragma omp という形になるが、オプションを入れない場合、通常処理されない。

コンパイル・コンパイルオプションでOpenMPを有効にする
gfortran -fopenmp test.f90

ifort -qopenmp test.f90
・オプションを指定しない場合はOpenMPの指示行は無視される。
!$OMP PARALLEL DO
do i = 1, 100a(i) = b(i) + c
enddo!$OMP END PARALLEL DO

指示行はFortranの場合 !$OMP から始まる。
行頭の ! は通常はコメントを意味する。
Cの場合は #pragma omp という形になるが、オプションを入れない場合、通常処理されない。

OpenMPの基本関数
OpenMPモジュール/ヘッダをロード[C] #include <omp.h>[F] use omp_lib*OpenMP関連の関数を使用するためのおまじない!$ use omp_lib
integer :: myid, nthreads#include <omp.h>
int myid, nthreads;

nthreads = omp_get_num_threads()
myid = omp_get_thread_num()

nthreads = omp_get_num_threads();
myid = omp_get_thread_num();

OpenMPの基本関数
最大スレッド数取得（Integer）[C][F] nthreads = omp_get_num_threads()
自スレッド番号取得（Integer）[C][F] myid = omp_get_thread_num()!$ use omp_lib
integer :: myid, nthreads#include <omp.h>
int myid, nthreads;

nthreads = omp_get_num_threads()
myid = omp_get_thread_num()

nthreads = omp_get_num_threads();
myid = omp_get_thread_num();

OpenMPの基本関数
時間を測る（倍精度型）[F][C] time = omp_get_wtime()!$ use omp_lib
real(8) :: dts, dte
dts = omp_get_wtime()
・・・ 処理 ・・・
dte = omp_get_wtime()
print *, dte-dts#include <omp.h>
double dts;
double dte;
dts = omp_get_wtime();・・・ 処理 ・・・
dte = omp_get_wtime();

なお、OpenMPモジュール（ヘッダ）のロードを忘れると、これらの関数を使用できずコンパイルエラーになるWorking Sharing構文〇 複数のスレッドで分担して実行する部分を指定〇 並列リージョン内で記述する#pragma omp parallel { } の括弧範囲内

指示文の基本形式は[C] #pragma omp xxx[F] !$omp xxx ~ !$omp end xxx◎for構文, do構文
ループを分割し各スレッドで実行◎section構文
各セクションを各スレッドで実行◎single構文
１スレッドのみ実行◎master構文
マスタースレッドのみ実行並列リージョンを指定#pragma omp parallel{#pragma omp forfor (i=0; i<100; i++) {a[i] = i}#pragma omp single{...
}#pragma omp forfor (...)
.....
}
}

スレッドの起動～終結[C] #pragma omp parallel { }括弧 { } 内が複数スレッドで処理される。
複数スレッドで処理（並列リージョン）

for構文#pragma omp parallel{#pragma omp forfor (i=0; i<100; i++) {a[i] = i}forループをスレッドで分割し、並列処理を行う#pragma omp single{
output(a);}・forループの前に指示行 #pragma omp for を入れる

#pragma omp forfor (i=0; i<100; i++) {b[i] = i}#pragma omp for を入れることでループ自体が分割され、各スレッドに処理が割り当てられる}[C] #pragma omp for#pragma omp parallel でスレッドを生成しただけでは、全てのスレッドが全ループを計算してしまう１スレッドのみで処理#pragma omp parallel{#pragma omp forfor (i=0; i<100; i++) {a[i] = i;}#pragma omp single{
output(a);}#pragma omp forfor (i=0; i<100; i++) {b[i] = i;}}[C] #pragma omp single { }逐次処理やデータの出力のような処理が入る場合、全スレッドで行う必要はなく、1スレッドのみで処理を行えばよい#pragma omp single { }を用いることで、{ } 内の記述は1スレッドのみで処理される

並列リージョンを指定!$omp parallel!$omp dodo i = 1, 100a(i) = i
enddo!$omp end do!$omp single
call output(a)!$omp end single!$omp dodo i = 1, 100b(i) = i
enddo!$omp end do!$omp end parallel

スレッドの起動[F] !$omp parallel
スレッドの終結[F] !$omp end parallel
複数スレッドで処理（並列リージョン）

do構文!$omp parallel!$omp dodo i = 1, 100a(i) = i
enddo!$omp end do

doループをスレッドで分割し、並列処理を行う!$omp single
call output(a)!$omp end single

・do の直前に指示行 !$omp do を入れる!$omp dodo i = 1, 100b(i) = i
enddo!$omp end do!$omp end parallel!$omp parallel でスレッドを生成しただけでは、全てのスレッドが全ループを計算してしまう[F] !$omp do ～ !$omp end do

・enddo の直後に指示行 !$omp end do を入れる!$omp do を入れることでループ自体が分割され、各スレッドに処理が割り当てられる

１スレッドのみで処理!$omp parallel!$omp dodo i = 1, 100a(i) = i
enddo!$omp end do!$omp single
call output(a)!$omp end single!$omp dodo i = 1, 100b(i) = i
enddo!$omp end do!$omp end parallel[F] !$omp single ~!$omp end single
逐次処理やデータの出力のような処理が入る場合、全スレッドで行う必要はなく、1スレッドのみで処理を行えばよい!$omp singleを用いることで、{ } 内の記述は1スレッドのみで処理される

OpenMPによるスレッド並列#pragma omp parallel{for (i=0; i<100; i++) {a[i] = i;}}スレッドを生成しただけでは、全スレッドが全ての処理を行ってしまい負荷分散にならない

スレッド０

スレッド1

スレッド2

スレッド3

for (i=0; i<100; i++)for (i=0; i<100; i++)for (i=0; i<100; i++)for (i=0; i<100; i++)

OpenMPによるスレッド並列#pragma omp parallel{#pragma omp forfor (i=0; i<100; i++) {a[i] = i;}}

ワークシェアリング構文を入れることにより、処理が分割され、正しく並列処理される。
#pragma omp for、 !$omp do はループを自動的にスレッド数で均等に分割する

スレッド０

スレッド1

スレッド2

スレッド3

for (i=0; i<25; i++)for (i=25; i<50; i++)for (i=50; i<75; i++)for (i=75; i<100; i++)

OpenMPの基本命令
スレッド生成とループ並列を１行で記述[C言語]#pragma omp parallel { }#pragma omp for→#pragma omp parallel for と書ける
#pragma omp parallel{#pragma omp forfor (i=0; i<100; i++) {a[i] = i;}}#pragma omp parallel forfor (i=0; i<100; i++) {a[i] = i;}OpenMPの基本命令
スレッド生成とループ並列を１行で記述[Fortran]!$omp parallel!$omp do→!$omp parallel do と書ける!$omp parallel!$omp dodo i = 1, 100a(i) = i
enddo!$omp end do!$omp end parallel!$omp parallel dodo i = 1, 100a(i) = i
enddo!$omp end parallel doプライベート変数について・OpenMPにおいて変数は基本的には共有（shared）であり、どのスレッドからもアクセス可能である。プライベート変数に指定した変数は各スレッドごとに値を保有し、他のスレッドからアクセスされない。
・並列化したループ演算の内部にある一時変数などは、プライベート変数に指定する必要がある。
・例外的に[C]#pragma omp for[F] !$omp parallel doの直後のループ変数はプライベート変数になる

プライベート変数についてプライベート変数を指定[C] #pragma omp parallel for private(a, b, ...)
[C] #pragma omp for private(a, b, ...)
#pragma omp parallel{#pragma omp for private(j, k)for (i=0; i<nx; i++) {for (j=0; j<ny; j++) {for (k=0; k<nz; k++) {f[i][j][k] = (double)(i * j * k);}}}}

ループ変数の扱いに関して並列化したループ変数は自動的にprivate変数になる。しかし多重ループの場合、内側のループに関しては共有変数のままである。

左の例の場合、i は自動的にprivateになるため必要ないが、j, k についてはprivate宣言が必要となる。

プライベート変数についてプライベート変数を指定[F] !$omp parallel private(a, b, ...)
[F] !$omp do private(a, b, ...)!$omp parallel!$omp do private(j, k)
do i = 1, nx
do j = 1, ny
do k = 1, nz
f(k, j, i) = dble(i * j * k)
enddo
enddo
enddo
!$omp end do!$omp end parallel

ループ変数の扱いに関して並列化したループ変数は自動的にprivate変数になる。しかし多重ループの場合、内側のループに関しては共有変数のままである。

左の例の場合、i は自動的にprivateになるため必要ないが、j, k についてはprivate宣言が必要となる。

プライベート変数について起こりがちなミス#pragma omp forfor (i=0; i<100; i++) {tmp = myfunc(i);a[i] = tmp;}#pragma omp for private(tmp)for (i=0; i<100; i++) {tmp = myfunc(i);a[i] = tmp;}

tmpを上書きしてしまい、正しい結果にならない

private宣言を入れる

並列化したループ内で値を設定・更新する場合は要注意→privateにすべきではないか確認する必要ありプライベート変数について
スレッド０
共有変数 tmp に 0 を代入

tmp = 0tmp = 25

共有変数 tmp は 25 を代入
a[0] には 25 が代入される

スレッド1

a[0] = tmpa[25] = tmp

private宣言なし#pragma omp forfor (i=0; i<100; i++) {tmp = myfunc(i);a[i] = tmp;}

処理順

プライベート変数について
スレッド０
スレッド0のプライベート変数 tmp に0 を代入

tmp = 0

スレッド1のプライベート変数 tmp に25 を代入

a[0] には 0 が代入される

スレッド1

tmp = 25a[0] = tmpa[25] = tmp

private宣言あり#pragma omp for private(tmp)for (i=0; i<100; i++) {tmp = myfunc(i);a[i] = tmp;}

処理順

多重ループに関して良くない例

改善案

for (i=0; i<nx; i++) {for (j=0; j<ny; j++) {#pragma omp parallel for private(i, j, k)for (k=0; k<nz; k++) {f[i][j][k] = (double)(i * j * k);}}}#pragma omp parallel{for (i=0; i<nx; i++) {for (j=0; j<ny; j++) {#pragma omp for private(i, j, k)for (k=0; k<nz; k++) {f[i][j][k] = (double)(i * j * k);}}}}

OpenMPを用いた並列化では、内側ループ、外側ループのどちらを並列化しても動作はするが、内側ループを並列化すると毎回スレッドの生成を行うため遅くなる。（上記の例では nx * ny 回のスレッド生成）

なお、並列化するループを変えたり、ループの計算順序を変更する可能性があるため、private宣言にはループ変数も書いた方が無難。

共有変数について共有（shared）変数を指定[C] #pragma omp parallel shared(a, b, ...)
[C] #pragma omp for shared(a, b, ...)
[F] !$omp parallel shared(a, b, ...)
[F] !$omp do shared(a, b, ...)
・指定しなければ基本的に共有変数であるため、省略可能。

スレッドの同期
nowait を明示しない限り、ワークシェアリング構文の終わりに自動的に同期処理が発生
スレッドの同期待ちをしない[C] #pragma omp for nowait[F] !$omp do ~ !$omp end do nowait
スレッドの同期をとる[C] #pragma omp barrier[F] !$omp barrier注意点等
・自動並列化と違い、並列化できるかどうかの判断はプログラマが行う・依存関係などにより並列化できないループであっても、明示してしまえば並列化されてしまう。

・スレッド内でのグローバル変数、プライベート変数を間違えるとRunごとに結果が変わってしまう。
→数回実行し、結果が変わらないことを確認・OpenMPは手軽だが、デバッグには注意が必要

演習問題2-1
OpenMPを用いて次のコードを完成させよ
Fortran： 2019spring/code/f90/openmp/1d_adv_omp.f90

C言語： 2019spring/code/c/openmp/1d_adv_omp.c*演習１にて使用したコードを用いても良い移流が出来たら流体にも挑戦！
1d_fluid_rk: １次元流体

2d_fluid: ２次元流体補足：マシン上の操作
コンパイル・ジョブの投入（C言語）> icc -qopenmp 1d_adv_omp.c> qsub run.sh#!/bin/bash#PBS -q S#PBS -l select=1:ncpus=8#PBS -N OpenMP#PBS -o output#PBS -j oesource /etc/profile.d/modules.shmodule load intel
export OMP_NUM_THREADS=8
export KMP_AFFINITY=disabled
cd ${PBS_O_WORKDIR}
dplace -x2 ./a.out

run.sh
使用ノード数、CPU数
ジョブ名
標準出力の出力先ファイル

Intelコンパイラ環境のロード
スレッド数を指定

実行補足：マシン上の操作
コンパイル・ジョブの投入（Fortran）> ifort -qopenmp 1d_adv_omp.f90> qsub run.sh#!/bin/bash#PBS -q S#PBS -l select=1:ncpus=8#PBS -N OpenMP#PBS -o output#PBS -j oesource /etc/profile.d/modules.shmodule load intel
export OMP_NUM_THREADS=8
export KMP_AFFINITY=disabled
cd ${PBS_O_WORKDIR}
dplace -x2 ./a.out

run.sh
使用ノード数、CPU数ジョブ名
標準出力の出力先ファイル

Intelコンパイラ環境のロード
スレッド数を指定

実行

演習問題2-1
解答例
Fortran：2019spring/code/f90/openmp/sample/1d_adv_omp_sample.f90
2019spring/code/f90/openmp/sample/1d_fluid_rk_omp_sample.f90
2019spring/code/f90/openmp/sample/2d_fluid_omp_sample.f90

C言語：
2019spring/code/c/openmp/sample/1d_adv_omp_sample.c
2019spring/code/c/openmp/sample/1d_fluid_rk_omp_sample.c

2019spring/code/c/openmp/sample/2d_fluid_omp_sample.c

演習問題2-2
OpenMPを用いて並列化したプログラムを、スレッド数を変えて実行し処理時間を計測せよ
スレッド数の指定方法：スクリプト run.sh の
export OMP_NUM_THREADS=xx の数値を変更
処理時間の計測
dts = omp_get_wtime()
・・・ 処理 ・・・
dte = omp_get_wtime()
print *, dte-dts