KOBE HPC サマースクール（初級）2023 講義4

Sep. 18-22, 2023 OpenMPによるスレッド並列計算 今村 俊幸 （理化学研究所 計算科学研究センター） imamura.toshiyuki@riken.jp 1 KOBE HPC サマースクール（初級） 2023年9月18-22日

Sep. 18-22, 2023 講義内容 1. スレッド並列とは 2. OpenMPによるループ処理の並列化 3. 差分化された偏微分方程式の並列化 2 4. アムダールの法則と並列化効率の評価

Sep. 18-22, 2023 計算機とムーアの法則 • スカラー計算機 • 単一の実行ユニットが命令を１つ１つ、逐次処理を行う計算機。ただし近 年では、複数のスカラー命令を同時に処理できるスー パースカラーが主流。 一般的なパソコンやスマートフォンなど は、このスカラー計算機に分類さ れる。 • Intel系、AMD系、ARM系 (FUJITSU の富岳、FXシリーズ）など • ベクトル計算機 3 • ベクトル演算（同じ計算式を各配列要素に対してそれぞれ計算するといっ た処理）を一括して行うことができる計算機。一時期は『スーパーコン ピュータ ≒ ベクトル機』と言われていたこともあった。 • NEC-SX（地球シミュレータ[初代]）、GPUアクセラレータ等

Sep. 18-22, 2023 計算機とムーアの法則 • 半導体の集積密度は１年半～２年ごとに倍になる • エンドユーザーからすれば、ただ待っているだけで計算機の性能 が向上し続け、大規模計算が可能になるという”ありがたい話” しかし 近年、ムーアの法則の限界が 囁かれるようになる →そもそも原子サイズより小 さな回路は実現不可能 →リーク電流、放熱の問題等 4 https://ja.wikipedia.org/wiki/%E3%83%A0%E3%83%BC%E3%82%A2%E3%8 1%AE%E6%B3%95%E5%89%87

Sep. 18-22, 2023 計算機とムーアの法則 大型計算機の主流はベクトルからスカラーへと変化してきた。また、 そのスカラー計算機も、製造プロセスの微細化によるクロックの高 速 化が頭打ちになっており、マルチコア化が図られるなど方針の転 換を 図られてきた。近年ではGPU（アクセラレータ）を計算利用するなど、 ハードウェアも多様化している。 →計算機の特性にあったプログラミングが必要 共有メモリ型 ・複数のコアを持つCPU（マルチコア）、SGI UVシリーズ ⇒ スレッド並列計算（OpenMP、自動並列など） 5 分散メモリ型 ・複数の計算機を通信システムで繋ぐ（クラスター） ⇒ プロセス並列計算（MPI、XcalableMPなど） アクセラレータ ・メインのCPUに加え、GPUやXeon Phiなどを計算に用いる ⇒ OpenACC、OpenMP 4.0、OpenCL、CUDAなど

Sep. 18-22, 2023 並列計算機に関して • 共有メモリ型並列計算機 - 各演算ユニット間でメモリ空間が共有されている計算機 - 一般的なマルチコアCPU搭載のPCや、スパコンのノード １つ１つも、共有メモ リ型の並列計算機と言える • - OpenMPを用いたスレッド並列計算が可能 プロセス スレッド スレッド スレッド スレッド Private Private Private Private Global 6 メモリ空間 各スレッド間で メモリ空間は共有

Sep. 18-22, 2023 並列計算機に関して • 分散メモリ型並列計算機 - 独立したメモリを持つ計算ノード間を通信しながら並列動作させる計算機 - 富岳に代表される近代的なスーパーコンピュータの主流 7 - MPIなどを用いたプロセス並列計算が必要となる ノード ノード プロセス プロセス メモリ メモリ ノード ノード プロセス間の データ通信が必要！ プロセス プロセス メモリ メモリ

Sep. 18-22, 2023 並列計算機のイメージ 非並列処理 スレッド並列 （共有メモリ） ・各CPU（スレッド）は全てのデータにアクセス可能 ・故にお互いが作業の邪魔してしまうこともある CPU0 MPI通信 CPU1 8 プロセス並列 （分散メモリ） ・各CPU（プロセス）は自身のメモリ空間上のデータのみにアクセス可能 ・他のプロセスのデータはMPI通信でアクセス

Sep. 18-22, 2023 OpenMPとは • Open Multi-Processing の略 ・ • 共有メモリ型計算機用の並列計算API（仕様） →ノード内のスレッド並列（ノード間は不可） • ユーザーが明示的に並列のための指示を与える →コンパイラの自動並列とは異なる • 標準化された規格であり、広く使われている • 指示行の挿入を行うことで並列化できる →既存の非並列プログラムに対し、元のコードの構造を大きく変える ことなく並列化できるため、比較的手軽 9 • ちなみに “OpenMPI” というライブラリが存在するが、こちらはMPI （Message Passing Interface）の実装の１つであり、OpenMPとは全 くの別物

Sep. 18-22, 2023 OpenMPによるスレッド並列 Fork-Joinモデル スレッド０ Fork スレッド０ スレッド1 スレッド2 スレッド3 並列リージョン Join スレッド０ 10 ... 非並列 !$omp parallel ... 並列リージョン !$omp end parallel ... 非並列 Fork Join ... 非並列 #pragma omp parallel { ... 並列リージョン } Join ... 非並列処理 Fork

Sep. 18-22, 2023

OpenMPの基本関数
• OpenMPモジュール/ヘッダをロード

[C] #include <omp.h>
[F] use omp_lib

11

*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();

Sep. 18-22, 2023

12

OpenMPの基本関数
・最大スレッド数取得
[C][F] nthreads = omp_get_max_threads()
・並列リージョン内のスレッド数取得
[C][F] nthreads = omp_get_num_threads()
・自スレッド番号取得
[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();

Sep. 18-22, 2023

OpenMPの基本関数
• 時間を測る（倍精度型）

• ある時点からの経過時間を取得できるので、測定したい区間の前後に以下の関数を挟
み込んで差を求めます。

[F][C] time = omp_get_wtime()

13

!$ 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モジュール（ヘッダ）のロードを忘れると、
これらの関数を使用できずコンパイルエラーになる

Sep. 18-22, 2023

並列リージョンを指定 (C言語)
#pragma omp parallel
{
#pragma omp for
for (i=0; i<100; i++) {
a[i] = i
}
#pragma omp single
{
...
}

14

#pragma omp for
for (...)
.....
}
}

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

Sep. 18-22, 2023 並列リージョンを指定 (Fortran) !$omp parallel !$omp do do i = 1, 100 a(i) = i enddo !$omp end do 15 !$omp single call output(a) !$omp end single !$omp do do i = 1, 100 ..... enddo !$omp end do !$omp end parallel スレッドの起動 [F] !$omp parallel スレッドの終結 [F] !$omp end parallel 複数スレッドで処理 （並列リージョン）

Sep. 18-22, 2023 演習準備 演習用ファイル（OpenMP） rokko:/home/guest59/openmp 上記のファイルをホームディレクトリにコピーする 16 mkdir ~/openmp cd ~/openmp cp -r /home/guest59/openmp/ ./

Sep. 18-22, 2023

演習1
演習1-1：omp_ex01_1.c
スレッドを４つ生成し、それぞれのスレッドで "Hello, World!" を出力せよ

演習1-2：omp_ex01_2.c
スレッド４つで実行、それぞれ自スレッド番号を取得し出力せよ

omp_ex01_1.c

#include <stdio.h>

17

int main(void){
printf("Hello, World!");
}

#include <stdio.h>
int main(void){
#pragma omp parallel
{
printf("Hello, World!");
}
}

Sep. 18-22, 2023 演習1 演習1-1：omp_ex01_1.f90 スレッドを４つ生成し、それぞれのスレッドで "Hello, World!" を出力せよ 演習1-2：omp_ex01_2.f90 スレッド４つで実行、それぞれ自スレッド番号を取得し出力せよ omp_ex01_1.f90 program omp_ex01_1 print*, "Hello, World!" 18 end program omp_ex01_1 !$omp parallel print*, "Hello, World!" !$omp end parallel end

Sep. 18-22, 2023 環境の構築 • 本実験はインテルコンパイラを使用します (GNUコンパイラでも可能ですが、今回はインテルコンパイラで統一) • 適切なモジュールのロード • module load intel コマンド成功時には何もメッセージが返ってこないので、以下で確認 • module list Currently Loaded Modulefiles: 1) intel/19.1.3 19 上記のようにモジュールが読み込まれていればOK. 問題があれば、講師やRAに相談し てください。

