Sep. 18-22, 2023

OpenMPによるスレッド並列計算
今村 俊幸（理化学研究所 計算科学研究センター）
imamura.toshiyuki@riken.jp

1

KOBE HPC サマースクール（初級）2023年9月18-22日

Sep. 18-22, 2023

講義内容
1. スレッド並列とは
2. OpenMPによるループ処理の並列化
3. 差分化された偏微分方程式の並列化

24. アムダールの法則と並列化効率の評価

Sep. 18-22, 2023

計算機とムーアの法則• スカラー計算機• 単一の実行ユニットが命令を１つ１つ、逐次処理を行う計算機。ただし近年では、複数のスカラー命令を同時に処理できるスー パースカラーが主流。
一般的なパソコンやスマートフォンなど は、このスカラー計算機に分類される。
• Intel系、AMD系、ARM系 (FUJITSU の富岳、FXシリーズ）など• ベクトル計算機

3• ベクトル演算（同じ計算式を各配列要素に対してそれぞれ計算するといった処理）を一括して行うことができる計算機。一時期は『スーパーコン
ピュータ ≒ ベクトル機』と言われていたこともあった。
• NEC-SX（地球シミュレータ[初代]）、GPUアクセラレータ等

Sep. 18-22, 2023

計算機とムーアの法則• 半導体の集積密度は１年半～２年ごとに倍になる• エンドユーザーからすれば、ただ待っているだけで計算機の性能が向上し続け、大規模計算が可能になるという”ありがたい話”
しかし
近年、ムーアの法則の限界が囁かれるようになる→そもそも原子サイズより小さな回路は実現不可能→リーク電流、放熱の問題等

4https://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を用いたスレッド並列計算が可能

プロセス
スレッド

スレッド

スレッド

スレッド

PrivatePrivatePrivatePrivateGlobal6

メモリ空間

各スレッド間でメモリ空間は共有

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

... 非並列Fork!$omp parallel... 並列リージョン!$omp end parallel... 非並列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 forfor (i=0; i<100; i++) {a[i] = i}#pragma omp single
{
...
}

14

#pragma omp forfor (...)
.....
}
}

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

Sep. 18-22, 2023

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

