第3回 配信講義　計算科学技術特論A （2023）

内容に関する質問は [email protected] まで 第３回 OpenMPの基礎 名古屋大学情報基盤センター 片桐孝洋 1 2023年度 計算科学技術特論A

講義日程と内容について  ２０２３年度 計算科学技術特論Ａ（木曜：１３：００－１４：３０ ）      2 第１回：プログラム高速化の基 礎、2023年4月13日  イントロダクション、ループアンローリング、キャッシュブロック化、 数値計算ライブラリの利用、その他 第２回：MPIの基礎、2023年4月20日  並列処理の基礎、MPIインターフェース、MPI通信の種類、その他 第３回：OpenMPの基礎、2023年4月27日  OpenMPの基礎、利用方法、その他 第４回：Hybrid並列化技法(MPIとOpenMPの応用)、2023年5月11日  背景、Hybrid並列化の適用事例、利用上の注意、その他 第５回：プログラム高速化実例と大規模学習への展開、2023年5月18日  プログラムの性能ボトルネック に関する考えかた（I/O、単体性能 (演算機ネック、メモリネック)、並列性能(バランス))、性能プロファイル、 機械学習におけるHPC、ほか 2023年度 計算科学技術特論A

参考書  「計算科学のためのHPC技術1 」 下司雅章 (編集), 片桐孝洋 , 中田真秀, 渡辺宙志, 山 本有作, 吉井範行, Jaewoon Jung, 杉田 有治, 石村和 也, 大石進一, 関根晃太, 森倉悠介, 黒田久泰，著 出版社: 大阪大学出版会 (2017/4/3) ISBN-10: 4872595866, ISBN-13: 978-4872595864 発売日： 2017/4/3 【本書の特徴】  計算科学に必要なＨＰＣ技術について、基礎的な事 項を解説している  片桐担当（１章～５章）  プログラム高速化の基礎、MPIの基礎、OpenMP の基礎、Hybrid並列化技法(MPIとOpenMPの応 用)、プログラム高速化の応用      3 2023年度 計算科学技術特論A

参考書  The Art of High Performance Computing for Computational Science, Vol. 1 Editor: Masaaki Geshi 出版社: Springer, Singapore Hardcover ISBN 978-981-13-6193-7 発売日： 2019/5 【本書の特徴】  「計算科学のためのHPC技術1」の英語版  片桐担当（１章～５章）  High-Performance Computing Basics(pp. 1-25), Basics of MPI Programming (pp. 27-44), Basics of OpenMP Programming (pp. 45-59), Hybrid Parallelization Techniques (pp. 61-68), Application of Techniques for High-Performance Computing (pp. 69-81)      4 2023年度 計算科学技術特論A

参考書（演習書）  「並列プログラミング入門： サンプルプログラムで学ぶOpenMPとOpenACC」  片桐 孝洋 著  東大出版会、ISBN-10: 4130624563、 ISBN-13: 978-4130624565、発売日： 2015年5月25日  【本書の特徴】  C言語、Fortran90言語で解説  C言語、Fortran90言語の複数のサンプルプログラムが入手可能 （ダウンロード形式）  本講義の内容を全てカバー  Windows PC演習可能(Cygwin利用)。スパコンでも演習可能。  内容は初級。初めて並列プログラミングを学ぶ人向けの 入門書 5 2023年度 計算科学技術特論A

OpenMP 超入門 指示文による簡単並列化 6 2023年度 計算科学技術特論A

OpenMPの概要 7 2023年度 計算科学技術特論A

OpenMPの対象計算機  OpenMPは共有メモリ計算機のためのプログラム言語 OpenMP 実行可能 コード OpenMP 実行可能 コード OpenMP 実行可能 コード OpenMP 実行可能 コード 共有 配列 A[ ] 同時に複数のPEが共有配列にアクセス ⇒並列処理で適切に制御をしないと、逐次計算の結果と一致しない 8 2023年度 計算科学技術特論A

OpenMPとは  OpenMP (OpenMP C and C++ Application Program Interface Version １.0)とは、共有メモリ型並列計算機用 にプログラムを並列化する以下： 指示文 2. ライブラリ 3. 環境変数 を規格化したものです。 1.   9 ユーザが、並列プログラムの実行させるための指示 を与えるものです。 コンパイラによる自動並列化ではありません。 分散メモリ型並列化（ＭＰＩなど）に比べて、データ分散 の処理の手間が無い分、実装が簡単です。 2023年度 計算科学技術特論A

OpenMPとマルチコア計算機（その１）   スレッド並列化を行うプログラミングモデル 近年のマルチコア計算機に適合   経験的な性能： 8～16スレッド並列以下の実行に向く 8～16スレッドを超えるスレッド実行で高い並列化効率を出す には、プログラミングの工夫が必要 1. 2. 3.  メインメモリ‐キャッシュ間のデータ転送能力が演算性能に比べ低い OpenMPで並列性を抽出できないプログラムになっている（後述） ccNUMAの影響（ソケットを超える実行時） ノード間の並列化はOpenMPではできない   ノード間の並列化はMPIを用いる 自動並列化コンパイラも、スレッド並列化のみ  10 HPF、 XcalableMP（筑波大、理研AICS） などのコンパイラではノード間の 並列化が可能だが、まだ完全に普及していない 2023年度 計算科学技術特論A

