2021年度
はじめての並列化

２０２１年 6月 ２日
東北大学サイバーサイエンスセンター
日本電気株式会社

本資料は，東北大学サイバーサイエンスセンターとＮＥＣの共同により作成された．
無断転載等は，ご遠慮下さい．

Page 2

 並列化概要編
 OpenMPプログラミング編

Page 3

演習問題の準備▌演算問題の環境を自分のホームディレクトリ配下にコピーします。
/mnt/stfs/ap/lecture/parallel/|-- practice_1|-- practice_2|-- practice_3|-- practice_4$ cd <環境をコピーしたいディレクトリ>$ cp –r /mnt/stfs/ap/lecture/parallel .

Page 4

演習問題１
演習問題2
演習問題3
演習問題4

並列化概要Page 5

1. 並列化とは▌並列処理・並列実行仕事（処理）を複数のコアに分割し，同時に実行すること▌並列化

1CPUの経過時間

処理A

処理B

処理A

処理A

処理A

処理A

処理B

処理B

処理B

処理B

処理C

処理C

処理C

処理C

並列化により短縮する時間
処理C

Page 6

4CPUの経過時間

並列処理を可能とするために，処理の分割を行うこと

処理Bを複数のコアで並列に実行

2. 並列化の効果▌姫野ベンチマーク・ポアッソン方程式解法をヤコビの反復法で解く場合に主要なループの処理
 AOBA-A 1コアの実行時間は約 29.2 秒
速度を測るもの
プログラムの一部抜粋
gosa=0.0
DO K=2,kmax-1
DO J=2,jmax-1
DO I=2,imax-1
S0=a(I,J,K,1)*p(I+1,J,K)+a(I,J,K,2)*p(I,J+1,K)
1+a(I,J,K,3)*p(I,J,K+1)2+b(I,J,K,1)*(p(I+1,J+1,K)-p(I+1,J-1,K)
3-p(I-1,J+1,K)+p(I-1,J-1,K))
4+b(I,J,K,2)*(p(I,J+1,K+1)-p(I,J-1,K+1)5-p(I,J+1,K-1)+p(I,J-1,K-1))
6+b(I,J,K,3)*(p(I+1,J,K+1)-p(I-1,J,K+1)
7-p(I+1,J,K-1)+p(I-1,J,K-1))
8+c(I,J,K,1)*p(I-1,J,K)+c(I,J,K,2)*p(I,J-1,K)
9+c(I,J,K,3)*p(I,J,K-1)+wrk1(I,J,K)
SS=(S0*a(I,J,K,4)-p(I,J,K))*bnd(I,J,K)
GOSA=GOSA+SS*SS
wrk2(I,J,K)=p(I,J,K)+OMEGA *SS
enddo
enddo
enddo

******** Program Information ********
Real Time (sec):User Time (sec):
Vector Time (sec):Inst. Count:
V. Inst. Count:
V. Element Count:
V. Load Element Count:FLOP Count:MOPS:
MOPS (Real):MFLOPS:MFLOPS (Real):
A. V. Length:
V. Op. Ratio (%):
L1 Cache Miss (sec):
CPU Port Conf. (sec):
V. Arith. Exec. (sec):
V. Load Exec. (sec):
VLD LLC Hit Element Ratio (%):
FMA Element Count:Power Throttling (sec):
Thermal Throttling (sec):Memory Size Used (MB):Non Swappable Memory Size Used (MB):29.197789
29.195756
29.140397
48890015000
9526901371
2379556386510
1290401280363
1411535186639
98044.806708
98038.074875
48347.227656
48343.908098
249.772334
98.624867
0.043942
0.000000
14.230997
14.907927
57.370385
443575440000
0.000000
0.000000
708.000000
98.000000

 複数のコアを用いることで実行時間の短縮が可能に
Page 7

3．演習問題1（practice_1)▌姫野ベンチマークの1コア実行
項目

対象

備考

作業ディレクトリ

practice_1

使用ソースファイル

sample1.f

編集 不要

ジョブファイル

run.sh

そのまま投入 手順①：作業ディレクトリを移動してください．
% cd parallel/practice_1

 手順②：コンパイルします．
% nfort sample1.f

 手順③：ジョブを投入します．
% qsub run.sh

Page 8

3．演習問題1（practice_1)

 手順④：結果を確認します．結果はp1-practice.o.XXXX(XXXXにはジョブIDが入ります)として格納されます．
% cat p1-practice.o.XXXX

Page 9

4. 並列化の効果（自動並列）▌自動並列化により複数のコアを利用し，実行時間を短縮
AOBA-A 1VE(8コア)
30

29.19

実行時間（秒）

25
20

15.15

