<- ホーム
ループタイリングをちゃんと理解したい
2026/07/18
ループタイリングの話します。
ループタイリングは割と基本的な最適化テクニックのような顔をして、実際はかなり理解するのが難しいと思っている。
調べるといきなり最適化前と最適化後のコードが出てくるので、それを読んでなんとなく理解したつもりになってしまいがちだ。
自分もそうであるという自覚を持ってきたので、ちゃんと理解しようと試みる。
まず、ループタイリングの定義としては、狭義と広義があると思っており、狭義の定義としては単なるループ分割のこと。
そして、広義の定義としてはループ分割を利用してCPUキャッシュを最適化する手法のことだ。
ただ、タイリングの「タイル」とは、あれは風呂の壁のタイルとかと同じ意味のタイルなんだが、文字通りメモリをタイルの集合と捉え、タイルごとにアクセスすることでCPUのキャッシュラインを可能な限り再利用し続ける、という意味がこめられていると思われるので、本来広義の意味のほうが正しいと思う。
ただいくつかのソフトウェア (Halideとか) はループ分割とループの順序変更というプリミティブ操作の組み合わせのことをタイルと呼んでいたりするので、そういう意味で狭義の意味合いもあるのかな、という感じ。
なので、まずはループ分割について。
本当なら先にループタイリングが解決する問題の話からすべきではあるのだが、ループタイリングは難しいので、いきなり問題の話からしてもかえってわかりにくいというか、問題と解決策を結びつけて考えるのが個人的には難しいと感じているので、手法の話からする。
ループ分割とは、文字通りループを分割することだ。ここでいう分割というのは実際にはネストのことで、同レベルで分割するわけではない。
例えば、32回ループしたいコードがあって、分割前のコードが以下だとすると、
for (x = 0; x < 32; x++) {
f(x);
}
分割するとこうなる。
factor = 8;
for (xo = 0; xo < 32; xo += factor) {
for (xi = 0; xi < factor; xi++) {
x = xo + xi
f(x)
}
}
これは32回ループする x について、4回ループする外側ループ xo と、8回ループする内側ループ xi に分割している。
xoは xo += factor で増えていくので、0, 8, 16, 24と増えていく。内側のループは0, 1, 2, 3, 4, 5, 6, 7と回る。
感覚的には、普通の8回ループ (内側ループ) を4周 (外側ループ) やっているから、合計回数は32で変わらない、というイメージだ。
xはxoとxiを足せば復元できるのもわかるだろう。
ループ分割とはこれだけのことで、これだとなんの意味が?という感じになるだろう。それはそのとおりで、つまりループ分割自体は最適化ではないのだ。最適化のための準備でしかない。例えば、32個の要素をSIMDで処理したいが、SIMDレーンは8個しかないから8個ずつに分けたいというときに、32個のループを4 * 8の2重ループにして、内側ループをSIMDで書いたりする。同じようにすればマルチスレッド化できたり、内側ループだけをアンロールしたりできる。大量データをマルチスレッドで処理したい時、まずデータをコア数分のチャンクに分けたりすることがあるが、あれと同じである。
さて、ループ分割がわかったので、ループタイリングの話ができる。
ループタイリングとは、ループ分割を巧みに使って、プログラムの挙動を変えずに、メモリにアクセスする範囲を限定しながらループを進めることで、CPUキャッシュをより効率よく使う手法である、という説明が最も私の中でしっくり来ている。
コードを出す前にイメージを作ったほうが良いのだが、例えば二重ループがあり、この二重ループによってメモリを飛び飛びにアクセスしてしまうような状況があるとする。
飛び飛びにアクセスするというのがどういうことか説明する前に、例えば普通の以下のようなループがあったとき、
v = 0
for (x = 0; x < 32; x++) {
v = memory[x]
}
これは、メモリの0番地、1番地、2番地、、とシーケンシャルにアクセスしていく。
このとき、例えば1回目のループでは memory[0] にアクセスするわけだが、CPUというのは本当にメモリの0番地だけを取りに行くのではなく、キャッシュラインという単位でまとめてデータを取ってきて、それをCPU側のキャッシュに載せておく。なので、2周目以降のループでは、データがキャッシュにある限りメモリアクセス不要である。ループを回していくといずれ、キャッシュラインにないデータアクセスが、キャッシュラインの大きさの限界を超えたタイミングで発生するがそれはしょうがない。
CPU側のキャッシュは計算機システムにとって、メモリアクセスよりも数桁速いという性質があるため、高速化には必要不可欠である。
このような単純なシーケンシャルループはCPUキャッシュ的に非常に有利なプログラムで、実際世の中にはこういうプログラム、つまり隣接したメモリに次々にアクセスしていくようなもののほうが多いので、キャッシュはプログラマが意識しなくても有効に効いてくれる。
ただ世の中にはこういう感じでもない問題があり、その典型例が行列積だ。飛び飛びアクセスが発生する。どういうプログラムかというと、まず2次元配列AとBがあるとして、
A = [
[a00, a01, a02, a03]
[a10, a11, a12, a13]
[a20, a21, a22, a23]
[a30, a31, a32, a33]
]
B = [
[b00, b01, b02, b03]
[b10, b11, b12, b13]
[b20, b21, b22, b23]
[b30, b31, b32, b33]
]
AとBの中身は何でも良いのだが、行列積Cは、
C = [0] * 16
for i = 0; i < 4; i++ {
for j = 0; j < 4; j++ {
for k = 0; k < 4; k++ {
C[i][j] += A[i][k] * B[k][j]
}
}
}
こういう計算でできる。行列積というのはどういう計算かというと、「Cのi行j列は、Aのi行目とBのj列目をベクトルの掛け算してsumを取ったやつ」というもので、これは行列積の定義とも一致する。最内のkループはそれをやっていて、これはもっとわかりやすくすると
import numpy as np
A = np.arange(16).reshape(4, 4)
B = np.arange(16).reshape(4, 4)
C = np.zeros((4, 4))
for i in range(4):
for j in range(4):
C[i][j] = (A[i, :] * B[:, j]).sum()
こういうPythonでも同じである。 C[i][j] = (A[i, :] * B[:, j]).sum() がまさに、「Cのi行j列が、Aのi行目とBのj列目をベクトルの掛け算してsumを取る」をやっている。
つまり、先程の三重ループは、外側のiとjは行と列のループであり、最内のkループは、あくまでベクトルの掛け算 -> sumをしたいのだが、ベクトルの掛け算はnumpyなしではループが必要になってしまうのでしかたなくkループを用意しているという話であり、本質的には二重ループのプログラムに近い。kループはベクトルに対して要素ごとにアクセスして掛け算し、 += によってsumを実現している。
で、この三重ループの方のプログラムに着目した場合、Aについてはシーケンシャルアクセスだが、Bについては飛び飛びのアクセスが起こっているのがわかるだろうか?
最内のkループに着目した時、Aについては A[i][k] とアクセスしているが、Bについては B[k][j] とアクセスしている。
ここをわかりやすくするために、i, j, kが実際にはどのような値を取りながらループするのかを明確にすると、次のようになる。
- ループ1 : i = 0, j = 0, k = 0
- ループ2 : i = 0, j = 0, k = 1
- ループ3 : i = 0, j = 0, k = 2
- ループ4 : i = 0, j = 0, k = 3
- ループ5 : i = 0, j = 1, k = 0
- ループ6 : i = 0, j = 1, k = 1
- ループ7 : i = 0, j = 1, k = 2
- ループ8 : i = 0, j = 1, k = 3
- ループ9 : i = 0, j = 2, k = 0
- ループ10: i = 0, j = 2, k = 1
- ループ11: i = 0, j = 2, k = 2
- ループ12: i = 0, j = 2, k = 3
- ループ13: i = 0, j = 3, k = 0
- ループ14: i = 0, j = 3, k = 1
- ループ15: i = 0, j = 3, k = 2
- ループ16: i = 0, j = 3, k = 3
- ループ17: i = 1, j = 0, k = 0
- ループ18: i = 1, j = 0, k = 1
- ループ19: i = 1, j = 0, k = 2
- ループ20: i = 1, j = 0, k = 3
- ループ21: i = 1, j = 1, k = 0
- ループ22: i = 1, j = 1, k = 1
- ループ23: i = 1, j = 1, k = 2
- ループ24: i = 1, j = 1, k = 3
- ループ25: i = 1, j = 2, k = 0
- ループ26: i = 1, j = 2, k = 1
- ループ27: i = 1, j = 2, k = 2
- ループ28: i = 1, j = 2, k = 3
- ループ29: i = 1, j = 3, k = 0
- ループ30: i = 1, j = 3, k = 1
- ループ31: i = 1, j = 3, k = 2
- ループ32: i = 1, j = 3, k = 3
- ループ33: i = 2, j = 0, k = 0
- ループ34: i = 2, j = 0, k = 1
- ループ35: i = 2, j = 0, k = 2
- ループ36: i = 2, j = 0, k = 3
- ループ37: i = 2, j = 1, k = 0
- ループ38: i = 2, j = 1, k = 1
- ループ39: i = 2, j = 1, k = 2
- ループ40: i = 2, j = 1, k = 3
- ループ41: i = 2, j = 2, k = 0
- ループ42: i = 2, j = 2, k = 1
- ループ43: i = 2, j = 2, k = 2
- ループ44: i = 2, j = 2, k = 3
- ループ45: i = 2, j = 3, k = 0
- ループ46: i = 2, j = 3, k = 1
- ループ47: i = 2, j = 3, k = 2
- ループ48: i = 2, j = 3, k = 3
- ループ49: i = 3, j = 0, k = 0
- ループ50: i = 3, j = 0, k = 1
- ループ51: i = 3, j = 0, k = 2
- ループ52: i = 3, j = 0, k = 3
- ループ53: i = 3, j = 1, k = 0
- ループ54: i = 3, j = 1, k = 1
- ループ55: i = 3, j = 1, k = 2
- ループ56: i = 3, j = 1, k = 3
- ループ57: i = 3, j = 2, k = 0
- ループ58: i = 3, j = 2, k = 1
- ループ59: i = 3, j = 2, k = 2
- ループ60: i = 3, j = 2, k = 3
- ループ61: i = 3, j = 3, k = 0
- ループ62: i = 3, j = 3, k = 1
- ループ63: i = 3, j = 3, k = 2
- ループ64: i = 3, j = 3, k = 3
となるのだが、ここでAへの A[i][k] アクセスについて考えると、ループ1からループ16まで、 A[0][0] A[0][1] A[0][2] A[0][3] へのアクセスを繰り返しているのがわかるだろうか?
そして、ループ17から32までは、A[1][0] A[1][1] A[1][2] A[1][3] へのアクセスを繰り返している。それ以降も同様だ。
これはシーケンシャルアクセスである。なぜかと言うと、Cやnumpyなどは多次元配列を次のようにメモリに配置するからだ。
プログラム上ではこうだが、
[
[a, b, c, d]
[e, f, g, h]
[i, j, k, l]
[m, n, o, p]
]
メモリ上ではこう
[a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p]
これは「Row-Major」と呼ばれる配置方法だ。これは多次元配列を線形のストレージ (メモリなど) に配置する方法を示している。
「Column-Major」という列優先の配置方法もあり、その場合は次のようにメモリ上に配置される。
[a, e, i, m, b, f, j, n, c, g, k, o, d, h, l, p]
縦にというか、1列目、2列目、、と配置されている。
CやC++、numpyのndarray 、PyTorchのTensorなどはRow Majorで、FortranやMATLAB、RなどはColumn Majorを採用している。ただしnumpyなどは配列を作るときに order 引数を指定すればColumn Majorにもできるのだが、実際やることはあまりないと思われる。
今更ながら、この文章はRow-Majorを前提としている。
で、Row-Majorを採用しているとした場合、 A[0][0] A[0][1] A[0][2] A[0][3] に繰り返しアクセスするのは、Aの0行目に連続アクセスしており、0行目はメモリ上連続に並んでいるので、この連続したエリアはキャッシュラインに乗りやすく、したがってAへのアクセスはキャッシュを有効活用しやすい。
17ループ目からはAの1行目にアクセスすることになるのだが、0行目の最後の要素と1行目の最初の要素は連続しているので、つまりAへのアクセスはずっとシーケンシャルアクセスになり、これはかなり嬉しい。
Bの場合話が変わってくる。Bだと B[k][j] とアクセスするので、実際のアクセスは、ループ1から8を見ると、次のようになっている。
- ループ1 : B[0][0]
- ループ2 : B[1][0]
- ループ3 : B[2][0]
- ループ4 : B[3][0]
- ループ5 : B[0][1]
- ループ6 : B[1][1]
- ループ7 : B[2][1]
- ループ8 : B[3][1]
0行目、1行目、2行目、3行目、0行目、1行目、、と進んでいるのがわかる。これ以降も同様である。これだと、Row-Majorにおいては、メモリの0番地 -> 4番地 -> 8番地 -> 12番地 -> 1番地 -> 5番地、、と進んでいくことになり、シーケンシャルでない飛び飛びのアクセス (ストライドアクセスなどという) になる。
キャッシュラインのサイズはCPUによる。ここからの話は、次のようなシステムを前提としよう。
- このシステムのキャッシュラインは、 1つの行列につき 8つの要素を格納できる
- キャッシュラインは1つしかない
- 1つの行列がキャッシュにない要素にアクセスすると、メモリアクセスが発生し、キャッシュラインの中身が入れ替えられる
最近の計算機なら通常キャッシュラインは複数あるので、これはかなり単純化したモデルになっているが、それ以外は特に間違ってはいない。
で、この場合、先のBへのアクセスは次のようになる。まず、最初に0番地にアクセスしたので、0, 1, 2, 3, 4, 5, 6, 7番地までをメモリから取ってきてキャッシュに乗せる。次に4番地にアクセスが来たが、4番地はキャッシュに乗っているのでok。しかし、次の8番地はメモリに乗っていないので、キャッシュを破棄し、8, 9, 10, 11, 12, 13, 14, 15番地をメモリから取得、キャッシュに乗せ、、となり、キャッシュが全く使えていないわけではないものの、もっと使いたい、困った。というのが問題意識である。
この前提で考えると、Bへのメモリアクセスは何回発生するかと言うと、連続する2回のループならキャッシュに載る、つまり奇数回目のループでキャッシュミスするので、32回メモリアクセスしている。
これをループ分割で解決できるのだが、考え方としては、行列をタイルとして分割して、そのタイルに集中してアクセスすることで、アクセス範囲を減らして、タイル内のデータが確実にキャッシュラインに載せられるようにするというものだ。
もっかいおさらいすると、行列Bは次のようになっている。
B = [
[b00, b01, b02, b03]
[b10, b11, b12, b13]
[b20, b21, b22, b23]
[b30, b31, b32, b33]
]
で、先のプログラムでは b00 , b10 , b20 , b30 、、とアクセスしていくのだが、 b00 にアクセスしたときにキャッシュに載るのは b00 b01 b02 b03 b10 b11 b12 b13 なので、 b20 のタイミングでメモリアクセスが発生する。
ここでタイリングとはどういうものかというと、4x4の行列を2x2の行列に分ける。左上、右上、左下、右下の4つの「タイル」に分け、タイルごとに処理をするのだ。このへんが最高にややこしいのでちゃんと考える。
例えば、Cの左上タイルをまず考えると、Cの左上を計算するにはAの左上・右上とBの左上・左下が必要だ。そこで、a00, b00にまずアクセスして、AとBの左上・右上タイルをキャッシュする。
すると、必要なタイルのうちキャッシュに足りないのはBの左下だけということになる。なので、まずAの左上、Bの左上を使って、Cの左上を途中まで計算する。これだけだと計算が途中で終わっているので、次にb20にアクセスする。すると、Bの左下、右下がキャッシュされる。これによって、Aの右上、Bの左下を使い、Cの左上の計算が終えられる。
次はCの右上だ。Cの右上を計算するにはAの左上・右上と、Bの右上・右下が必要だ。Aの上は既にキャッシュされているので、b00にアクセスし、Bの上をキャッシュする。そして、Aの左上とBの右上からCの右上を途中まで計算する。次に、b20にアクセスしてBの下をキャッシュし、Aの右上・Bの右下からCの右上の計算を終える。
このように、Cの左下、Cの右下、、と計算していくのだが、これ、タイリングしない場合と考え方は同じである。c00, c01, c02, c03、ここまで来たら次の行に移って、、とするのを、タイルごとにまとめてやっているだけである。これでちゃんとメモリアクセスが減るので、ようやく実際のコードを見てみよう。
T = 2; // タイルサイズ
// 外側ループ
for io = 0; io < 4; io+=T {
for jo = 0; jo < 4; jo+=T {
for ko = 0; ko < 4; ko+=T {
// 内側ループ
for ii = io; ii < min(4, io+T); ii++ {
for ji = jo; ji < min(4, jo+T); ji++ {
for ki = ko; ki < min(4, ko+T); ki++ {
C[ii][ji] += A[ii][ki] * B[ki][ji]
}
}
}
}
}
}
こうである。外側ループを分割して、io, jo, koは[0, 2]を取る。内側ループのii, ji, kiは外側が0のとき[0, 1]、2のとき[2, 3]を取る。ループ変数を全て列挙すると次のようになる。
- ループ1 : io = 0, jo = 0, ko = 0, ii = 0, ji = 0, ki = 0
- ループ2 : io = 0, jo = 0, ko = 0, ii = 0, ji = 0, ki = 1
- ループ3 : io = 0, jo = 0, ko = 0, ii = 0, ji = 1, ki = 0
- ループ4 : io = 0, jo = 0, ko = 0, ii = 0, ji = 1, ki = 1
- ループ5 : io = 0, jo = 0, ko = 0, ii = 1, ji = 0, ki = 0
- ループ6 : io = 0, jo = 0, ko = 0, ii = 1, ji = 0, ki = 1
- ループ7 : io = 0, jo = 0, ko = 0, ii = 1, ji = 1, ki = 0
- ループ8 : io = 0, jo = 0, ko = 0, ii = 1, ji = 1, ki = 1
- ループ9 : io = 0, jo = 0, ko = 2, ii = 0, ji = 0, ki = 2
- ループ10: io = 0, jo = 0, ko = 2, ii = 0, ji = 0, ki = 3
- ループ11: io = 0, jo = 0, ko = 2, ii = 0, ji = 1, ki = 2
- ループ12: io = 0, jo = 0, ko = 2, ii = 0, ji = 1, ki = 3
- ループ13: io = 0, jo = 0, ko = 2, ii = 1, ji = 0, ki = 2
- ループ14: io = 0, jo = 0, ko = 2, ii = 1, ji = 0, ki = 3
- ループ15: io = 0, jo = 0, ko = 2, ii = 1, ji = 1, ki = 2
- ループ16: io = 0, jo = 0, ko = 2, ii = 1, ji = 1, ki = 3
- ループ17: io = 0, jo = 2, ko = 0, ii = 0, ji = 2, ki = 0
- ループ18: io = 0, jo = 2, ko = 0, ii = 0, ji = 2, ki = 1
- ループ19: io = 0, jo = 2, ko = 0, ii = 0, ji = 3, ki = 0
- ループ20: io = 0, jo = 2, ko = 0, ii = 0, ji = 3, ki = 1
- ループ21: io = 0, jo = 2, ko = 0, ii = 1, ji = 2, ki = 0
- ループ22: io = 0, jo = 2, ko = 0, ii = 1, ji = 2, ki = 1
- ループ23: io = 0, jo = 2, ko = 0, ii = 1, ji = 3, ki = 0
- ループ24: io = 0, jo = 2, ko = 0, ii = 1, ji = 3, ki = 1
- ループ25: io = 0, jo = 2, ko = 2, ii = 0, ji = 2, ki = 2
- ループ26: io = 0, jo = 2, ko = 2, ii = 0, ji = 2, ki = 3
- ループ27: io = 0, jo = 2, ko = 2, ii = 0, ji = 3, ki = 2
- ループ28: io = 0, jo = 2, ko = 2, ii = 0, ji = 3, ki = 3
- ループ29: io = 0, jo = 2, ko = 2, ii = 1, ji = 2, ki = 2
- ループ30: io = 0, jo = 2, ko = 2, ii = 1, ji = 2, ki = 3
- ループ31: io = 0, jo = 2, ko = 2, ii = 1, ji = 3, ki = 2
- ループ32: io = 0, jo = 2, ko = 2, ii = 1, ji = 3, ki = 3
- ループ33: io = 2, jo = 0, ko = 0, ii = 2, ji = 0, ki = 0
- ループ34: io = 2, jo = 0, ko = 0, ii = 2, ji = 0, ki = 1
- ループ35: io = 2, jo = 0, ko = 0, ii = 2, ji = 1, ki = 0
- ループ36: io = 2, jo = 0, ko = 0, ii = 2, ji = 1, ki = 1
- ループ37: io = 2, jo = 0, ko = 0, ii = 3, ji = 0, ki = 0
- ループ38: io = 2, jo = 0, ko = 0, ii = 3, ji = 0, ki = 1
- ループ39: io = 2, jo = 0, ko = 0, ii = 3, ji = 1, ki = 0
- ループ40: io = 2, jo = 0, ko = 0, ii = 3, ji = 1, ki = 1
- ループ41: io = 2, jo = 0, ko = 2, ii = 2, ji = 0, ki = 2
- ループ42: io = 2, jo = 0, ko = 2, ii = 2, ji = 0, ki = 3
- ループ43: io = 2, jo = 0, ko = 2, ii = 2, ji = 1, ki = 2
- ループ44: io = 2, jo = 0, ko = 2, ii = 2, ji = 1, ki = 3
- ループ45: io = 2, jo = 0, ko = 2, ii = 3, ji = 0, ki = 2
- ループ46: io = 2, jo = 0, ko = 2, ii = 3, ji = 0, ki = 3
- ループ47: io = 2, jo = 0, ko = 2, ii = 3, ji = 1, ki = 2
- ループ48: io = 2, jo = 0, ko = 2, ii = 3, ji = 1, ki = 3
- ループ49: io = 2, jo = 2, ko = 0, ii = 2, ji = 2, ki = 0
- ループ50: io = 2, jo = 2, ko = 0, ii = 2, ji = 2, ki = 1
- ループ51: io = 2, jo = 2, ko = 0, ii = 2, ji = 3, ki = 0
- ループ52: io = 2, jo = 2, ko = 0, ii = 2, ji = 3, ki = 1
- ループ53: io = 2, jo = 2, ko = 0, ii = 3, ji = 2, ki = 0
- ループ54: io = 2, jo = 2, ko = 0, ii = 3, ji = 2, ki = 1
- ループ55: io = 2, jo = 2, ko = 0, ii = 3, ji = 3, ki = 0
- ループ56: io = 2, jo = 2, ko = 0, ii = 3, ji = 3, ki = 1
- ループ57: io = 2, jo = 2, ko = 2, ii = 2, ji = 2, ki = 2
- ループ58: io = 2, jo = 2, ko = 2, ii = 2, ji = 2, ki = 3
- ループ59: io = 2, jo = 2, ko = 2, ii = 2, ji = 3, ki = 2
- ループ60: io = 2, jo = 2, ko = 2, ii = 2, ji = 3, ki = 3
- ループ61: io = 2, jo = 2, ko = 2, ii = 3, ji = 2, ki = 2
- ループ62: io = 2, jo = 2, ko = 2, ii = 3, ji = 2, ki = 3
- ループ63: io = 2, jo = 2, ko = 2, ii = 3, ji = 3, ki = 2
- ループ64: io = 2, jo = 2, ko = 2, ii = 3, ji = 3, ki = 3
ループ1からループ8に着目してほしいのだが、実際の計算は次のようになっている。
- ループ1 : C[0][0] += A[0][0] * B[0][0]
- ループ2 : C[0][0] += A[0][1] * B[1][0]
- ループ3 : C[0][1] += A[0][0] * B[0][1]
- ループ4 : C[0][1] += A[0][1] * B[1][1]
- ループ5 : C[1][0] += A[1][0] * B[0][0]
- ループ6 : C[1][0] += A[1][1] * B[1][0]
- ループ7 : C[1][1] += A[1][0] * B[0][1]
- ループ8 : C[1][1] += A[1][1] * B[1][1]
これを見るとわかるのだが、A、B、C全て0行0列、0行1列、1行0列、1行1列の4つ、つまり左上タイルにしか触っていない。ループ1のタイミングで行列ごとに2行目まではキャッシュできるため、
この8回のループでは全ての要素がキャッシュに載っており、メモリアクセスは最初の1回しか発生しない。これは、「AとBの左上からCの左上を途中まで計算している」状態だ。どの計算でどこをどう計算しているかを全て列挙すると、
- Cの左上 (前半) : Aの左上 * Bの左上
- ループ1 : C[0][0] += A[0][0] * B[0][0]
- ループ2 : C[0][0] += A[0][1] * B[1][0]
- ループ3 : C[0][1] += A[0][0] * B[0][1]
- ループ4 : C[0][1] += A[0][1] * B[1][1]
- ループ5 : C[1][0] += A[1][0] * B[0][0]
- ループ6 : C[1][0] += A[1][1] * B[1][0]
- ループ7 : C[1][1] += A[1][0] * B[0][1]
- ループ8 : C[1][1] += A[1][1] * B[1][1]
- Cの左上 (後半) : Aの右上 * Bの左下
- ループ9 : C[0][0] += A[0][2] * B[2][0]
- ループ10: C[0][0] += A[0][3] * B[3][0]
- ループ11: C[0][1] += A[0][2] * B[2][1]
- ループ12: C[0][1] += A[0][3] * B[3][1]
- ループ13: C[1][0] += A[1][2] * B[2][0]
- ループ14: C[1][0] += A[1][3] * B[3][0]
- ループ15: C[1][1] += A[1][2] * B[2][1]
- ループ16: C[1][1] += A[1][3] * B[3][1]
- Cの右上 (前半) : Aの左上 * Bの右上
- ループ17: C[0][2] += A[0][0] * B[0][2]
- ループ18: C[0][2] += A[0][1] * B[1][2]
- ループ19: C[0][3] += A[0][0] * B[0][3]
- ループ20: C[0][3] += A[0][1] * B[1][3]
- ループ21: C[1][2] += A[1][0] * B[0][2]
- ループ22: C[1][2] += A[1][1] * B[1][2]
- ループ23: C[1][3] += A[1][0] * B[0][3]
- ループ24: C[1][3] v= A[1][1] * B[1][3]
- Cの右上 (後半) : Aの右上 * Bの右下
- ループ25: C[0][2] += A[0][2] * B[2][2]
- ループ26: C[0][2] += A[0][3] * B[3][2]
- ループ27: C[0][3] += A[0][2] * B[2][3]
- ループ28: C[0][3] += A[0][3] * B[3][3]
- ループ29: C[1][2] += A[1][2] * B[2][2]
- ループ30: C[1][2] += A[1][3] * B[3][2]
- ループ31: C[1][3] += A[1][2] * B[2][3]
- ループ32: C[1][3] += A[1][3] * B[3][3]
- Cの左下 (前半) : Aの左下 * Bの左上
- ループ33: C[2][0] += A[2][0] * B[0][0]
- ループ34: C[2][0] += A[2][1] * B[1][0]
- ループ35: C[2][1] += A[2][0] * B[0][1]
- ループ36: C[2][1] += A[2][1] * B[1][1]
- ループ37: C[3][0] += A[3][0] * B[0][0]
- ループ38: C[3][0] += A[3][1] * B[1][0]
- ループ39: C[3][1] += A[3][0] * B[0][1]
- ループ40: C[3][1] += A[3][1] * B[1][1]
- Cの左下 (後半) : Aの右下 * Bの左下
- ループ41: C[2][0] += A[2][2] * B[2][0]
- ループ42: C[2][0] += A[2][3] * B[3][0]
- ループ43: C[2][1] += A[2][2] * B[2][1]
- ループ44: C[2][1] += A[2][3] * B[3][1]
- ループ45: C[3][0] += A[3][2] * B[2][0]
- ループ46: C[3][0] += A[3][3] * B[3][0]
- ループ47: C[3][1] += A[3][2] * B[2][1]
- ループ48: C[3][1] += A[3][3] * B[3][1]
- Cの右下 (前半) : Aの左下 * Bの右上
- ループ49: C[2][2] += A[2][0] * B[0][2]
- ループ50: C[2][2] += A[2][1] * B[1][2]
- ループ51: C[2][3] += A[2][0] * B[0][3]
- ループ52: C[2][3] += A[2][1] * B[1][3]
- ループ53: C[3][2] += A[3][0] * B[0][2]
- ループ54: C[3][2] += A[3][1] * B[1][2]
- ループ55: C[3][3] += A[3][0] * B[0][3]
- ループ56: C[3][3] += A[3][1] * B[1][3]
- Cの右下 (後半) : Aの右下 * Bの右下
- ループ57: C[2][2] += A[2][2] * B[2][2]
- ループ58: C[2][2] += A[2][3] * B[3][2]
- ループ59: C[2][3] += A[2][2] * B[2][3]
- ループ60: C[2][3] += A[2][3] * B[3][3]
- ループ61: C[3][2] += A[3][2] * B[2][2]
- ループ62: C[3][2] += A[3][3] * B[3][2]
- ループ63: C[3][3] += A[3][2] * B[2][3]
- ループ64: C[3][3] += A[3][3] * B[3][3]
となる。
Cに着目した場合、左上、右上、左下、右下とアクセスしている1回のキャッシュには上半分もしくは下半分が載せられるので、Cに関するメモリアクセスは2回しか起こらない。Aも同様である。
Bに着目したすると、左上、左下、右上、右下、左上、左下、右上、右下というアクセスをしている。つまり8回のメモリアクセスが発生する。最適化前よりメモリアクセスが24回減っていて、嬉しい、というのがループタイリングの効果である。これがもっと大きな行列になると全く無視できないようなオーバーヘッドがかかってくるという話。
さて、ここまで読んで勘のいい方は気づいたと思うのだが、Bの上半分をキャッシュしているが、Bのアクセスは左上の次に左下に行ってしまっている。つまり右上のキャッシュは無駄になっている。これがもったいないなと思うのだが、これは実は内側ループの順序を変えると解決できる。このようにする。
T = 2; // タイルサイズ
// 外側ループ
for io = 0; io < 4; io+=T {
for ko = 0; ko < 4; ko+=T { // これと
for jo = 0; jo < 4; jo+=T { // これを入れ替え
// 内側ループ
for ii = io; ii < min(4, io+T); ii++ {
for ki = ko; ki < min(4, ko+T); ki++ { // これと
for ji = jo; ji < min(4, jo+T); ji++ { // これも入れ替え
C[ii][ji] += A[ii][ki] * B[ki][ji]
}
}
}
}
}
}
こうすると、まずループ変数の挙動は次のようになる。
- ループ1 : io = 0, ko = 0, jo = 0, ii = 0, ki = 0, ji = 0
- ループ2 : io = 0, ko = 0, jo = 0, ii = 0, ki = 0, ji = 1
- ループ3 : io = 0, ko = 0, jo = 0, ii = 0, ki = 1, ji = 0
- ループ4 : io = 0, ko = 0, jo = 0, ii = 0, ki = 1, ji = 1
- ループ5 : io = 0, ko = 0, jo = 0, ii = 1, ki = 0, ji = 0
- ループ6 : io = 0, ko = 0, jo = 0, ii = 1, ki = 0, ji = 1
- ループ7 : io = 0, ko = 0, jo = 0, ii = 1, ki = 1, ji = 0
- ループ8 : io = 0, ko = 0, jo = 0, ii = 1, ki = 1, ji = 1
- ループ9 : io = 0, ko = 0, jo = 2, ii = 0, ki = 0, ji = 2
- ループ10: io = 0, ko = 0, jo = 2, ii = 0, ki = 0, ji = 3
- ループ11: io = 0, ko = 0, jo = 2, ii = 0, ki = 1, ji = 2
- ループ12: io = 0, ko = 0, jo = 2, ii = 0, ki = 1, ji = 3
- ループ13: io = 0, ko = 0, jo = 2, ii = 1, ki = 0, ji = 2
- ループ14: io = 0, ko = 0, jo = 2, ii = 1, ki = 0, ji = 3
- ループ15: io = 0, ko = 0, jo = 2, ii = 1, ki = 1, ji = 2
- ループ16: io = 0, ko = 0, jo = 2, ii = 1, ki = 1, ji = 3
- ループ17: io = 0, ko = 2, jo = 0, ii = 0, ki = 2, ji = 0
- ループ18: io = 0, ko = 2, jo = 0, ii = 0, ki = 2, ji = 1
- ループ19: io = 0, ko = 2, jo = 0, ii = 0, ki = 3, ji = 0
- ループ20: io = 0, ko = 2, jo = 0, ii = 0, ki = 3, ji = 1
- ループ21: io = 0, ko = 2, jo = 0, ii = 1, ki = 2, ji = 0
- ループ22: io = 0, ko = 2, jo = 0, ii = 1, ki = 2, ji = 1
- ループ23: io = 0, ko = 2, jo = 0, ii = 1, ki = 3, ji = 0
- ループ24: io = 0, ko = 2, jo = 0, ii = 1, ki = 3, ji = 1
- ループ25: io = 0, ko = 2, jo = 2, ii = 0, ki = 2, ji = 2
- ループ26: io = 0, ko = 2, jo = 2, ii = 0, ki = 2, ji = 3
- ループ27: io = 0, ko = 2, jo = 2, ii = 0, ki = 3, ji = 2
- ループ28: io = 0, ko = 2, jo = 2, ii = 0, ki = 3, ji = 3
- ループ29: io = 0, ko = 2, jo = 2, ii = 1, ki = 2, ji = 2
- ループ30: io = 0, ko = 2, jo = 2, ii = 1, ki = 2, ji = 3
- ループ31: io = 0, ko = 2, jo = 2, ii = 1, ki = 3, ji = 2
- ループ32: io = 0, ko = 2, jo = 2, ii = 1, ki = 3, ji = 3
- ループ33: io = 2, ko = 0, jo = 0, ii = 2, ki = 0, ji = 0
- ループ34: io = 2, ko = 0, jo = 0, ii = 2, ki = 0, ji = 1
- ループ35: io = 2, ko = 0, jo = 0, ii = 2, ki = 1, ji = 0
- ループ36: io = 2, ko = 0, jo = 0, ii = 2, ki = 1, ji = 1
- ループ37: io = 2, ko = 0, jo = 0, ii = 3, ki = 0, ji = 0
- ループ38: io = 2, ko = 0, jo = 0, ii = 3, ki = 0, ji = 1
- ループ39: io = 2, ko = 0, jo = 0, ii = 3, ki = 1, ji = 0
- ループ40: io = 2, ko = 0, jo = 0, ii = 3, ki = 1, ji = 1
- ループ41: io = 2, ko = 0, jo = 2, ii = 2, ki = 0, ji = 2
- ループ42: io = 2, ko = 0, jo = 2, ii = 2, ki = 0, ji = 3
- ループ43: io = 2, ko = 0, jo = 2, ii = 2, ki = 1, ji = 2
- ループ44: io = 2, ko = 0, jo = 2, ii = 2, ki = 1, ji = 3
- ループ45: io = 2, ko = 0, jo = 2, ii = 3, ki = 0, ji = 2
- ループ46: io = 2, ko = 0, jo = 2, ii = 3, ki = 0, ji = 3
- ループ47: io = 2, ko = 0, jo = 2, ii = 3, ki = 1, ji = 2
- ループ48: io = 2, ko = 0, jo = 2, ii = 3, ki = 1, ji = 3
- ループ49: io = 2, ko = 2, jo = 0, ii = 2, ki = 2, ji = 0
- ループ50: io = 2, ko = 2, jo = 0, ii = 2, ki = 2, ji = 1
- ループ51: io = 2, ko = 2, jo = 0, ii = 2, ki = 3, ji = 0
- ループ52: io = 2, ko = 2, jo = 0, ii = 2, ki = 3, ji = 1
- ループ53: io = 2, ko = 2, jo = 0, ii = 3, ki = 2, ji = 0
- ループ54: io = 2, ko = 2, jo = 0, ii = 3, ki = 2, ji = 1
- ループ55: io = 2, ko = 2, jo = 0, ii = 3, ki = 3, ji = 0
- ループ56: io = 2, ko = 2, jo = 0, ii = 3, ki = 3, ji = 1
- ループ57: io = 2, ko = 2, jo = 2, ii = 2, ki = 2, ji = 2
- ループ58: io = 2, ko = 2, jo = 2, ii = 2, ki = 2, ji = 3
- ループ59: io = 2, ko = 2, jo = 2, ii = 2, ki = 3, ji = 2
- ループ60: io = 2, ko = 2, jo = 2, ii = 2, ki = 3, ji = 3
- ループ61: io = 2, ko = 2, jo = 2, ii = 3, ki = 2, ji = 2
- ループ62: io = 2, ko = 2, jo = 2, ii = 3, ki = 2, ji = 3
- ループ63: io = 2, ko = 2, jo = 2, ii = 3, ki = 3, ji = 2
- ループ64: io = 2, ko = 2, jo = 2, ii = 3, ki = 3, ji = 3
計算としては次のようになる。
- Cの左上 (前半) : Aの左上 * Bの左上
- ループ1 : C[0][0] += A[0][0] * B[0][0]
- ループ2 : C[0][1] += A[0][0] * B[0][1]
- ループ3 : C[0][0] += A[0][1] * B[1][0]
- ループ4 : C[0][1] += A[0][1] * B[1][1]
- ループ5 : C[1][0] += A[1][0] * B[0][0]
- ループ6 : C[1][1] += A[1][0] * B[0][1]
- ループ7 : C[1][0] += A[1][1] * B[1][0]
- ループ8 : C[1][1] += A[1][1] * B[1][1]
- Cの左上 (後半) : Aの左上 * Bの右上
- ループ9 : C[0][2] += A[0][0] * B[0][2]
- ループ10: C[0][3] += A[0][0] * B[0][3]
- ループ11: C[0][2] += A[0][1] * B[1][2]
- ループ12: C[0][3] += A[0][1] * B[1][3]
- ループ13: C[1][2] += A[1][0] * B[0][2]
- ループ14: C[1][3] += A[1][0] * B[0][3]
- ループ15: C[1][2] += A[1][1] * B[1][2]
- ループ16: C[1][3] += A[1][1] * B[1][3]
- Cの右上 (前半) : Aの右上 * Bの左下
- ループ17: C[0][0] += A[0][2] * B[2][0]
- ループ18: C[0][1] += A[0][2] * B[2][1]
- ループ19: C[0][0] += A[0][3] * B[3][0]
- ループ20: C[0][1] += A[0][3] * B[3][1]
- ループ21: C[1][0] += A[1][2] * B[2][0]
- ループ22: C[1][1] += A[1][2] * B[2][1]
- ループ23: C[1][0] += A[1][3] * B[3][0]
- ループ24: C[1][1] += A[1][3] * B[3][1]
- Cの右上 (後半) : Aの右上 * Bの右下
- ループ25: C[0][2] += A[0][2] * B[2][2]
- ループ26: C[0][3] += A[0][2] * B[2][3]
- ループ27: C[0][2] += A[0][3] * B[3][2]
- ループ28: C[0][3] += A[0][3] * B[3][3]
- ループ29: C[1][2] += A[1][2] * B[2][2]
- ループ30: C[1][3] += A[1][2] * B[2][3]
- ループ31: C[1][2] += A[1][3] * B[3][2]
- ループ32: C[1][3] += A[1][3] * B[3][3]
- Cの左下 (前半) : Aの左下 * Bの左上
- ループ33: C[2][0] += A[2][0] * B[0][0]
- ループ34: C[2][1] += A[2][0] * B[0][1]
- ループ35: C[2][0] += A[2][1] * B[1][0]
- ループ36: C[2][1] += A[2][1] * B[1][1]
- ループ37: C[3][0] += A[3][0] * B[0][0]
- ループ38: C[3][1] += A[3][0] * B[0][1]
- ループ39: C[3][0] += A[3][1] * B[1][0]
- ループ40: C[3][1] += A[3][1] * B[1][1]
- Cの左下 (後半) : Aの左下 * Bの右上
- ループ41: C[2][2] += A[2][0] * B[0][2]
- ループ42: C[2][3] += A[2][0] * B[0][3]
- ループ43: C[2][2] += A[2][1] * B[1][2]
- ループ44: C[2][3] += A[2][1] * B[1][3]
- ループ45: C[3][2] += A[3][0] * B[0][2]
- ループ46: C[3][3] += A[3][0] * B[0][3]
- ループ47: C[3][2] += A[3][1] * B[1][2]
- ループ48: C[3][3] += A[3][1] * B[1][3]
- Cの右下 (前半) : Aの右下 * Bの左下
- ループ49: C[2][0] += A[2][2] * B[2][0]
- ループ50: C[2][1] += A[2][2] * B[2][1]
- ループ51: C[2][0] += A[2][3] * B[3][0]
- ループ52: C[2][1] += A[2][3] * B[3][1]
- ループ53: C[3][0] += A[3][2] * B[2][0]
- ループ54: C[3][1] += A[3][2] * B[2][1]
- ループ55: C[3][0] += A[3][3] * B[3][0]
- ループ56: C[3][1] += A[3][3] * B[3][1]
- Cの右下 (後半) : Aの右下 * Bの右下
- ループ57: C[2][2] += A[2][2] * B[2][2]
- ループ58: C[2][3] += A[2][2] * B[2][3]
- ループ59: C[2][2] += A[2][3] * B[3][2]
- ループ60: C[2][3] += A[2][3] * B[3][3]
- ループ61: C[3][2] += A[3][2] * B[2][2]
- ループ62: C[3][3] += A[3][2] * B[2][3]
- ループ63: C[3][2] += A[3][3] * B[3][2]
- ループ64: C[3][3] += A[3][3] * B[3][3]
こうすると、CとAについては左上、右上、左下、右下というアクセスになるのでメモリアクセスは2回で変わらず。
Bについては、左上、右上、左下、右下、左上、右上、左下、右下というアクセスになるので、メモリアクセスがさらに4回に減らせる。
行列積の最適化ではこのjループとkループの入れ替えはかなり定番である。
このようなことがなぜできるかというと、AとCはiiによって行が決まり、iiはji・kiよりも外側にあるため、ji・kiが変化してもiiは変化しないためである。Bはkiでの行アクセスになるため順序に敏感である。
私がループタイリングについてわかっているのは以上のようなことである。
発展させていくと、さらに複数のキャッシュラインの活用やL1/L2/L3キャッシュの階層ごとの活用など、本当はもっと奥が深いのだと思うが、入門にはこのくらいで十分だろう。