OpenMPとマルチコア計算機（その２）  典型的なスレッド数  マルチコアＣＰＵ： 20～80スレッド／ノード  Fujitsu FX1000 (ARM A64FX)  48物理コア  Intel Xeon Gold 6230×4ソケット （「不老」クラウドシステムの１ノード）  80物理コア  メニーコアＣＰＵ： 60～280スレッド／ノード  Intel Xeon Phi (Intel MIC(Many Integrated Core) 、Knights Landing)  68物理コア、136～272論理コア（HT利用時）   ☞生産終了・・・ 2023年現在、80スレッド超のOpenMP実行環境が手元に！  11 性能を出すためには、相当のプログラム上の工夫が必要 2023年度 計算科学技術特論A

OpenMPコードの書き方の原則  Ｃ言語の場合  #pragma omp で始まるコメント行  Ｆｏｒｔｒａｎ言語の場合  !$omp で始まるコメント行 12 2023年度 計算科学技術特論A

OpenMPのコンパイルの仕方  逐次コンパイラのコンパイルオプションに、OpenMP用の オプションを付ける    例）富士通Fotran90コンパイラ frt –Kfast,openmp foo.f 例）富士通Cコンパイラ fcc –Kfast,openmp foo.c 注意  OpenMPの指示がないループは逐次実行  コンパイラにより、自動並列化によるスレッド並列化との 併用ができる場合があるが、できない場合もある  OpenMPの指示行がある行はOpenMPによるスレッド並列化、 指示がないところはコンパイラによる自動並列化  13 例）富士通Fortran90コンパイラ frt –Kfast,parallel,openmp foo.f 2023年度 計算科学技術特論A

OpenMPの実行可能ファイルの実行     OpenMPのプログラムをコンパイルして生成した実行可能 ファイルの実行は、そのファイルを指定することで行う スレッド数を、環境変数OMP_NUM_THREADSで指定 例）OpenMPによる実行可能ファイルがa.outの場合 $ export OMP_NUM_THREADS=１6 $ ./a.out 注意  逐次コンパイルのプログラムと、OpenMPによるプログラムの 実行速度が、OMP_NUM_THREADS=１にしても、異なることがある （後述）    14 この原因は、OpenMP化による処理の増加（オーバーヘッド） 高スレッド実行で、このオーバーヘッドによる速度低下が顕著化 プログラミングの工夫で改善可能 2023年度 計算科学技術特論A

OpenMPの実行モデル 15 2023年度 計算科学技術特論A

OpenMPの実行モデル（C言語） OpenMP指示文 ブロックＡ #pragma omp parallel ｛ ブロックＢ ｝ ブロックＣ ブロックＡ スレッドの起動 スレッド０ （マスタースレッド） スレッド１ ブロックＢ ブロックＢ ※スレッド数ｐは、 環境変数 OMP_NUM_THREADS で指定する。 スレッドの終結 ブロックＣ 16 2023年度 計算科学技術特論A スレッドｐ … ブロックＢ

OpenMPの実行モデル（Fortran言語） OpenMP指示文 ブロックＡ !$omp parallel ブロックＢ !$omp end parallel ブロックＣ ブロックＡ スレッドの起動 スレッド０ （マスタースレッド） スレッド１ ブロックＢ ブロックＢ ※スレッド数ｐは、 環境変数 OMP_NUM_THREADS で指定する。 スレッドの終結 ブロックＣ 17 2023年度 計算科学技術特論A スレッドｐ … ブロックＢ

Work sharing構文    parallel指示文のように、複数のスレッドで実行する場合 において、OpenMPで並列を記載する処理（ブロックB）の 部分を並列領域 (Parallel Region)と呼ぶ。 並列領域を指定して、スレッド間で並列実行する処理を 記述するOpenMPの構文をWork Sharing構文と呼ぶ。 Work Sharing構文は、以下の２種がある。 並列領域内で記載するもの 1.    parallel指示文と組み合わせるもの 2.   18 for構文（do構文） sections構文 single構文 (master構文)、など parallel for 構文 (parallel do構文) parallel sections構文、など 2023年度 計算科学技術特論A

代表的な指示文 19 2023年度 計算科学技術特論A

※Fortran言語の場合は For構文（do構文） #pragma omp parallel for for (i=0; i<100; i++){ a[i] = a[i] * b[i]; } !$omp parallel do ～ !$omp end parallel do 上位の処理 スレッドの起動 スレッド０ スレッド１ スレッド2 スレッド3 for (i=0; i<25; i++){ a[i] = a[i] * b[i]; } for (i=25; i<50; i++){ a[i] = a[i] * b[i]; } for (i=50; i<75; i++){ a[i] = a[i] * b[i]; } for (i=75; i<100; i++){ a[i] = a[i] * b[i]; } スレッドの終結 ※指示文を書くループが 並列化をしても、 正しい結果になることを ユーザが保障する。 20 下位の処理 2023年度 計算科学技術特論A

For構文の指定ができない例 for (i=0; i<100; i++) { •ループ並列化指示すると、 逐次と結果が異なる a[i] = a[i] +1; （a[i-1]が更新されていない b[i] = a[i-1]+a[i+1]; 場合がある） } for (i=0; i<100; i++) { a[i] = a[ ind[i] ]; } •ind[i]の内容により、 ループ並列化できるか どうか決まる •a[ind[i]]が既に更新された 値でないとき、 ループ並列化できる 21 計算科学技術特論A 2023年度