15
10

8.03
4.46

5
0

1

2

コア数

Page 1048

5．演習問題2（practice_2)▌姫野ベンチマークの自動並列実行
項目

対象

備考

作業ディレクトリ

practice_2

使用ソースファイル

sample2.f

編集 不要

ジョブファイル

run.sh

そのまま投入 手順①：作業ディレクトリを移動してください．
% cd parallel/practice_2

 手順②：コンパイルします．
% nfort -mparallel sample2.f

 手順③：ジョブを投入します．
% qsub run.sh

Page 11

5．演習問題2（practice_2)

 手順④：結果を確認します．結果はp2-practice.o.XXXX(XXXXにはジョブIDが入ります)として格納されます．
% cat p2-practice.o.XXXX

ジョブファイル（run.sh）の OMP_NUM_THREADS で設定した数値が，実行時の並列数となります．
OMP_NUM_THREADSの数値を1～8（AOBA-Aは8コア）に変更することで，実行時間が変化することを確認してください．

Page 12

6. 並列化率▌並列に実行可能(あるいは効果のある)部分と並列に実行不可(あるいは効果のない)部分を見つけ，並列に実行可能部分を複数のコア(CPU)に割り当てる．
▌できるだけ多くの部分を並列化の対象としなければ，コア数に応じた効果が得られない．

並列化率α

並列化対象部分の処理時間＝
全体の処理時間(並列化の対象部分と非対象部分の時間の合計)

Page 13

7. 並列処理の限界（アムダールの法則）
N=4

N=8

N=16

N=32

並列化の効果(性能倍率S)

32.0

1
24.0

S=

1 – α+

α
n

16.0

•

Nはコア(CPU)数•

並列化率が100%から下がるにしたがって性能倍率は急速に
低下する

8.0

0.0

並列化率α
並列化率と並列化の効果の関係(アムダールの法則)

並列化率100%はあり得ない(データの入出力で必ず逐次処
理発生)が，可能な限り100%に近づかなければ並列化の大きな効果は得られない
Page 14

8. 並列処理モデル
コンピュータアーキテクチャに応じた処理の分担(分割)のさせ方によって幾つかの並列処理がある
1. 分散メモリ並列処理––

AOBA-A(マルチノード)
AOBA-B(マルチノード)

2. 共有メモリ並列処理––

CPU(core)

CPU(core)

CPU(core)

CPU(core)

CPU(core)

CPU(core)

CPU(core)

CPU(core)

AOBA-A(シングルノード)AOBA-B(シングルノード)

 MPI(Message Passing Interface)は分散メモリ並列処理のための並列手法
 OpenMPは共有メモリ並列処理のための並列手法
Page 15

8. 並列処理モデル▌並列処理モデルの特徴
共有メモリ並列処理

分散メモリ並列処理

メモリ空間

すべてのコアが同じメモリ空間
をアクセス可能

ネットワーク上のすべてのメモ
リ空間をデータ通信によりアクセス可能

利用可能な並列化手法

自動並列化
OpenMP並列化
MPI並列化

MPI並列化(※1)

利用可能なAOBA-Aのコア数

8コア(1VE)

2,048コア(256VE)

利用可能なAOBA-Aの
メモリ容量

48GByte

約12TByte※1：共有メモリ内は共有メモリ並列処理，分散メモリ間は分散メモリ並列処理のハイ
ブリッド並列処理が可能

Page 16

9. 並列化の対象
物理現象

Page 17

プログラム化プログラム

空間的並列性

空間による並列化（領域分割法）現象の並列性処理による並列化

9. 並列化の対象空間による並列化（イメージ）
領域分割法
処理１
Ｄ
Ｃ
Ｂ
Ａ

各々，コア(CPU)に割り当てる

Page 18

処理２
．
．
．

処理ｎ

空
空
空
間
間
間
Ａ
Ｂ
Ｃ
プロセス０ プロセス１ プロセス２

…

処理１

処理１

処理１

処理２

処理２

処理２

．
．

．
．

．
．

処理 ｎ

処理 ｎ

処理 ｎ

時
間
短
縮

9. 並列化の対象空間による並列化の例▌DOループ（FORTRAN）単位での並列処理例）領域分割法

do iz=1,100
do iy=1,100
do ix=1,100
処理１
：
処理ｎ
enddo
enddo
enddo

do iz=1,25
do iy=1,100
do ix=1,100
処理１
do iz=26,50：
do iy=1,100
処理ｎ
do ix=1,100
enddo
処理１
do iz=51,75
enddo：
do iy=1,100
enddo
処理ｎ
do ix=1,100
do iz=76,100
enddo
処理１
do ：
iy=1,100
enddo
do ix=1,100
enddo
処理ｎ
処理１
enddo
：
enddo
処理ｎ
enddo
enddo
enddo
enddo