Sep. 18-22, 2023

コンパイル
・コンパイルオプションでOpenMPを有効にする
icc -qopenmp -o omp_ex01_1 omp_ex01_1.c
gcc -fopenmp -o omp_ex01_1 omp_ex01_1.c
コンパイルオプションを指定しない場合はOpenMPの指示行はコメントとして認
識される。

20

#pragma omp parallel for
{
for (i=0; i<100; i++) {
a[i] = b[i] + c;
}

指示行はCの場合は #pragma omp ...とい
う形式で記述する。
オプションを付けない場合、指示行は
無視される

Sep. 18-22, 2023 コンパイル ・コンパイルオプションでOpenMPを有効にする ifort -qopenmp -o omp_ex01_1 omp_ex01_1.f90 gfortran -fopenmp -o omp_ex01_1 omp_ex01_1.f90 コンパイルオプションを指定しない場合はOpenMPの指示行はコメントとして認 識される。 21 !$omp parallel do do i = 1, 100 a(i) = b(i) + c enddo !$omp end parallel do OpenMPで用いる指示行は、Fortranの 場合 !$OMP から始まる。行頭に! があ る行は通常、コメントとして処理され る。

Sep. 18-22, 2023

スレッド数の指定
・シェルの環境変数で与える（推奨）
export OMP_NUM_THREADS=4 （bashの場合)
setenv OMP_NUM_THREADS 4 （tcshの場合）
・プログラム内部で設定することも可能
!$ use omp_lib
call omp_set_num_threads(4)