※Fortran言語の場合は !$omp parallel sections ～ !$omp end parallel sections Sections構文 #pragma omp parallel sections { スレッド数が３の場合 #pragma omp section sub1(); スレッド０ スレッド１ #pragma omp section sub2(); sub1(); sub2(); #pragma omp section sub3(); sub4(); #pragma omp section sub4(); } スレッド2 sub3(); スレッド数が４の場合 スレッド０ sub1(); 22 スレッド１ スレッド2 sub2(); 2023年度 sub3(); 計算科学技術特論A スレッド3 sub4();

※Fortran言語の場合は !$omp critical ～ !$omp end critical Critical補助構文 #pragma omp critical { s = s + x; } スレッド０ スレッド１ スレッド2 スレッド3 s= s+x s= s+x s= s+x s= s+x 23 2023年度 計算科学技術特論A

Private補助構文 #pragma omp parallel for private(c) for (i=0; i<100; i++){ a[i] = a[i] + c * b[i]; } ただし、c にループに入る前の値を 代入して使う場合は firstprivate(c) スレッド０ for (i=0; i<25; i++){ a[i] = a[i] + c0*b[i]; } ※変数ｃが各スレッドで 別の変数を確保して実行 →高速化される 上位の処理 スレッド１ スレッドの起動 スレッド2 for (i=25; i<50; i++){ a[i] = a[i] + c1*b[i]; } スレッド3 for (i=50; i<75; i++){ a[i] = a[i] + c2*b[i]; } for (i=75; i<100; i++){ a[i] = a[i] + c3* b[i]; } スレッドの終結 下位の処理 24 2023年度 計算科学技術特論A

Private補助構文の注意（Ｃ言語） #pragma omp parallel for private( j ) for (i=0; i<100; i++) { for (j=0; j<100; j++) { a[ i ] = a[ i ] + amat[ i ][ j ]* b[ j ]; } •ループ変数 ｊ が、各スレッドで別の変数を確保して実行される。 •private( j ) がない場合、各スレッドで 共有変数の j のカウントを独立で行ってしまい、逐次と加算結果が異なる。 →演算結果が逐次と異なり、エラーとなる。 25 2023年度 計算科学技術特論A

Private補助構文の注意（Fortran言語） !$omp parallel do private( j ) do i=1, 100 do j=1, 100 a( i ) = a( i ) + amat( i , j ) * b( j ) enddo enddo !$omp end parallel do •ループ変数 ｊ が、各スレッドで別の変数を確保して実行される。 •private( j ) がない場合、各スレッドで 共有変数の j のカウントを独立で行ってしまい、逐次と加算結果が異なる。 →演算結果が逐次と異なり、エラーとなる。 26 2023年度 計算科学技術特論A

リダクション補助構文 （C言語）  内積値など、スレッド並列の結果を足しこみ、１つの結果を 得たい場合に利用する   上記の足しこみはスレッド毎に非同期になされる reduction補助構文が無いと、ddotは共有変数になるため、 並列実行で逐次の結果と合わなくなくなる #pragma omp parallel for reduction(+: ddot ) for (i=1; i<=100; i++) { ddot += a[ i ] * b[ i ] } ddotの場所はスカラ変数のみ記載可能（配列は記載できません） 27 2023年度 計算科学技術特論A

リダクション補助構文 （Fortran言語）  内積値など、スレッド並列の結果を足しこみ、１つの結果を 得たい場合に利用する   上記の足しこみはスレッド毎に非同期になされる reduction補助構文が無いと、ddotは共有変数になるため、 並列実行で逐次の結果と合わなくなくなる !$omp parallel do reduction(+: ddot ) do i=1, 100 ddot = ddot + a(i) * b(i) enddo !$omp end parallel do ddotの場所はスカラ変数のみ記載可能（配列は記載できません） 28 2023年度 計算科学技術特論A

reduction補助構文の注意  reduction補助構文は、排他的に加算が行われるので、 一般的に性能が悪い    経験的に、8~16スレッド並列を超える場合、性能劣化が激しい 以下のように、ddot用の配列を確保して逐次で加算する ほうが高速な場合もある（ただし、問題サイズ、ハードウェア依存） ） !$omp parallel do private ( i ) スレッド数分のループを作成：最大ｐスレッド利用 do j=0, p-１ 各スレッドでアクセスするインデックス範囲を事前に設定 do i=istart( j ), iend( j ) ddot_t( j ) = ddot_t( j ) + a(i) * b(i) 各スレッドで用いる、ローカルなddot用の enddo 配列ddot_t()を確保し、0に初期化しておく enddo !$omp end parallel do ddot = 0.0d0 do j=0, p-１ 逐次で足しこみ ddot = ddot + ddot_t( j ) enddo 29 2023年度 計算科学技術特論A

その他、よく使うOpenMPの関数 30 2023年度 計算科学技術特論A

最大スレッド数取得関数   最大スレッド数取得には、omp_get_num_threads()関数 を利用する 型はinteger (Fortran言語)、int (C言語)  Fortran90言語の例  C言語の例 use omp_lib Integer nthreads #include <omp.h> int nthreads; nthreads = omp_get_num_threads() nthreads = omp_get_num_threads(); 31 2023年度 計算科学技術特論A