より外側のループで並列化することが重要

Page 19

OpenMPプログラミング編

Page 20

１．OpenMP概要
OpenMPとは 共有メモリマシン環境における並列化方法の１つ．
 指示行によって並列化可能な場所を指示．
 プログラム実行時に，並列数を指定して実行．
 Fortranおよび，C言語で使用可能．
共有メモリシステム

メモリ

Page 21

CPU

CPU

CPU

(core)(core)(core)

・・・

CPU(core)

１．OpenMP概要
OpenMP並列の基本構造
 !$OMP parallel と !$OMP end parallel指示文に囲まれた範囲の処理を複数のスレッドで実行する．
 スレッド数は実行時の環境変数で指定する．

処理A

逐次実行

!$OMP PARALLEL

処理B

並列実行
・・・

!$OMP END PARALLEL

処理C

Page 22

逐次実行

環境変数で指定したスレッド数で並列に実行

2．OpenMP指示文(指示行)の書き方
 Fortranプログラムの場合「!$OMP [指示文]」を1カラム目から指定する．DO文を並列化する場合
はDO文の前に「parallel do」を指示し，ENDDO文の後ろに「end parallel
do」を指示する．
（例）!$OMP parallel dodo k=1,nz
do j=1,ny
:
enddo
enddo
!$OMP end parallel do Cプログラムの場合「#pragma omp [指示文]」で指定する．for文を並列化する場合は「parallel for」と指示する．
（例）#pragma omp parallel for
for(i=0;i<100;i++){:Page 23

3. OpenMP化の例（総和計算プログラム）
 1から1000の総和を求める(逐次実行プログラム)
program sample1
parameter(n=1000)
integer sum
sum=0
do i=1,nsum=sum+i
enddo
print *,"Total = ",sum
stop
end

 並列化の対象はDOループ

 1から1000の総和を計算するループが対象．
 ループ内の変数sumは総和演算を行っており，各スレッドが個々に変数sumの内容を更新すると正しい結果が得られない．

Page 24

3. OpenMP化の例（総和計算プログラム）
 DOループの並列化は parallel do ～ end parallel do指示文でス
レッド並列化を指示!$OMP parallel dodo i=1,nsum=sum+i
enddo!$OMP end parallel do

スレッド０
do i=1,125
sum=sum+i
enddo

Page 25

 !$OMP parallel do～!$OMP end parallel do
で囲んだDOループはスレッド数で分割され並列に実行される
 この例では，8スレッドで実行する場合，ループは1～125，126～250，251～375，376～
500，501～625，626～750，751～875，
876～1000とスレッドに割り当てられる

スレッド１
do i=126,250
sum=sum+i
enddo

スレッド２
do i=251,625
sum=sum+i
enddo

スレッド7

・・・

do i=876,1000
sum=sum+i
enddo3. OpenMP化の例（総和計算プログラム）
複数のスレッドが同時に同じ変数に書き込む可能性がある場合には注意が必要
 各スレッドの計算結果は変数sumに格納するが，全スレッドが同じ領域の値を更新すると正しい結果は得られないため，特別な指定が必要．
 スレッドごとに部分和を格納する領域を用意し，各スレッドの計算結果の部分和を変数sumに足し込む必要がある．
 このような演算を「リダクション演算」と呼ぶ．
 変数sumに対する演算はリダクション演算であることを明示する．
!$OMP parallel do reduction(+:sum)
do i=1,nsum=sum+i
enddo!$OMP end parallel do reduction(オペレータ:変数名リスト） オプションを記述する．オペレータの種類として和(+)，差(-)，積(＊)のほかに最大・最小値を求めるMAX，MINやビット操作を行うIAND，IOR，IEOR，論理演算のAND.，.OR.，.EQV.，.NEQV.を使用することができる.
 複数の変数を指定する場合は，変数名リストをカンマ(,)で区切って記述する．
Page 26

4．OpenMPの記述法
スレッドごとに独立した変数が必要になる場合には注意が必要

 PRIVATE指定
スレッドごとに独立した変数が必要な場合には，PRIVATE句が必要である．
①ループ中で作業配列・変数として使用されている②並列化の対象となるループ変数の次元を持たない配列（例）
do k=1,nz
do j=1,ny
do i=1,nxa(i,j) = b(i,j,k)*c(i,k)
d(i,j,k) = d(i,j,k) + a(i,j)*e(j,k):!$OMP parallel do private(a)do k=1,nz
do j=1,ny
do i=1,nx:

Page 27

 配列aはi,jの次元しか持たず，ループ中で定義・参照されている．
 このような配列はスレッドごとに領域を確保しなければ，正しい結果を得ることができない．
 配列aに対して!$OMP parallel do指示文
でprivate指定を行うことにより，配列aはスレッドごとに独立した領域に確保される．
 private(a,b,c)のように複数の配列，変数を指定することが可能．

4．OpenMPの記述法
 do～end do指示文(parallel～end parallel) 連続する複数のループが並列化の対象となる場合，parallel do～end
parallel do指示文を複数書くのではなく，do～end do指示文を組合せる．
!$OMP parallel!$OMPdo
並列化対象のループ!$OMP end do!$OMPdo
並列化対象のループ!$OMP end do!$OMP end parallel parallel～end parallelの範囲を広げることにより，OpenMPによるオー
バーヘッドを軽減することが可能になる．

Page 28

4．OpenMPの記述法
 実行時の総スレッド数や自スレッド番号は関数呼出しで取得可能．
omp_get_num_threads( ):総スレッド数の取得
omp_get_thread_num( ) :自スレッド番号の取得

 スレッド番号で処理を制御する際に使用．
 OpenMP実行時ルーチンを使用する場合は，moduleのomp_libを使用することを宣言する．

program sample2
use omp_lib!$OMP parallel
print *,"Hello World",omp_get_num_threads(),omp_get_thread_num()!$OMP end parallel
stop
end

総スレッド数と自スレッド番号を取得しprintする

Page 29# setenv OMP_NUM_THREADS 4# ./a.out
Hello World
40
Hello World
4
1
Hello World
4
2
Hello World
43

5. コンパイル・実行スクリプト OpenMPプログラムをコンパイルするには，OpenMP専用のコンパイルオプションを指定してコンパイルする必要がある．
 並列数は，実行時環境変数（OMP_NUM_THREADS）で指定す
る．
○ AOBA-Aの場合
nfort -fopenmp [オプション] ソースファイル名※OpenMPを利用する場合は-fopenmpオプションが必須※自動並列化-mparallelとの併用が可能

実行のスクリプト例#!/bin/sh#PBS -q ①#PBS --venode ②#PBS -l elapstim_req=③#PBS -v OMP_NUM_THREADS=④
cd $PBS_O_WORKDIR./a.out (実行文)
Page 30① AOBA-Aのキューを指定．
通常利用の場合は「sx」と指定．※占有利
用の方は別途お知らせします．(必須)② 使用VE数を指定．(必須)③ 使用計算時間（経過時間）を設定．
「hh:mm:ss」のように指定．
例えば1時間半の場合には#PBS -l elapstim_req=01:30:00と指定．(設定を強く推奨)④ スレッド数を指定（最大は8） ．
(共有並列の場合必須)

5. コンパイル・実行スクリプト○ AOBA-Bの場合
サイバーサイエンスセンターのAOBA-BでもOpenMPの使用は可能．
OpenMPの使い方はAOBA-Aと同じ．
フロントエンド(front.cc.tohoku.ac.jp)上で，以下のようにコンパイルする．
・AOCCコンパイラ
flang –fopenmp [オプション] ソースファイル名
・GCCコンパイラ
gfortran –fopenmp [オプション] ソースファイル
・Intelコンパイラ
ifort -qopenmp [オプション] ソースファイル名

Page 31

5. コンパイル・実行スクリプト
コンパイルオプションについては，以下を参照
AOBA-A「AOBA-A 利用方法」
https://www.ss.cc.tohoku.ac.jp/sx-aurora/#toc4「Fortran Compiler ユーザーズガイド」
https://www.hpc.nec /documentation
AOBA-B

Page 32「AOBA-B 利用方法」
https://www.ss.cc.tohoku.ac.jp/lx406rz-2/#toc4「AMD Optimizing C/C++ Compiler」
https://developer.amd.com/amd-aocc/5. コンパイル・実行スクリプト AOBA-Aのジョブクラス
利用形態

キュー名

VE数

実行形態（※）

最大経過時間
既定値/最大値

メモリサイズ

無料

sxf

1

1VE

1時間/1時間

48GB×VE数

共有

sx

1

1VE

72時間/720時間

2～256

8VE単位で確保(VHを共用しない)

2～8

1VE単位で確保(VHを共用する)

sxmix※OpenMPのみの並列では複数VEは利用できません．
※MPIとのハイブリッド実行については2021年度後期の講習会にて詳しく説明します．

 AOBA-Bのジョブクラス
利用形態

キュー名

ノード数
(-bオプション)

最大経過時間
既定値/最大値

メモリサイズ

共有

lx

1～16

72時間/720時間

256GB×ノード数※OpenMPのみの並列では複数ノードは利用できません．
※MPIとのハイブリッド実行については2021年度後期の講習会にて詳しく説明します．
Page 33

6．演習問題３（practice_3)▌姫野ベンチマークコードをOpenMPで並列化し，逐次実行（1core実行），
自動並列実行と実行性能を比較する．
項目

対象

備考

作業ディレクトリ

practice_3

使用ソースファイル

sample3.f

編集 必要

ジョブファイル

run.sh

そのまま投入 手順①：作業ディレクトリを移動してください．
% cd parallel/practice_3

 手順②：エディタでソースファイルを編集し，OpenMP並列化を実施してください．

Page 34

6．演習問題３（practice_3)

 ヒント：サブルーチンjacobiの以下のループで202行目のループは並列化不可(ただし、204行目、224行目からの3重ループは並列化可能）
202: +------>
203: |
204: |+----->
205: ||+---->
206: |||+--->
207: ||||
208: ||||
209: ||||
210: ||||
211: ||||
212: ||||
213: ||||
214: ||||
215: ||||
216: ||||
217: ||||
218: ||||
219: ||||
220: |||+--221: ||+---222: |+----223: |224: |+----->
225: ||+---->
226: |||+--->
227: ||||
228: |||+--229: ||+---230: |+----231: |
232: +-----Page 35

DO loop=1,nn
gosa=0.0
DO K=2,kmax-1
DO J=2,jmax-1
DO I=2,imax-1
S0=a(I,J,K,1)*p(I+1,J,K)+a(I,J,K,2)*p(I,J+1,K)
1+a(I,J,K,3)*p(I,J,K+1)2+b(I,J,K,1)*(p(I+1,J+1,K)-p(I+1,J-1,K)
S0
3-p(I-1,J+1,K)+p(I-1,J-1,K))（エスゼロ）
4+b(I,J,K,2)*(p(I,J+1,K+1)-p(I,J-1,K+1)5-p(I,J+1,K-1)+p(I,J-1,K-1))
6+b(I,J,K,3)*(p(I+1,J,K+1)-p(I-1,J,K+1)
7-p(I+1,J,K-1)+p(I-1,J,K-1))
8+c(I,J,K,1)*p(I-1,J,K)+c(I,J,K,2)*p(I,J-1,K)
9+c(I,J,K,3)*p(I,J,K-1)+wrk1(I,J,K)
SS=(S0*a(I,J,K,4)-p(I,J,K))*bnd(I,J,K)
GOSA=GOSA+SS*SS
wrk2(I,J,K)=p(I,J,K)+OMEGA *SS
enddo
enddo
enddo
C
DO K=2,kmax-1
DO J=2,jmax-1
DO I=2,imax-1
p(I,J,K)=wrk2(I,J,K)
enddo
enddo
enddo
C
enddo

6．演習問題３（practice_3)

 手順③：コンパイルします．
