1

サンプルソースコードはITOの /home/tmp/gpu 以下に置いてあります※実質的に「演習の答え」となるものもあるので注意※PGI compiler 19.4で実行した場合の出力に準拠
18.10でも19.4でも基本的には同じであるが、18.10では出力されていた
Accelerator kernel generated（並列実行カーネルが作成できた旨）が出力されなくなったことを反映

九州大学情報基盤研究開発センター
2019.06.04 13:30-17:30
2019.06.13 13:30-17:30

2• 単純なベクトル計算を題材としてOpenACCの基本を学ぶ– コンパイルの仕方、実行の仕方– CPU-GPU間のデータ転送について• 行列積を題材としてOpenACCの基本を学ぶ– 並列化ループの指定方法について• CG法を題材としてOpenACCの基本を学ぶ– 少し複雑なコードのOpenACC化について• その他、OpenACCに関する話題– 最適化のための一般的なヒントなど

講習会の時間設定的に最後まで到達できません。
自主学習などにご利用ください。

3• GPUプログラムを簡単に記述するために開発された指示文規格– GPU向けのOpenMP、のようなもの• マルチGPU・マルチノードについてはMPIなどと組み合わせて利用– 幾つかの会社が独自に開発していたものが共通規格として集約された• 初登場が2011年、まだ10年経っていない– CUDAでしか書けない処理も多いが、「とにかく高い並列度で一気に計算すれば良い」という典型的なGPU向けプログラムでは十分高性能– CUDAを使うべきプログラム• GPU上の高速共有メモリやシャッフル命令を意識したアルゴリズム• インスタンスIDを意識したアルゴリズム– CUDAはThreadIDなど「ID」を意識して並列処理を記述する– OpenACCは「ID」の概念そのものがない• その他、最新のハードウェア機能をフル活用したい場合– Tensor core、RT core、半精度演算– C(C/C++)とFortranに対応

4

計算ノード

CPU

メインメモリ（DDRメモリ）GPU

デバイスメモリ（GDDR*、HBM2）• CPUとGPUが搭載された計算ノード• CPUとGPUはそれぞれ個別のメモリを持ち、相手側のメモリに直接は触れない（データ転送が必要、ある程度の時間がかかる）

5• C (vector0.c)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

#include <stdio.h>#include <stdlib.h>#include <string.h>
int main(int argc, char **argv){
int i, n=10;
double v1[10], v2[10];for(i=0; i<n; i++){v1[i] = (double)(i+1);v2[i] = 0.0;}
#pragma acc kernels
for(i=0; i<n; i++){v2[i] = v1[i] * 2.0;}
for(i=0;i<n;i++){ printf(" %.2f", v1[i]); }printf("¥n");for(i=0;i<n;i++){ printf(" %.2f", v2[i]); }printf("¥n");return 0;• Fortran (vector0.f90)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

program main
implicit none
integer :: i, n=10
double precision :: v1(10), v2(10)do i=1, nv1(i) = dble(i)v2(i) = 0.0d0
enddo

OpenACC並列化対象＝GPU上で実行される!$acc kernels
do i=1, nv2(i) = v1(i) * 2.0d0
enddo!$acc end kernels

do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)enddowrite(*,*)""
do i=1,nwrite(*,'(1H F8.2)',advance="NO")v2(i)enddowrite(*,*)""
end program main}

➢ 単純なプログラムであればkernels指示文で対象を指定するだけでGPU化が可能

6• 指定した部分だけがGPU上で実行される
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

#include <stdio.h>#include <stdlib.h>#include <string.h>
int main(int argc, char **argv){
int i, n=10;
double v1[10], v2[10];for(i=0; i<n; i++){v1[i] = (double)(i+1);v2[i] = 0.0;}

• 全体としてCPUが「主」、GPUが「従」の関係• 指定されていない部分はCPU上で実行される
C
P
U
G
P
U

#pragma acc kernels
for(i=0; i<n; i++){v2[i] = v1[i] * 2.0;}
for(i=0;i<n;i++){ printf(" %.2f", v1[i]); }printf("¥n");Cfor(i=0;i<n;i++){ printf(" %.2f", v2[i]); }P
printf("¥n");return 0;}

U• CPUからGPUに対して計算指示が行われる• GPU上のメモリを確保、CPUからGPUへ必要なデータを転送
GPU上の計算コアにより並列計算される• CPUはGPUの計算終了を待つ• GPUからCPUへ結果データを転送、GPU上のメモリを解放7• pgccまたはpgfortranでコンパイルする– -acc OpenACC指示文を有効化– -ta OpenACCの対象ハードウェア（対象GPUの種類）を指定– -Minfo=accel OpenACC化に関する情報を出力– -tp 対象CPUを指定
pgcc -Minfo=accel -acc -ta=tesla:cc60 -tp=skylake -o vector0_c_acc vector0.c
main:
16, Generating implicit copyout(v2[:]) データ転送やGPU上の計算コアの使い方はGenerating implicit copyin(v1[:])
コンパイラが適当に判断してくれた
17, Loop is parallelizable※最適でない（問題が起きる）こともある
Generating Tesla code
17, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */pgfortran -Minfo=accel -acc -ta=tesla:cc60 -tp=skylake -cpp -o vector0_f_acc vector0.f90
main:
13, Generating implicit copyout(v2(:))
Generating implicit copyin(v1(:))
14, Loop is parallelizable※出力情報の具体的な読み方は後述する
Generating Tesla code
14, !$acc loop gang, vector(32) ! blockidx%x threadidx%x8• ./a.out– CPU向けの実行可能ファイルと同様にそのまま実行できる
C版の場合$ ./a.out
1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00
2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 20.00
Fortran版の場合$ ./a.out
1.00 2.00
2.00 4.00

3.00
6.00

4.00
8.00

5.00
10.00

6.00
12.00

7.00
14.00

8.00
16.00

9.00
18.00

10.00
20.00

• LD_LIBRARY_PATHなどの対応は必要– module loadにより解決（envコマンドにより環境変数をみれば確認可能）– GPU向けにコンパイルしたものはGPUが使えないと実行時エラー• ITOで実行するにはバッチジョブの使い方を知る必要がある

9• ITOのように多人数で利用するシステムでは、計算ノードに直接ログインするのではなく「バッチジョブシステム」を利用すると便利– 計算機の稼動率を高める– 特定ユーザに占有させない– セキュリティを高める• 利用手順

ログインノード
フロントエンド

ITOポータル

1. ジョブスクリプトの作成
2. ジョブの投入
3. システムがスケジューリング
して実行
4. 結果の確認

バックエンド

計算結果

