1

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

2

画像処理用の GPUを、多数の演算器で構成された並列計算機として使う技術が注目されています。そのための並列プログラミ
ング手法について、記述が簡単な
OpenACCという並列プログラミングモデ
ルを中心に、紹介します。一部、OpenMPや MPIの知識を前提とした内容を含みます
が、知らない方でも理解できます。

3• 第一部（本資料を用いる）– ITOの説明– 実際にログインしてみる、PGIコンパイラをインストールする– GPUと並列計算についての説明• 第二部（別資料を用いる）– OpenACCの説明と実習– ～17:30：質疑・自由演習• 配付したITOのアカウントは1週間有効です– 受講後にも色々と試してみてください4• ITOへのログイン（別紙）• PGIコンパイラ（Community Edition）の導入– OpenACCのために必要なソフトウェアを入手する• インストール済のPGIコンパイラは古い、利用者に制限がある– ログインノードにSSHログインし、以下のコマンドでWebアクセス• elinks pgroup.com– pgilinux-2018-1810-x86-64.tar.gz を入手• エラーメッセージが出て終了してしまうが恐らくダウンロード自体はできている、md5sum で正しくDLできたかを確認• 正しいチェックサムは 9d2254112147d634e5fa422defd9e66b– インストール先を作成：mkdir opt– 作業場所を作成：mkdir tmp– 入手したファイルを展開：
tar zxvf ./pgilinux-2018-1810-x86-64.tar.gz –C tmp– 実際のインストール作業：cd tmp → ./install5•./install•

問題が起きたらCtrl+Cでキャンセルし、~/optを消して（名前を変えて）やり直す– Press enter to continue... enter– ライセンス情報を読む→ Do you accept these terms? (accept,decline) accept– Please choose install option: 1– Installation directory? [/opt/pgi] ~/opt/pgi– ライセンス情報を読む → accept– Do you wish to update/create links in the 2018 directory? (y/n) y– Press enter to continue... enter– Do you want to install Open MPI onto you system? (y/n) y– Do you want to enable NVIDIA GPU support in Open MPI? (y/n) y– Do you wish to obtain permanent license key or configure license service?
(y/n) n
– Do you want the files in the install directory to be read-only? (y/n) n• y にしても構わないのだが、あとで更新・追加する際に引っかかる可能性あり

6• ファイルの確認– ls ~/opt/pgi以下にlinux86-64, linux86-64-llvm, modulefilesディレク
トリとlicense関係のファイルが存在することを確認• modulefileの使用– module load ~/opt/pgi/modulefiles/pgi/18.10– pgccなどのコマンドが利用可能となることを確認– pgaccelinfoコマンドによりOpenACC対応GPUの情報などが得られる、が、ログインノードで実行してもCUDAドライバのバージョン情報しか
得られない

7

• 以下の内容を適当なファイルに保存する（仮にtest.shとする）#!/bin/bash#PJM -L "rscunit=ito-b"#PJM -L "rscgrp=ito-b-lecture-4"#PJM -L "vnode=1"#PJM -L "vnode-core=9"•#PJM -L "elapse=00:05:00"#PJM -o out.txt•#PJM -j#PJM -S

コピペする場合は記号が全角に化けるなどしないか気を付けること書かれている内容については後で説明する

module load ~/opt/pgi/modulefiles/pgi/18.10
pgaccelinfo• pjsub test.shを実行する• しばし待つ（バックエンドでジョブが実行される）• out.txtファイルが作られたら中身を確認する– （作られてすぐは中身が書き込まれていないことあり、もう少し待つ）

8• この建物の地下1Fに設置されているスーパーコンピュータ

バックエンド サブシステムA

バックエンド サブシステムB

フロントエンド サブシステム• スーパーコンピュータ＝「パソコンのようなもの」が大量に並べてつなげてある、ようなイメージ– どの程度「パソコンのようなもの」であるかはシステム次第– ITOは「普通のパソコンに近い」部分が多い• 普通のパソコンで使っているものより「ちょっとスゴイ」ものを多数並べて全体として「とてもスゴイ」性能を得ている

