OpenMPによるスレッド並列計算

寺尾 剛史（理化学研究所 計算科学研究センター）takeshi.terao@riken.jp
KOBE HPC サマースクール（初級）2022年8月30-31日

講義内容・スレッド並列とは・OpenMPによるループ処理の並列化
・差分化された偏微分方程式の並列化
・アムダールの法則と並列化効率の評価

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

ベクトル計算機
ベクトル演算（同じ計算式を各配列要素に対してそれぞれ計算するといった処理）を一括して行うことができる計算機。一時期は『スーパーコンピュータ ≒ ベクトル機』と言われていた
こともあった。

NEC-SX（地球シミュレータ[初代]）、GPUアクセラレータ等

計算機とムーアの法則
半導体の集積密度は１年半～２年ごとに倍になる
エンドユーザーからすれば、ただ待っているだけで計算機の性能が向上し続け、大規模計算が可能になるという”ありがたい話”

しかし
近年、ムーアの法則の限界が囁かれるようになる→そもそも原子サイズより小さな回路は実現不可能→リーク電流、放熱の問題等

https://ja.wikipedia.org/wiki/%E3%83%A0%E3%83%BC%E3%82%A2%E3%8
1%AE%E6%B3%95%E5%89%87

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

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

スレッド

スレッド

スレッド

PrivatePrivatePrivatePrivateGlobalメモリ空間

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

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

ノード

プロセス

プロセス

メモリ

メモリ

ノード

ノード

プロセス

プロセス

メモリ

メモリ

プロセス間のデータ通信が必要！

並列計算機のイメージ
非並列処理

スレッド並列（共有メモリ）・各CPU（スレッド）は全てのデータにアクセス可能・故にお互いが作業の邪魔してしまうこともある
CPU0

MPI通信

CPU1

プロセス並列（分散メモリ）・各CPU（プロセス）は自身のメモリ空間上のデータのみにアクセス可能・他のプロセスのデータはMPI通信でアクセス

OpenMPとは・Open Multi-Processing の略

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

・指示行の挿入を行うことで並列化できる→既存の非並列プログラムに対し、元のコードの構造を大きく変えることなく並列化できるため、比較的手軽

ちなみに “OpenMPI” というライブラリが存在するが、こち
らはMPI（Message Passing Interface）の実装の１つであり、OpenMPとは全くの別物

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

スレッド０

スレッド1

スレッド2

スレッド3

並列リージョン

Join

スレッド０

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

Fork

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の基本関数
・最大スレッド数取得[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();

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モジュール（ヘッダ）のロードを忘れると、これらの関数を使用できずコンパイルエラーになる並列リージョンを指定#pragma omp parallel{#pragma omp forfor (i=0; i<100; i++) {a[i] = i}#pragma omp single{...
}#pragma omp forfor (...)
.....
}
}

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

並列リージョンを指定!$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, 100
.....
enddo
!$omp end do!$omp end parallel

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

演習準備
演習用ファイル（OpenMP）
rokko:/home/guest59/openmp

上記のファイルをホームディレクトリにコピーする
mkdir ~/openmp
cd ~/openmp
cp -r /home/guest59/openmp/ ./

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

演習1-2：omp_ex01_2.c
スレッド４つで実行、それぞれ自スレッド番号を取得し出力せよ
omp_ex01_1.c#include <stdio.h>
int main(void){
printf("Hello, World!");}#include <stdio.h>
int main(void){#pragma omp parallel{
printf("Hello, World!");
}
}

演習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!"
end

program omp_ex01_1!$omp parallel
print*, "Hello, World!"!$omp end parallelendスレッド数の指定
・シェルの環境変数で与える（推奨）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);ただしスレッド数を変えて実行する時など、毎回コンパイルが必要となってしまうため、今回は環境変数による指定を推奨する。

コンパイル・コンパイルオプションでOpenMPを有効にする
icc -qopenmp omp_ex01_1.c
gcc -fopenmp omp_ex01_1.c

オプションを指定しない場合はOpenMPの指示行はコメントとして認識される。

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

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

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

オプションを指定しない場合はOpenMPの指示行はコメントとして認識される。

!$omp parallel dodo i = 1, 100a(i) = b(i) + c
enddo!$omp end parallel do

OpenMPで用いる指示行は、Fortranの
場合 !$OMP から始まる。行頭に ! がある行は通常、コメントとして処理され
る。

操作補足
ジョブスクリプト
omp_run.sh

ジョブ投入方法はマシンの環境によって変わるため、実際にスパコンを使う時にはユーザーマニュアルを参考にすること#!/bin/bash#PBS –q S#PBS -l select=1:ncpus=4#PBS -N omp_JOB#PBS –o omp_out#PBS -j oe

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

source /etc/profile.d/modules.shmoduleコマンドのための環境設定
module load intel
Intelコンパイラ環境の読み込み
export KMP_AFFINITY=disabledexport OMP_NUM_THREADS=4
cd ${PBS_O_WORKDIR}./a.out
実行

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

