一度キャッシュしたデータを何度も再利用することで、 メモリアクセスの平均時間を短縮し、 プログラムを高速化することができる。
このために、データをある固まりとして扱うように、 プログラムを改良することをブロック化という。
ブロック化では、点(1要素)や線(ベクトル)に幅を持たせる。 これをブロック幅と呼ぶ。
行列積相当の3重ループにブロック化は有効だが、 内積や外積相当の2重ループにはうまく適用できないことが多い。
行列の場合、小行列に分けて考えるアプローチがある。
例えば n×n 行列 Aを a_ij (i=1, 2, ..., n; j=1, 2, ..., n)
と1要素単位で考えるのではなく、ブロック幅 w を考えて N = n/w とし
(例えば w = 16 などとする)、
w×w の小行列 A_IJ (I=1, 2, ..., N; J=1, 2, ..., N)に分けて考える。
ただし、n は w で割り切れるとする。
(一般には割り切れないので実際のプログラミングでは辻褄あわせが必要となる。)
つまり
[ A_11 A_12 A_13 ]
A = [ A_21 A_22 A_23 ]
[ A_31 A_32 A_33 ]
などとなる。
すると、要素 a_ijについての定式化が
小行列 A_IJ についての定式化にうまく応用できることがある。
n×n 行列 Aを ブロック幅 w を考えて N = n/w とし、
w×w の小行列A_IJ (I=1, 2, ..., N; J=1, 2, ..., N)に分けて考える。
ただし、n は w で割り切れるとする。
(一般には割り切れないので実際のプログラミングでは辻褄あわせが必要となる。)
[ A_11 A_12 A_13 ]
A = [ A_21 A_22 A_23 ]
[ A_31 A_32 A_33 ]
下三角行列 L, 上三角行列 U も同様に、L_IJ,
U_IJ で考える。
例えば、小行列 L_II は下三角行列、
I<J のとき L_IJ は、零行列である。
[ L_11 0 0 ]
L = [ L_21 L_22 0 ]
[ L_31 L_32 L_33 ]
すると
j
a_ij = Σ l_ik u_kj (i ≧ j)
k=1
i
a_ij = Σ l_ik u_kj (i < j)
k=1
と同様に、
J
A_IJ = Σ L_IK U_KJ (I ≧ J)
K=1
I
A_IJ = Σ L_IK U_KJ (I ≦ J)
K=1
であるから、総和の最後の部分（小行列積）を分けると
J-1
A_IJ = Σ L_IK U_KJ + L_IJ U_JJ (I ≧ J)
K=1
I-1
A_IJ = Σ L_IK U_KJ + L_II U_IJ (I ≦ J)
K=1
となり、これから、関係式★:
J-1
L_IJ = (A_IJ - Σ L_IK U_KJ) U_JJ^(-1) (I > J)
K=1
I-1
U_IJ = L_II^(-1) (A_IJ - Σ L_IK U_KJ) (I < J)
K=1
あるいは関係式★2:
J-1
A_IJ - Σ L_IK U_KJ = L_IJ U_JJ (I > J)
K=1
I-1
A_IJ - Σ L_IK U_KJ = L_II U_IJ (I < J)
K=1
I-1
A_II - Σ L_IK U_KI = L_II U_II (I = J)
K=1
が得られる。
これをLU分解の関係式☆
j-1
l_ij = a_ij - Σ l_ik u_kj (i ≧ j)
k=1
i-1
u_ij = (a_ij - Σ l_ik u_kj) / l_ii (i < j)
k=1
と比較すると、
以下の違いがある。
u_ii=1)で考えたので、
u_ii での除算がなかったが、上記関係式★では、
U_JJ の逆行列が用いられている。
（U_JJ は対角成分が1であっても、単位行列とは限らないため）
U_IJ に関する関係式は類似している。
L_II や
U_JJ については述べていない。
これは、関係式★2のほうでは3つ目の等式として述べている。
L_IJ、U_IJ
を求めるために使う手順には次のものがある。
L_II や U_JJ を求める。
L_II を含む 2つ目の等式を用いて、
左辺から U_IJ を求める。
これは、関係式★の2つ目の等式で右辺から左辺を求めることと等価である。
U_JJ を含む 1つ目の等式を用いて、
左辺から L_IJ を求める。
これは、関係式★の1つ目の等式で右辺から左辺を求めることと等価である。
LU分解の関係式☆と、 小行列に基づく関係式★(★2)の類似から、 ブロック化されたLU分解の3種類のアルゴリズムを考えることができる。 つまり、 ブロック化right-looking法、 ブロック化left-looking法、 ブロック化Crout法である。
特に、ブロック化right-looking法は、 M=1, 2, ... N についての以下の関係式に基づき、
A_II^(M) = L_II U_II (I=J=M) A_IJ^(M) = L_II U_IJ (I=M, J=M+1,...,N) A_IJ^(M) = L_IJ U_JJ (J=M, I=M+1,...,N) A_IJ^(M+1) = A_IJ^(M) - L_IK U_KJ (I=M+1,...,N, J=M+1,...,N, K=M)ブロック形式ガウス法と呼ばれる。 アルゴリズムとしては M として K を用いて、K ループを構成し、 I ループや J ループを用いればよい。
A_II^(M) を L_II と U_II
に LU分解するのに用いる。
U_IJ を求めるのに用いる。
L_IJ を求めるのに用いる。
A_IJ を M 段階から、M+1 段階へ更新するために用いる。
L_IK U_KJ が用いられており、これが小行列であることから、
一度キャッシュした小行列のデータを何度も再利用することで、
メモリアクセスの平均時間を短縮し、
プログラムを高速化することができると期待される。
また後述するように、
行列積に対してさらにブロック化を適用することも考えられる。
なお、行列積 L_IK U_KJ の結果を一度得てから
A_IJ を更新するのではなく、
一般化された行列積計算として、A_IJ を更新先
（内積による総和の部分和で更新する先）
として直接用いればよい。
ブロック化right-looking法のアルゴリズムでは、K ループ中に、
更新用の IJループが含まれる、3重ループが使われている。
さらに、行列積 L_IK U_KJ に基づく更新の計算には、
小行列についての要素単位の ijk ループが含まれている。
[
補足:
改めて、行列 A, B, C を用いて行列積 C = A B を考えると、
(行列 C の i 行 j 列目の要素を c_ij と書けば)
n
c_ij = Σ a_ik b_kj
k=1
であるからである。]
よって全体としてはループは6重(KIJijk)に入れ子であり、
Kループ中では5重(IJijk)の入れ子となる。この5重ループの
入れ子の交換(Jとi)を行うと、
IiJjkとなり、さらにIiをまとめてi'、
Jjをまとめてj'とすることができ、
i'j'k の3重ループとなる。
(Kループを含めると4重の入れ子)
i'j'k 3重ループは、
L_IK (I=K+1,...,N)を縦に連結した (w*(N-K))×w行列と、
U_KJ (J=K+1,...,N)を横に連結した w×(w*(N-K))行列との、
行列積で、
A_IJ (I=K+1,...,N, J=K+1,...,N)を縦横に連結した
(w*(N-K))×(w*(N-K))行列を更新することに相当する。
見方を変えると、 n×n 行列 A を最初に、多数の w×w の小行列に分割するのではなく、 左上に 1つの w×w 正方行列の小行列、 左下に 1つの (n-w)×w の小行列（ブロック幅w）、 右上に 1つの w×(n-w) の小行列（ブロック幅w）、 右下に 1つの (n-w)×(n-w) 正方行列の小行列、 に4分割し、 (n-w)×(n-w) の小行列に対して、 同じ分割を繰り返し（再帰的に）適用することに相当する。 このような分割に基づいてもLU分解のブロック化（ブロック幅w）を行うことができる。
i'j'kの3重ループに対して、
後述の行列積のブロック化を行ったり、
多重ループの入れ子の交換（例えばj'とk）や
多段ループ・アンローリングを含む高速化手法を適用したりすることで、
高速化が望める。
(本演習ではそうしないが、現実には、
高性能ライブラリが提供する行列積ルーチンが使うべきところ)
i'j'kの3重ループの行列積に関するブロック化の一つには、
見方を変える前の5重(IJijk)の入れ子に戻す方式がある。
ブロック幅wで，right-looking 法を用いる場合という、さらに別の考え方として、 右下の(w*(N-K))×(w*(N-K))行列を行列積で更新することが主要部分であり、 それ以外の部分はブロック化前のアルゴリズムと同様にすることも考えられる。 効果もある上に、最も簡便なブロック化方式ともいえる。
例: right-looking (down-looking?) アルゴリズムの場合 n=8のとき: w=2 で，K=0 (k=0, 1) の場合の処理後，K=1 (k=2, 3)のとき: k=2: ........ ........ ..*!!!!! ........ ........ ........ ........ ........ (*から読み出し．各!からの読み出しとそこへの書き込み) k=2: a[i][j]のwを超える部分の更新は後回し: ........ ........ ...***** ..*!!!!! ..*!.... ..*!.... ..*!.... ..*!.... (*から読み出し．各!からの読み出しとそこへの書き込み) k=3: ........ ........ ........ ...*!!!! ........ ........ ........ ........ (*から読み出し．各!からの読み出しとそこへの書き込み) k=3: k=2で後回しにしたa[i][j]の更新も合わせて(幅wだけまとめて) a[i][j]の更新 ........ ........ ....**** ....**** ..**!!!! ..**!!!! ..**!!!! ..**!!!! (*から読み出し．各!からの読み出しとそこへの書き込み) このとき、単純なkijループとしたのではブロック化の意味がないことに注意。例えば、簡便なブロック化方式のアルゴリズムは以下のようにして得られる。
手順の簡易メモ:: (各ステップで(動かして)正しさを確認しながら変形していくとよい)
ブロック幅wとする
k を kk,w,k と2重ループに
つまり k:[0,n-1) を kk:[w:0,n-1) k:[kk,min(kk+w,n-1)) に
# k:[n,m)は k が n ≦ k ＜ m を動くループであることを表す。
# kk:[w:n,m)は kk が w の倍数として n ≦ kk ＜ m
# を動くループであることを表す。
kk:[w:0,n-1) を kk:[w:0,n-w),kk:[w:n-w,n-1) に分割
kk:[w:0,n-w) k:[kk,min(kk+w,n-1)) を kk:[w:0,n-w) k:[kk,kk+w) に
kk:[w:n-w,n-1) k:[kk,min(kk+w,n-1)) を k:[kk,n-1) に (kkは上の最後の値のまま)
kk:[w:0,n-w) k:[kk,kk+w)にて
i:[k+1,n) を i:[k+1,kk+w),i:[kk+w,n) に分割
kk:[w:0,n-w) k:[kk,kk+w) i:[kk+w,n) にて
j:[k+1,n) を j:[k+1,kk+w),j:[kk+w,n) に分割
kk:[w:0,n-w) にて
k:[kk,kk+w) i:[kk+w,n) j:[kk+w,n) を 主要部分として最後に移動 (まとめる)
主要部分 k i j ループ => i k j ループ (あるいは i j k ループ)
n×n 行列 A, B, C を ブロック幅 w を考えて N = n/w とし、
w×w の小行列A_IJ, B_IJ, C_IJ
(I=1, 2, ..., N; J=1, 2, ..., N)に分けて考える。
ただし、n は w で割り切れるとする。
[ A_11 A_12 A_13 ]
A = [ A_21 A_22 A_23 ]
[ A_31 A_32 A_33 ]
[ B_11 B_12 B_13 ]
B = [ B_21 B_22 B_23 ]
[ B_31 B_32 B_33 ]
[ C_11 C_12 C_13 ]
C = [ C_21 C_22 C_23 ]
[ C_31 C_32 C_33 ]
行列積 C = A B を考える。
n
c_ij = Σ a_ik b_kj
k=1
と同様にして、
N
C_IJ = Σ A_IK B_KJ
K=1
よって、C_IJ を IJK ループで、
その中の行列積 A_IK B_KJ を
ijk ループで計算する場合を考えると、
IJKijkの6重ループとなる。
多重ループの入れ子の交換をして、
IKJikjの6重ループとした場合について、
http://en.wikipedia.org/wiki/Locality_of_referenceに例がある。
多重ループの入れ子の交換をして Iiのようになれば、
まとめて i' などとできるので、必ずしも6重ループである必要はない。例えば、
交換後が JKIikj ループだった場合、
Iiをi'にまとめると
JKi'kjの5重ループとなる。
ここでは正方行列を例にとって示したが、 正方行列の行列積以外の行列積についても同様にブロック化は適用できる。