自スレッド番号取得関数   自スレッド番号取得には、omp_get_thread_num()関数を 利用する 型はinteger (Fortran言語)、int (C言語)  Fortran90言語の例  C言語の例 use omp_lib Integer myid #include <omp.h> int myid; myid = omp_get_thread_num() myid = omp_get_thread_num(); 32 2023年度 計算科学技術特論A

時間計測関数   時間計測には、omp_get_wtime()関数を利用する 型はdouble precision (Fortran言語)、double (C言語)  Fortran90言語の例  C言語の例 use omp_lib double precision dts, dte #include <omp.h> double dts, dte; dts = omp_get_wtime() 対象の処理 dte = omp_get_wtime() print *, “Elapse time [sec.] =”,dte-dts dts = omp_get_wtime(); 対象の処理 dte = omp_get_wtime(); printf(“Elapse time [sec.] = %lf ¥n”, dte-dts); 33 2023年度 計算科学技術特論A

その他の構文 34 2023年度 計算科学技術特論A

※Fortran言語の場合は Single構文    !$omp single ～ !$omp end single Single補助構文で指定されたブロックを、 どれか１つのスレッドに割り当てる どのスレッドに割り当てられるかは予測できない nowait補助構文を入れない限り、同期が入る プログラムの開始 #pragma omp parallel for ｛ ブロックA #pragma omp single { ブロックB } … } 35 スレッドの起動 スレッド０ スレッド１ （マスタースレッド） ブロックA ブロックA ブロックＢ 同期処理 2023年度 計算科学技術特論A スレッドｐ … ブロックA

Master構文  使い方は、single補助構文文と同じ  ただし、master補助構文で指定した処理 （先ほどの例の「ブロックB」の処理）は、 必ずマスタースレッドに割り当てる  終了後の同期処理が入らない  36 そのため、場合により高速化される 2023年度 計算科学技術特論A

Flush構文   物理メモリとの一貫性を取る Flush構文で指定されている変数のみ、その場所で一貫性を取る。 それ以外の共有変数の値は、メモリ上の値との一貫性は無い。     演算結果はレジスタ上に保存されるだけ。 メモリに計算結果を書き込んでいない つまり、flush補助指定文を書かないと、スレッド間で同時に足しこんだ結果 が、実行ごとに異なる。 barrier補助構文、critical補助構文の出入口、parallel構文の出口、 for、sections、single構文の出口では、暗黙的にflushされている。 Flushを使うと性能は悪くなる。できるだけ用いない。 #pragma omp flush (対象となる変数名の並び) 省略すると、 全ての変数が対象 37 2023年度 計算科学技術特論A

Threadprivate構文



スレッドごとにプライベート変数にするが、スレッド内で大域アクセスできる
変数を宣言する。
スレッドごとに異なる値をもつ大域変数の定義に向く。

たとえば、スレッドごとに異なるループの開始値と終了値の設定
…
#include <omp.h>
void main() {
int myid, nthreds, istart, iend;
…
#pragma omp threadprivate (istart,
#pragma omp parallel private (myid, nthreds,
iend)
istart, iend) {
…
nthreds = omp_num_threds();
void kernel() {
myid = omp_get_thread_num();
int i;
istart = myid * (n/nthreads);
for (i=istart; i<iend; i++) {
iend = (myid+１)*(n/nthreads);
for (j=0; j<n; j++) {
if (myid == (nthreads-１)) {
a[ i ] = a[ i ] + amat[ i ][ j ] * b[ j ];
nend = n;
}
スレッド毎に異なる値を持つ
}
}
大域変数を、Parallel構文中
kernel();
}
で定義する
}
…


38

2023年度

計算科学技術特論A

スケジューリング 39 2023年度 計算科学技術特論A

スケジューリングとは（その１）  Parallel do構文（Parallel for構文）では、対象ループの 範囲（例えば１～ｎの長さ）を、単純にスレッド個数分に 分割（連続するように分割）して、並列処理をする。 1  スレッド0 スレッド1 スレッド2 スレッド3 スレッド4 n ループ変数の流れ （反復空間） このとき、各スレッドで担当したループに対する計算負荷 が均等でないと、スレッド実行時の台数効果が悪くなる 1 スレッド0 スレッド1 スレッド2 スレッド3 スレッド4 n 計算負荷 40 2023年度 計算科学技術特論A

スケジューリングとは（その２）  負荷分散を改善するには、割り当て間隔を短くし、かつ、 循環するように割り当てればよい。 1 n 計算負荷   41 最適な、割り当て間隔（チャンクサイズとよぶ）は、計算機 ハードウェアと、対象となる処理に依存する。 以上の割り当てを行う補助構文が用意されている。 2023年度 計算科学技術特論A

ループスケジューリングの補助構文 （その１）  schedule (static, n)   1 ループ長をチャンクサイズで分割し、スレッド0番から順番に （スレッド０、スレッド１、・・・というように、ラウンドロビン方式と 呼ぶ）、循環するように割り当てる。 nにチャンクサイズを指定できる。 Schedule補助指定文を記載しないときのデフォルトはstaticで、 かつチャンクサイズは、ループ長/スレッド数。 スレッド0 42 スレッド1 スレッド2 スレッド3 2023年度 計算科学技術特論A