10• 処理してほしい内容を記述したファイル（ジョブスクリプト）を作成しシステムに投入– ジョブとして受付• 資源の空き状況に応じて順に処理される– 要求内容や空き状況によっては先を越される– 適切な資源量の指定が重要
リソースグループ（待ち行列）
ジョブの使用資源量（CPU,メモリ）で選択

ジョブ

小規模ジョブ用

中規模ジョブ用

大規模ジョブ用#!/bin/sh#PJM -L "vnode=4"#PJM -L "vnode-core=36"#PJM -L "rscunit=ito-b"#PJM -L "rscgrp=ito-g-16"#PJM -L "elapse=10:00"
mpiexec -np 16 ./a.out

ジョブスクリプトの例
計算機群

11

• バッチジョブの投入– pjsub

ジョブスクリプトファイル名• バッチジョブの状況確認– pjstat• 終了済のジョブの情報を見る場合はpjstat -H• 特定ジョブの詳細情報をみる場合はpjstat -S ジョブID• -Hと-Sの組み合わせも可能（終了済の特定ジョブの詳細を確認）• バッチジョブのキャンセル– pjdelジョブID– 各コマンドに追加できるオプションは--helpオプションを付けることで
確認可能

12

• コマンド$ pjsub オプション

ジョブスクリプトファイル名– オプション：使用する資源等に関する指定• スクリプトファイル内に書いても良いし、オプションで与えても良い• オプションで与えた情報が優先される– ジョブスクリプトファイル：依頼する処理内容• シェルスクリプトとして記述– 例）ジョブスクリプトファイル test.sh を投入$ pjsub test.sh[INFO] PJM 0000 pjsub Job 28246 submitted.

ジョブ ID13#!/bin/bash#PJM -L "rscunit=ito-b"#PJM -L "rscgrp=ito-b-lecture"#PJM -L "vnode=1"#PJM -L "vnode-core=9"#PJM -L "elapse=00:05:00"#PJM -j#PJM -S（スクリプト記述にbashを使用）
サブシステムBを利用
講習会用のリソースグループを指定
1GPU使えれば良いので1/4ノード
対応するコア数は95分で打ち切り標準出力とエラー出力を統合
ノード利用情報を出力

modulefileを用いて環境設定
module load ~/opt/pgi/modulefiles/pgi/19.4
OpenACC対応デバイス情報の確認
pgaccelinfo
プログラムの実行./a.out「動作確認」で書いていた#PJM -o out.txt
は、標準出力の内容を指定のファイルへ書き出すというオプション。
エラー出力についても同様に-eで指定可能。
-o –e -j 全て指定すると-oで指定したファイルに標準出力とエラー出力の両方が書き出される。

14• サンプルプログラム（vector0.c または vector.f90）をコンパイル
し、バッチジョブとして実行してみる• 実行結果が確認できたら、次は環境変数PGI_ACC_TIMEをセットして実行してみる– export PGI_ACC_TIME=1 をプログラム実行よりも前の行に書く– GPUがどのような仕事をしたのかが確認できるようになる#!/bin/bash#PJM 省略
module load ~/opt/pgi/modulefiles/pgi/19.4
pgaccelinfo
export PGI_ACC_TIME=1./a.out

15• export PGI_ACC_TIME=1 によって出力される情報の例
Accelerator Kernel Timing data/home/usr0/m70000a/work/gitprojects/testprograms/openacc/vector0/vector0.c
main NVIDIA devicenum=0time(us): 36
16: compute region reached 1 time
←GPU上での計算
17: kernel launched 1 time
grid: [1] block: [32]device time(us): total=4 max=4 min=4 avg=4
elapsed time(us): total=459 max=459 min=459 avg=459
16: data region reached 2 times←CPU-GPU間の通信
16: data copyin transfers: 1
device time(us): total=13 max=13 min=13 avg=13
21: data copyout transfers: 1
device time(us): total=19 max=19 min=19 avg=19– 最適化やデバッグに活用できることがある• 何回実行されたか？
• どれだけ時間がかかったか？
• GPUカーネルの実行形状は？

16• もっと細かく確認したい場合には環境変数 PGI_ACC_NOTIFY
も有効– 1,2,4,8のビット組み合わせで指定、以下の情報が出力される• 1: GPUカーネル起動• 2: データ転送• 4: regionのentry/exit• 8: wait/sync– 例：export PGI_ACC_NOTIFY=3• 3=1と2の論理和、GPUカーネル起動情報とデータ転送情報が出
力される

17

• 指示文（directive）により全てを記述する– 指示文：コンパイラに対して指示を行う特殊なコメント• C/C++：#pragma acc ~• Fortran：!$acc ~• 基本的には「無視してしまっても問題が起きない文」– コンパイラが対応していない場合もコンパイルと実行自体は可能、もちろんGPUは使えない• 具体的な指示文の例– 並列計算の方法を指示するもの• kernels, parallel• loop, seq, collapse• gang/num_gangs, worker/num_workers, vector/vector_length– データの移動について指示するもの• data, enter/exit data, copy{,in,out}, present, update, create, delete

18• 「この範囲内をGPU上で並列実行したい」ことを示す– kernelsとparallelではコンパイラによる解釈の仕方が異なる• parallel：基本的に利用者が細かく指定する• kernels：ある程度コンパイラが判断、手動で調整（上書き）可能• 最適化をしていくと結局同じようなコードになる、はずである• 範囲の途中で離脱するような構造は不可（forループのbreakなど）– 利用する指示節にも違いが生じる
kernelsと共に利用するもの• async / wait
• device_type• if• default• copy系

parallelと共に利用するもの• async / wait
• device_type• if• default• copy系• num_gangs / num_workers / vector_length• reduction• private• どちらを用いても良いが、本講習会ではkernelsを用いる

19• C (vector0.c)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

• Fortran (vector0.f90)#include <stdio.h>#include <stdlib.h>#include <string.h>

1 program main2
implicit none
3
integer :: i, n=10
4
double precision :: v1(10), v2(10)int main(int argc, char **argv)
5{6
do i=1, n
int i, n=10;7v1(i) = dble(i)double v1[10], v2[10];8v2(i) = 0.0d0
9
enddo
for(i=0; i<n; i++){10v1[i] = (double)(i+1);
11 !$acc kernels
v2[i] = 0.0;12
do i=1, n}
13v2(i) = v1(i) * 2.0d0
14
enddo#pragma acc kernels
15 !$acc end kernels
for(i=0; i<n; i++){16v2[i] = v1[i] * 2.0;17
do i=1,n}18
write(*,'(1H F8.2)',advance="NO")v1(i)19
enddo
for(i=0;i<n;i++){ printf(" %.2f", v1[i]); }20
write(*,*)""
printf("¥n");21
do i=1,nfor(i=0;i<n;i++){ printf(" %.2f", v2[i]); }22
write(*,'(1H F8.2)',advance="NO")v2(i)（OpenMPと同様に）printf("¥n");
23
enddo
24
write(*,*)""➢
C/C++では指示文直後のループや{} で括った部分（構造化ブロック）がreturn 0;
25 end program main
指示文の適用対象となる}
26➢ Fortranではendで閉じる必要がある