共有ストレージ

ノード間高速ネットワーク

9より詳細な情報はWebを参照
https://www.cc.kyushu-u.ac.jp/scp/system/ITO/01_intro.html• 主に計算に使う部分：バックエンドとフロントエンド– バックエンド サブシステムA• バックエンド1台の性能はPC• 1CPUあたり18コア×3GHzと比べて「ちょっとスゴイ」• （2CPU + 192GBメモリ）×2000• 多数並べると全体として「とてもスゴイ」• 合計 6.91 PFLOPS– バックエンド サブシステムB• 1CPUあたり18コア×2.3GHz
今回の利用対象• （2CPU + 384GBメモリ + 4GPU）×128• 合計 3.05 PFLOPS– フロントエンド• 1CPUあたり18コア×2.3GHz、（2CPU+384GBメモリ+1GPU）×160• 1CPUあたり22コア×2.2GHz、（16CPU+12TBメモリ+1GPU）×4
 1PFLOPS＝ 1秒間に1000兆回の実数計算ができる• 2018年のハイエンドAndroid端末は500GFLOPS程度• 1000=1K, 1000K=1M, 1000M=1G, 1000G=T, 1000T=1P

10• 画像処理用のハードウェア– CPUやマザーボードに組み込まれたチップ、または拡張スロットに搭載するビデオカードとして普及– 本来の役割：高速・高解像度描画、3D描画処理（透視変換、陰影・照明）、画面出力3次元画像描画の手順① (2, 2)② (8, 3)③ (5, 7)
オブジェクト単位、頂点単位、ピクセル単位で同時（並列）処理が可能– 同時にできることが多いため、それにあわせたハードウェアへと進化– 数値シミュレーション（科学技術計算）を高速に行えるハードウェアと両立することがわかり、GPUの重要な市場の一つとなってきた– 現在では特にビッグデータ、機械学習、AI処理などで性能を発揮

11• GPUはCPUと比べて単純な処理を（順序を問わず）たくさん行うことが求められるため、それに適したハードウェアへと進化してきた– HW構成バランスのイメージ
マルチコアCPU

メニーコアプロセッサ

GPU

計算ユニット
メモリ
キャッシュなど
制御部など– GPUの特徴：たくさんの計算ユニット、高速なメモリ、OSレス• OSを動かし様々な仕事をせねばならないCPUとは大きく異なる• 単体で利用できない、CPUによるサポートが必要• 現代の計算加速装置（アクセラレータ）の代表格• CPUよりスゴイ、というよりも、得意とする仕事が違う12※名古屋大版

Tesla P100

Xeon Gold 6140
（サブシステムB）

Xeon Gold 6154
（サブシステムA）

SPARC64 XIfx（「京」後継機）

アーキテクチャ
名PascalSkylake-SP

SPARC64 IXfx

コア数3584 FP32/1792 FP64
(56 SMs)

18（HTにより36スレッド実行可能だが、HTは無効化している）

32（ + 2アシスタン
トコア）

動作周波数

1.328 GHz – 1.480
GHz

2.3 GHz - 3.7 GHz

3.0 GHz – 3.7 GHz

2.2 GHz

搭載メモリ

HBM2
16 GB

DDR4
192 GB/socket

DDR4
96 GB/socket

HMC
32GB

HPL
ピボット選択付
きLU分解

4 TFLOPS程度※NVIDIA提供バイナ
リ実測値

0.9 TFLOPS程度※MKL実測値

1.1 TFLOPS程度※MKL実測値

0.97 TFLOPS※名古屋大提供値

STREAM Triad
A[i]=B[i]+C[i]550 GB/s

95 GB/s

95 GB/s

210 GB/s↑非常に高い性能を持つ、使わなければもったいない！