ループスケジューリングの補助構文 （その２）  schedule(dynamic, n)  ループ長をチャンクサイズで分割し、処理が終了したスレッド から早い者勝ちで、処理を割り当てる。 nにチャンクサイズを指定できる。 1 スレッド0 43 スレッド1 スレッド2 スレッド3 2023年度 計算科学技術特論A

ループスケジューリングの補助構文 （その３）  schedule(guided, n)  ループ長をチャンクサイズで分割し、徐々にチャンクサイズを 小さくしながら、処理が終了したスレッドから早い者勝ちで、 処理を割り当てる。nにチャンクサイズを指定できる。  チャンクサイズの指定が１の場合、残りの反復処理をスレッド数で割った おおよその値が各チャンクのサイズになる。 チャンクサイズは １ に向かって指数的に小さくなる。 チャンクサイズに １より大きい k を指定した場合、チャンク サイズは指数 的に k まで小さくなるが、最後のチャンクは k より小さくなる場合がある。 チャンクサイズが指定されていない場合、デフォルトは １ になる。    1 スレッド0 44 スレッド1 スレッド2 スレッド3 2023年度 計算科学技術特論A

ループスケジューリングの補助構文
の使い方
 Fortran90言語の例
!$omp parallel do private( j ) schedule(dynamic,１０)
do i=１, n
do j=indj(i), indj (i+１)-１
y( i ) = amat( j ) * x( indx( j ) )
enddo
enddo
!$omp end parallel do

 C言語の例

45

#pragma omp parallel for private( j ) schedule(dynamic,１０)
for (i=0; i<n; i++) {
for ( j=indj(i); j<indj (i+１); j++) {
y[ i ] = amat[ j ] * x[ indx[ j ]];
}
}
2023年度

計算科学技術特論A

ループスケジューリングにおける プログラミング上の注意  dynamic、guidedのチャンクサイズは性能に大きく影響      実行時のチャンクサイズのチューニングが必須で、チューニングコスト が増える。 staticのみで高速実装ができる（場合がある）    46 チャンクサイズが小さすぎると負荷バランスは良くなるが反面、 システムのオーバヘッドが大きくなる。 一方、チャンクサイズが大きすぎと負荷バランスが悪くなる半面、 システムのオーバヘッドが小さくなる。 両者のトレードオフがある。 dynamicなどの実行時スケジューリングは、システムのオーバーヘッド が入るが、staticはオーバーヘッドは（ほとんど）無い。 事前に負荷分散が均衡となるループ範囲を調べた上で、 staticスケジューリングを使うと、最も効率が良い可能性がある。 ただし、プログラミングのコストは増大する 2023年度 計算科学技術特論A

Staticスケジューリングのみで負荷バランス を均衡化させる実装例  疎行列‐ベクトル積へ適用した例（詳細は後述） スレッド個数分のループ （スレッドごとのループ担当範囲 を知るために必要） !$omp parallel do private(S, J_PTR,I) DO K=１, NUM_SMP DO I=KBORDER(K-１)+１, KBORDER(K) 事前に調べて設定しておいた、 S=0.0D0 負荷分散が均衡となる DO J_PTR=IRP(I), IRP(I+１)-１ スレッドごとのループ範囲 S=S + VAL( J_PTR ) * X(ICOL( J_PTR)) （各スレッドは、連続しているが、 END DO 不均衡なループ範囲を設定） Y(I)=S END DO 実行前に、各スレッドが担当するループ範囲について、 END DO 連続する割り当てで、かつ、それで負荷が均衡する !$omp end parallel do 問題に適用できる。 ※実行時に負荷が動的に変わっていく場合は適用できない 47 2023年度 計算科学技術特論A

OpenMPのプログラミング上の注意 （全般） 48 2023年度 計算科学技術特論A

OpenMPによるプログラミング上の注意点 OpenMP並列化は、  parallel構文を用いた単純なforループ並列化 が主になることが多い。  複雑なOpenMP並列化はプログラミングコストがかかるので、 OpenMPのプログラミング上の利点が失われる parallel構文による並列化は  private補助構文の正しい使い方 を理解しないと、バグが生じる！ 49 2023年度 計算科学技術特論A

Private補助構文に関する注意（その１）  OpenMPでは、対象となる直近のループ変数以外は、 private変数で指定しない限り、全て共有変数になる。  デフォルトの変数は、スレッド間で個別に確保した変数でない  ループ変数に関する共有変数の例 !$omp parallel do 宣言なしにプライベート変数として確保されるのは、 このi-ループ変数のみ do i=1, 100 このj-ループ変数は、private宣言なしでは共有変数になる do j=1, 100 ←スレッド間で早い者勝ちで加算 ←並列実行時にバグ tmp = b(i) + c(i) a( i ) = a( i ) + tmp enddo この変数tmpは、private宣言なしでは共有変数になる ←スレッド間で早い者勝ちで値が代入 ←並列実行時にバグ enddo !$omp end parallel do 50 2023年度 計算科学技術特論A