20

• C

コンパイラの判断によって

pgcc -Minfo=accel -acc -ta=tesla:cc60 -tp=skylake 以下省略
copyout, copyinという命令が生成された
main:
16, Generating implicit copyout(v2[:])
Generating implicit copyin(v1[:])17, Loop is parallelizable※implicit：暗黙的な（プログラムには書かれていなかったが、コンパイラが判断しました）
Generating Tesla code
17, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */• Fortran

アクセラレータ(GPU)向けのカーネルが生成された
Tesla(NVIDIA GPU)向けのコードが生成されたループの並列化が行われた

pgfortran -Minfo=accel -acc -ta=tesla:cc60 以下省略
main:
13, Generating implicit copyout(v2(:))
Generating implicit copyin(v1(:))
14, Loop is parallelizable
Generating Tesla code
14, !$acc loop gang, vector(32) ! blockidx%x threadidx%x• どのように判断・処理されたのかを確認することは非常に重要※blockidx や threadidx については後述21• C (vector0.c)
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

• C (vector1.c)

int main(int argc, char **argv){
int i, n=10;
double v1[10], v2[10];4 int main(int argc, char **argv)
5 {6
int i, n=10;
7
double v1[10], v2[10];8for(i=0; i<n; i++){9for(i=0; i<n; i++){v1[i] = (double)(i+1);10v1[i] = (double)(i+1);v2[i] = 0.0;11v2[i] = 0.0;
}
12
}
13
#pragma acc kernels
14 #pragma acc kernels copyout(v2[:]) copyin(v1[:])for(i=0; i<n; i++){15 #pragma acc loop gang, vector(32)v2[i] = v1[i] * 2.0;
16
for(i=0; i<n; i++){}17v2[i] = v1[i] * 2.0;
18}
for(i=0;i<n;i++){ printf(" %.2f", v1[i]); }19
printf("¥n");20
for(i=0;i<n;i++){ printf(" %.2f", v1[i]); }for(i=0;i<n;i++){ printf(" %.2f", v2[i]); }21
printf("¥n");printf("¥n");22
for(i=0;i<n;i++){ printf(" %.2f", v2[i]); }23
printf("¥n");return 0;
24}
25
return 0;
26 }コンパイラの判断により、「copyout」「copyin」「loop gang, vector(32)」が自動的に挿入されていたと思えば良い22• Fortran (vector0.f90)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

program main
implicit none
integer :: i, n=10
double precision :: v1(10), v2(10)• Fortran (vector1.f90)

1 program main2
implicit none
3
integer :: i, n=10
4
double precision :: v1(10), v2(10)5do i=1, n6
do i=1, nv1(i) = dble(i)7v1(i) = dble(i)v2(i) = 0.0d08v2(i) = 0.0d0
enddo
9
enddo
10
!$acc kernels
11 !$acc kernels copyout(v2(:)) copyin(v1(:))
do i=1, n
12 !$acc loop gang, vector(32)v2(i) = v1(i) * 2.0d013
do i=1, n
enddo
14
v2(i) = v1(i) * 2.0d0※end loopは不要!$acc end kernels
15
enddo（書いてもエラーには
16 !$acc end kernels
ならないようだ）
do i=1,n17
write(*,'(1H F8.2)',advance="NO")v1(i)18do i=1,nenddo
19
write(*,'(1H F8.2)',advance="NO")v1(i)write(*,*)""
20
enddo
do i=1,n21
write(*,*)""
write(*,'(1H F8.2)',advance="NO")v2(i)22do i=1,n
enddo
23
write(*,'(1H F8.2)',advance="NO")v2(i)コンパイラの判断により、「copyout」「copyin」
write(*,*)""
24
enddo「loop gang, vector(32)」が自動的に挿入されていたと思えば良い
end program main25
write(*,*)""
26 end program main

23

• 指示文行を次の行に継続させることも可能– 長くなってしまったときなどに• C/C++とFortranで少し違うので注意#pragma acc kernels copyout(v2[:]) copyin(v1[:])#pragma acc loop gang, vector(32)for(i=0; i<n; i++){v2[i] = v1[i] * 2.0;}!$acc kernels copyout(v2(:)) copyin(v1(:))!$acc loop gang, vector(32)do i=1, nv2(i) = v1(i) * 2.0d0
enddo!$acc end kernels#pragma acc kernels ¥copyout(v2[:]) copyin(v1[:])#pragma acc loop gang, vector(32)for(i=0; i<n; i++){v2[i] = v1[i] * 2.0;}!$acc kernels &!$acc copyout(v2(:)) copyin(v1(:))!$acc loop gang, vector(32)do i=1, nv2(i) = v1(i) * 2.0d0
enddo!$acc end kernels

24• CPUとGPUは個別のメモリを持っており、直接相手側のメモリにアクセスできない– 例外あり、詳しくは後述• 適切なデータ送受信を行わねば正しい計算が行えない– GPUカーネル起動時：メインメモリからデバイスメモリへのデータ転送– GPUカーネル終了後：デバイスメモリからメインメモリへのデータ転送• 単純なプログラムでは自動的にデータ転送を行ってくれるが、ある程度複雑な場合には明示する必要がある– 特にC/C++ではコンパイラが長さを認識できない配列を扱うことが多いため注意が必要– GPUカーネルが生成されなかったり、実行時にエラーしたりする原因となるCPU

メインメモリ（DDRメモリ）GPU

デバイスメモリ（HBM2）

25• kernels指示文に追加して配列のデータ転送を明示する– GPUカーネル実行前に、デバイスメモリを確保し、ホストからデバイス
へコピーする• copyin– GPUカーネル終了後に、デバイスからホストへ書き戻し、デバイスメモ
リを破棄する※OpenACC2.5からは常に• copyout「present_or_*」の挙動となり、存在していれば使い回してくれる。
– copyin + copyout
実際にどう扱われるかはコンパイル• copy時のメッセージやPGI_ACC_NOTIFY– デバイスメモリを確保するのみを用いて確認すること。
• create– 既にデバイスメモリに存在していることをコンパイラに伝える※データを使い回す、該当するものが無ければ実行時エラー• present– 存在していない場合のみcopy{,in,out}する• present_or_copy{,in,out}26• 配列全体ではなく一部のみを送受信することも可能• 注意：C/C++とFortranでは部分配列の指定方法が異なる– C/C++：先頭と長さを指定する#pragma acc kernels copy(A[head:length])– Fortran：開始点と終了点を指定する!$acc kernels copy(A(begin:end))– A[:N], A(:N) のような省略表記も可能（先頭からN要素が送られる）