13• 56のStreaming Multiprocessor (SM)• SMあたり32のFP64コアと64のFP32コア、キャッシュなど• コアのスケジューリングはWARP単位（SIMD長32のようなもの）

14• TOP500 1~10位中、5システムがNVIDIAのGPUを活用– 上位50中の13、上位100中の24、全500中の128
Total Cores

ACC Cores

Rmax[TFlop/s]

Rpeak[TFlop/s]

Name

Site

Country

1Summit

ORNL

US

23978242196480

143500200794.88Volta GV100

2Sierra

LLNL

US

1572480

1382400

94640

125712Volta GV100

3Sunway TaihuLight Wuxi

China

10649600

4Tianhe-2A

Guangzhou

China4981760455475261444.5

100678.66 Matrix-2000

5Piz Daint

CSCS

Switzerland

38787231942421230

27154.3Tesla P100

6Trinity

LANL/SNL

US

97907293014.594

ACC

125435.9 None20158.7

41461.15 None

7ABCI
8SuperMUC-NG

AIST
LeibnizJapan
Germany

39168030585634816019880
19476.6

Tesla V100
32576.635 SXM226873.856 None

9Titan

ORNL

US56064026163217590

27112.55Tesla K20x

10Sequoia

LLNL

US

1572864

17173.224

20132.659 None

15

• 何故GPUを活用する必要があるのか？
– ITOにGPUが搭載されているから……というのも事実ではあるが– GPUは非常に高い性能をもつハードウェアであり、うまく活用できれば大変強力な研究の道具となる：成果を得るための加速装置– GPUは国内外の様々な計算環境に導入されているため、利用スキルを得ておくことはきっとプラスになる• どうやって使う？
– 並列化されたソフトウェア（アプリケーション）を使う– （並列化されていないプログラムから）並列化されたライブラリやフレームワークを使う– 自分で並列計算を実装する• 自分が扱いたい任意の問題（アプリ）をGPU化できる• 扱いたい問題が既にGPU化されているならそれを使っても良いが、GPUのことを理解しているかどうかでさらに高い性能が得られるかど
うかが決まる、こともある

16• GPU（などのアクセラレータ）向けのプログラムを簡単に書くことのできる並列化プログラミング言語– NVIDIAのGPUに限定されない• C/C++やFortranで書かれたプログラムに簡単な「指示文」を追加するだけでGPU対応が可能– コンパイラ向けのコメントを記述する– CPU向けの並列化に広く用いられているOpenMPを踏襲、理解も利用も容易– 最適化を行う上ではある程度GPUの知識があった方が良い• 「どんなプログラムでもGPU化できる」「どんなプログラムでも高速化できる」わけではないが、多くのプログラムに対して有用– 自作コードのGPU化を行う上で最も簡単でコストパフォーマンスの高い方法（だと思う）– コストパフォーマンスが高い＝労力に対して十分な性能が得られる

17• 規格などの情報– https://www.openacc.org/• 主な対応コンパイラ– 商用：PGI（無償版あり）、Cray、国家超级计算无锡中心– 研究：Omni、OpenARC、OpenUH、ROSEACC– OSS：GCC• プロファイラやデバッガなど– allinea MAP/DDT、PGI pgprof、NVIDIA nvprof/nvvp• 日本語で書かれた資料– ソフテック社（PGIコンパイラの代理店）の資料• OpenACC ディレクティブによるプログラミング by PGI Compilers
https://www.softek.co.jp/SPG/Pgi/OpenACC/

18• 並列計算とは何か？並列プログラミングとは何か？
– 並列計算：何らかの計算処理を同時並行的に行うこと• 並列：同じ処理を同時に行う、ある処理を幾つかのサブ処理に分けて同時並行的に行う• 平行：何らかの処理を同時並行的に行う– 並列プログラミング（並列化プログラミング）：並列化・並列計算を行
うためのプログラミング• なぜ並列計算を行うのか？
– 短時間で終わらせたい計算があるから、計算機（CPUやパソコンなど計算を行うHW）の持つ能力をフル活用したいから– 現代のプロセッサは複数の計算コアを搭載– PC用CPUもスマートフォン向けCPUも4～8程度のコア数が主流– 大規模なスーパーコンピュータの総コア数は100万以上– → HWの性能を十分に引き出すには並列計算が必須