% nfort -fopenmp sample3.f

 手順④：ジョブを投入します．
% qsub run.sh 手順⑤：結果を確認します．結果はp3-practice.o.XXXX(XXXXにはジョブIDが入ります)として格納されます．
% cat p3-practice.o.XXXX

ジョブファイル（run.sh）の OMP_NUM_THREADS で設定した数値が，実行時の並列数となります．
OMP_NUM_THREADSの数値を1～8（AOBA-Aは8コア）に変更することで，実行時間が変化することを確認してください．
また，OpenMPの実行結果と，演習問題1，演習問題2の実行結果を比較してください．
Page 36

６．演習問題３（practice_3)

つづき 時間計測は以下のように修正ずみ．
38 C ------------------39 C "use portlib" statement on the next line is for Visual fortran40 C to use UNIX libraries. Please remove it if your system is UNIX.
41 C ------------------42 ! use portlib
43
use omp_lib
44
IMPLICIT REAL*4(a-h,o-z)45real*8 t1,t2

行追加

(省略)
81 !
Cpu0=dtime(time0)82
t1=omp_get_wtime( )83 C
84 C Jacobi iteration
85
call jacobi(nn,gosa)
86 C
87 !
cpu1= dtime(time1)
88
t2=omp_get_wtime( )89 !
cpu = cpu190
cpu = t2-t1
91
flop=real(kmax-2)*real(jmax-2)*real(imax-2)*34.0*real(nn)