27• GPUカーネルを何度も実行する場合はどうなるだろうか？
int main(int argc, char **argv){
int i, j, n=10;
double v1[10];for(i=0; i<n; i++){v1[i] = (double)(i+1);}

for(j=0; j<10; j++){#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}}

for(i=0; i<n; i++){printf(" %.2f", v1[i]);}printf("¥n");return 0;}

program main
implicit none
integer :: i, j, n=10
double precision :: v1(10)do i=1, nv1(i) = dble(i)enddo
do j=1, 10!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels
enddo

do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)enddowrite(*,*)""
デバイスメモリの生成とend program main

(vector2.c)

繰り返し
デバイスメモリの生成
計算
デバイスメモリの破棄(vector2.f90)

破棄を繰り返してしまう（余計な時間がかかる）

28• メインメモリやデバイスメモリの転送速度に対してCPU-GPU間のデータ転送速度はずっと低速、頻繁な通信は避けたい– CPU - メインメモリ• ITOの場合はDDR4で100GB/s/socket弱（STREAM Triad実測）– GPU - デバイスメモリ• HBM2、550GB/s程度（STREAM Triad実測）– 後継機のV100では800GB/sを超える– CPU - GPU• PCI Express Gen.3 x16~200 GB/s
メインメモリ（DDRメモリ）CPU– 理論性能でも最大16GB/s（×双方向）– GPU - GPU• NVLink、20GB/s*1or2（×双方向）– 1GPUあたり4本のNVLink– Power系CPUではCPU-GPU間でもNVLinkが使える

PCIe、~32GB/s
デバイスメモリ（HBM2）GPU~1 TB/s29• data指示文– ループ並列化のタイミング以外で、データのみを操作できる#pragma acc data copyin(A) copyout(B)
構造化ブロック!$acc data copyin(A) copyout(B))
構造化ブロック!$acc end data

• enter/exit data指示文– 構造化ブロックを囲まずに自由な位置で送受信を行うことも可能#pragma acc enter data copyin(A)!$acc enter data copyin(A)#pragma acc exit data copyout(B)!$acc exit data copyout(A)– どちらを使っても良い• enter/exit data指示文の方が便利だが、プログラムの見通しが悪くならないように注意が必要– GPU化範囲の前でとにかく全部送信したいとき？
– 複数ソースコードにプログラムが分割されているとき？

30• GPUカーネルを何度も実行する場合などにdata指示文が有効
int main(int argc, char **argv){
int i, j, n=10;
double v1[10];for(i=0; i<n; i++){v1[i] = (double)(i+1);}#pragma acc data copy(v1[:])for(j=0; j<10; j++){#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}}

for(i=0; i<n; i++){printf(" %.2f", v1[i]);}printf("¥n");return 0;}(vector3.c)

program main
implicit none
integer :: i, j, n=10
double precision :: v1(10)do i=1, nv1(i) = dble(i)enddo!$acc data copy(v1(:))
do j=1, 10!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels
enddo!$acc end data

※コンパイラが判断に失敗する場合はkernelsに
present節を加えると良い
acc kernels present(v1)（v1のアクセスに間接参照がある場合などに効果的）
デバイスメモリの生成繰り返し
計算
デバイスメモリの破棄

do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)enddo
生成と破棄は最初と最後にwrite(*,*)""(vector3.f90)
end program main

のみ行われ、無駄がない31• 多次元配列の送受信も可能– ただし、連続したメモリの範囲しか扱えない• 低次元分は全て転送する必要がある#pragma acc data copyin(A[head:length][0:N])– C/C++• C/C++は右側の次元が低次元– Fortran !$acc data copyin(A(1:N, begin:end))• Fortranは左側の次元が低次元

32• 部分転送時の範囲には変数を用いても良い• 配列長については注意が必要– 動的に確保した配列などコンパイル時に配列の長さが分からないものは長さを明示しておく必要がある• 範囲は変数による指定で良い、間接参照をしているときなどに注意
double *v1, *v2;
int *index;
v1 = (double*)malloc(sizeof(double)*n));v2 = (double*)malloc(sizeof(double)*n));index = (int*)malloc(sizeof(int)*n);#pragma acc kernels
for(i=0; i<n; i++){v2[i] = v1[index[i]] * 2.0;}

double precision, allocatable :: v1(:), v2(:)
integer,allocatable :: index(:)
allocate(v1(n), v2(n), index(n))!$acc kernels
do i=1, nv2(i) = v1(index(i)) * 2.0d0
enddo!$acc end kernels➢ v1の範囲（長さ）がうまく認識できず、コンパイルはできるが実行時エラー➢ copyin(v1[n]) および copyin(v1(n)) を加えると正しく動作する➢ 具体例は vector12.c/f90 および vector13.c/f90 を参照（コンパイル時に出力される情報を比較してみよう）

33• 完全なDeep copyができない問題– 動的な要素を持つ集合的な要素をまとめて転送できない問題– C/C++：動的に確保された配列をメンバとして持つ構造体やクラスをDeep copyできない– Fortran：allocatable属性やpointer属性を持つメンバを含む派生型をDeep copyできない• 解決方法– 必要な分だけ手動でcopyする– コンパイル時の -ta オプションに deepcopy を追加• PGI Fortranのみ対応、完全なDeep copyができる（はず）– Unified memoryを使う• メインメモリとデバイスメモリを同一に扱う技術• コンパイル時の -ta オプションに managed を追加、色々と制限があるため注意が必要

34• サンプルプログラム（vector2.c, vector2.f90）をコンパイルし、
バッチジョブとして実行してみる• さらに、data指示文を挿入したプログラム（vector3.c,
vector3.f90）も実行し、比較してみる– 環境変数PGI_ACC_TIMEをセットして実行すると容易に比較が可能
Accelerator Kernel Timing data/home/usr0/m70000a/work/gitprojects/testprograms/openacc/lecture201811/vector2.cmain NVIDIA devicenum=0time(us): 522
回数を示す値が変化すること
14: compute region reached 10 times
15: kernel launched 10 timesに気がつくはずである
grid: [1] block: [32]device time(us): total=22 max=4 min=2 avg=2
elapsed time(us): total=734 max=532 min=21 avg=73
14: data region reached 20 times
14: data copyin transfers: 10
device time(us): total=93 max=34 min=6 avg=9
18: data copyout transfers: 10
device time(us): total=407 max=351 min=6 avg=40• 余裕がある人は、配列の長さや本数を増やすなどして比較してみる