#include <omp.h>
omp_set_num_threads(4);

22

ただしスレッド数を変えて実行する時など、毎回コンパイルが必要となっ
てしまうため、今回は環境変数による指定を推奨する。

Sep. 18-22, 2023

操作補足

ジョブスクリプト
run.sh

ジョブ投入方法はマシンの環境によって変わるため、実際に
スパコンを使う時にはユーザーマニュアルを参考にすること

#!/bin/bash
#PBS –q S
#PBS -l select=1:ncpus=4
#PBS -N omp_JOB
#PBS -j oe

キューを指定
リソース確保(4 コア)
ジョブ名
出力ファイル
標準エラー出力と結合

source /etc/profile.d/modules.sh
moduleコマンドのための環境設定
module load intel
Intelコンパイラ環境の読み込み

23

export KMP_AFFINITY=disabled
export OMP_NUM_THREADS=4
cd ${PBS_O_WORKDIR}
./a.out
実行

AFFINITYをdisabledにする
スレッド並列数の設定
作業ディレクトリへ移動

Sep. 18-22, 2023 Working Sharing構文 〇複数のスレッドで分担して実行する部分を指定 〇 並 列 リージョン内で記述する #pragma omp parallel { } の括弧範囲内 24 指示文の基本形式は [C] #pragma omp xxx [F] !$omp xxx ~ !$omp end xxx ◎for構文, do構文 ループを分割し各スレッドで実行 ◎section構文 各セクションを各スレッドで実行 ◎single構文 １スレッドのみ実行 ◎master構文 マスタースレッドのみ実行

Sep. 18-22, 2023

for構文 (C言語)
#pragma omp parallel
{
#pragma omp for
for (i=0; i<100; i++) {
a[i] = i
}
#pragma omp for
for (i=0; i<100; i++) {
b[i] = i
}

25

}

forループをスレッドで分割し、
並列処理を行う
[F] #pragma omp for
・forループの前に指示行 #pragma omp for を入れる
#pragma omp parallel でスレッドを生成しただけで
は、全てのスレッドが全ループを計算してしまう
#pragma omp for を入れることでループ自体が分割
され、各スレッドに処理が割り当てられる

Sep. 18-22, 2023 do構文 (Fortran) !$omp parallel !$omp do do i = 1, 100 a(i) = i enddo !$omp end do !$omp do do i = 1, 100 b(i) = i enddo !$omp end do !$omp end parallel doループをスレッドで分割 し、並列処理を行う [F] !$omp do ～ !$omp end do ・do の直前に指示行 !$omp do を入れる ・enddo の直後に指示行 !$omp end do を入れる !$omp parallel でスレッドを生成しただけでは、全 てのスレッドが全ループを計算してしまう 26 !$omp do を入れることでループ自体が分割され、 各スレッドに処理が割り当てられる

Sep. 18-22, 2023

OpenMPによるスレッド並列

27

#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++)

Sep. 18-22, 2023

OpenMPによるスレッド並列

28

#pragma omp parallel
{
#pragma omp for
for (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++)

Sep. 18-22, 2023

OpenMPの基本命令(C言語)
スレッド生成とループ並列を１行で記述

29

[C言語]
#pragma omp parallel { }
#pragma omp for
→#pragma omp parallel for と書ける
#pragma omp parallel
{
#pragma omp for
for (i=0; i<100; i++) {
a[i] = i;
}
}

#pragma omp parallel for
for (i=0; i<100; i++) {
a[i] = i;
}

Sep. 18-22, 2023 OpenMPの基本命令(Fortran) スレッド生成とループ並列を１行で記述 30 [Fortran] !$omp parallel !$omp do →!$omp parallel do と書ける !$omp parallel !$omp do do i = 1, 100 a(i) = i enddo !$omp end do !$omp end parallel !$omp parallel do do i = 1, 100 a(i) = i enddo !$omp end parallel do

Sep. 18-22, 2023

演習2
演習2-1：omp_ex02_1.c
サンプルのプログラムはループがスレッドで分割されていない。
指示文を挿入（もしくは修正）し、ループを正しく並列化せよ。
まずは omp_ex02_1.c をそのまま動かして挙動を確認しよう
omp_ex02_1.c

31

#include <stdio.h>
...
#pragma omp parallel
{
for (i=0; i<10; i++) {
printf("myid=%d, i=%d", omp_get_thread_num(), i);
}
}
}