修正(一つ上の行と置換)

Page 37

６．演習問題３（practice_3)
 サブルーチンinitmtの以下の148 !$OMP
149 !$OMP
150
do k=1,mkmax
151
do j=1,mjmax
152
do i=1,mimax153
a(i,j,k,1)=0.0154
a(i,j,k,2)=0.0155
a(i,j,k,3)=0.0156
a(i,j,k,4)=0.0157
b(i,j,k,1)=0.0
158
b(i,j,k,2)=0.0159b(i,j,k,3)=0.0160
c(i,j,k,1)=0.0
161
c(i,j,k,2)=0.0
162
c(i,j,k,3)=0.0163
p(i,j,k) =0.0
164
wrk1(i,j,k)=0.0165bnd(i,j,k)=0.0
166
enddo
167
enddo
168
enddo
169 !$OMP①

つづきの部分を埋めてくだい。

③

171 !$OMP172
do k=1,kmax
173
do j=1,jmax
174
do i=1,imax
175
a(i,j,k,1)=1.0
176
a(i,j,k,2)=1.0
177
a(i,j,k,3)=1.0
178
a(i,j,k,4)=1.0/6.0179
b(i,j,k,1)=0.0180b(i,j,k,2)=0.0
181b(i,j,k,3)=0.0③②
182
c(i,j,k,1)=1.0
183
c(i,j,k,2)=1.0
184
c(i,j,k,3)=1.0
185p(i,j,k) =float((k-1)*(k-1))/float((kmax-1)*(kmax-1))
186
wrk1(i,j,k)=0.0
187bnd(i,j,k)=1.0
188
enddo
189
enddo
190
enddo
191 !$OMP
192 !$OMP① 150行目から168行目までのループが並列化可能か検討します．
② 172行目から190行目までのループが並列化可能か検討します．
③ ①および②ともに並列化可能であれば，2つの連続するループが並列化の対象となるため，OpenMP並列化のオーバーヘッドを軽減するよう， parallel～end parallel 指示文およびdo～end do指示文を組み合わせて記述します．
Page 38