35• data指示文内のGPUカーネル間でデータの確認や更新を行いたい
場合はどうすれば良いだろうか？
#pragma acc data copy(v1){#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}!$acc data copy(v1)!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels

並列化ループの途中で値を確認したい・更新したい
たとえばここでv1の値を出力したら何が見えるのだろうか？
#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}}

計算前の値が見える!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels!$acc end data

data範囲のあとであれば計算結果が全てメインメモリに書き戻されているのだが……？
for(i=0; i<n; i++){printf(" %.2f", v1[i]);}printf("¥n");(vector20.c)

do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)enddowrite(*,*)""(vector20.f90)

36• デバイスメモリの生成・破棄を伴わないデータ転送（更新）にはupdateを用いる– ホストからデバイス：update device– デバイスからホスト：update host または update self
– 一部のみの更新も可能、範囲の指定方法はdata指示文と同様#pragma acc data copy(v1){#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}!$acc data copy(v1)!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo!$acc end kernels#pragma acc update host(v1)for(i=0; i<n; i++){printf(" %.2f", v1[i]);}printf("¥n");!$acc update host(v1)do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)enddowrite(*,*)""#pragma acc kernels
for(i=0; i<n; i++){v1[i] = v1[i] * 2.0;}
(vector21.c)}!$acc kernels
do i=1, nv1(i) = v1(i) * 2.0d0
enddo(vector21.f90)!$acc end kernels!$acc end data

37

• スカラ変数はfirstprivateとなる– 並列化範囲外の値を引き継ぎ、互いに干渉しあわない• 配列はデバイスメモリにて共有される– 互いに干渉する– private指示節により変更することも可能• 参考：OpenMPの場合– private/shared指示節で指定– 何も指定しないとsharedとなりスレッド間で干渉する– Fortranのみ並列化範囲内の逐次ループのループカウンタはprivate扱い

38• 単純なベクトル計算を題材としてOpenACCの基本を学ぶ– コンパイルの仕方、実行の仕方– CPU-GPU間のデータ転送について• 行列積を題材としてOpenACCの基本を学ぶ– 並列化ループの指定方法について• CG法を題材としてOpenACCの基本を学ぶ– 少し複雑なコードのOpenACC化について• その他、OpenACCに関する話題– 最適化のための一般的なヒントなど

39• 現実のプログラムではコンパイラが全てのループの並列化の判断を行うのは難しい、プログラマが並列化の判断をする必要がある• loop指示文– 並列化対象ループを指定する– さらに以下の指示節と組み合わせることで動作の制御が可能➢ independent指示節とseq指示節• 対象ループを並列実行するか逐次実行するかを明示する• 強制力があり、コンパイラによる判断は行われなくなる➢ collapse(n)指示節：collapse(2), collapse(3)など• 多重ループをまとめて並列化する• 並列度の低い階層ループの並列化などに極めて重要➢ reduction指示節：reduction(+:a）など• 計算結果の集約などを行う• 多くの場合はコンパイラが正しく判断してくれるため書く必要はない40• C (matmul0.c)• Fortran (matmul0.f90)

double **a=NULL, **b=NULL, **c=NULL;// mallocでa,b,cを確保#pragma acc kernels
for(i=0; i<n; i++){for(j=0; j<n; j++){for(k=0; k<n; k++){c[i][j] += a[i][k] * b[k][j];}}}n=10で実行してみたC
37: kernel launched 1 time

double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n), b(n,n), c(n,n)!$acc kernels
do i=1, n
do j=1, n
do k=1, nc(j,i) = c(j,i) + a(k,i) * b(j,k)
enddo
enddo
enddo
!$acc end kernels

grid: [1] block: [1]device time(us): total=108 max=108 min=108 avg=108
elapsed time(us): total=130 max=130 min=130 avg=130

Fortran

40: kernel launched 1 time
grid: [1x10] block: [128]device time(us): total=5 max=5 min=5 avg=5
elapsed time(us): total=418 max=418 min=418 avg=418

どうやら実行時間に大きな差があるようだ、何故だろう？

41• Fortran (matmul0.f90)
double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n), b(n,n), c(n,n)!$acc kernels
do i=1, n
do j=1, n
do k=1, nc(j,i) = c(j,i) + a(k,i) * b(j,k)
enddo
enddo
enddo
!$acc end kernels

pgfortran -Minfo=accel -acc -ta=tesla:cc60 ¥-tp=skylake -o m0f_f_acc matmul0.f90
main:
37, Generating implicit copyin(a(1:n,1:n))Generating implicit copy(c(1:n,1:n))Generating implicit copyin(b(1:n,1:n))38, Loop is parallelizable➢ copy関係はコンパイラの判断で特に問題はない
39, Loop is parallelizable
40, Complex loop carried dependence of c prevents parallelization
Loop carried dependence of c prevents parallelization
Loop carried backward dependence of c prevents vectorization
Inner sequential loop scheduled on accelerator
Generating Tesla code
38, !$acc loop gang ! blockidx%y
39, !$acc loop gang, vector(128) ! blockidx%x threadidx%x
40, !$acc loop seq
sequential
ループが逐次実行されることを意味する➢ 配列cの依存関係に関するメッセージは出ているが、3重ループの外側2つが並列化された42• C (matmul0.c)
double **a=NULL, **b=NULL, **c=NULL;// mallocでa,b,cを確保#pragma acc kernels
for(i=0; i<n; i++){for(j=0; j<n; j++){for(k=0; k<n; k++){c[i][j] += a[i][k] * b[k][j];}}pgcc -Minfo=accel -acc -ta=tesla:cc60 -tp=skylake -o m0c_c_acc matmul0.c}
main:

36, Generating implicit copyin(b[:n][:n])➢ copy関係はコンパイラの判断でGenerating implicit copy(c[:n][:n])特に問題はない
Generating implicit copyin(a[:n][:n])37, Complex loop carried dependence of a->->,c->->,b->-> prevents parallelization
Accelerator serial kernel generated
Generating Tesla code
➢ 依存関係があり並列化できず、逐次実行コー
37, #pragma acc loop seq
ドが生成された旨が出力されている
38, #pragma acc loop seq➢ 正しく実行はできるが、逐次実行のため低速
39, #pragma acc loop seq
38, Complex loop carried dependence of a->->,c->->,b->-> prevents parallelization
39, Complex loop carried dependence of a->->,c->->,b->-> prevents parallelization
Loop carried dependence due to exposed use of c[i1][i2] prevents parallelization