Sep. 18-22, 2023 演習2 演習2-1：omp_ex02_1.f90 サンプルのプログラムはループがスレッドで分割されていない。 指示文を挿入（もしくは修正）し、ループを正しく並列化せよ。 まずは omp_ex02.f90 をそのまま動かして挙動を確認しよう omp_ex02_1.f90 32 program omp_ex02_1 ... !$omp parallel do i=1, 10 print*, 'myid =', omp_get_thread_num(), 'i =', i enddo !$omp end parallel end

Sep. 18-22, 2023

演習2
演習2-2：omp_ex02_2.c
演習2-1と同様に指示行を挿入し、ループを並列化せよ。また、
計算結果（u）が、スレッド数を1,2,4と変えても変わらないこと
を確認せよ。
omp_ex02_2.c

33

#include <stdio.h>
...
#pragma omp parallel
{
for (i=0; i<100; i++) {
u[i] = sin(2.0*pi*(double)(i+1)/100.0);
//
printf("myid=%d, i=%d", omp_get_thread_num(), i);
}
}

Sep. 18-22, 2023 演習2 演習2-2：omp_ex02_2.f90 演習2-1と同様に指示行を挿入し、ループを並列化せよ。また、 計算結果（u）が、スレッド数を1,2,4と変えても変わらないこと を確認せよ。 omp_ex02_2.f90 34 program omp_ex02_2 ... !$omp parallel do i=1, 100 u(i) = sin(2.0*pi*dble(i)/100.0) ! print*, 'myid =', omp_get_thread_num(), 'i =', i enddo !$omp end parallel end

Sep. 18-22, 2023 プライベート変数について ・OpenMPにおいて変数は基本的には共有（shared）であり、 どのスレッドからもアクセス可能である。プライベート変数に 指定した変数は各スレッドごとに値を保有し、他のスレッドか らアクセスされない。 ・並列化したループ演算の内部にある一時変数などは、プライ ベート変数に指定する必要がある。 35 ・例外的に [C]#pragma omp for [F] !$omp parallel do の直後のループ変数はプライベート変数になる

Sep. 18-22, 2023

36

プライベート変数について
プライベート変数を指定
[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宣言が必要となる。

Sep. 18-22, 2023 37 プライベート変数について プライベート変数を指定 [F] !$omp parallel do 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宣言が必要となる。

Sep. 18-22, 2023

プライベート変数について
起こりがちなミス
#pragma omp for
for (i=0; i<100; i++) {
tmp = myfunc(i);
a[i] = tmp;
}
tmpを上書きしてしまい、
正しい結果にならない

#pragma omp for private(tmp)
for (i=0; i<100; i++) {
tmp = myfunc(i);
a[i] = tmp;
}
private宣言を入れる

38

並列化したループ内で値を設定・更新する場合は要注意
→privateにすべきではないか確認する必要あり

Sep. 18-22, 2023

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

tmp = 0

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

スレッド1

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

private宣言なし

39

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

処理順

Sep. 18-22, 2023

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

tmp = 0

スレッド１のプライベート変数
tmp に 25 を代入
a[0] には 0 が代入される

スレッド1

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

private宣言あり

40

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

処理順

Sep. 18-22, 2023 OpenMP外で定義されたプライベート変数 • firstprivate/lastprivate tmp=100 スレッド０ スレッド０のプライベート変数 tmp に 1加算 (=101) tmp = tmp+1 スレッド１のプライベート変数 tmp に 25 を加算 (値は不定) a[0] には 101 が代入される スレッド1 tmp = tmp+25 a[0] = tmp private変数はあくまでも 各スレッドのfork時に生成された変数であり、 初期値は不定。マスタースレッドのみ OpenMP外側で定義されたものと同等 a[25] = tmp 41 全てのスレッドでfork時にOpenMP外側での 値とする場合はfirstprivateの属性が必要 Join時、マスタースレッドの値が保持されるが、 シーケンシャルな動作と同等にしたい場合は lastprivate属性の指定が必要 処理順 tmp=???

Sep. 18-22, 2023

多重ループに関して
良くない例
for (i=0; i<nx; i++) {
for (j=0; j<ny; j++) {
#pragma omp parallel for private(k)
for (k=0; k<nz; k++) {
f[i][j][k] = (double)(i * j * k);
}
}
}

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

42

OpenMPを用いた並列化では、内側ループ、外側ループのどちらを並
列化しても良い。ただし、内側ループを並列化する場合、毎回forkjoinしてしまうとスレッド生成回数がものすごいことになる。
（上記の例では１つの３次元ループで nx * ny 回）
なお、並列化するループを変えたり、ループの計算順序を変更する
可能性があるため、private宣言にはループ変数も書いた方が無難。

Sep. 18-22, 2023 共有変数について 共有（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, ...) 43 ・指定しなければ基本的に共有変数であるため、 省略可能。