Working Sharing構文〇 複数のスレッドで分担して実行する部分を指定〇 並列リージョン内で記述する#pragma omp parallel { } の括弧範囲内

指示文の基本形式は[C] #pragma omp xxx[F] !$omp xxx ~ !$omp end xxx◎for構文, do構文
ループを分割し各スレッドで実行◎section構文
各セクションを各スレッドで実行◎single構文
１スレッドのみ実行◎master構文
マスタースレッドのみ実行

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

[F] #pragma omp for#pragma omp parallel でスレッドを生成しただけでは、全てのスレッドが全ループを計算してしまう#pragma omp for を入れることでループ自体が分割され、各スレッドに処理が割り当てられるdo構文!$omp parallel!$omp dodo i = 1, 100a(i) = i
enddo!$omp end do

doループをスレッドで分割し、並列処理を行う!$omp dodo i = 1, 100b(i) = i
enddo!$omp end do!$omp end parallel・do の直前に指示行 !$omp do を入れる

[F] !$omp do ～ !$omp end do

・enddo の直後に指示行 !$omp end do を入れる!$omp parallel でスレッドを生成しただけでは、全てのスレッドが全ループを計算してしまう!$omp do を入れることでループ自体が分割され、各スレッドに処理が割り当てられる

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

演習2
演習2-1：omp_ex02_1.c
サンプルのプログラムはループがスレッドで分割されていない。
指示文を挿入（もしくは修正）し、ループを正しく並列化せよ。
まずは omp_ex02_1.c をそのまま動かして挙動を確認しよう
omp_ex02_1.c#include <stdio.h>...
#pragma omp parallel{for (i=0; i<10; i++) {printf("myid=%d, i=%d", omp_get_thread_num(), i);
}
}
}

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

program omp_ex02_1...
!$omp parallel
do i=1, 10
print*, 'myid =', omp_get_thread_num(), 'i =', ienddo!$omp end parallel
end

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

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

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

プライベート変数について起こりがちなミス#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宣言を入れる

並列化したループ内で値を設定・更新する場合は要注意→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;}

処理順

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

tmp = 0tmp = 25

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

スレッド1

a[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(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);}}}}

OpenMPを用いた並列化では、内側ループ、外側ループのどちらを並列化しても良い。ただし、内側ループを並列化する場合、毎回forkjoinしてしまうとスレッド生成回数がものすごいことになる。
（上記の例では１つの３次元ループで 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

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!$omp end sections!$omp end parallel#pragma omp section
処理B}}

処理の分岐Section構文
スレッドごとに処理を分岐させる#pragma omp parallel sections

スレッド０#pragma omp section
処理0

スレッド1
#pragma omp section
処理1

待ち時間
待ち時間

スレッド2
#pragma omp section
処理2

スレッド3
#pragma omp section
処理3

待ち時間

同期

・各スレッドに割り当てられた処理の負荷が異なると、無駄な待ちが発生する
・ロードバランスに注意演習3
演習3-1：omp_ex03.c
演習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*…);}

演習3
演習3-1：omp_ex03.f90
演習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

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

１スレッドのみで処理!$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, 100
.....
enddo
!$omp end do!$omp end parallel[F] !$omp single ~!$omp end single
一般的に、スレッドの立ち上げ回数は極力減らすほうがオーバーヘッドが少なくなるため良いとされる。
逐次処理やデータの出力のような処理が入る場合、!$omp singleとすると１スレッドのみで処理される

リダクション変数を指定
reduction(演算子:変数）・並列計算時はそれぞれのスレッドで別々の値を持ち、並列リージョン終了時に各スレッドの値が足し合わされる（総和）変数・総和の他、積などを求めることも可能
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;}

演習4
演習4-1：omp_ex04.f90 / omp_ex04.c
1から100までを足し合わせるプログラムをOpenMPで並列化せよ。
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);並列化できないプログラム
・どのループを並列化するか（可能か）という判断は全てプログラマが行う
並列化できないプログラム

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

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

演習5
演習5-1：
熱伝導問題（Laplace方程式）のプログラム（laplace.c / laplace.f90）をOpenMPでスレッド並列化して計算せよ基礎方程式
𝜕2𝑇 𝜕2𝑇
+ 2=0
2
𝜕𝑥
𝜕𝑦
中心差分（2次精度）

1
𝑇𝑖,𝑗 = (𝑇𝑖−1,𝑗 + 𝑇𝑖+1,𝑗 + 𝑇𝑖,𝑗−1 + 𝑇𝑖,𝑗+1 )4

演習5
演習5-1：
熱伝導問題のプログラム（laplace.c / laplace.f90）をOpenMPでスレッド並列化して計算せよ

SOR反復法
1 𝑘+1
𝑘
𝑘+1
𝑘
𝜌𝑖,𝑗 = 𝑇𝑖−1,𝑗 + 𝑇𝑖+1,𝑗+ 𝑇𝑖,𝑗−1+ 𝑇𝑖,𝑗+1
4
𝑘+1
𝑘
𝑘
𝑇𝑖,𝑗
= 𝑇𝑖,𝑗+ 𝜔 𝜌𝑖,𝑗 − 𝑇𝑖,𝑗
境界条件