43• C (matmul1.c)
double **a=NULL, **b=NULL, **c=NULL;// mallocでa,b,cを確保#pragma acc kernels#pragma acc loop independentfor(i=0; i<n; i++){#pragma acc loop independentfor(j=0; j<n; j++){#pragma acc loop seq
for(k=0; k<n; k++){c[i][j] += a[i][k] * b[k][j];}}}• Fortran (matmul1.f90)
double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n), b(n,n), c(n,n))!$acc kernels!$acc loop independent
do i=1, n!$acc loop independent
do j=1, n!$acc loop seq
do k=1, nc(j,i) = c(j,i) + a(k,i) * b(j,k)
enddo
enddo
enddo
!$acc end kernels• 一般的に、Fortranプログラムの方がコンパイラによる並列化判断が適切に働く• C/C++はポインタ参照の都合で不具合が起きないよう保守的な判断がされやすい• loop指示文で指定すればコンパイラの判断を上書きできる44• C言語、loop independent指定なし
37, Complex loop carried dependence of a->->,c->->,b->-> prevents parallelization
Accelerator serial kernel generated
Accelerator kernel generated
Generating Tesla code
37, #pragma acc loop seq
38, #pragma acc loop seq
39, #pragma acc loop seq
38, Complex loop carried dependence of a->->,c->->,b->-> prevents parallelization
39, Complex loop carried dependence of a->->,c->->,b->-> prevents parallelization
Loop carried dependence due to exposed use of c[i1][i2] prevents parallelization• C言語、loop independent指定あり
38, Loop is parallelizable
40, Loop is parallelizable
42, Complex loop carried dependence of a->->,c->->,b->-> prevents parallelization
Loop carried dependence of c->-> prevents parallelization
Loop carried backward dependence of c->-> prevents vectorization
Generating Tesla code
38, #pragma acc loop gang /* blockIdx.y */40, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */42, #pragma acc loop seq

45• C (matmul2.c)
double **a=NULL, **b=NULL, **c=NULL;// mallocでa,b,cを確保#pragma acc kernels#pragma acc loop independent collapse(2)for(i=0; i<n; i++){for(j=0; j<n; j++){#pragma acc loop seq
for(k=0; k<n; k++){c[i][j] += a[i][k] * b[k][j];}}}• Fortran (matmul2.f90)
double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n), b(n,n), c(n,n))!$acc kernels!$acc loop independent collapse(2)do i=1, n
do j=1, n!$acc loop seq
do k=1, nc(j,i) = c(j,i) + a(k,i) * b(j,k)
enddo
enddo
enddo
!$acc end kernels

C38, #pragma acc loop gang, vector(128) collapse(2) /* blockIdx.x threadIdx.x */39, /* blockIdx.x threadIdx.x collapsed */41, #pragma acc loop seq
Fortran 39, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
40, ! blockidx%x threadidx%x collapsed
42, !$acc loop seq• collapseを指定すると、ループを融合してから並列化する• 短いループがネストしている際に有用• 実行時間的にはそれぞれのループを並列化した場合と変わらないことが多い• 並列実行形状（後述）を細かく指定する場合などにうまく使い分けると良い46• C (matmul1.c)
double **a=NULL, **b=NULL, **c=NULL;// mallocでa,b,cを確保#pragma acc kernels#pragma acc loop independentfor(i=0; i<n; i++){#pragma acc loop independentfor(j=0; j<n; j++){#pragma acc loop seq
for(k=0; k<n; k++){c[i][j] += a[i][k] * b[k][j];}}}• Fortran (matmul1.f90)
double precision, allocatable :: a(:,:), b(:,:), c(:,:)
allocate(a(n,n), b(n,n), c(n,n))!$acc kernels!$acc loop independent
do i=1, n!$acc loop independent
do j=1, n!$acc loop seq
do k=1, nc(j,i) = c(j,i) + a(k,i) * b(j,k)
enddo
enddo
enddo
!$acc end kernels

実際にコンパイルして実行してみよう➢ 並列化されたか？
➢ 実行時間は短くなったか？
➢ independentやseqを変更してみるとどうか？
➢ collapse指定を変更するとどうか？

47• ここまで、ループをどのようにGPU上の計算コアに割り当てるかは「おまかせ」だった• 簡単なプログラムでは特に問題ないことが多いが、明示的に調整したい場合もある– ネストしたループ（多重ループ）はどのように計算コアに割り当たっている？
– WARP長(32)にあわせたループ構造にしたが、コンパイラはそれに合わせた実行をしてくれているのか？

• 実はコンパイル時のメッセージにどのように割り当てるかが出力されていた– gang, vector と blockIdx, threadIdx という概念が存在するようだ
pgcc -Minfo=accel -acc -ta=tesla:cc60 -tp=skylake -o vector0_c_acc vector0.c
17, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */pgfortran -Minfo=accel -acc -ta=tesla:cc60 -tp=skylake -o vector0_f_acc vector0.f90
14, !$acc loop gang, vector(32) ! blockidx%x threadidx%x48• OpenACCではハードウェアに階層的な並列性があることを想定しており、上位階層から順にgang, worker, vectorとなっている• これらをどのように組み合わせて実行するかを指定できる• CUDAの並列実行モデルが階層的になっているため、それにあわせて言語設計された、という方が正しい– grid, threadblock, threadの三階層構造– OpenACCのおおまかな並列実行モデルの対応付け
gangvectorblockId 0
threadId 0,1,2…

blockId 1
threadId 0,1,2…

blockId 2
threadId 0,1,2…

49Streaming Multiprocessor(SM)内の並列性はvector、
SM単位の並列性はgang• 基本的には• 連続メモリアクセスする最内側のループはvector• より外側のループはgangくらいのイメージ• HWの制約上実際には32コア単位で動作していることを覚えておくと良い性能が出ることもある• ()で数字を与えた場合はその数単位で割り当てられる•

gang : block
vector : thread