Sep. 18-22, 2023 スレッドの同期 nowait を明示しない限り、ワークシェアリング 構文の終わりに自動的に同期処理が発生 スレッドの同期待ちをしない [C] #pragma omp for nowait [F] !$omp do ~ !$omp end do nowait スレッドの同期をとる [C] #pragma omp barrier 44 [F] !$omp barrier

Sep. 18-22, 2023 Section構文 スレッドごとに処理を分岐させる [C] #pragma omp sections {#pragma omp section, ...} [F] !$omp sections, !$omp section, ... , !$omp end sections !$omp parallel !$omp sections !$omp section 処理A !$omp section 処理B #pragma omp parallel { #pragma omp sections { #pragma omp section 処理A #pragma omp section 処理B !$omp end sections !$omp end parallel 45 } }

Sep. 18-22, 2023 処理の分岐 Section構文 スレッドごとに処理を分岐させる #pragma omp parallel sections スレッド０ #pragma omp section 処理0 46 待ち時間 スレッド1 #pragma omp section 処理1 待ち時間 同期 スレッド2 #pragma omp section 処理2 スレッド3 #pragma omp section 処理3 待ち時間 ・各スレッドに割り当てられた処理の負荷が異なると、無駄な待ちが発生する ・ロードバランスに注意

Sep. 18-22, 2023

演習3
演習3-1：omp_ex03.c

47

演習2で行ったループ演算の並列化は、SECTION構文を用いて分割
することも可能である。未完成のサンプルコード（omp_ex03.c）を
完成させ、同じ処理を行っていることを確認せよ。なおスレッド数
は４にすること。
omp_ex02_2.c

omp_ex03.c

#include <stdio.h>
...
for (i=0; i<100; i++) {
u[i] = sin(2.0*pi*…);
}

#include <stdio.h>
for (i=0; i<25; i++) {
u[i] = sin(2.0*pi*…);
}
for (i=25; i<50; i++) {
u[i] = sin(2.0*pi*…);
}

Sep. 18-22, 2023 演習3 演習3-1：omp_ex03.f90 48 演習2で行ったループ演算の並列化は、SECTION構文を用いて分割 することも可能である。未完成のサンプルコード（omp_ex03.f90） を完成させ、同じ処理を行っていることを確認せよ。なおスレッド 数は４にすること。 omp_ex02.f90 omp_ex03.f90 program ex02 ... do i=1, 100 u(i) = sin(2.0*pi*...) enddo program ex03 ... do i=1,25 u(i) = sin(2.0*pi*...) enddo end do i=26,50 u(i) = sin(2.0*pi*...) enddo

Sep. 18-22, 2023

１スレッドのみで処理
#pragma omp parallel
{
#pragma omp for
for (i=0; i<100; i++) {
a[i] = i;
}
#pragma omp single
{
output(a);
}

49

#pragma omp for
for (i=0; i<100; i++) {
b[i] = i;
}
}

[C] #pragma omp single { }

一般的に、スレッドの立ち上げ回数は
極力減らすほうがオーバーヘッドが少
なくなるため良いとされる。
逐次処理やデータの出力のような処理
が入る場合
#pragma omp single { }
とすると１スレッドのみで処理される

Sep. 18-22, 2023 １スレッドのみで処理 !$omp parallel !$omp do do i = 1, 100 a(i) = i enddo !$omp end do 50 !$omp single call output(a) !$omp end single !$omp do do i = 1, 100 ..... enddo !$omp end do !$omp end parallel [F] !$omp single ~ !$omp end single 一般的に、スレッドの立ち上げ回数は 極力減らすほうがオーバーヘッドが少 なくなるため良いとされる。 逐次処理やデータの出力のような処理 が入る場合、 !$omp single とすると１スレッドのみで処理される

Sep. 18-22, 2023

１スレッドのみで処理
#pragma omp parallel
{
#pragma omp for
for (i=0; i<100; i++) {
a[i] = i;
}
#pragma omp master
{
output(a);
}

51

#pragma omp for
for (i=0; i<100; i++) {
b[i] = i;
}
}

[C] #pragma omp master { }