Private補助構文に関する注意（その２） Private補助構文に記載する変数を減らすため、 対象部分を関数化し、かつ、その関数の引数を増やすと、 関数呼び出し時間が増加し、スレッド並列化の効果を 相殺することがある   呼び出し関数の引数が多い例 !$omp parallel do do i=1, 100 call foo(i,arg1,arg2,arg3, arg4,arg5, ….., arg100) enddo !$omp end parallel do 51 関数引数は自動的にプライベート変数に なるため、private補助構文に記載する 変数を削減できる ← しかし、関数呼び出し時のオーバーヘッド が増加する ← スレッド実行時においても、関数呼び出し のオーバーヘッドが無視できなくなり、 台数効果が制限される ※解決法：大域変数で引き渡して引数を削減 2023年度 計算科学技術特論A

Private補助構文に関する注意のまとめ   OpenMPでは、宣言せずに利用する変数は、 すべて共有変数（shared variable）になる C言語の大域変数、Fotran90言語のcommon変数は、 そのままでは共有変数になる   プライベート変数にしたい場合は、Threadprivate宣言が必要 外側ループをParallel構文などで並列化する場合   52 ループ内で呼ばれる関数（手続き）内で宣言される 変数はPrivateになる C言語で、ループ内で明示的に宣言される変数 (例：int a;) は Privateになる。 2023年度 計算科学技術特論A

Parallel構文の入れ子に関する注意（その１）   Parallel構文は、do構文で分離して記載できる １ループが対象の場合、分離するとdo構文の場所で ループごとにforkするコードを生成するコンパイラがあり、 速度が低下する場合がある  この逆の場合もある。双方を確認する必要あり。 !$omp parallel !$omp do private(j,tmp) do i=1, 100 do j=1, 100 tmp = b( ｊ ) + c( ｊ ) a( i ) = a( i ) + tmp enddo enddo !$omp end do !$omp end parallel 53 Parallel構文の 対象が１ループ なら parallel do で指定 2023年度 !$omp parallel do private(j,tmp) do i=1, 100 do j=1, 100 tmp = b( ｊ ) + c( ｊ ) a( i ) = a( i ) + tmp enddo enddo !$omp end parallel do 計算科学技術特論A

Parallel構文の入れ子に関する注意（その２）   Parallel構文は、do構文で分離して記載できる 複数ループの内側を並列化したい場合は、分離した 方が高速になる   ただし、外側ループを並列化できる時はその方が性能が良い 外側ループにデータ依存があり、並列化できない場合 !$omp parallel do i=1, n !$omp do do j=1, n <並列化できる式> enddo !$omp end do enddo !$omp end parallel do i=1, n !$omp parallel do do j=1, n <並列化できる式> enddo !$omp end parallel do enddo 54 2023年度 計算科学技術特論A

データ依存関係を壊しバグになる例  間接参照があるインデックスに対して加算する例  間接参照のパターン、および、スレッド実行のタイミング次第で、 逐次処理と結果が一致し、正常動作だと勘違いする場合がある   理論的には間違っている OpenMPの共有変数は、データ一貫性の保証はしない  データの一貫性の保証には、critical補助指定文などによる相互排除が必要  バグになるプログラム例 !$omp parallel do private( j ) do i=1, n j = indx( i ) a( j ) = a( j ) + 1 enddo !$omp end parallel do 55 2023年度 !$omp parallel do private( j ) do i=１, n j = indx( i ) !$omp critical a( j ) = a( j ) + １ !$omp end critical enddo !$omp end parallel do 計算科学技術特論A

Critical補助構文による速度低下  先述のように、critical補助構文を入れないといけない場合、 特に、高スレッド数での実行で性能が低下する    1. 高性能化するには、基本的にはアルゴリズムを変更するしかない。 この場合、以下の３つのアプローチがある。 スレッド内アクセスのみに限定し、critical補助構文をはずす  間接参照されるデータについて、理論的に、割り当てられたスレッド内の データしかアクセスしないように、アルゴリズムを変更する スレッド間アクセスを最小化 2.  3. ハードウェアサポートがある場合、atomic節による実装のほうが高速であ ることが多いが、スレッド数が増えると性能が低下する criticalの並列領域に同時に入るスレッド数が減るように、間接参照するデータを 事前に調べ、間接参照するデータの順番を変更する。 スレッド間アクセス部分をループから分離し、逐次処理にする  例）内積演算におけるリダクション補助指定文 56 2023年度 計算科学技術特論A

OpenMPを用いた並列化の欠点 （その１）  OpenMPは単純なループを並列化することに向く  実用アプリケーションにおける複雑なループは、そのままでは OpenMP化に向いていないことがある。 private補助構文中に書かれる変数名の数が膨大になる  外側ループからOpenMP並列化する場合、内部で使っている 変数の数が多いことがある  private変数リストに変数を書き忘れても、コンパイラによる エラーは出ない。（並列化の責任はユーザにあるため）  実行すると、タイミングに依存し計算結果が逐次と異なる。 どこが間違っているかわからないので、デバックが大変になる。  解決策：コンパイラによっては、最適化情報を出力することがで きる。その情報から、ちゃんとprivate化されているか確認する。 1. 57 2023年度 計算科学技術特論A

OpenMPを用いた並列化の欠点 （その２） 2.   3.   高スレッド実行時に性能が出ない場合のチューニングが困難 一般に、8~16スレッド未満では性能が出るが、8~16スレッド以上 で性能が劣化する。 1. 近年のハードウェアはメモリアクセスの性能が低い 2. ループそのものに並列性がない（ループ長が短い） 3. ccNUMAの影響（ソケットを跨ぐ実行） 解決するには、アルゴリズムの変更、実装の変更、が必要 になり、OpenMPの利点である容易なプログラミングを損なう 複雑なスレッドプログラミングには向かない 単純な数値計算のカーネルループを、parallel for構文で記載 する方針で仕様が作られている（と思われる） 複雑な処理は、PthreadなどのnativeなスレッドAPIで書くほうが やりやすい 58 2023年度 計算科学技術特論A

