コンピュータ・システムの勉強の一環としてC++を用いた行列-行列積演算の高速化を行ったうちの「スレッド並列化編」としてここで簡単にまとめておきます。
勉強したことをまとめがてらアウトプットしているので、間違い等ございましたらご指摘いただけると助かります。
※参考文献は本ページ末尾にあります
並列化概要
まず、並列化の目的は処理の高速化です。以下の図のように複数の演算器を用いてタスクをうまく並列化することができれば、時間が大幅に短縮できるはずです。
並列化と一口に言っても、今回のOpenMPを用いたスレッド並列にあたるノード内の並列化と、MPIを用いるノード外での並列化などがあります。
一般的なパソコンの中に搭載されているような、マルチコア・プロセッサやメニーコア・プロセッサを用いた計算では、OpenMPで行うようなスレッド並列化のようなプログラムモデルが適しています。
一方、スーパーコンピュータのような、複数のプロセッサを備えた環境ではMPIを用いた分散メモリ型並列化によって、ノード間並列化を行うプログラムモデルが必要とされたりします。
ほかにも、OpenACCやCUDAを用いてGPUを利用するものや、OpenCLを用いてその他のプロセッサと協調計算を行うプログラムモデルなんかもあったりします。
OpenMPを用いたスレッド並列
今回用いるマルチスレッドによる並列化について少しだけ踏み込んでまとめておきます。
今回評価に用いているCPUは「Intel Core i7-9700K」です。これは８コア８スレッドを持っています。コアは演算処理を行うやつのことです。コア＝CPUという表現がされていることもあります。
パタヘネによると、1つのチップ（CPU）上に載ったプロセッサのことをメーカーがコアと呼び、これらを複数載せたマイクロプロセッサののことを「マルチコア・マイクロプロセッサ」と呼んでいるそうです。スレッドは、そのようなプロセッサがあるわけでなく1種のプロセスのことを言います。
共有記憶型マルチプロセッサでは、各プロセッサ内のスレッドが共有メモリ上の単一物理アドレス空間にアクセスしながら処理を行っています。
また、共有メモリへのアクセスにおいては、異なるアドレスへのアクセスであれば複数のスレッドが同時にメモリへ読み出し処理を行った場合でも、メモリコントローラーによって同時にメモリアクセスを処理することができる（メモリインターリーブ）ため、共有メモリを用いる場合でもメモリアクセスがボトルネックになることはないそうです。ただし、メモリ帯域幅には限界があります。
評価に用いる行列-行列積（GEMM）
全ての要素に値を入れた密な行列AとBを入力として行列の積（General Matrix Multiply）を計算し、結果として行列Cを出力する処理を考えます。例を示すと以下のような処理です。
前提条件として「行列Ａの列数＝行列Ｂの行数」となっていることが挙げられます。得られる結果の行列Ｃの大きさは、「行×列」＝「行列Ａの行数×行列Ｂの列数」です。
行列積演算では、複数の配列により行列の要素を格納するので、多数のメモリアクセスが発生し、キャッシュメモリによる処理速度の差が表れやすい処理だと考えられます。
また、複数の配列を並列に処理することで、並列化による高速化が期待できると思われます。
このような行列のデータ構造は、画像処理にも応用できるはずです。
並列化のためのコード詳細
データ構造と、並列化の具体的な手法について見ていきます。
コードの全文は、以下Githubに置いています。
https://github.com/sukima-log/parallel_test
データ構造
メモリの確保を行って、そこに行列データを用意していきます。
/* メモリ確保 */
int** a; // 入力行列A InputMatrixA
int** b; // 入力行列B InputMatrixB
long** c; // 出力行列C OutputMatrixC
// 行先頭位置
a = (int**) malloc(sizeof(int*) * ROW_A);
b = (int**) malloc(sizeof(int*) * ROW_B);
c = (long**) malloc(sizeof(long*) * ROW_A);
// 行配列
for(size_t i=0; i<ROW_A; i++) {
a[i] = (int*) malloc(sizeof(int) * COL_A);
c[i] = (long*) malloc(sizeof(long) * COL_B);
}
for(size_t i=0; i<ROW_B; i++) {
b[i] = (int*) malloc(sizeof(int) * COL_B);
}09～16行目で各行の要素を格納する行数個の配列をそれぞれの行列ぶん確保し、06～08行目で各行の配列先頭位置を格納する配列を確保しています。図にすると以下です。
疎行列などの特別な行列であれば、値の格納形式にも様々な方式があるようですが、今回は密行列なので特に工夫もなくメモリを確保しています。しいて意識するなら、メモリアクセスがランダムにならないようにメモリアライメントに気を使っておくべきなのかもしれませんが、今回はmallocの仕様に任せています。
計算精度
桁あふれを一応考慮して、Int型（4Byte）で入力行列を用意して
→ 環境によっても変わりますがどちらも4Byteになっていました。
OpenMPのプラグマ指示子
各行列へのアクセスの方向としては、「キャッシュメモリ編」で用いたキャッシュメモリを意識してループを回していきます。
size_t i, j, k;
#pragma omp parallel for private(j, k)
for(i=0; i<ROW_A; i++) {
for(j=0; j<ROW_B; j++) {
for(k=0; k<COL_B; k++) {
c[i][k] += (long)a[i][j] * (long)b[j][k];
}
}
}内側ループのjとkは、各スレッド間で共有の変数としたくないので、privateを用いてスレッド間で個別に確保する変数に指定しています。ここは地味に忘れやすいポイントなので気を付ける必要があります。
行列AとBはデータの読み出しのみ行われますが、行列Cは値の加算が行われるため読み出しと書き込みの両方が行われます。
行列Cでは行方向に並列化を行った場合、書き込みについてはデータ依存がないです。読み出しについても、i方向に並列化されているので各スレッド間で別々の配列にアクセスしておりデータ依存は無いはずです。
行列Aも同様に、i方向に並列化された行列では複数のスレッド間で別々の配列を用いて処理しているはずなので、各スレッド間の読み出しにデータ依存はないはずです。
問題はBなのですが、すべてのスレッドで同じメモリアドレスの値が読み出しされているはずです。この辺りは、書籍でもぼやかされた表現がされていることが多い気がしますが、複数のスレッドが同時に同じメモリアクセスを行う「共有読み出し」と呼ばれる読み出し方式がサポートされてるらしいです。
この辺りは、物理的に全く同じメモリアドレスに複数のスレッドが全く同時にアクセスすることは無理な気がするので、メモリインターリーブという仕組みによって複数のスレッドが1つずつ同じアドレスにアクセスするタイミングをずらしながら同時にメモリにアクセスしているものなのかなと思ってます。
少なくとも同じメモリアドレスに書き込みは行っていないので、データの整合性上の問題が生じる競合は発生しません。よって、OpenMPにより自動的に必要な同期機構が提供されているくらいの理解でいいような気はします。
速度評価
上記の並列化によりどれだけ計算速度に差が出るのか評価していきます。
実装環境
main: main.o
g++ -o main -O2 -mcmodel=medium -std=c++11 -Ofast -fopenmp main.o
main.o: main.cpp
g++ -c -O2 -mcmodel=medium -std=c++11 -Ofast -fopenmp main.cpp
時間計測
図 行列サイズ変化による行列-行列積時間
プロセスを監視するコマンドである、htopコマンドからも8スレッドすべてを用いてプログラムを実行している様子を確認することができました。
並列版は8倍の並列数を持つので純粋に8倍の処理速度を持っているのかというと、それほど都合良くはいきませんが、かなり高速化できていることが分かります。具体的には4000×4000の行列サイズだと4.8倍くらい高速化できています。
現在のハードウェアでは、メモリからのデータ転送時間に比べ、演算速度の方がかなり高速と言われているので、スレッド数に比例した処理性能は得られないことがほとんどだと言われています。いわゆるメモリ律速な処理になっているというやつだと思います。
まとめ
並列化について調べてるとメモリについて理解が浅いことに気が付いたり、いろんな思想があったりと割と楽しいので今後も勉強していきたいと思います。
参考文献：
- "並列プログラミング入門：サンプルプログラムで学ぶOpenMPとOpenACC ",片桐 孝洋
- "コンピュータの構成と設計 第5版 下",デイビッド・A・パターソン and ジョン・L・ヘネシー