19• 元々の逐次計算– 単純な繰り返しループ計算
for(i=0; i<N; i++){A[i] = B[i] + C[i];D[i] = E[i] + F[i];}• ループ内の処理を分割し、同
時に計算
PE1

for(i=0; i<N; i++){A[i] = B[i] + C[i];}PE2for(i=0; i<N; i++){D[i] = E[i] + F[i];}PE: Processor Element、CPUやGPUなど計算をする「何か」• ループそのものを分割し、同
時に計算
PE1

for(i=0; i<N/2; i++){A[i] = B[i] + C[i];D[i] = E[i] + F[i];}PE2for(i=N/2; i<N; i++){A[i] = B[i] + C[i];D[i] = E[i] + F[i];}• 本講習会で扱うOpenACCによる並列化はこちらのイメージ– OpenMPもこちら

20• 単純なループによる並列化とは異なり、やや大きい粒度の並列計算などでよく活用• OpenMPにおいて近年サポートが活発• GPU（OpenACC）にはあまり適していないため今回は扱わない

21• 並列化できないプログラムも少なくない– 計算順序に制約がある場合は困難• 工夫により可能となることもある

for(i=1; i<N; i++){A[i] = A[i-1] + B[i];C[i] = A[i] + D[i];}• （プログラム高速化全般に言えることだが）速くした部分しか速くならない
逐次処理A

一部並列 処理A

処理B
B

time
頑張って処理Bを並列高速化したが、全体の実行時間的には……

– 探索問題などでは分割数（並列度）よりずっと速くなることもある★解• 探索対象★
• 逐次探索★
• 並列探索– プログラムのどこをGPUに担当させるかをちゃんと考える必要がある22• 多数のコア全てを満たすに十分な仕事を与える– 数千のコア全てをフル稼働させる並列度が必要– コア数の数倍以上の仕事があることが望ましい• メモリアクセスを待つ間に別の仕事ができる• WARPを意識する– 内部的には32演算器単位で動作しており、分岐の際などに注意が必要• Voltaから変わってしまった• 連続メモリアクセスについて考える– コアレスなメモリアクセスが高速• 詳細は終盤（OpenACCプログラムの最適化）にて解説する– 基本的には、並列度が高くて分岐の少ない単純なプログラムが良い23• 並列計算環境が普及した今日では様々な「並列化のための道具」が使われている– バラバラだと提供側・利用側ともに不便なため共通化されることが多い• ある環境向けに書いたものが別の環境でも利用できる（互換性）• 性能まで互換性があるとは限らない点には注意が必要（性能可搬性）– ある環境では有効であった最適化が別の環境では性能低下要因になることも– 環境が変わってもその環境ごとに常に最大の性能が得られるようにするための研究も行われている→自動チューニング

• 自動化はできないのか？
– 全く不可能なわけではない、ある程度はコンパイラ等が行ってくれる– 「コンパイラ様のご機嫌を伺う」コードを書く必要がある、コンパイラ
によって差が大きい、簡単なコードでないとうまくいかないことが多い• 単純な行列積などコードの形状と最適化のパターンが決まっているものは人が書くよりコンパイラやライブラリの方が得意– コード記述者が明示することで適切な最適化を行うのが現在のトレンド

24• MPI：プロセス間の通信規格、特に複数ノード利用時に必須– MPI_Send, MPI_Recv, MPI_Gatherなどの関数を使う• pthread：スレッド並列処理のための関数群– pthread_create, pthread_……
– 現在ではアプリケーションコードで直接利用することはあまりない• OpenMP：スレッド並列処理、主にループ並列化向けの指示文規格– #pragma omp parallel for 、!$omp parallel do– 最近の規格ではタスク処理やGPU対応などを拡充– 現状では利用できる環境が少ないが、将来的にはOpenACCを吸収？