プログラム実例 59 2023年度 計算科学技術特論A

行列-行列積のコードのOpenMP化の例 （Ｃ言語） 以下のようなコードになる  #pragma omp parallel for private (j, k) 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]; } } } 60 2023年度 計算科学技術特論A

行列-行列積のコードのOpenMP化の例 （Fortran言語） 以下のようなコードになる  !$omp parallel do private (j, k) do i=1, n do j=1, n do k=1, n C(i, j) = C(i, j) + A(i, k) * B(k, j) enddo enddo enddo !$omp end parallel do 61 2023年度 計算科学技術特論A

OpenMPの高速化技法： ファーストタッチ 62 2023年度 計算科学技術特論A

ファーストタッチとは    ファーストタッチとは、マルチコア計算機の中でも、 ccNUMA (Cache Coherent Non-Uniform Memory Access) のハードウェア向けの、メモリ最適化の方法 OpenMPによる並列プログラミングでも重要な技法 ccNUMAのメモリ構造の特性を利用する アクセス遅い CPU0 CPU1 メモリ0 メモリ1 メモリ2 メモリ3 CPU2 CPU3 アクセス速い 63 2023年度 ccNUMAの ハードウェア 計算科学技術特論A

ファーストタッチの原理    64 ccNUMA型のハードウェアでは、確保した配列は、 各コアで、その配列に初めてアクセスした時、 各コアに最も近いメモリに配列が置かれる この原理を利用し、本計算と同じデータ・アクセスパターン （＝ループ構造）で、プログラム上最も先に、 OpenMP指示文を用いて配列を初期化すると、 CPUに近いメモリに配列データがセットされる 本計算と同じループ構造で、確保した配列の初期化 （例えば0クリア、もしくが、データのセット）をするだけで、 ファーストタッチが実現できる 2023年度 計算科学技術特論A

ファーストタッチの例 (C言語の例) #pragma omp parallel for private( j ) for (i=0; i<100; i++) { for (j=0; j<100; j++) { a[ i ] = 0.0; amat[ i ][ j ] = 0.0; } …. #pragma omp parallel for private( j ) for (i=0; i<100; i++) { for (j=0; j<100; j++) { a[ i ] = a[ i ] + amat[ i ][ j ]* b[ j ]; } 65 2023年度 ファーストタッチの ための初期化 （プログラムの 一番最初に 実行すること） ファーストタッチ データを利用した 本計算 計算科学技術特論A

ファーストタッチの例 (Fortran言語の例) !$omp parallel do private( j ) do i=1, 100 do j=1, 100 a( i ) = 0.0d0 amat( i , j ) =0.0d0 enddo enddo !$omp end parallel do …. !$omp parallel do private( j ) do i=1, 100 do j=1, 100 a( i ) = a( i ) + amat( i , j ) * b( j ) enddo enddo !$omp end parallel do 66 2023年度 ファーストタッチの ための初期化 （プログラムの 一番最初に 実行すること） ファーストタッチ データを利用した 本計算 計算科学技術特論A

ファーストタッチの効果 うまくいく場合で、経験的に約２～５倍高速化  効果的な例    必ずしも、効果があるとは限らない   67 キャッシュに載るようなサイズやアクセス領域 での実行時 キャッシュミスヒットが大きい場合 そもそも、ファーストタッチの実装が困難 ☞次のスライド 2023年度 計算科学技術特論A

ファーストタッチの実装上の注意  ccNUMAのアーキテクチャでないと効果がない   スーパーコンピュータ「富岳」、スーパーコンピュータ「不老」では、 ハードウェア的に効果が期待できる 対象となる配列を自ら確保し、演算も自ら行う 「手製の」プログラムでないと効果がない  数値計算ライブラリを使う場合     68 配列データはユーザが用意する。 一般的に、配列データの値を設定するプログラムが先に 動いて、その後、数値計算ライブラリを呼ぶ。 このとき、数値計算ライブラリ内でのアクセスパターンがわからないので、 配列データを設定するプログラムのアクセスパターンが数値計算ライブラリ 内のデータアクセスパターンと異なる。 以上の理由から、ファーストタッチできない。 2023年度 計算科学技術特論A

task構文 69 2023年度 計算科学技術特論A

task構文 (OpenMP 3.1)  今までは、スレッド並列で実行   70 全ての実行は、対象実行の前に並列性（スレッド数）を OMP_NUM_THREADSで設定 タスク構文  タスク並列で行う  対象の場所の並列性を柔軟に増やしたり、 実行するコア割り当てを実行時に決めることがで きる ⇒記述の柔軟性が高い  欠点：タスク生成コストが高く、実行時間で スレッド並列に劣ることがある 2023年度 計算科学技術特論A

task構文の例 !$omp parallel num_threads(3) !$omp single !$omp task タスクA !$omp end task !$omp task タスクB !$omp end task !$omp task タスクC !$omp end task !$omp end single !$omp end parallel 71 スレッドの起動 スレッド０ スレッド１ スレッド2 タスクA生成 タスクB生成 タスクA実行 タスクC生成 タスクC実行 2023年度 計算科学技術特論A タスクB実行