一般的に、スレッドの立ち上げ回数は
極力減らすほうがオーバーヘッドが少
なくなるため良いとされる。
逐次処理やデータの出力のような処理
が入る場合
#pragma omp master { }
とするとマスタースレッドのみで処理さ
れる。
終了時にスレッド間同期は行われない。

Sep. 18-22, 2023 １スレッドのみで処理 !$omp parallel !$omp do do i = 1, 100 a(i) = i enddo !$omp end do 52 !$omp master call output(a) !$omp end master !$omp do do i = 1, 100 ..... enddo !$omp end do !$omp end parallel [F] !$omp master ~ !$omp end master 一般的に、スレッドの立ち上げ回数は 極力減らすほうがオーバーヘッドが少 なくなるため良いとされる。 逐次処理やデータの出力のような処理 が入る場合、 !$omp master とするとマスタースレッドのみで処理さ れる 終了時にスレッド間同期は行われない。

Sep. 18-22, 2023

排他処理
#pragma omp parallel
{
#pragma omp for
for (i=0; i<100; i++) {
a[i] = i;
}
#pragma omp critical
{
a[0]=func(a,b,100);
}

53

#pragma omp for
for (i=0; i<100; i++) {
b[i] = i;
}
}

[C] #pragma omp critical { }

共有メモリではshared変数の書き込
みと読み込みのタイミングによる意図
しない動作がおこります(競合状態)。
ブロック内の処理を行うスレッドが１
つだけになるよう排他処理する場合
#pragma omp critical { }
とすると排他処理される

Sep. 18-22, 2023 排他処理 !$omp parallel !$omp do do i = 1, 100 a(i) = i enddo !$omp end do 54 !$omp critical a(1)=func(a,b,100) !$omp end critical !$omp do do i = 1, 100 ..... enddo !$omp end do !$omp end parallel [F] !$omp critical ~ !$omp end critical 共有メモリではshared変数の書き込 みと読み込みのタイミングによる意図 しない動作がおこります(競合状態)。 ブロック内の処理を行うスレッドが１ つだけになるよう排他処理する場合 !$omp critical とすると排他処理される

Sep. 18-22, 2023