もう少し詳細に言えば……
• gangはHWレベルで同期できない単位の粗粒度並列性（CUDAではSM単位）• workerはHWレベルで同期できる単位の細粒度並列性（CUDAではSM内のWARP群）• vectorはworker内部でのSIMDやベクトル並
列処理（CUDAではWARP）※外側のループほど上位（gang側）でなければならない：OpenACC 2.0以降で厳密化• とりあえずgangとvectorを意識しておくと高い性能が出やすくなることもあるかも知れない、程度に考えておこう50• 単純なベクトル計算を題材としてOpenACCの基本を学ぶ– コンパイルの仕方、実行の仕方– CPU-GPU間のデータ転送について• 行列積を題材としてOpenACCの基本を学ぶ– 並列化ループの指定方法について• CG法を題材としてOpenACCの基本を学ぶ– 少し複雑なコードのOpenACC化について• その他、OpenACCに関する話題– 最適化のための一般的なヒントなど

51• 単純なCG法の計算カーネル（反復計算部）を題材として、プログ
ラムのOpenACC化を考えてみる– 簡単にするため、行列は密行列、前処理は対角スケーリング• CG法（共役勾配法、Conjugate Gradient Method）– 対称正定値行列を係数とする連立一次方程式Ax=bを解く手法• 行列A、既知のベクトルb、未知のベクトルx– 基本アルゴリズム• Wikipediaから引用、前処理なし• 実際のコードは計算順序が変更されている版• 単純な行列Aとランダム行列xxからbを求めておき
Ax=bを解いてxとxxが（ほぼ）一致することを確認する、という構造にしてある• 時間測定を簡単に書くためOpenMP関数を利用– コンパイル時に-mpオプションを加える必要あり52

Compute r (0) = b-[A]x (0)for i= 1, 2, …
solve [M]z (i-1) = r (i-1)ρi-1 = r (i-1) z (i-1)
if i=1p (1) = z (0)else
βi-1 = ρi-1 /ρi-2
p (i) = z (i-1) + βi-1 z (i)endif
q (i) = [A]p (i)αi = ρi-1 /p (i) q (i)x (i) = x (i-1) + αi p (i)r (i) = r (i-1) - αi q (i)check convergence |r|
end

•
•
•
•

初期値、x(0)は適当な値（今回は0ベクトル）
収束するまで繰り返す
前処理、今回はrを対角要素で割るだけ
リダクション• コピー• リダクション• ベクトル積和•••••

行列ベクトル積
リダクション
ベクトル積和ベクトル積和収束判定（中身はリダクションと平方根）※上付き文字は反復回数に対応• Ax=b：Aは行列、xとbはベクトル• z,r,p,qはベクトル• α・β・ρはスカラー（ベクトルのリダクション結果）

for(iter=1; iter<=maxiter; iter++){printf("iter %d ", iter);// {z} = [Minv]{r}for(i=0;i<N;i++){z[i] = dd[i]*r[i]; }// {rho} = {r}{z}※対角要素分の1だけのベクトルddを用意済
rho = 0.0;
for(i=0;i<N;i++){rho += r[i]*z[i]; }// {p} = {z} if iter=1// beta = rho/rho1 otherwise
if(iter==1){
for(i=0;i<N;i++){p[i] = z[i];}}else{
beta = rho/rho1;
for(i=0;i<N;i++){p[i] = z[i] + beta*p[i];}}// {q} = [A]{p}for(i=0;i<N;i++){★行列ベクトル積
q[i] = 0.0;
行ごとの計算を並列に行える
for(j=0;j<N;j++){q[i] += A[i*N+j]*p[j];}}// alpha = rho / {p}{q}pq = 0.0;
for(i=0;i<N;i++){pq += p[i]*q[i]; }alpha = rho / pq;// {x} = {x} + alpha*{p}// {r} = {r} - alpha*{q}for(i=0;i<N;i++){x[i] += + alpha*p[i];r[i] += - alpha*q[i];}// check converged
dnrm = 0.0;
for(i=0;i<N;i++){dnrm += r[i]*r[i]; }resid = sqrt(dnrm/bnrm);
if(resid <= cond){break;}if(iter == maxiter){break;}rho1 = rho;}★

★

• 行列やベクトルに対する単純な計算ばかりで構成されているため並列化・OpenACC化は容易– ★：行列要素同士のコピーや四則演算– ★：集約演算– reduction、dot product– 複雑な前処理を適用する
場合は難易度が上がる– 具体的にはどのような手順でOpenACC化すれば良いだろうか？

★

★

★

★

★

53

54

• 並列化できることがわかっているループに指示文を挿入– サンプルコードcg2.cおよびcg2.f90に指示文を挿入する• Fortran版ではコンパイル時に-cppオプションを加える必要あり– 出力部の調整のために #if を使っているため– はじめはループ1つだけに指示文を挿入して実行してみよう• step by stepで少しずつ並列化できるのはOpenACCの強みの一つ• 注意点（コンパイラのメッセージも参考に）– 特にCの場合はindependent節も活用する必要あり• 主に、配列の参照先アドレスが重複する可能性を考慮するためか– （既に例示したように）配列長が認識できずにデータ転送に躓くことがあるので適切にcopy指示節を追加する• C言語版の行列Aに要注意• もちろん、in/outを意識してcopyin/copyoutとしても良い– reductionが適切に生成されているかを確認• （問題なく生成されると思って良いが、念のため）

55• C言語版
CPU向け（OpenMP並列化）
pgcc -tp=skylake -mp -o cg2c_c cg2.c
GPU向け（OpenACC並列化）
pgcc -Minfo=accel -acc -ta=tesla:cc60 -tp=skylake -mp -cpp -o cg2c_c_acc cg2.c• Fortran版
CPU向け（OpenMP並列化）
pgfortran -tp=skylake -mp -cpp -o cg2f_f cg2.f90
GPU向け（OpenACC並列化）
pgfortran -Minfo=accel -acc -ta=tesla:cc60 -tp=skylake -mp -cpp -o cg2f_f_acc cg2.f90• ジョブスクリプト#!/bin/bash#PJM 省略
module load ~/opt/pgi/modulefiles/pgi/19.4./cg2c_c./cg2c_c_acc• cg2.cとcg2.f90はOpenACC指示文を挿入していないソースコード••

kernelsとloopを適切に指定したのがcg3.cとcg3.f90
cg3.cではコンパイラが行列Aの大きさに迷わないようcopyinも追加している56