６．演習問題３（practice_3)
 サブルーチン jacobi の以下のの部分を埋めてくだい.

218
DO loop=1,nn
219
gosa=0.0
220 !$OMP③
221 !$OMP
222
DO K=2,kmax-1①
223
DO J=2,jmax-1
224
DO I=2,imax-1
225
S0=a(I,J,K,1)*p(I+1,J,K)+a(I,J,K,2)*p(I,J+1,K)226
1+a(I,J,K,3)*p(I,J,K+1)
227
2+b(I,J,K,1)*(p(I+1,J+1,K)-p(I+1,J-1,K)
228
3-p(I-1,J+1,K)+p(I-1,J-1,K))229
4+b(I,J,K,2)*(p(I,J+1,K+1)-p(I,J-1,K+1)
230
5-p(I,J+1,K-1)+p(I,J-1,K-1))
231
6+b(I,J,K,3)*(p(I+1,J,K+1)-p(I-1,J,K+1)232
7-p(I+1,J,K-1)+p(I-1,J,K-1))
233
8+c(I,J,K,1)*p(I-1,J,K)+c(I,J,K,2)*p(I,J-1,K)
234
9+c(I,J,K,3)*p(I,J,K-1)+wrk1(I,J,K)
235
SS=(S0*a(I,J,K,4)-p(I,J,K))*bnd(I,J,K)
236
GOSA=GOSA+SS*SS
237
wrk2(I,J,K)=p(I,J,K)+OMEGA *SS
238
enddo
239
enddo
240
enddo
241 !$OMPPage 39

つづき243 !$OMP
244
DO K=2,kmax-1
245
DO J=2,jmax-1
246
DO I=2,imax-1
247
p(I,J,K)=wrk2(I,J,K)
248
enddo
249
enddo
250
enddo
251 !$OMP
252 !$OMP
253 C
254
enddo

②

③

① 222行目から240行目までのループが並列化可能か検討します．並列化に必要なプライベート
指定とリダクション演算指定を検討します．
② 244行目から250行目までのループが並列化可能か検討します．
③ ①および②ともに並列化可能であれば，2つの連続するループが並列化の対象となるため，OpenMP並列化のオーバーヘッドを軽減するよう， parallel～end parallel 指示文およびdo～end do指示文を組み合わせて記述します．

7．演習問題４（practice_4)▌姫野ベンチマークコードをOpenMPで並列化し，AOBA-Bで実行し，AOBA-Aと比較する．
項目

対象

備考

作業ディレクトリ

practice_4

使用ソースファイル

sample3.f

演習問題3の結果を使用

ジョブファイル

run.sh

そのまま投入 手順①：作業ディレクトリを移動してください．
% cd parallel/practice_4

 手順②：ソースファイルをコピーします．
% cp ../practice_3/sample3.f .

 手順③：コンパイルします．
% flang –fopenmp sample3.f
Page 40

7．演習問題４（practice_4)

 手順③：ジョブを投入します．
