密行列積の高速化 OpenBLAS編

# はじめに
どうもアドベントカレンダーとしてはかなりの分野違いをしてしまっているようです。
皆さんが複雑な数式を紹介する中で、本記事は実数を要素に持つ行列`A,B,C`に対して

```math
C = A \times B
```

これしか計算していない上に各行列が何を表しているのかについては不問です。
その代わりこの式をどのように計算機上で高速に計算をするのかについて紹介しています。
「なんか泥臭いことやっとるな」くらいの認識でOKです。

# OpenBLASに関するTips的な
早速OpenBLASを使用する際のTipsについて簡単に説明したいと思います。
「え、多次元配列の処理はjuliaでやるからOpenBLAS使わないよ」と言う方、
[公式githubのページ](https://github.com/JuliaLang/julia)の"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なので以下のような結果が得られました。

```sh
 % 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における行列積の高速化手段について調べようと思います。
内容としては[去年の自分の記事](https://qiita.com/higucheese/items/062591dbad8d15cdef06)で記述したものの答え合わせ的なものになっています。
なのでマルチノードでの行列積の高速化については本記事では触れていません。

MPIを使ってマルチノードでの行列積を行う場合にはデータ通信にかかる時間がボトルネックとなるので、
どのノードでどこの部分行列積を行うかについての調整がパフォーマンスに効いてくると思われます。
[こちらの記事](https://qiita.com/kanji_/items/495ea6fead99973cc0c4)でSUMMAのアルゴリズムが紹介されている他、foxのアルゴリズムなど興味のある方は検索してみてください。

## OpenBLASの歴史
[OpenBLAS](https://www.openblas.net/)とは、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`としたものを、

```shell
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命令までサポート
- 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`関数を呼び出します。

```cpp
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$に対して以下のような計算を行います。

```math
C = \alpha A \times B + \beta C
```

> `CblasRowMajor`: 配列の内容がメモリにどのように並んでいるか。これは> 通常のC言語の並びであることを示す。
> `CblasNoTrans`: 入力する行列$A,B$が転置されていないことを示す。逆は> `CblasTrans`
> `LDA,LDB,LDC`: $M,N,K$のいずれかを表す。詳細は[Intel MKLの公式ページ](https://software.intel.com/en-us/mkl-developer-reference-c-cblas-gemm)に記載されています。

すると`OpenBLAS/interface/gemm.c:219`に飛ばされました。
関数の先頭を見てみましょう。

```c
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`は以降の関数に渡される構造体です。

```c
  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`では...

```c
  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文内では他に転置に関しても入れ替えていて、
元々の行列積を転置して計算していることがわかります。

```math
C^T = \alpha B^T \times A^T + \beta C^T
```

なぜこのようなことをするのかについてですが、
元々の実装がFortran系の`CblasColMajor`用に設計されており、
C系の`CblasRowMajor`は上記の方法を使ってその実装を使いまわしているようです。

次に部分行列に格納するためのメモリ確保が入ります(`blas_memory_alloc`)。
`driver/others/memory.c:938`に飛ばされます。
ここでは何やら部分行列を載せるためのメモリ領域に加えて、
複数回の使い回しを考慮した制御ビットと、
マルチスレッド用に読み書き制御ビットを付与した構造体を確保しているようです。

```c
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倍確保するのかですが、これは自分の予想ですが片方を使っている間にもう片方へメモリの移し替えをやっているのではないでしょうか。つまりデータ移動と計算を同時に行うためということを考えています。

もう一つこの関数で興味深いところがあります。
それがここです。

```c
#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への対応などもされていますが、
ここではこれ以上踏み込まないことにします。

さて、メモリ確保を終えてさらに進みます。

```c
#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$の積がある一定の水準を下回ったときにはマルチスレッド化しないという点です。
なんというか..結構テキトーに決められているような気がしますね。
とはいえ動的に、もしくは自動チューニングによって決定されるべきなのかは、
利便性や実際の性能の変化も含めて考えるべきことなのかもしれません。

ようやくインターフェースも終わりが見えてきました。

```c
    (gemm[(transb << 2) | transa])(&args, NULL, NULL, sa, sb, 0);
```

これは使用するスレッド数が1であるときに呼ばれる関数です。
これに対してスレッド数が1よりも大きい場合はこうなります。

```c
	(gemm[16 | (transb << 2) | transa])(&args, NULL, NULL, sa, sb, 0);
```

`gemm`関数ポインタのアドレスが16増えていることがわかります。
この`gemm`関数は以下のように定義されています。

```c
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$の処理を先にやってしまいます。
その処理を行うのがこの部分です。見やすくなるように少し修正しています。

```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`)はアセンブリで記述されています。

```asm
.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ではいよいよ行列積をやるわけですが、
元の大きなサイズの行列積をやっても性能はメモリのバンド幅に律速されます。
計算の都合上行列の端から端までデータを読み込んで使うわけですが、
行列のサイズが大きいとこのデータも大きくなりキャッシュメモリに乗り切らなくなるためです。
こうなるとプロセッサはメモリからデータがやってくるのを延々と待つことになります。

そこで行列のタイリングをやります。これは行列を複数の部分行列に分割して計算することで、
一度にロードするメモリの量を減らしてキャッシュメモリに載っているデータの再利用を促します。
数式を使うとこんな感じになります。

```math
A = [a_{ij}]_{m \times k},B = [b_{ij}]_{k \times n}
```

であるとき、

```math
  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`)を見てみます。
ここも見やすくなるようにかなり手を加えました。
完全に正しく動くコードではなく疑似コードに近いものになっています。

```c
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`)を見てみるとこんな感じになっています。

```asm
.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`はマクロになっていて内部はこのようになっています。

```asm
.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などを使っていた場合はあれらの環境変数を設定しても意味はないはずです。

以上ありがとうございました。