リダクション変数を指定
reduction(演算子:変数）

・並列計算時はそれぞれのスレッドで別々の値を
持ち、並列リージョン終了時に各スレッドの値が
足し合わされる（総和）shared変数
・総和の他、積などを求めることも可能

55

integer :: i, sum
sum = 0
!$omp parallel do reduction(+:sum)
do i=1, 10000
sum = sum + 1
enddo
!$omp end parallel do

int i, sum;
sum = 0;
#pragma omp parallel for reduction(+:sum)
for (i=0; i<10000; i++) {
sum = sum + 1;
}

Sep. 18-22, 2023

演習4
演習4-1：omp_ex04.f90 / omp_ex04.c
1から100までを足し合わせるプログラムをOpenMPで並列化せよ。

56

reduction変数の指定を忘れると正しく動かないため注意。
omp_ex04.f90

omp_ex04.c

program ex04
...

#include <stdio.h>
...

do i=1,100
a=a+i
enddo

for (i=1; i<=100; i++) {
a = a + i;
}

print*, a

printf("%d¥n", a);

Sep. 18-22, 2023

並列化できないプログラム
・どのループを並列化するか（可能か）
という判断は全てプログラマが行う
並列化できないプログラム

57

a[0] = 0
for (i=1; i<=100; i++) {
a[i] = a[i-1] + 1
}

このプログラムはi 番目の計算
を行うためにはi-1 番目の計算
結果が必要であり（データの依
存関係）、並列化できないプロ
グラムである。
しかしOpenMPの指示行を入れ
れば、コンピュータは無理やり
このプログラムを並列化し、間
違った計算を行う。

Sep. 18-22, 2023 演習5 演習5-1： 熱伝導問題のプログラム（laplace.c / laplace.f90） を OpenMPでスレッド並列化して計算せよ（ただし 並列化可能かも調べよ） 基礎方程式 𝜕𝜕2𝑇𝑇 𝜕𝜕2𝑇𝑇 + 2 =0 2 𝜕𝜕𝑥𝑥 𝜕𝜕𝑦𝑦 58 中心差分（2次精度） 𝑇𝑇𝑖𝑖,𝑗𝑗 1 = (𝑇𝑇𝑖𝑖−1,𝑗𝑗 + 𝑇𝑇𝑖𝑖+1,𝑗𝑗 + 𝑇𝑇𝑖𝑖,𝑗𝑗−1 + 𝑇𝑇𝑖𝑖,𝑗𝑗+1) 4

Sep. 18-22, 2023 演習5 演習5-1： 熱伝導問題のプログラム（laplace.c / laplace.f90） を OpenMPでスレッド並列化して計算せよ（ただし 正しく並列化されているかも調べよ） 59 SOR反復法 1 𝑘𝑘+1 𝑘𝑘+1 𝑘𝑘 𝑘𝑘 + 𝑇𝑇𝑖𝑖,𝑗𝑗−1 𝜌𝜌𝑖𝑖,𝑗𝑗 = 𝑇𝑇𝑖𝑖−1,𝑗𝑗 + 𝑇𝑇𝑖𝑖+1,𝑗𝑗 + 𝑇𝑇𝑖𝑖,𝑗𝑗+1 4 𝑘𝑘+1 = 𝑇𝑇 𝑘𝑘 𝑘𝑘 𝑇𝑇𝑖𝑖,𝑗𝑗 + 𝜔𝜔 𝜌𝜌 − 𝑇𝑇 𝑖𝑖,𝑗𝑗 𝑖𝑖,𝑗𝑗 𝑖𝑖,𝑗𝑗 境界条件 1 𝑇𝑇𝑖𝑖,0 = 𝑇𝑇𝑁𝑁+1,𝑗𝑗 = 𝑇𝑇𝑖𝑖,𝑁𝑁+1 = 0 𝑇𝑇0,𝑗𝑗 = sin(𝜋𝜋 ∗ 𝑖𝑖/𝑛𝑛) 0 1

Sep. 18-22, 2023 SOR法の補足 ・基本は、ガウス・ザイデル法（𝜔𝜔 = 1 のとき同じ式になる） 1 𝑘𝑘+1 𝑘𝑘+1 𝑘𝑘 𝑘𝑘 + 𝑇𝑇𝑖𝑖,𝑗𝑗−1 + 𝑇𝑇𝑖𝑖+1,𝑗𝑗 + 𝑇𝑇𝑖𝑖,𝑗𝑗+1 𝜌𝜌𝑖𝑖,𝑗𝑗 = 𝑇𝑇𝑖𝑖−1,𝑗𝑗 4 𝑘𝑘+1 = 𝑇𝑇 𝑘𝑘 𝑘𝑘 𝑇𝑇𝑖𝑖,𝑗𝑗 𝑖𝑖,𝑗𝑗 + 𝜔𝜔 𝜌𝜌𝑖𝑖,𝑗𝑗 − 𝑇𝑇𝑖𝑖,𝑗𝑗 ・注意せずに単純にOpenMPでパラレルforなどの指定をすると「競合状態」 に陥り正しく計算ができないことや、結果が実行毎に異なる場合があります。 ・一般に、オーダリングを用いて並列化（ハイパープレーン、２色など） ・２色オーダリングの場合： ブロックを先に計算 1 𝑘𝑘 𝑘𝑘 𝑘𝑘 𝑘𝑘 + 𝑇𝑇 + 𝑇𝑇 𝑇𝑇 + 𝑇𝑇 𝜌𝜌𝑖𝑖,𝑗𝑗 = 𝑖𝑖,𝑗𝑗+1 𝑖𝑖−1,𝑗𝑗 𝑖𝑖+1,𝑗𝑗 𝑖𝑖,𝑗𝑗−1 4 60 １． ２．次に、 ブロックを計算 1 𝑘𝑘+1 𝑘𝑘+1 + 𝑇𝑇𝑘𝑘+1 𝑘𝑘+1 𝑇𝑇 + 𝑇𝑇 + 𝑇𝑇 𝜌𝜌𝑖𝑖,𝑗𝑗 = 𝑖𝑖,𝑗𝑗+1 𝑖𝑖+1,𝑗𝑗 𝑖𝑖,𝑗𝑗−1 4 𝑖𝑖−1,𝑗𝑗

Sep. 18-22, 2023 演習5 演習5-1： 熱伝導問題のプログラム（ laplace.c / laplace.f90） をOpenMPでスレッド並列化して計算せよ 補足（laplace.c / laplace.f90) n：縦、横それぞれのグリッド数 ITMAX : 最大反復回数 eps：反復ベクトルの差のノルムの閾値 61 （どこまで値を収束させるか）

Sep. 18-22, 2023 演習5 gnuplotを使って結果をプロット username@rokko:~/khpc2019/omp> gnuplot gnuplot> set pm3d gnuplot> set ticslevel 0 gnuplot> set cbrange[0:1] カラーバーの値の範囲指定 gnuplot> set palette defined (0 “blue”, 1 "red") gnuplot> splot “data.d" with pm3d カラーバーの値の範囲等は 計算のパラメータによって 62 適宜変更すること カラーバーの色設定

