はじめに
どうもアドベントカレンダーとしてはかなりの分野違いをしてしまっているようです。
皆さんが複雑な数式を紹介する中で、本記事は実数を要素に持つ行列A,B,Cに対して
C = A \times B
これしか計算していない上に各行列が何を表しているのかについては不問です。
その代わりこの式をどのように計算機上で高速に計算をするのかについて紹介しています。
「なんか泥臭いことやっとるな」くらいの認識でOKです。
OpenBLASに関するTips的な
早速OpenBLASを使用する際のTipsについて簡単に説明したいと思います。
「え、多次元配列の処理はjuliaでやるからOpenBLAS使わないよ」と言う方、
公式githubのページの"Platform-Specific Build Notes"を見ると、
実はjulia内部でOpenBLASを使っている(らしい)ので間接的に使っていることになるんですね。
NumPy/SciPyを使うような方も以下の方法は(多分)有効です。
使い方と言っても普段は何も考えずに使用していても問題ありませんが、
例えばかなり重たい行列積をやるときにはプログラムを起動する際の環境変数に、
OMP_NUM_THREADS=(計算機の物理コア数)
OMP_NUM_THREADS=spread
を追加すると嬉しいことがあるかもしれません。(https://github.com/xianyi/OpenBLAS/issues/1653)
juliaをjupyterから使う場合にはjupyter notebookを起動する際に設定する必要があります。
ここで言う計算機の物理コア数は論理コアの数ではないのに注意です。
使用するOSがLinuxであればLANG=C lscpuを実行してCore(s) per socketの項を確認すればOKです。
筆者の使用する環境ではCPUがIntel(R) Core(TM) i7-7600U CPU @ 2.80GHzなので以下のような結果が得られました。
% LANG=C lscpu | grep "Core(s)"
Core(s) per socket: 2
もし強い計算機(Xeon搭載レベル)を使っていてソケット数が1よりも大きい場合は、
Core(s) per socketで得られた数字にソケット数をかけたものが物理コア数になります。
Macの場合はsystem_profiler SPHardwareDataTypeからTotal Number of Coresで、
Windowsの場合はTask Manager→Performance→Coresで確かめられるようです。
Mac/Linux: http://mikamisan.hatenablog.com/entry/2017/01/09/101739
Windows: https://wiki.windward.net/Returning_Users/Cross-Product_Articles/Determining_the_Number_of_Cores_for_Your_System
本記事のこれ以降の内容を読むと、
なぜこのような環境変数を設定すると良いのかがわかるかもしれません。
OpenBLASを探るだけ
Intel Math Kernel Library(MKL)と遜色ないパフォーマンスを出すOpenBLASのコードを探り、
IntelのCPUにおける行列積の高速化手段について調べようと思います。
内容としては去年の自分の記事で記述したものの答え合わせ的なものになっています。
なのでマルチノードでの行列積の高速化については本記事では触れていません。
MPIを使ってマルチノードでの行列積を行う場合にはデータ通信にかかる時間がボトルネックとなるので、
どのノードでどこの部分行列積を行うかについての調整がパフォーマンスに効いてくると思われます。
こちらの記事でSUMMAのアルゴリズムが紹介されている他、foxのアルゴリズムなど興味のある方は検索してみてください。
OpenBLASの歴史
OpenBLASとは、Basic Linear Algebra Subprograms(BLAS)と呼ばれる科学計算ライブラリが起源のオープンソースプロジェクトです。元々のBLASはFORTRANで記述された、ベクトル積・行列ベクトル積・行列積を行うためのライブラリだったようです。
そこから時代は流れ、これら以外にも逆行列計算からLU分解など様々な機能が追加されたGotoBLAS(後藤BLAS)が現れ、GotoBLAS2へと発展しました。これが後藤さんの異動により開発が中止となりました。
それをZhang Xianyiさんらが名前を変えてメンテナンスをしているものがOpenBLAS、という歴史のようです。
元々のGotoBLAS2は商用不可・再配布禁止などライセンス的には厳しかったようですが、
OpenBLASはBSDライセンスで開発されているため、何かと使いやすいのではないかと思います。
(2018/12/15 追記)
@NakataMaho さんからご指摘をいただきました。ありがとうございます。
GotoBLAS2は最終的にはオープンソースになりました。で、後藤さんのIntelへの異動につき開発は中止、OpenBLASとしてZhang Xianyiによって引き継ぎ、が正しいと思います。
探り環境
OpenBLAS
OpenBLASバージョンは現在の最新(v0.2.20)を使用し、
一部のコンパイル時のデバッグオプションを-g→-g3としたものを、
make USE_OPENMP=1 CC=gcc FC=gfortran ONLY_CBLAS=1 DEBUG=1
このオプションでビルドしたものを使います。
マシン
- CPU: Intel(R) Core(TM) i7-7600U CPU @ 2.80GHz
- コア数: 2
- Hyper Threadingは有効(1コアあたり2スレッド 合計4スレッド)
- AVX2命令までサポート
- コア数: 2
- Mem: 16GB(DDR4)
- OS:
Linux (my name) 4.9.0-6-amd64 #1 SMP Debian 4.9.88-1+deb9u1 (2018-05-07) x86_64 GNU/Linux- hugepage機能はoff(pagesizeは4096k)
- gcc, gfortran: v6.3.0
測定コード
単にcblas_sgemmを呼んでいるだけですが、興味があれば。
https://github.com/higucheese/cblas_sgemm
探りについて
扱うもの
C言語用の単精度密行列積(cblas_sgemm)の挙動を調べます。GEMMはGEneral Matrix Multiplyの略です。
- マルチスレッディングの有無
- 行列のサイズによる挙動の変化
- 転置操作の有無
扱わないもの
- 倍精度用(dgemm), 単精度複素数(cgemm,cgemm3m), 倍精度複素数(zgemm,zgemm3m)など他の型用行列積
(そういえばなぜ整数用の関数はないのでしょうか...?)
探り内容
インターフェース
C++で書かれた測定用コードからcblas_sgemm関数を呼び出します。
void MulMat::multiply() {
cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,M,N,K,alpha,A,LDA,B,LDB,beta,C,LDC);
}
この式で単精度浮動小数点数の値を持った行列$A(M,K),B(K,N),C(M,N)$と、
同じ型の係数$\alpha, \beta$に対して以下のような計算を行います。
C = \alpha A \times B + \beta C
CblasRowMajor: 配列の内容がメモリにどのように並んでいるか。これは> 通常のC言語の並びであることを示す。
CblasNoTrans: 入力する行列$A,B$が転置されていないことを示す。逆は>CblasTrans
LDA,LDB,LDC: $M,N,K$のいずれかを表す。詳細はIntel MKLの公式ページに記載されています。
するとOpenBLAS/interface/gemm.c:219に飛ばされました。
関数の先頭を見てみましょう。
void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANSPOSE TransB,
blasint m, blasint n, blasint k,
#ifndef COMPLEX
FLOAT alpha,
#else
FLOAT *alpha,
#endif
FLOAT *a, blasint lda,
FLOAT *b, blasint ldb,
#ifndef COMPLEX
FLOAT beta,
#else
FLOAT *beta,
#endif
FLOAT *c, blasint ldc) {
どうやらcblas_sgemmはCNAMEのエイリアスだったようです。
同じインターフェースCNAMEをsgemmだけではなくdgemm,c/zgemmにも使っているのでしょう。
そこから変数宣言や何の計算をするのか(ここでは単精度浮動小数点数の行列積)のモードが設定されていきます。
そして早速CblasRowMajorかCblasColMajorで分岐があります。
まずはCblasColMajor(FORTRANと同じメモリの並び順)での記述を見てみます。
argsは以降の関数に渡される構造体です。
if (order == CblasColMajor) {
args.m = m;
args.n = n;
args.k = k;
args.a = (void *)a;
args.b = (void *)b;
args.c = (void *)c;
これに対してCblasRowMajorでは...
if (order == CblasRowMajor) {
args.m = n;
args.n = m;
args.k = k;
args.a = (void *)b;
args.b = (void *)a;
args.c = (void *)c;
とこのように行列$A$と$B$を入れ替えています。
このif文内では他に転置に関しても入れ替えていて、
元々の行列積を転置して計算していることがわかります。
C^T = \alpha B^T \times A^T + \beta C^T
なぜこのようなことをするのかについてですが、
元々の実装がFortran系のCblasColMajor用に設計されており、
C系のCblasRowMajorは上記の方法を使ってその実装を使いまわしているようです。
次に部分行列に格納するためのメモリ確保が入ります(blas_memory_alloc)。
driver/others/memory.c:938に飛ばされます。
ここでは何やら部分行列を載せるためのメモリ領域に加えて、
複数回の使い回しを考慮した制御ビットと、
マルチスレッド用に読み書き制御ビットを付与した構造体を確保しているようです。
static void *alloc_mmap(void *address){
void *map_address;
if (address){
map_address = mmap(address,
BUFFER_SIZE,
MMAP_ACCESS, MMAP_POLICY | MAP_FIXED, -1, 0);
} else {
map_address = mmap(address,
BUFFER_SIZE,
MMAP_ACCESS, MMAP_POLICY, -1, 0);
}
if (map_address != (void *)-1) {
LOCK_COMMAND(&alloc_lock);
release_info[release_pos].address = map_address;
release_info[release_pos].func = alloc_mmap_free;
release_pos ++;
UNLOCK_COMMAND(&alloc_lock);
}
#ifdef OS_LINUX
my_mbind(map_address, BUFFER_SIZE, MPOL_PREFERRED, NULL, 0, 0);
#endif
return map_address;
}
common_x86_64.h:#define BUFFER_SIZE (32 << 20)
common.h:#define MMAP_POLICY (MAP_PRIVATE | MAP_ANONYMOUS)
common.h:#define MMAP_ACCESS (PROT_READ | PROT_WRITE)
去年の自分の記事では元の行列と全く同じサイズの領域を確保していました。
それに対してOpenBLASでは複数のアドレスに対してmmapを使い、
BUFFER_SIZE分のメモリを確保しています。
このBUFFER_SIZEですが上にあるとおりx86_64プロセッサだと32MiB分のサイズがあることがわかります。
更にこのmem_alloc関数は計NUM_BUFFERS回呼ばれます。
common.h:#define NUM_BUFFERS (MAX_CPU_NUMBER * 2)
NUM_BUFFERSはこのように定義されているので、どうやらOpenBLASが一回の行列積に使うメモリのサイズは、
コアの数×2×32MiBと表せそうです。
ちなみになぜコアの数の2倍確保するのかですが、これは自分の予想ですが片方を使っている間にもう片方へメモリの移し替えをやっているのではないでしょうか。つまりデータ移動と計算を同時に行うためということを考えています。
もう一つこの関数で興味深いところがあります。
それがここです。
#ifdef OS_LINUX
my_mbind(map_address, BUFFER_SIZE, MPOL_PREFERRED, NULL, 0, 0);
#endif
このmy_mbindですが、
近年流行りのNon-Uniform Memory Access (NUMA)システムを想定した上で、
プロセッサにより近い位置にあるメモリに使用するバッファをおいていると思われます。
#ifdef OS_LINUXからこの機能はLinux上でしか動きませんが、
実際にNUMAシステムを動かすときの知見が得られて面白いですね。
他にもHUGEPAGEへの対応などもされていますが、
ここではこれ以上踏み込まないことにします。
さて、メモリ確保を終えてさらに進みます。
#ifdef SMP
mode |= (transa << BLAS_TRANSA_SHIFT);
mode |= (transb << BLAS_TRANSB_SHIFT);
nthreads_max = num_cpu_avail(3);
nthreads_avail = nthreads_max;
#ifndef COMPLEX
MNK = (double) args.m * (double) args.n * (double) args.k;
if ( MNK <= (65536.0 * (double) GEMM_MULTITHREAD_THRESHOLD) )
nthreads_max = 1;
#else
MNK = (double) args.m * (double) args.n * (double) args.k;
if ( MNK <= (8192.0 * (double) GEMM_MULTITHREAD_THRESHOLD) )
nthreads_max = 1;
#endif
Makefile.system:GEMM_MULTITHREAD_THRESHOLD=4
まずSMPって何？という気持ちになりますが、
これはおそらく対称型マルチプロセッシング(Symmetric Multi Processing)の略称を表していると思われます。
今時の普通のCPU(良い表現がわかりません)であれば大抵の場合でSMPは有効と思われます。
ここでは使うスレッド数(処理をいくつに分割するか)を決めており、
ポイントは$M,N,K$の積がある一定の水準を下回ったときにはマルチスレッド化しないという点です。
なんというか..結構テキトーに決められているような気がしますね。
とはいえ動的に、もしくは自動チューニングによって決定されるべきなのかは、
利便性や実際の性能の変化も含めて考えるべきことなのかもしれません。
ようやくインターフェースも終わりが見えてきました。
(gemm[(transb << 2) | transa])(&args, NULL, NULL, sa, sb, 0);
これは使用するスレッド数が1であるときに呼ばれる関数です。
これに対してスレッド数が1よりも大きい場合はこうなります。
(gemm[16 | (transb << 2) | transa])(&args, NULL, NULL, sa, sb, 0);
gemm関数ポインタのアドレスが16増えていることがわかります。
このgemm関数は以下のように定義されています。
static int (*gemm[])(blas_arg_t *, BLASLONG *, BLASLONG *, FLOAT *, FLOAT *, BLASLONG) = {
#ifndef GEMM3M
GEMM_NN, GEMM_TN, GEMM_RN, GEMM_CN,
GEMM_NT, GEMM_TT, GEMM_RT, GEMM_CT,
GEMM_NR, GEMM_TR, GEMM_RR, GEMM_CR,
GEMM_NC, GEMM_TC, GEMM_RC, GEMM_CC,
#if defined(SMP) && !defined(USE_SIMPLE_THREADED_LEVEL3)
GEMM_THREAD_NN, GEMM_THREAD_TN, GEMM_THREAD_RN, GEMM_THREAD_CN,
GEMM_THREAD_NT, GEMM_THREAD_TT, GEMM_THREAD_RT, GEMM_THREAD_CT,
GEMM_THREAD_NR, GEMM_THREAD_TR, GEMM_THREAD_RR, GEMM_THREAD_CR,
GEMM_THREAD_NC, GEMM_THREAD_TC, GEMM_THREAD_RC, GEMM_THREAD_CC,
#endif
#else
GEMM3M_NN, GEMM3M_TN, GEMM3M_RN, GEMM3M_CN,
GEMM3M_NT, GEMM3M_TT, GEMM3M_RT, GEMM3M_CT,
GEMM3M_NR, GEMM3M_TR, GEMM3M_RR, GEMM3M_CR,
GEMM3M_NC, GEMM3M_TC, GEMM3M_RC, GEMM3M_CC,
#if defined(SMP) && !defined(USE_SIMPLE_THREADED_LEVEL3)
GEMM3M_THREAD_NN, GEMM3M_THREAD_TN, GEMM3M_THREAD_RN, GEMM3M_THREAD_CN,
GEMM3M_THREAD_NT, GEMM3M_THREAD_TT, GEMM3M_THREAD_RT, GEMM3M_THREAD_CT,
GEMM3M_THREAD_NR, GEMM3M_THREAD_TR, GEMM3M_THREAD_RR, GEMM3M_THREAD_CR,
GEMM3M_THREAD_NC, GEMM3M_THREAD_TC, GEMM3M_THREAD_RC, GEMM3M_THREAD_CC,
#endif
#endif
};
当然っちゃ当然ですがインターフェースから先に進むときに呼ばれる関数はそれぞれ異なりますね。
使用するスレッド数が1である場合はGEMM_NNに入っていくところでインターフェース部分は終了です。
スレッド数が1よりも大きいときはGEMM_THREAD_NNに入っていきます。
GEMM本体
それではGEMM本体を見ていきたいと思います。
driver/level3/level3.cにやってきました。
これ以降は転置されてこの関数に渡されている$B^T,A^T,C^T$を改めて$A,B,C$と表記したいと思います。
理由は単純にソースコード内部の表記と合わせるためです。
ちょっとここでlevel3が何を指しているのかを確認しますと、
- Level1: ベクトルの処理またはベクトル同士の演算
- Level2: ベクトルと行列同士の演算
- Level3: 行列同士の演算
となります。
以前はLevel1を組み合わせてLevel2、Level2を組み合わせてLevel3を構成するのかなと思っていたのですが、
実際にはそれぞれのLevelで専用の実装がされているようです。
ここではまず最初に上の式で出てきた$\beta C$の処理を先にやってしまいます。
その処理を行うのがこの部分です。見やすくなるように少し修正しています。
if (beta) {
#if !defined(XDOUBLE) || !defined(QUAD_PRECISION)
#ifndef COMPLEX
if (beta[0] != ONE
#else
if ((beta[0] != ONE) || (beta[1] != ZERO)
#endif
#else
if (((beta[0].x[1] != 0x3fff000000000000UL) || beta[0].x[0] != 0)
#endif
) {
BETA_OPERATION(0, m, 0, n, beta, c, ldc);
}
}
$\beta = 1$の場合は飛ばされますがそうでない場合はBETA_OPERATIONが行われます。
ここ(kernel/x86_64/gemm_beta.S)はアセンブリで記述されています。
.L202:
#ifdef OPTERON
prefetchw 32 * SIZE(C1)
#endif
MOVSD 0 * SIZE(C1), %xmm8
MOVSD 1 * SIZE(C1), %xmm9
MOVSD 2 * SIZE(C1), %xmm10
MOVSD 3 * SIZE(C1), %xmm11
MOVSD 4 * SIZE(C1), %xmm12
MOVSD 5 * SIZE(C1), %xmm13
MOVSD 6 * SIZE(C1), %xmm14
MOVSD 7 * SIZE(C1), %xmm15
MULSD %xmm0, %xmm8
MULSD %xmm0, %xmm9
MULSD %xmm0, %xmm10
MULSD %xmm0, %xmm11
MULSD %xmm0, %xmm12
MULSD %xmm0, %xmm13
MULSD %xmm0, %xmm14
MULSD %xmm0, %xmm15
MOVSD %xmm8, 0 * SIZE(C1)
MOVSD %xmm9, 1 * SIZE(C1)
MOVSD %xmm10, 2 * SIZE(C1)
MOVSD %xmm11, 3 * SIZE(C1)
MOVSD %xmm12, 4 * SIZE(C1)
MOVSD %xmm13, 5 * SIZE(C1)
MOVSD %xmm14, 6 * SIZE(C1)
MOVSD %xmm15, 7 * SIZE(C1)
addq $8 * SIZE, C1
decq I
jg .L202
ALIGN_4
やっていることはAVXのベクトル演算を使って$C$の要素すべてに$\beta$をかけているだけです。
仮に$\beta = 0$の時はBETA_OPERATIONをやらずに直接$\alpha A \times B$の結果を$C$に書き込めば良いはずですが、
元々Fused-Multiply-Add (FMA)という掛け算と足し算を同時に行う命令がAVX2にはあるので、
そのケースのためにアセンブリのパターンをもう一つ用意する必要はない(性能に大きな差はない)ということかもしれません。
ここLevel3ではいよいよ行列積をやるわけですが、
元の大きなサイズの行列積をやっても性能はメモリのバンド幅に律速されます。
計算の都合上行列の端から端までデータを読み込んで使うわけですが、
行列のサイズが大きいとこのデータも大きくなりキャッシュメモリに乗り切らなくなるためです。
こうなるとプロセッサはメモリからデータがやってくるのを延々と待つことになります。
そこで行列のタイリングをやります。これは行列を複数の部分行列に分割して計算することで、
一度にロードするメモリの量を減らしてキャッシュメモリに載っているデータの再利用を促します。
数式を使うとこんな感じになります。
A = [a_{ij}]_{m \times k},B = [b_{ij}]_{k \times n}
であるとき、
A = \left[
\begin{array}{rrrr}
A_{11} & A_{12} & \dots & A_{1k_t}\\
A_{21} & A_{22} & & A_{2k_t}\\
\vdots & & \ddots & \vdots\\
A_{m_t1} & A_{m_t2} & \dots & A_{m_tk_t}\\
\end{array}
\right]
,
B = \left[
\begin{array}{rrrr}
B_{11} & B_{12} & \dots & B_{1n_t}\\
B_{21} & B_{22} & & B_{2n_t}\\
\vdots & & \ddots & \vdots\\
B_{k_t1} & B_{k_t2} & \dots & B_{k_tn_t}\\
\end{array}
\right]
,
C_{ij} = \sum^{k_t}_k A_{i k} \times B_{k j}
では実際にコード(driver/level3/level3.c:287)を見てみます。
ここも見やすくなるようにかなり手を加えました。
完全に正しく動くコードではなく疑似コードに近いものになっています。
for(js = 0; js < n; js += GEMM_R){
min_j = GEMM_R;
for(ls = 0; ls < k; ls += min_l){
min_l = GEMM_Q;
/* First, we have to move data A to L2 cache */
min_i = GEMM_P;
ICOPY_OPERATION(min_l, min_i, a, lda, ls, 0, sa);
for(jjs = js; jjs < js + min_j; jjs += min_jj){
min_jj = GEMM_UNROLL_N;
OCOPY_OPERATION(min_l, min_jj, b, ldb, ls, jjs,
sb + min_l * (jjs - js));
KERNEL_OPERATION(min_i, min_jj, min_l, alpha,
sa, sb + min_l * (jjs - js), c, ldc, m, jjs);
}
for(is = 0 + min_i; is < m; is += min_i){
min_i = GEMM_P;
ICOPY_OPERATION(min_l, min_i, a, lda, ls, is, sa);
KERNEL_OPERATION(min_i, min_j, min_l, alpha, sa, sb, c, ldc, is, js);
} /* end of is */
} /* end of js */
} /* end of ls */
GEMM_R: 21056
GEMM_Q: 384
GEMM_P: 768
GEMM_UNROLL_N: 4
sa: Aをバッファするためのメモリの先頭アドレス
sb: Bをバッファするためのメモリの先頭アドレス
タイリングとしては$N$方向にGEMM_Rごとに、そして$K$方向にはGEMM_Qごとに、$M$方向にはGEMM_Pごとに分割しています。
そしてICOPY_OPERATIONで$A$を、OCOPY_OPERATIONで$B$をそれぞれsa,sbにロードしてきます。
KERNEL_OPERATIONはsa,sbにロードされた$A,B$を使って行列積を計算する関数になります。
最初にある程度の大きさの$A$の部分行列を持ってきたあと、
$B$の部分行列をちまちま持ってきながらKERNEL_OPERATIONを実行します。
これは$B$のロードと計算を同時に行うという意図があります。
そして$B$の部分行列がある程度溜まってきたところで、
今度は$A$の部分行列を持ってきながらKERNEL_OPERATIONを実行します。
キャッシュメモリのサイズが許す限り行列積のサイズは大きい方が、
ロードする値の個数に対して計算回数が多くなるので有利になるためです。
そうして$M$方向に一通り計算が終わったあとは$K$方向にGEMM_Qだけインデックスをずらします。
ここはlsループに対応しています。
最終的に$N$方向に計算が終われば行列積は完了となります。
試しにKERNEL_OPERATIONの中(kernel/x86_64/sgemm_kernel_16x4_haswell.S)を見てみるとこんな感じになっています。
.L4_12:
prefetcht0 A_PR1(AO, %rax, SIZE)
prefetcht0 B_PR1(BO, BI , SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
prefetcht0 B_PR1(BO, BI , SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
je .L4_16
prefetcht0 A_PR1(AO, %rax, SIZE)
prefetcht0 B_PR1(BO, BI , SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
prefetcht0 B_PR1(BO, BI , SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
prefetcht0 A_PR1(AO, %rax, SIZE)
KERNEL16x4_SUB
je .L4_16
jmp .L4_12
ALIGN_4
データのプリフェッチをしながらKERNEL16x4_SUBを実行することで、
ここでもしっかりデータのロードと計算を同時に行っています。
KERNEL16x4_SUBはマクロになっていて内部はこのようになっています。
.macro KERNEL16x4_SUB
vmovups -16 * SIZE(AO, %rax, SIZE), %ymm0
vmovups -8 * SIZE(AO, %rax, SIZE), %ymm1
vbroadcastss -4 * SIZE(BO, BI, SIZE), %ymm2
vbroadcastss -3 * SIZE(BO, BI, SIZE), %ymm3
VFMADD231PS_( %ymm4,%ymm2,%ymm0 )
VFMADD231PS_( %ymm5,%ymm2,%ymm1 )
VFMADD231PS_( %ymm6,%ymm3,%ymm0 )
VFMADD231PS_( %ymm7,%ymm3,%ymm1 )
vbroadcastss -2 * SIZE(BO, BI, SIZE), %ymm2
vbroadcastss -1 * SIZE(BO, BI, SIZE), %ymm3
VFMADD231PS_( %ymm8,%ymm2,%ymm0 )
VFMADD231PS_( %ymm9,%ymm2,%ymm1 )
VFMADD231PS_( %ymm10,%ymm3,%ymm0 )
VFMADD231PS_( %ymm11,%ymm3,%ymm1 )
addq $ 4 , BI
addq $ 16, %rax
.endm
この計算は例えるなら8要素のベクトル$A$に、1要素の$B$をそれぞれにかけ合わせて8要素のベクトル$C$に足し合わせる、といったようなことをしています。
AVX2は256bit幅のベクトル演算をサポートしていて、単精度浮動小数点は32bit(=4Byte)なので$256/32=8$として8という数字が出ています。
このように行列積の奥の奥ではベクトル演算器をフルに使うような実装をしていますが、
全体的にはどのようにメモリのロードと計算を同時に行うかもポイントになっています。
マルチスレッドでは
上記のシングルスレッド用の実装では$M$は0からmまでとなっていましたが、
これを使えるスレッド数で分割します。
例えばスレッドの数が4でmが4で割り切れる場合は、
それぞれのスレッドで$M$方向のレンジが[0,m/4],[m/4,m/2],[m/2,3m/4],[3m/4,m]となります。
おわりに
結構端折って説明したのでどれほどわかっていただけるか不安ですが、
全部の関数をたどるのでなければOpenBLASはかなり読みやすい(アセンブリに抵抗が少なければ)ので、
自分で探ってみるのはいかがでしょうか？(やらないか)
ハードウェアレベルのゴリゴリのチューニングが好きな人に特におすすめです。
最初の方に書いたTipsの理由付けですが、
十分なサイズの行列積の場合は常に物理コア(=実際に計算するところ)が忙しいため、
Intelのハイパースレッディング技術で論理コアを増やしても意味がないか、
もしくはパフォーマンスに悪影響がでるためです。
ちなみにここまでの話はすべてOpenMPを用いて並列化を行っていることを前提にしていました。
勝手に並列化してくれるような科学計算ライブラリであれば、
大抵の場合でOpenMPを使ってくれているとは思いますが、
Pthreadなどを使っていた場合はあれらの環境変数を設定しても意味はないはずです。
以上ありがとうございました。