% qsub run.sh 手順④：結果を確認します．結果はp4-practice.o.XXXX(XXXXにはジョブIDが入ります)として格納されます．
% cat p4-practice.o.XXXX

ジョブファイル（run.sh）の OMP_NUM_THREADS で設定した数値が，実行時の並列数となります．
OMP_NUM_THREADSの数値を8以上（AOBA-Bは最大128コア）に変更することで，実行時間が変化することを確認してください．
また，OpenMPの実行結果と，演習問題3の実行結果を比較してください．

Page 41

8．参考(参考文献)
 OpenMPの仕様はhttp://www.openmp.org/で公開されている．
 現在日本語で出版されている参考書籍は「OpenMPによる並列プログラ
ミングと数値計算法」(牛島省著,丸善株式会社)や「OpenMP並列プログ
ラミング」(菅原清文著，株式会社カットシステム)などがある．

Page 42

演習問題解答例

Page 43

演習問題3

回答例

C*********************************************************************
C（省略）
C ------------------!
use portlib
use omp_lib
IMPLICIT REAL(4)(a-h,o-z)real(8) :: t1,t2
C
C
PARAMETER (mimax=513,mjmax=257,mkmax=257)
PARAMETER (mimax=257,mjmax=129,mkmax=129)
C
PARAMETER (mimax=129,mjmax=65,mkmax=65)
C
PARAMETER (mimax=65,mjmax=33,mkmax=33)
CC
ttarget specifys the measuring period in sec
PARAMETER (ttarget=60.0)
CC Arreycommon /pres/ p(mimax,mjmax,mkmax)
common /mtrx/ a(mimax,mjmax,mkmax,4),+b(mimax,mjmax,mkmax,3),c(mimax,mjmax,mkmax,3)common /bound/ bnd(mimax,mjmax,mkmax)
common /work/ wrk1(mimax,mjmax,mkmax),wrk2(mimax,mjmax,mkmax)
CC Other constantscommon /others/ imax,jmax,kmax,omegaC!
dimension time0(2),time1(2)!

Page 44integer ts, te, tr

演習問題3

回答例

C
omega=0.8
imax=mimax-1
jmax=mjmax-1
kmax=mkmax-1
CC Initializing matrixes
call initmt
write(*,*) ' mimax=',mimax,' mjmax=',mjmax,' mkmax=',mkmax
write(*,*) ' imax=',imax,' jmax=',jmax,' kmax=',kmax
CC Start measuringC
nn=10000
write(*,*) ' Start rehearsal measurement process.'
write(*,*) ' Measure the performance in 10000 times.'C!
cpu0=dtime(time0)
t1=omp_get_wtime()C
C Jacobi iteration
call jacobi(nn,gosa)
C!
cpu1= dtime(time1)
t2=omp_get_wtime()!
cpu = cpu1
cpu=t2-t1
flop=real(kmax-2)*real(jmax-2)*real(imax-2)*34.0*real(nn)xmflops2=flop/cpu*1.0e-6
write(*,*) ' MFLOPS:',xmflops2!
write(*,'(a10,f10.3)') ' time(s):',(te-ts)/dble(tr)write(*,'(a10,f10.3)') ' time(s):',cpu
write(*,*) 'gosa:',gosa
Page 45

演習問題3

回答例

CC
end the test loop!
nn=ifix(ttarget/(cpu/3.0))!
write(*,*) 'Now, start the actual measurement process.'!
write(*,*) 'The loop will be excuted in',nn,' times.'!
write(*,*) 'This will take about one minute.'!
write(*,*) 'Wait for a while.'C
C Jacobi iteration!
cpu0=dtime(time0)!
call jacobi(nn,gosa)
C!
cpu1= dtime(time1)!
cpu = cpu1!
flop=real(kmax-2)*real(jmax-2)*real(imax-2)*34.0*real(nn)!
xmflops2=flop*1.0e-6/cpu
CCCC
xmflops2=nflop/cpu*1.0e-6*float(nn)C!
write(*,*) ' Loop executed for ',nn,' times'!
write(*,*) ' Gosa :',gosa!
write(*,*) ' MFLOPS:',xmflops2, ' time(s):',cpu!
score=xmflops2/82.84!
write(*,*) ' Score based on Pentium III 600MHz :',score
C
!
pause
stop
END

Page 46

演習問題3

回答例

C
C
C**************************************************************
subroutine initmt
C**************************************************************IMPLICIT REAL(4) (a-h,o-z)
C
C
PARAMETER (mimax=513,mjmax=257,mkmax=257)
PARAMETER (mimax=257,mjmax=129,mkmax=129)
C
PARAMETER (mimax=129,mjmax=65,mkmax=65)
C
PARAMETER (mimax=65,mjmax=33,mkmax=33)
C

