subroutine subblock-matmul({ ,}, { }, { }):
do = 1: /
do = 1: /
kernel8x24( ,, , ) end
end
図6.4 計算カーネルを⽤いた部分⾏列積の擬似プログラム．
前節の計算カーネルを⽤いると，部分⾏列積は図 6.4のようなプログラムとなる．実際にはこの 2 重 ループのうち内側ループは上三⾓の構造に合わせて開始位置を調整し，計算量を抑える．このプログ ラムにおいて𝑢 は内側ループで繰り返し使われる．𝑋̄ のパッキングを⾏わない場合は，この𝑢 への アクセスは⾏列𝐴へ直接アクセスすることになる．𝑣 は外側ループで繰り返し使われるため，キャッ シュに乗せる．前者の格納に必要なデータ量は64𝑏 Byteであり，後者は8𝑏 𝑏 Byteである．よって データ量だけを考えて後者を L2 キャッシュに収めるようなブロック幅を計算すると，𝑏 ,𝑏 を約 181 にすればよいが，計算カーネルのブロック幅や，𝑏 と𝑏 の性能への影響の違い，キャッシュ制御機構 の影響を考慮したうえで最適な𝑏,𝑏 を理論的に求めることは困難である．そこで，⼿動でチューニ ングした結果，性能の良かった𝑏 = 144,𝑏 = 200を⽤いる．この場合，𝑣 のデータ量は 225KB，𝑢 は 12.5KB となり，𝑢 は⼗分 L1D キャッシュに収まる⼤きさとなるが，次に⽰すコア内での並列化の ため，𝑢 も L2 キャッシュに収めるようにする．
このように部分⾏列積ではキャッシュの使⽤量を考えるため，キャッシュを共有するコア内部での 並列化を同時に取り扱うと都合がよい．subblock-matmulにおいて，外側ループを並列化すると，
𝑢 は共有できないが，𝑣 は共有でき，L2 キャッシュのデータの共有ができる．この並列化において は，内側ループの⻑さの違いや，スレッド間の多少の実⾏時間の違いを吸収するため，外側ループ番 号の𝑖をアトミック変数として動的な並列化を⾏う．
6.3.1 内積⽅向分割
内積⽅向分割は，𝐴 = 𝐴 𝐴 ⋯ 𝐴 と𝑁 個に分割しておき，より⼩さな DSYRK である
𝐶 = 𝐴 𝐴 をすべて並列に⾏った後，最後に総和𝐶 = ∑ 𝐶 を計算するものである．これはプログラ
ムblocked-matmul2において，中央の𝑘に関するループを並列化したことに相当する．
この並列化⼿法の利点は，演算量や計算時間を均等化しやすい点である．𝐴の分割を 1 列単位で⾏
えば，各𝐶 間の演算量の違いは 1 列分以下となるため，𝑛が𝑁 と⽐べて⼗分⼤きければ，演算量の違 いは⼩さくなる．また，各𝐶 の計算領域は𝐴 の列数の違いを除けば同じとなるため，計算時間の違 いも⼩さくなると予想される．
⼀⽅，問題点は，分割によって内積⽅向のデータ再利⽤性が低下することと，総和の計算のための コストが増えること，そして𝐶 のためのメモリ領域が必要になることである．𝐴 の列数がキャッシュ ブロックの⼤きさ𝑏 よりも⼩さくなると，その分だけデータ再利⽤性が低下するため，性能が劣化す る．また，𝐴 の列数が𝑏 よりも多少⼤きい場合でも，端数処理のために𝑏 よりも⼩さい列数のブロ ックを計算するため，同じ問題が起きる．総和の計算は，⾏列積の内積のうちいくつかの⾜し算を最 後に処理することであるため，⾏列積の演算量を変化させない．しかし積和命令を持つような CPU で は，積と和を同時に⾏えないため演算量が増えることと同然となる．また，総和の計算は演算量に対 してデータ移動量が⼤きいため，データ移動時間を演算時間の中に隠蔽することができない．このた め総和の計算の実⾏効率は⾏列積⾃体と⽐べて⼩さくなり，相対的に総和のコストが⼤きく⾒えやす い．以上 2 つの問題点は𝑛が𝑁 と⽐べて⼗分⼤きければ無視できるようになる．
𝐶 を格納するために必要なメモリ量は，ユーザーから与えられた𝐶を活⽤すると 1 つ分だけ省略で きるので，8𝑚 (𝑁 − 1) Byteとなる．ただし，今回は実装していないが，𝐶 の対称性を利⽤すればメ モリ量を約半分にすることも可能である．この式を⽤いて実際に必要なメモリ量を計算すると，例え
ば𝑚 = 5, 000，𝑁 = 56の場合，約 10GB のメモリ領域を必要とすることになり，現実的でない．また
データパッキングを⾏う場合，𝑋̄や ̄𝑌もコアごとに⽤意しなければならない．とくに𝑋̄はもとから⼤
きな領域となっているため，問題である．
そこで，内積⽅向の分割は単独では利⽤せず，別の並列化⼿法（以下で説明する 1 次元，2 次元分割）
と組み合わせて利⽤する．例えば𝑁 = 10としておいて，各𝐶 を別の並列化⼿法を⽤いて，𝑁 = 5ま
たは𝑁 = 6並列で実⾏する．こうすれば，必要なメモリ量を⼤幅に⼩さくでき，なおかつ，以下の並
列化⼿法においても現れる，並列度が⼤きい場合のデメリットを緩和できる．
ここで，𝑁 の決定基準が問題となる．𝑁 の値は，メモリ量や計算時間にも影響するが，ここでは，
使⽤するメモリ量の最⼤値を決めて，その範囲で最も⼤きな𝑁 を選択する，固定量メモリ基準を使う．
いま全体のコア数を𝑁 ，𝐶 のために⽤意するメモリ量の最⼤値を要素数として𝑀 = 25 × 10 のよう に決めたとき，
̄
𝑚 = 8⌈𝑚/8⌉ (6.8)
とおき，
𝑁 = min 𝑀
̄
𝑚 + 1, 𝑁 (6.9)
と決める．𝑚̄ はアライメントを調整した𝑚であり，𝑚 × ̄̄ 𝑚の⼤きさの⾏列が𝑀の中に何個⼊るか計 算して，それに 1 を⾜した値を𝑁 としている．ただし𝑁 より⼤きくならないように調整する．これ
は𝑁 = 56だとすると，𝑚 = 672のときの最⼤値56をとり，𝑚が⼤きくなるにつれ段階的に減少し，
𝑚 ≥ 5, 000のとき最⼩値1をとる．このとき𝑋̄ を格納するために必要な領域は，8 ̄𝑚𝑁 𝑏 Byteである．
̄
𝑚𝑁 の最⼤値は𝑀によって決まるため，𝑏 ≥ max( ̄𝑚𝑁 )として𝑋̄を 1 つ⽤意しておき，𝑁 > 1の場合 は𝑋̄のメモリ領域を𝑁 個に分割して⽤いるようにすれば，𝑋̄のためのメモリ領域も固定サイズになる．
固定量メモリ基準によって𝑁 を𝑚 に反⽐例するように決めることは，使⽤するメモリ量を⼀定に できるだけでなく，次の観点で合理性がある．まず，𝑚が⼗分⼤きい場合は，以下の並列化⼿法を⽤
いても⼗分⾼い性能が実現できるため，⼩さな𝑁 でよいが，逆に𝑚が⼩さい場合は，以下の並列化
⼿法の性能が落ちるため，𝑁 を⼤きくすべきである．ただし，この基準は𝑛を計算に使⽤していない ため，「𝑛が⼩さく，総和のコストが相対的に⼤きくなる状態では，𝑚が⼩さくても𝑁 を⼩さくする」
などの調整は⾏えない．このような計算時間の最⼩化を⾏うためには，計算時間のモデル化が必要で あるため，今後の課題とする．
6.3.2 1 次元分割
1 次元静的分割⼿法は内積⽅向分割を⾏った後の⾏列𝐶 を𝑁 個の列ブロックに分割し，それぞれ の列ブロックに対する計算を並列に⾏う⼿法である．（内積⽅向分割と組み合わせない場合は𝐶 = 𝐶 をコア数だけ分割する．）すなわち𝐶 を：
𝐶 = 𝐷 𝐷 ⋯ 𝐷 , (6.10)
のように分割する．これはプログラムblocked-matmul2において最内側ループの𝑗に関するループ を並列化したものである．ただし𝐶 は上三⾓部分のみを計算するため，列ブロック幅を均等に分割す れば演算量のばらつきが⼤きくなってしまう．そこで列ブロック幅をうまく選択することで実⾏時間 のばらつきを⼩さくしたいが，正確な実⾏時間を予測するためには性能モデルが必要であるため，こ こではコアごとの演算量を均等化することを⽬標とする．
いまそれぞれの列ブロックの列幅を𝑐 , 𝑐 , … , 𝑐 とおき，列ブロックの開始位置𝑗 = ∑ 𝑐（ただ
し𝑗 = 0）とおく．𝑗 は本来，整数だが，実数の範囲で考えると最適値は簡単に計算でき，
𝑗 − 𝑗
2 = 𝐷 (6.11)
と漸化式で求められる．ただし，𝐷 = であり，上三⾓部分の⾯積をコア数で割った値である．
GotoBLAS や BLIS で実装されている⼿法は，単純にこれを整数で丸めたものであり，次の漸化式を 使っている．
𝑗 = round 𝑗 + 2𝐷 . (6.12)
ただし，終端は⾏列サイズと⼀致するよう𝑗 = 𝑚とする．roundは丸め関数で，計算カーネルのブ ロック幅 8 へ丸める．このように求めた𝑐 は𝑏 より⼤きくなりうるが，その場合は𝑏 ごとに分割し て，同じコアで計算させればよい．
𝑗 を整数の範囲で考えた場合，最適値を求めることは整数最適化問題となる．すなわち min . max ∑ 𝑖 𝑠 = 1, 2, … , 𝑁 ,
sub. to. 0 = 𝑖 ≤ 𝑖 ≤ 𝑖 ≤ ⋯ ≤ 𝑖 = 𝑚, (6.13)
の最適解を求めることになる．これは⼀⾒複雑な問題だが，実際には，逆問題を考えることで簡単に 解くことができる．つまり，「コア数（の上限）が与えられたとき，演算量の最⼤値を最⼩化する問題」
の代わりに，逆問題「1 コアあたりの演算量の上限が与えられたときに，必要となるコア数の最⼩値を 求める問題」を考える．この逆問題は，1 番⽬のコアから順に上限を満たす中で最も⼤きな列幅を選択 する貪欲解法によって計算可能である．なぜならば，残りの列が多くなればなるほど必要となるコア 数が増⼤するからである．また，この逆問題を，演算量の上限が与えられたときのコア数の最⼩値を 返す関数だと考えると，これは単調⾮増加関数である．そこで演算量の上限を⼆分探索することでコ ア数が𝑁 以下のものの内，もっとも⼩さな演算量の上限を得ることができ，これは元の整数最適化 問題の最適値となる．このアルゴリズムは，列幅を単語の⻑さと⾒ることで，word wrapping の逆問 題として理解することができ，その解説がインターネット上に存在する [136]．
1 次元分割の問題点は，列ブロックの幅が⼩さくなり，データの再利⽤性が低くなりやすいことで ある．式（6.11）で概算すると，𝑁 = 𝑁 = 56の場合𝑚 = 5000のとき𝑐 は約 45 となる．これはキ ャッシュブロック幅𝑏 の約 3 分の 1 である．この問題は内積⽅向分割との組み合わせで緩和できる．
6.3.3 2 次元分割
1 次元分割では，1 列ブロック当たり 1 コアを割り当てるので，右端の列ブロックの幅が⼩さくなっ た．2 次元分割では，右の列ブロックに⾏くほど演算量を増やす代わりに多くのコアを割り当てるこ とで，列ブロックの幅を均等にする．1 つの列ブロックをさらに⾏ブロックに分割することで，列ブロ ック内での並列化を⾏う．これはプログラムblocked-matmul2において，最内側処理を並列化した ことに相当する．例として，図6.5に分割の様⼦を⽰す．
まず，列ブロックの個数を𝑁 ，𝑠番⽬の列ブロックに割り当てるコアの数を𝑚 とおく．このとき，
1 コアあたりの演算量がなるべく等しくなるように分割するため，1 次元静的分割の問題を拡張した次 の問題の最適解を，列ブロックの分割位置とする．
min . max ∑ 𝑠 = 1, 2, … , 𝑁 , sub. to. 0 = 𝑗 ≤ 𝑗 ≤ 𝑗 ≤ ⋯ ≤ 𝑗 = 𝑚,
(6.14) この問題に対しても 1 次元分割の時と同じアルゴリズムが使えるので効率的に計算できる．𝑁 や𝑚 の設定の仕⽅が問題であるが，分割した領域が正⽅形となることを⽬指して，𝑁 ≈ 2 𝑁 ，𝑚 ≈ 𝑠を 整数に丸めた値とする．