OpenMP 4.0 72 2023年度 計算科学技術特論A

OpenMP 4.0  ２０１３年７月仕様公開   デバイス（GPU等）へのOpenMP演算のオフロード指定   SIMD構文 スレッドとコアへの割り当て指定（NUMAアフィニティ）   Terms構文 SIMD指定   Target構文 複数の並列デバイスを指定   http://www.openmp.org/mp-documents/OpenMP4.0.0.pdf Proc_bind節 GPU利用について、後述のOpenACCと同等の機能 73 2023年度 計算科学技術特論A

• http://www.openmp.org/mp-documents/OpenMP4.0.0.pdf

OpenACCへの展開 74 2023年度 計算科学技術特論A

OpenACCへの展開  GPUを、OpenMPのように、ディレクティブで指定して使う OpenACCが普及しつつある。   OpenMP化されたプログラムは、比較的簡単に、 OpenACCに変換できる   OpenMP の Parallel構文 → OpenACC の Kernel構文 か Parallel構文 に書き換え 注意する点は：   75 OpenMP 4.0でもGPUを扱える。どちらが普及するかわからない。 CPU→GPU、および、GPU→CPUのデータ移動の最小化 データ転送の対象となる配列を指定するData構文が重要 2023年度 計算科学技術特論A

Data構文の節 !$acc data … !$acc end data GPU copyin A データの転送 A 結果の書戻し create present copyout CPUメモリ 76 デバイスメモリ 2023年度 計算科学技術特論A

do iter = 1, MAX_ITER !$acc kernels do i=1, n do j=1, n b(i) = A(i, j) * … enddo enddo !$acc end kernels … !$acc kernels do i=1, n do j=1, n b(i) = b(i) + A(i, j) * … enddo enddo !$acc end kernels … enddo 77 CPUメモリ A(i, j) b(i) デバイスメモリ データの転送 結果の書戻し b(i) b(i) 2023年度 b(i) デバイスメモリ CPUメモリ A(i, j) A(i, j) データの転送 結果の書戻し 計算科学技術特論A A(i, j) b(i) b(i)

!$acc data copyin(A) create(b) do iter = 1, MAX_ITER !$acc data present(A, b) !$acc kernels do i=1, n do j=1, n b(i) = A(i, j) * … enddo enddo !$acc end kernels !$acc end data … !$acc data present(A, b) !$acc kernels do i=1, n do j=1, n b(i) = b(i) + A(i, j) * … enddo enddo !$acc end kernels !$acc end data … enddo !$acc78end data CPUメモリ A(i, j) デバイスメモリ データの転送 A(i, j) b(i) デバイスメモリ A(i, j) b(i) デバイスメモリ上にあるデータのみで演算 （ＣＰＵメモリからの転送、および、 ＣＰＵメモリへの書き戻しが無い） 2023年度 計算科学技術特論A

レポート課題（その１） 問題レベルを以下に設定  問題のレベルに関する記述： •L00: きわめて簡単な問題。 •L10： ちょっと考えればわかる問題。 •L20： 標準的な問題。 •L30： 数時間程度必要とする問題。 •L40： 数週間程度必要とする問題。複雑な実装を必要とする。 •L50： 数か月程度必要とする問題。未解決問題を含む。 ※L４０以上は、論文を出版するに値する問題。  教科書のサンプルプログラムは以下が利用可能 （ただし、MPIの部分をコメントアウトする必要あり）   79 Mat-Mat-noopt-fx.tar Mat-Vec-fx.tar 2023年度 計算科学技術特論A

レポート課題（その２） 1. 2. 3. [L10] 行列‐行列積のコードをOpenMPで並列化せよ。 また、１スレッド実行に対する台数効果を測定せよ。 [L10] 行列‐行列積のコードについて、ファーストタッチを 実装し、性能を評価せよ。 [L20]疎行列‐行列積のコードについて、OpenMPで並列化 せよ。また、１スレッド実行に対する台数効果を測定せよ。 80 2023年度 計算科学技術特論A

レポート課題（その３） 4. 5. 6. 7. [L10] データスコープ属性とは何か調べよ。また、 firstprivate, lastprivate補助構文の機能は何かを調べよ。 [L10] Barrier指示文、Nowait補助構文について調べよ。ま たどのように利用するか例を記載して説明せよ。 [L10] 本講義で取り上げていない、OpenMPの 実行時ライブラリ関数を調べ、その機能と利用方法を記せ。 [L10] OMP_NUM_THREADS以外のOpenMPで定義された 環境変数を調べ、その機能を説明せよ。 81 2023年度 計算科学技術特論A

レポート課題（その４） 8. 9. 10. 11. [L10] スケジューラの補助指示構文runtimeの機能調べよ。 また、OpenMPの環境変数との関係を説明せよ。 [L15] OpenMP version 3.0、もしくは、4.0の仕様を調べよ。 [L15] OpenACC version 1.0、もしくは2.0 の仕様を調べよ。 [L10～] 自分の持っている逐次コードを、OpenMPで並列化 せよ。スレッド数を変化させて、台数効果を調べよ。 82 2023年度 計算科学技術特論A