• CUDA：NVIDIA社のGPU向け、GPUを普及させた最大要因の一つ、
NVIDIA GPUの能力をフル活用したければ必須• OpenCL：「汎用版CUDA」のようなもの、FPGAなどでも利用可能• コンパイラによる自動並列化：SIMD並列化など一部の処理で有用

25• 単純な配列コピーの例
GPU上で行われる処理__global__ void gpukernel(int N, float* C, float* A){（GPUカーネル）
int tid = blockIdx.x*blockDim.x + threadIdx.x;
int nt = gridDimd.x * blockDim.x;
for(int i=tid; i<N; i+=nt){• 同時にたくさん実行される
C[i] = A[i];• 自分のIDを元に計算するべき範囲を特定する}}
int main(int argc, char **argv){
int N = 100000;
float *A, *C;
CPU上で行われる処理
float *d_A, *d_C;
A = (float*)malloc(sizeof(float)*N);C = (float*)malloc(sizeof(float)*N);• 専用の関数などを用いて様々な処理を行う
cudaMalloc((void**)&d_A, sizeof(float)*N);
cudaMalloc((void**)&d_C, sizeof(float)*N);
cudaMemcpy(d_A, A, sizeof(float)*N, cudaMemcpyHostToDevice);
cudaMemcpy(d_C, C, sizeof(float)*N, cudaMemcpyHostToDevice);gpukernel<<<4,4>>>(N, d_C, d_A);
cudaMemcpy(C, d_C, sizeof(float)*N, cudaMemcpyDeviceToHost);
cudaFree(d_A);
return 0;
慣れてしまえばそれほど難しいものではないが、記述量がそ}れなりにあり「カジュアルなGPU利用」にはあまり向かない26• GPUの性能を活用できる様々なライブラリが公開されている– https://developer.nvidia.com/gpu-accelerated-libraries– Deep Learningライブラリ• cuDNN, TensorRT, DeepStream SDK– 数値計算・数学ライブラリ• cuBLAS, cuSPARSE, cuRAND, cuFFT, ...
– 通信、C++クラス、etc.
• NCCL, Thrust, OpenCV, MAGMA, ...
– 行いたい処理が対応している場合は容易にGPUの性能を活用することができる1
サンプルソースコードはITOの /home/tmp/gpu 以下に置いてあります※実質的に「演習の答え」となるものもあるので注意

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

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

– 基本的には単体GPUを駆動するための手段、マルチGPU・マルチノード
についてはMPIなどと組み合わせて利用する3• 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化が可能

4• 指定した部分だけが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上のメモリを解放5• 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※最適でない（問題が起きる）こともある
Accelerator kernel generated
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
Accelerator kernel generated
Generating Tesla code
14, !$acc loop gang, vector(32) ! blockidx%x threadidx%x6• ./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 5.00 6.00 7.00 8.00 9.00 10.00
8.00 10.00 12.00 14.00 16.00 18.00 20.00• LD_LIBRARY_PATHなどの対応は必要– module loadにより解決可能– （envコマンドにより環境変数をみれば確認可能）– GPU向けにコンパイルしたものはGPUがないと実行時エラー

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

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

ITOポータル

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

バックエンド計算結果

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

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

9

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

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

10

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

ジョブ ID11#!/bin/bash#PJM -L "rscunit=ito-b"#PJM -L "rscgrp=ito-b-lecture-4"#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分で打ち切り標準出力とエラー出力を統合
ノード利用情報を出力

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

12• サンプルプログラム（vector0.c または vector.f90）をコンパイル
し、バッチジョブとして実行してみる• 実行結果が確認できたら、次は環境変数PGI_ACC_TIMEをセットして実行してみる– export PGI_ACC_TIME=1 をプログラム実行よりも前の行に書く– GPUがどのような仕事をしたのかが確認できるようになる• 例 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

