課題1 で完成させたプログラムの高速化を試みる。LU分解 を行うC言語の関数 に対して高速化を試みる。元の関数は消さずに残しておき、 速度比較ができるよう にする。速度比較は行列のサイズ n の変化に対し， LU分解に要する時間がどう変化するか調べること。
以下にはさまざまな高速化手法を示すが、 現在重要になっているのは キャッシュなどに関する局所性を高めることである。 TLBミス，キャッシュミスについて資料・説明 や学科の講義に基づきよく理解しておくことが望ましい。
すくなくともキャッシュミスを減らすような形で高速化を行うことを必須とするが、 すべての高速化手法を試す必要はない。 余力がある場合、他の実験履修者の速度を上回りたい場合などに参考にすること。
並列プログラムのサンプルの 最初の説明は高速化に関するものに なっているので参考に．
補足: ここの課題で行っているような密行列行列演算などの高速化は、 高速化自体が目的であれば「共通パターン」についてライブラリ化した BLAS ルーチン等を利用するのが現実的である。 しかし、本演習は各自、高速化を体験しようという趣旨なので、 BLASのような既存ライブラリを利用しないことになる。
本演習・実験で用いる計算機システムでは，Intel
386 プロセッサ(と浮動小数点コプロセッサ387)の命令セットに基づくものになって
いる。計算機室のプロセッサは何世代も後継の Intel Core 2 Duo で，
Core2 に基づくものになっているが，命令セットは同じと考えてもよい。
浮動小数点演算はPentium 4 で追加されている SSE2，
あるいは Prescott (Pentium 4 改良版)で追加されている SSE3
を用いたほうがよい可能性がある．SSE3を浮動小数点演算に用いるためには，
gcc [その他のオプション] -msse3 -mfpmath=sse foo.c -o foo
のようにする．
逆に、SSE2(SSE3)を用いない場合は
浮動小数点演算はレジスタスタックを用いているので、少し特別に考える
必要がある。
なお、計算機室の Intel Core 2 Duo は SSE4.1 もサポートしているので、
gcc [その他のオプション] -msse4.1 -mfpmath=sse foo.c -o foo
とすることもできる。
もちろん、gcc に Core 2 Duo向けの最適化を指示することも有効であろう。
計算機室のCore 2 Duo の場合は、 -mtune=core2 を指定する。
(なお、gcc バージョン3.4以降にのみ -mtune があり、
それ以前は -mcpu であった。また -march とすると指定したマシン(だけ)の
命令を生成する。)
準備で述べたように， 計算機室の gcc にはバージョンが3つある(gcc, gcc34, gcc44)。 新しいものを使うことで性能が向上すると期待できる。
for(j=0;j<1000;j++) for(i=0;i<1000;i++) a[i][j] += c;を
for(i=0;i<1000;i++) for(j=0;j<1000;j++) a[i][j] += c;のように入れ替える。 いつも単純に入れ替えられるとは限らない(なぜか?)。
for(i=1;i<1000;i++) a[i] += b[j]*c;は
double w = b[j]*c; for(i=1;i<1000;i++) a[i] += w;などとできる可能性がある。
LU分解に関するブロック化(blocking)を含めた詳細については、
こちらに示す。
なお、LU分解のブロック化などは検索でも見つかる。
しかし，検索でみつかるもののほとんどは，
LU分解の分類で示した，
Doolittle型(l_ii=1)について扱っている。
本演習では，クラウト(Crout)型(u_ii=1)について
扱っているので，注意すること。
(ヒント: 本演習では配列の添え字を除くと，CではなくFortran でのブロック化LU分解に近くなる)
for(i=0;i<2*n;i++) a[i] += 1.0をアンローリングして
for(i=0;i<2*n;i+=2){a[i] += 1.0 ; a[i+1] += 1.0}
とすることができる。この例では繰り返し回数が偶数だったが、
奇数にも対応するにはどうしたらよいか考えよ。
また、ループの本体に2倍の処理が含まれるように展開したが、
もっと大きくすることも考えられる。
for(i=0;i<n;i++) for(j=0;j<n;j++) a[i][j] = ... b[i][j] ...↓
for(i=0;i<n;i+=2) for(j=0;j<n;j++) a[i][j] = ... b[i][j] ...; a[i+1][j] = ... b[i+1][j] ...「...」の部分に共通部分式があれば、 計算やレジスタへの読み出しを共通化できる。
prefetcht0命令、一度しかアクセスしないなら
prefetchnta命令が使える。
これらは、gcc の拡張機能を用いると:
asm("prefetcht0 %0"::"m"(a[k+1][j]):"memory");
のようにCプログラム中に記述することができる。ただし、
__builtin_prefetch(&a[k+1][j], 0, 3)
を用いることで，他機種としても同じように書くことはできる．
for(i=0;i<n;i++) a[i] = b[i] + c[i];は、
for(i=0;i<n;i++) {
t2 = b[i];
t3 = c[i];
t1 = t2 + t3;
a[i] = t1
}
であるが、bからの読み込みとcからの読み込みが終わらないと加算ができず、
加算が終わらないとaへの書き出しができない。
これを:
for(...;...;i+=2) {
s2 = b[i+1];
s3 = c[i+1];
t1 = t2 + t3
a[i-1] = s1;
if(i...) goto Lout;
t2 = b[i+2];
t3 = c[i+2];
s1 = s2 + s3;
a[i] = t1
}
のような形にできれば、
bからの読み込みとcからの読み込みを一つ前の反復に移動、
aへの書き出しを一つ後のの反復に移動することができ、
遅延に耐えることができる。
もちろんプログラムの意味を変えないように、
入口と出口で補正コードが必要であるが、省略した。