A:
9.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 9.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 9.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 9.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 9.000000 1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000 9.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 9.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 9.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 9.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 9.000000
x:
8.400000
8.700000
7.800000
1.600000
9.400000
3.600000
8.700000
9.300000
5.000000
2.200000b:131.900000
134.300000
127.100000
77.500000
139.900000
93.500000
134.300000
139.100000
104.700000
82.300000
iter 1 1.089293e-01 1.000000e-16
iter 2 2.331098e-16 1.000000e-16
iter 3 5.450800e-18 1.000000e-16
time: 0.000010 sec , 0.000003 sec/iter
result(x):8.400000
8.700000
7.800000
1.600000
9.400000
3.600000
8.700000
9.300000
5.000000
2.200000
1 8.400000 8.400000: 3.552714e-15
2 8.700000 8.700000: 5.329071e-15
3 7.800000 7.800000: -8.881784e-16
4 1.600000 1.600000: -8.881784e-16
5 9.400000 9.400000: 0.000000e+00
6 3.600000 3.600000: -1.776357e-15
7 8.700000 8.700000: 0.000000e+00
8 9.300000 9.300000: 0.000000e+00
9 5.000000 5.000000: -8.881784e-16
10 2.200000 2.200000: -1.776357e-15←既知の行列A←ベクトルX（求める答え）←既知のベクトルb←反復計算の履歴←計算結果ベクトルX←計算結果ベクトルXと最初に設定したベクトルXの比較

57• 正しく計算できているはずだが、調べて見るとデータ転送が多い– PGI_ACC_TIMEなどで確認できる• 現在のデータ転送状況イメージ
Compute r (0) = b-[A]x (0)for i= 1, 2, …
solve [M]z (i-1) = r (i-1)ρi-1 = r (i-1) z (i-1)
if i=1p (1) = z (0)else
βi-1 = ρi-1 /ρi-2
p (i) = z (i-1) + βi-1 z (i)endif
q (i) = [A]p (i)αi = ρi-1 /p (i) q (i)x (i) = x (i-1) + αi p (i)r (i) = r (i-1) - αi q (i)check convergence |r|
end

CPU

GPU

計算をする度に、関係する配列全てを送受信している↓
更新があったデータだけを送受信すれば良いはずでは？

58• data節を用いてデータ転送を削減してみる– 送受信が必要なデータはどれだろうか？
#pragma acc data copyin(？) copyout(？)!$acc data copyin(？) copyout(？)

Compute r (0) = b-[A]x (0)for i= 1, 2, …
solve [M]z (i-1) = r (i-1)ρi-1 = r (i-1) z (i-1)
if i=1p (1) = z (0)else
βi-1 = ρi-1 /ρi-2
p (i) = z (i-1) + βi-1 z (i)endif
q (i) = [A]p (i)αi = ρi-1 /p (i) q (i)x (i) = x (i-1) + αi p (i)r (i) = r (i-1) - αi q (i)check convergence |r|end※リダクション結果のス
カラ変数は自動的にCPU側に送られるため気にしなくてよい– 余裕がある人は、update節を用いて途中の配列データを覗いてみよう• デバッグの際などに役立つことがある59• 実は最初に全てのデータを転送し、結果ベクトルのみ回収すれば良かった– もちろん、わかっていればいきなりデータ通信の最適化を行っても良い– 実際に実行時間レベルで意味がある（GPUを使うことでCPUより大幅な高速化が得られる）のは大規模な行列の場合• コードを書き換えたり実行時に引数を与えれば問題サイズを変更できるようになっている– データを出力している部分をなんとかしないと（コメントアウトなど）大量のデータが書き出されてしまうため注意• 問題が簡単過ぎてすぐに反復計算が終わってしまうため、もっと長時間実行したい場合は行列をいじったり収束条件を厳しくしたりする必
要がある
•
•

cg4.cとcg4.f90はcg3.cとcg3.f90にdata指示文を加えて最初に全てのデータをGPUへ転送したもの
cg5.cとcg5.f90はさらにkernels指示文にpresent節を加えてデータ転送が不要であることを明示したもの

60• 単純なベクトル計算を題材としてOpenACCの基本を学ぶ– コンパイルの仕方、実行の仕方– CPU-GPU間のデータ転送について• 行列積を題材としてOpenACCの基本を学ぶ– 並列化ループの指定方法について• CG法を題材としてOpenACCの基本を学ぶ– 少し複雑なコードのOpenACC化について• その他、OpenACCに関する話題– 最適化のための一般的なヒントなど

61• GPUの特徴にあわせた適切な実装を行うことで性能が向上する• 特に以下の点は性能への影響が大きいため気を付けたい– 並列度• GPUは多数の計算コアによる並列計算によって高性能を得ているため
並列化対象ループに十分な長さがあるようにする– vector(threadIdx)は32以上、gang(blockIdx)はSM数以上– 短いループはcollapseで結合させるなどする– 連続メモリアクセス• vectorループでは連続メモリアクセスを心がける– GPUの得意なメモリアクセス方式となる
for(i=0;i<N;i++){for(j=0;j<N;j++){○ a[i][j] = b[i][j] + c[i][j];× a[j][i] = b[j][i] + c[j][i];}}do i=1, N
do j=1, N× a(i,j) = b(i,j) + c(i,j);○ a(j,i) = b(j,i) + c(j,i);
}
}

62

• kernels/parallel内部（OpenACC並列化対象内部＝GPU上）で関数を実行する場合はroutine指示文が必要// プロトタイプ宣言にも指示文が必要#pragma acc routine
void calc(int n, double *v);// 呼び出し元#pragma acc data copy(v1){#pragma acc kernels{
calc(n, v1);}}// 呼び出し対象の関数#pragma acc routine
void calc(int n, double *v){
int i;
for(i=0; i<n; i++){v[i] = v[i] * 2.0;}}

! 呼び出し対象の関数
module modcontains
subroutine calc(n,v)!$acc routine
integer :: ndouble precision :: v(*)!$acc loop
do i=1, nv(i) = v(i) * 2.0d0
enddo
end subroutine calc
end module mod! 呼び出し元!$acc data copy(v1)!$acc kernels
call calc(n,v1)!$acc end kernels!$acc end data

➢ Cでは関数名の前に、Fortranでは関数名の次に指示文を挿入➢ 関数内でも並列計算を行わせるにはroutineの後にさらにgangなどの並列実行形状の指定も必要

63• デバッガ、プロファイラ– PGIコンパイラ付属のpgdbgやpgprof、NVIDIA社の提供するnvvpやnvprof、その他サードパーティ製のソフトの幾つかが利用可能– GUIが表示されるものはX転送が必要• CUDAやMPIとの連携– OpenACCと外部とでデータ（ポインタ）をやりとりする方法を提供– host_data, use_device, deviptrなどを活用• （Deep copyの話で少しだけ触れたが）CUDA Unified Memory– CPUとGPUが連続したメモリアドレス空間を利用する技術– データ転送を記述しなくても必要に応じてCPU-GPU間のデータ転送が勝手に行われるため、プログラミングが容易になる– 幾つかの条件（制限）があり、性能にもペナルティが生じる– コンパイル時に-ta=tesla:managedオプションをつけるだけで良い