Sep. 18-22, 2023 演習5 演習5-2： スレッド並列化した熱伝導問題のプログラムを、並列数を 1, 2, 4, 8, 16, … と変えて実行し、計算時間を計測せよ 並列数はジョブスクリプトにて（8スレッドの場合） CPU確保： #PBS -l select=1:ncpus=8 スレッド数： export OMP_NUM_THREADS=8 63 の２か所を書き換える

Sep. 18-22, 2023 演習5 演習5-2： 時間計測（簡易版）： ジョブスクリプトの実行で time を使用 time dplace ./a.out 時間計測（関数）： omp_get_wtime 関数を使用 計算開始時に dts = omp_get_wtime() 計算終了時に dte = omp_get_wtime() 64 とし、 dte - dts を出力する

Sep. 18-22, 2023 並列数と計算性能 スピードアップ 45 40 35 30 25 20 15 10 5 0 1 2 4 8 16 32 40 No. of threads 理論値 スピードアップ 65 並列数 1 2 4 8 16 32 40 実測値(秒) 2.56E+00 1.30E+00 6.61E-01 3.42E-01 1.77E-01 1.23E-01 1.16E-01 スピード アップ値 1.00E+00 1.97E+00 3.87E+00 7.49E+00 1.45E+01 2.08E+01 2.21E+01

Sep. 18-22, 2023 スケジューリング #pragma omp parallel for スレッド０ for (i=0; i<50; i++) スレッド1 for (i=50; i<100; i++) 待ち時間 66 ・各スレッドに割り当てられた処理の負荷が異なると、無駄な待ちが発生する ・ロードバランスに注意

Sep. 18-22, 2023 スケジューリング #pragma omp parallel for schedule(static, 25) スレッド０ for (i=0; i<25; i++) スレッド０ スレッド1 for (i=25; i<50; i++) for (i=50; i<75; i++) 待ち時間 スレッド1 67 for (i=75; i<100; i++) ・各スレッドに割り当てられた処理の負荷が異なると、無駄な待ちが発生する ・ロードバランスに注意

Sep. 18-22, 2023 スケジューリング #pragma omp parallel for schedule(dynamic, 25) スレッド０ for (i=0; i<25; i++) スレッド０ for (i=50; i<75; i++) スレッド０ 68 for (i=75; i<100; i++) スレッド1 for (i=25; i<50; i++) 待ち時間 ・各スレッドに割り当てられた処理の負荷が異なると、無駄な待ちが発生する ・ロードバランスに注意

Sep. 18-22, 2023

スケジューリング
• 例：（上三角行列）＊（ベクトル）
• 𝑦𝑦 = 𝐴𝐴 ∗ 𝑥𝑥, (𝑎𝑎𝑖𝑖𝑗𝑗 = 0, 𝑖𝑖 > 𝑗𝑗)
for (i=0; i<n; i++) {
for (j=i; j<n; j++) {
y[i] += a[i][j]*x[j];
}}

69

• チャンクサイズが性能に大きく影響する。
• 負荷バランスとシステムのオーバーヘッドの
トレードオフ。
• チューニングコストが増加する。

Sep. 18-22, 2023 並列化率に関して CPUをN個使って並列計算した時、計算速 度がN倍になるのが理想だが・・・ 並列化率の問題 プログラム内に並列化できない処理が含ま れていると、その部分が並列計算における ボトルネックになる（アムダールの法則） 70 並列計算を行うためのコスト 並列化を行うことで、逐次実行には不要 だ った処理が増えることがある（スレッド起 動、MPIのプロセス間通信等）

Sep. 18-22, 2023 並列化率に関して 依存関係により並列化 できない処理 逐次処理 ① ② ③ ④ ⑤ ⑥ 並列化可能な処理 ② ２並列 ④ ① ⑥ ③ ⑤ 2スレッド並列 計算時間は2/3（計算速 度は1.5倍） ② 71 ４並列 ① ④ ③ ⑤ 例：並列化できる部分が プログラム全体のうち 2/3の場合 ⑥ 4スレッド並列 計算時間は半分（計算速 度は2倍）

Sep. 18-22, 2023 並列化率に関して アムダールの法則 プログラムの並列化できる割合をP とし、プ ロセッサ数をn とすると、並列計算した時の 性能向上率は 1 で与えられる。 𝑃𝑃 1 − 𝑃𝑃 + 𝑛𝑛 72 これをアムダールの法則と呼ぶ。

Sep. 18-22, 2023 並列化率に関して • アムダールの法則 • 例えば、プログラム全体の９割は並列化できるが１割は逐次処理が残っ てしまうような場合、どれだけプロセッサを投入しても計算速度は１０倍 以上にはならない。 73 • 富岳のような大型計算機を用いる上では、如何にして並列化率を上げるかが 重要である。