CC Arreycommon /pres/ p(mimax,mjmax,mkmax)
common /mtrx/ a(mimax,mjmax,mkmax,4),+b(mimax,mjmax,mkmax,3),c(mimax,mjmax,mkmax,3)common /bound/ bnd(mimax,mjmax,mkmax)
common /work/ wrk1(mimax,mjmax,mkmax),wrk2(mimax,mjmax,mkmax)
CC other constantscommon /others/ imax,jmax,kmax,omega
C!$OMP parallel!$OMP dodo k=1,mkmax
do j=1,mjmax
do i=1,mimaxa(i,j,k,1)=0.0a(i,j,k,2)=0.0
Page 47

演習問題3

回答例

a(i,j,k,3)=0.0a(i,j,k,4)=0.0b(i,j,k,1)=0.0b(i,j,k,2)=0.0b(i,j,k,3)=0.0c(i,j,k,1)=0.0c(i,j,k,2)=0.0c(i,j,k,3)=0.0
p(i,j,k) =0.0
wrk1(i,j,k)=0.0
bnd(i,j,k)=0.0
enddo
enddo
enddo
!$OMP end do
C!$OMP dodo k=1,kmax
do j=1,jmax
do i=1,imax
a(i,j,k,1)=1.0
a(i,j,k,2)=1.0
a(i,j,k,3)=1.0
a(i,j,k,4)=1.0/6.0
b(i,j,k,1)=0.0b(i,j,k,2)=0.0b(i,j,k,3)=0.0c(i,j,k,1)=1.0c(i,j,k,2)=1.0

Page 48

演習問題3

回答例

c(i,j,k,3)=1.0
p(i,j,k) =float((k-1)*(k-1))/float((kmax-1)*(kmax-1))
wrk1(i,j,k)=0.0
bnd(i,j,k)=1.0
enddo
enddo
enddo
!$OMP end do!$OMP end parallelC
returnend
C
C*************************************************************
subroutine jacobi(nn,gosa)
C*************************************************************IMPLICIT REAL(4)(a-h,o-z)
C
C
PARAMETER (mimax=513,mjmax=257,mkmax=257)
PARAMETER (mimax=257,mjmax=129,mkmax=129)
C
PARAMETER (mimax=129,mjmax=65,mkmax=65)
C
PARAMETER (mimax=65,mjmax=33,mkmax=33)
C
CC Arreycommon /pres/ p(mimax,mjmax,mkmax)
common /mtrx/ a(mimax,mjmax,mkmax,4),+b(mimax,mjmax,mkmax,3),c(mimax,mjmax,mkmax,3)common /bound/ bnd(mimax,mjmax,mkmax)
common /work/ wrk1(mimax,mjmax,mkmax),wrk2(mimax,mjmax,mkmax)
Page 49

演習問題3

回答例

CC other constantscommon /others/ imax,jmax,kmax,omega
CC
DO loop=1,nn
gosa=0.0!$OMP parallel!$OMP do private(S0,SS) reduction(+:GOSA)
DO K=2,kmax-1
DO J=2,jmax-1
DO I=2,imax-1
S0=a(I,J,K,1)*p(I+1,J,K)+a(I,J,K,2)*p(I,J+1,K)
1+a(I,J,K,3)*p(I,J,K+1)2+b(I,J,K,1)*(p(I+1,J+1,K)-p(I+1,J-1,K)
3-p(I-1,J+1,K)+p(I-1,J-1,K))
4+b(I,J,K,2)*(p(I,J+1,K+1)-p(I,J-1,K+1)5-p(I,J+1,K-1)+p(I,J-1,K-1))
6+b(I,J,K,3)*(p(I+1,J,K+1)-p(I-1,J,K+1)
7-p(I+1,J,K-1)+p(I-1,J,K-1))
8+c(I,J,K,1)*p(I-1,J,K)+c(I,J,K,2)*p(I,J-1,K)
9+c(I,J,K,3)*p(I,J,K-1)+wrk1(I,J,K)
SS=(S0*a(I,J,K,4)-p(I,J,K))*bnd(I,J,K)
GOSA=GOSA+SS*SS
wrk2(I,J,K)=p(I,J,K)+OMEGA *SS
enddo
enddo
enddo
!$OMP end do

Page 50

演習問題3

回答例

C!$OMP do
DO K=2,kmax-1
DO J=2,jmax-1
DO I=2,imax-1
p(I,J,K)=wrk2(I,J,K)
enddo
enddo
enddo
!$OMP end do!$OMP end parallelC
enddo
CC End of iteration
returnend

Page 51