13• もっと細かく確認したい場合には環境変数 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カーネル起動情報とデータ転送情報が出
力される

14

• 指示文（directive）により全てを記述する– 指示文：コンパイラに対して指示を行う特殊なコメント• C/C++：#pragma acc ~• Fortran：!$acc ~• 基本的には「無視してしまっても問題が起きない文」– コンパイラが対応していない場合もコンパイルと実行自体は可能• 具体的な指示文の例– 並列計算の方法を指示するもの• 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

15• 「この範囲内を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を用いる

16• 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)

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

17

• 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
Accelerator kernel generated
Generating Tesla code
17, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */アクセラレータ(GPU)向けのカーネルが生成された
Tesla向けのコードが生成されたFortran
ループの並列化が行われた
pgfortran -Minfo=accel -acc
-ta=tesla:cc60 以下省略
main:
13, Generating implicit copyout(v2(:))
Generating implicit copyin(v1(:))
14, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
14, !$acc loop gang, vector(32) ! blockidx%x threadidx%x• どのように判断・処理されたのかを確認することは非常に重要

18• 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)

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

• Fortran (vector1.f90)

1 program main2
implicit none
3
integer :: i, n=10
4
double precision :: v1(10), v2(10)56do i=1, n
do i=1, n7v1(i) = dble(i)v1(i) = dble(i)8v2(i) = 0.0d0v2(i) = 0.0d0
9
enddo
enddo
10
11 !$acc kernels copyout(v2(:)) copyin(v1(:))!$acc kernels
12 !$acc loop gang, vector(32)do i=1, n13
do i=1, nv2(i) = v1(i) * 2.0d0
14v2(i) = v1(i) * 2.0d0
enddo※end loopは不要
15
enddo!$acc end kernels（書いてもエラーには
16 !$acc end kernels
ならないようだ）
17
do i=1,n18
do i=1,nwrite(*,'(1H F8.2)',advance="NO")v1(i)19write(*,'(1H F8.2)',advance="NO")v1(i)
enddo
20
enddo
write(*,*)""
21
write(*,*)""
do i=1,n22
do i=1,nwrite(*,'(1H F8.2)',advance="NO")v2(i)23write(*,'(1H F8.2)',advance="NO")v2(i)enddoコンパイラの判断により、「copyout」「copyin」
24
write(*,*)""
enddo
「loop
gang,
vector(32)」が自動的に挿入されていたと思えば良い25
write(*,*)""
end program main
26 end program main
program main
implicit none
integer :: i, n=10
double precision :: v1(10), v2(10)20• 指示文行を次の行に継続させることも可能– 長くなってしまったときなどに• 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

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

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

デバイスメモリ（HBM2）

22• 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}23• 配列全体ではなく一部のみを送受信することも可能• 注意：C/C++とFortranでは部分配列の指定方法が異なる– C/C++：先頭と長さを指定する#pragma acc kernels copy(A[head:length])– Fortran：開始点と終了点を指定する!$acc kernels copy(A(begin:end))– :N のような省略表記も可能（先頭からN要素が送られる）

24• 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)

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

25• メインメモリやデバイスメモリの転送速度に対して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26• 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)– data指示文とは好みで使い分けてよい• プログラム全体の見通しが悪くならないように注意– GPU化範囲の前でとにかく全部送信したいとき？
– 複数ソースコードにプログラムが分割されているとき？

27• 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)※コンパイラが判断に失敗する場合はkernelsに
present節を加えると良い
acc kernels present(v1)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

デバイスメモリの生成繰り返し
計算
デバイスメモリの破棄

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