15!$omp single
call output(a)!$omp end single!$omp dodo 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 parallelend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 dodo i = 1, 100a(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.shmoduleコマンドのための環境設定
module load intel
Intelコンパイラ環境の読み込み

23

export KMP_AFFINITY=disabledexport 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 forfor (i=0; i<100; i++) {a[i] = i}#pragma omp forfor (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 dodo i = 1, 100a(i) = i
enddo!$omp end do!$omp dodo i = 1, 100b(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 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++)Sep. 18-22, 2023

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

29

[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;}Sep. 18-22, 2023

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

30[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Sep. 18-22, 2023

演習2
演習2-1：omp_ex02_1.c
サンプルのプログラムはループがスレッドで分割されていない。
指示文を挿入（もしくは修正）し、ループを正しく並列化せよ。
まずは omp_ex02_1.c をそのまま動かして挙動を確認しよう
omp_ex02_1.c31#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 =', ienddo!$omp end parallelendSep. 18-22, 2023

演習2
演習2-2：omp_ex02_2.c
演習2-1と同様に指示行を挿入し、ループを並列化せよ。また、計算結果（u）が、スレッド数を1,2,4と変えても変わらないことを確認せよ。
omp_ex02_2.c33#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, 100u(i) = sin(2.0*pi*dble(i)/100.0)! print*, 'myid =', omp_get_thread_num(), 'i =', ienddo!$omp end parallelend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 forfor (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 = 25a[0] = tmpa[25] = tmp

private宣言なし

39#pragma omp forfor (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 = 25a[0] = tmpa[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 forfor (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 section
処理B#pragma omp parallel{#pragma omp sections{#pragma omp section
処理A!$omp end sections!$omp end parallel#pragma omp section
処理B!$omp parallel!$omp sections!$omp section
処理A

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, 100u(i) = sin(2.0*pi*...)
enddo

program ex03...
do i=1,25u(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 forfor (i=0; i<100; i++) {a[i] = i;}#pragma omp single{
output(a);}49#pragma omp forfor (i=0; i<100; i++) {b[i] = i;}}[C] #pragma omp single { }一般的に、スレッドの立ち上げ回数は極力減らすほうがオーバーヘッドが少なくなるため良いとされる。
逐次処理やデータの出力のような処理が入る場合#pragma omp single { }とすると１スレッドのみで処理される

Sep. 18-22, 2023

１スレッドのみで処理!$omp parallel!$omp dodo i = 1, 100a(i) = i
enddo!$omp end do

50!$omp single
call output(a)!$omp end single!$omp dodo 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 forfor (i=0; i<100; i++) {a[i] = i;}#pragma omp master{
output(a);}51#pragma omp forfor (i=0; i<100; i++) {b[i] = i;}}[C] #pragma omp master { }一般的に、スレッドの立ち上げ回数は極力減らすほうがオーバーヘッドが少なくなるため良いとされる。
逐次処理やデータの出力のような処理が入る場合#pragma omp master { }とするとマスタースレッドのみで処理される。
終了時にスレッド間同期は行われない。

Sep. 18-22, 2023

１スレッドのみで処理!$omp parallel!$omp dodo i = 1, 100a(i) = i
enddo!$omp end do

52!$omp mastercall output(a)!$omp end master!$omp dodo 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 forfor (i=0; i<100; i++) {a[i] = i;}#pragma omp critical{a[0]=func(a,b,100);}

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

Sep. 18-22, 2023

排他処理!$omp parallel!$omp dodo i = 1, 100a(i) = i
enddo!$omp end do

54!$omp criticala(1)=func(a,b,100)!$omp end critical!$omp dodo 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*, aprintf("%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(𝜋𝜋 ∗ 𝑖𝑖/𝑛𝑛)

01Sep. 18-22, 2023

SOR法の補足
・基本は、ガウス・ザイデル法（𝜔𝜔 = 1 のとき同じ式になる）
1
𝑘𝑘+1
𝑘𝑘+1
𝑘𝑘
𝑘𝑘
+ 𝑇𝑇𝑖𝑖,𝑗𝑗−1+ 𝑇𝑇𝑖𝑖+1,𝑗𝑗+ 𝑇𝑇𝑖𝑖,𝑗𝑗+1
𝜌𝜌𝑖𝑖,𝑗𝑗 = 𝑇𝑇𝑖𝑖−1,𝑗𝑗
4
𝑘𝑘+1 = 𝑇𝑇 𝑘𝑘
𝑘𝑘
𝑇𝑇𝑖𝑖,𝑗𝑗
𝑖𝑖,𝑗𝑗 + 𝜔𝜔 𝜌𝜌𝑖𝑖,𝑗𝑗 − 𝑇𝑇𝑖𝑖,𝑗𝑗

・注意せずに単純にOpenMPでパラレルforなどの指定をすると「競合状態」に陥り正しく計算ができないことや、結果が実行毎に異なる場合があります。
・一般に、オーダリングを用いて並列化（ハイパープレーン、２色など）
・２色オーダリングの場合：

ブロックを先に計算
1
𝑘𝑘
𝑘𝑘
𝑘𝑘
𝑘𝑘
+
𝑇𝑇
+
𝑇𝑇
𝑇𝑇
+
𝑇𝑇
𝜌𝜌𝑖𝑖,𝑗𝑗 =𝑖𝑖,𝑗𝑗+1
𝑖𝑖−1,𝑗𝑗
𝑖𝑖+1,𝑗𝑗
𝑖𝑖,𝑗𝑗−1460

１．

２．次に、 ブロックを計算
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=863の２か所を書き換える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++)スレッド1for (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++)スレッド1for (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• 富岳のような大型計算機を用いる上では、如何にして並列化率を上げるかが重要である。