1

𝑇𝑖,0 = 𝑇𝑁+1,𝑗 = 𝑇𝑖,𝑁+1 = 0
𝑇0,𝑗 = sin(𝜋 ∗ 𝑖/𝑛)
0

1

演習5
演習5-1：
熱伝導問題のプログラム（ laplace.c / laplace.f90）
をOpenMPでスレッド並列化して計算せよ
補足（laplace.c / laplace.f90)
n：縦、横それぞれのグリッド数
ITMAX : 最大反復回数

eps：反復ベクトルの差のノルムの閾値（どこまで値を収束させるか）演習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
カラーバーの値の範囲等は計算のパラメータによって適宜変更すること

カラーバーの色設定

演習5
演習5-2：
スレッド並列化した熱伝導問題のプログラムを、並列数を1, 2, 4, 8, 16, … と変えて実行し、計算時間を計測せよ

並列数はジョブスクリプトにてCPU確保：#PBS -l select=1:ncpus=8

スレッド数：

export OMP_NUM_THREADS=8の２か所を書き換える演習5
演習5-2：
時間計測（簡易版）：
ジョブスクリプトの実行で time を使用
time dplace ./a.out

時間計測（関数）：
omp_get_wtime 関数を使用
計算開始時に
dts = omp_get_wtime()
計算終了時に
dte = omp_get_wtime()とし、 dte - dts を出力する

SOR法の補足
・基本は、ガウス・ザイデル法（𝜔 = 1 のとき同じ式になる）
1 𝑘+1
𝑘
𝑘+1
𝑘
𝜌𝑖,𝑗 = 𝑇𝑖−1,𝑗 + 𝑇𝑖+1,𝑗+ 𝑇𝑖,𝑗−1+ 𝑇𝑖,𝑗+14
𝑘+1
𝑘𝑘𝑇𝑖,𝑗 = 𝑇𝑖,𝑗+ 𝜔 𝜌𝑖,𝑗 − 𝑇𝑖,𝑗・オーダリングを用いて並列化（ハイパープレーン、２色など）
・２色オーダリングの場合：

１．

ブロックを先に計算

1 𝑘
𝑘
𝑘
𝑘
𝜌𝑖,𝑗 = 𝑇𝑖−1,𝑗 + 𝑇𝑖+1,𝑗+ 𝑇𝑖,𝑗−1+ 𝑇𝑖,𝑗+1
4
２．次に、ブロックを計算

1 𝑘+1
𝑘+1
𝑘+1
𝑘+1
𝜌𝑖,𝑗 = 𝑇𝑖−1,𝑗 + 𝑇𝑖+1,𝑗+ 𝑇𝑖,𝑗−1+ 𝑇𝑖,𝑗+14

並列数と計算性能

並列数

1

2

4

8

16

32

40

実測値

5.18E+00

2.69E+00

1.44E+00

8.08E-01

6.25E-01

9.04E-01

1.12E+00

理論値

5.18E+00

2.59E+00

1.30E+00

6.48E-01

3.24E-01

1.62E-01

1.30E-01

スケジューリング#pragma omp parallel for

スレッド０for (i=0; i<50; i++)スレッド1for (i=50; i<100; i++)

待ち時間

・各スレッドに割り当てられた処理の負荷が異なると、無駄な待ちが発生する
・ロードバランスに注意スケジューリング#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
for (i=75; i<100; i++)・各スレッドに割り当てられた処理の負荷が異なると、無駄な待ちが発生する
・ロードバランスに注意スケジューリング#pragma omp parallel for schedule(dynamic, 25)

スレッド０for (i=0; i<25; i++)スレッド０
for (i=50; i<75; i++)スレッド1for (i=25; i<50; i++)スレッド０
for (i=75; i<100; i++)

待ち時間

・各スレッドに割り当てられた処理の負荷が異なると、無駄な待ちが発生する
・ロードバランスに注意スケジューリング• 例：（上三角行列）＊（ベクトル）• 𝑦 = 𝐴 ∗ 𝑥, (𝑎𝑖𝑗 = 0, 𝑖 > 𝑗)for (i=0; i<n; i++) {for (j=i; j<n; j++) {y[i] += a[i][j]*x[j];}}• チャンクサイズが性能に大きく影響する。
• 負荷バランスとシステムのオーバーヘッドのトレードオフ。
• チューニングコストが増加する。

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

並列化率に関して依存関係により並列化できない処理

逐次処理

①

②

③

④

⑤

⑥

並列化可能な処理

②
２並列

④

①

⑥
③

⑤

2スレッド並列計算時間は2/3（計算速度は1.5倍）

②
４並列

①

④
③
⑤

例：並列化できる部分がプログラム全体のうち2/3の場合⑥

4スレッド並列計算時間は半分（計算速度は2倍）

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

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