のみ行われ、無駄がない28• 多次元配列の送受信も可能– ただし、連続したメモリの範囲しか扱えない#pragma acc data copyin(A[head:length][0:N])– C/C++• C/C++は右側の次元が低次元、低次元分は全て転送する必要がある– Fortran !$acc data copyin(A(1:N, begin:end))• Fortranは左側の次元が低次元、低次元分は全て転送する必要がある29• 部分転送時の範囲には変数を用いても良い• 配列長については注意が必要– 動的に確保した配列などコンパイル時に配列の長さが分からないものは長さを明示しておく必要がある• 範囲は変数による指定で良い、間接参照をしているときなどに注意
double *v1, *v2;
int *index;
v1 = (double*)malloc(sizeof(double)*n));v2 = (double*)malloc(sizeof(double)*n));index = (int*)malloc8sizeof(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, n
具体例は
v2(i) = v1(index(i)) * 2.0d0 vector12.c/f90
enddo
vector13.c/f90!$acc end kernels
を参照

 v1の範囲がうまく認識できず、コンパイルはできるが実行時エラー
 copyin(v1[n]) および copyin(v1(n)) を加えると正しく動作する• Deep copyはできない（可能となりつつある）– C/C++：動的に確保された配列をメンバとして持つ構造体やクラス– Fortran：allocatable属性やpointer属性を持つメンバを含む派生型

30• サンプルプログラム（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• 余裕がある人は、配列の長さや本数を増やすなどして比較してみる31• 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)

32• デバイスメモリの生成・破棄を伴わないデータ転送（更新）には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

33

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

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

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

36• 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 配列cの依存関係に関するメッAccelerator kernel generated
セージは出ているが、3重ループ
Generating Tesla code
38, !$acc loop gang ! blockidx%yの外側2つが並列化された
39, !$acc loop gang, vector(128) ! blockidx%x threadidx%x
40, !$acc loop seq

37• 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
Accelerator kernel generated
 依存関係があり並列化できない、逐次コード
Generating Tesla codeが生成された旨が出力されている
34, #pragma acc loop seq
 正しく実行はできるが、逐次実行のため低速
35, #pragma acc loop seq
36, #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

38• 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指示文で指定すればコンパイラの判断を上書きできる39• loop independent指定なし
37, Complex loop carried dependence of a->->,c->->,b->-> prevents parallelization
Accelerator serial kernel generated
Accelerator kernel generated
Generating Tesla code
34, #pragma acc loop seq
35, #pragma acc loop seq
36, #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• 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
Accelerator kernel generated
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

40• 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

38, #pragma acc loop gang, vector(128) collapse(2) /* blockIdx.x threadIdx.x */39,/* blockIdx.x threadIdx.x collapsed */41, #pragma acc loop seq
39, !$acc loop gang, vector(128) collapse(2) ! blockidx%x threadidx%x
40,! blockidx%x threadidx%x collapsed
42, !$acc loop seq• collapseを指定すると、ループを融合してから並列化する• 短いループがネストしている際に有用• 実行時間的にはそれぞれのループを並列化した場合と変わらないことが多い• 並列実行形状（後述）を細かく指定する場合などにうまく使い分けると良い41• 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指定を変更するとどうか？

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

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化すれば良いだろうか？

★

★

★

★

★

★

44

45

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

46• 正しく計算できているはずだが、調べて見るとデータ転送が多い– 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

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

47• 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節を用いて途中の配列データを覗いてみよう• デバッグの際などに役立つことがある

48• 実は最初に全てのデータを転送し、結果ベクトルのみ回収すれば良かった– もちろん、わかっていればいきなりデータ通信の最適化を行っても良い– 実際に実行時間レベルで意味がある（GPUを使うことでCPUより大幅な高速化が得られる）のは大規模な行列の場合• コードを書き換えたり実行時に引数を与えれば問題サイズを変更できるようになっている– データを出力している部分をなんとかしないと（コメントアウトなど）大量のデータが書き出されてしまうため注意• 問題が簡単過ぎてすぐに反復計算が終わってしまうため、もっと長時間実行したい場合は行列をいじったり収束条件を厳しくしたりする必
要がある

49

• 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などの並列実行形状の指定（後述）も必要

50• ここまで、ループをどのように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

51• 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…

52Streaming 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以降で厳密化

53• 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);
}
}

54

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