LAPACKに無謀な勝負を挑むの巻 - プログラマーになりたい。

注意(追記)：以下の文章に重大なミスがあった可能性があります。

すっごいいまさらですが。おれ、もしかしてすっごい勘違いをしていたかもしれない。

以下、冒頭のほうに
>　→C-BLAS使用(after): 0.22秒
などというトンデモな数字がありますが、これは、やっぱおかしすぎる。
これ、もしかしたら計算結果をprintfとかしてなくてコンパイラが「どうせ計算したって使わないんだから、これんだろ」と判断、バッサリやられたのかもしれないという悪寒がしているのです。

つまり正味
int main(){
return 0;
}
ぐらいに空っぽだということ。
だとしたらこの数字も十分あり得ると思う。
ぜんぜん分かってなかった！

LAPACKに無謀な勝負を挑む

事実上の世界標準にして事実上の世界最速とうたわれる線形代数系の数値計算用ライブラリ、「LAPACK(レイパック)」。
これはそのLAPACKに無謀な勝負を挑んで見事に負け続けている男達の物語である——(例の番組風に)。

そうですよ、どーせ勝てない。とてもがんばってもたぶん「迫る」ぐらいが限度なの。それはわかってるの。

「行列の積(行列のかけ算)」だけに絞って、どうやったら速くできるかという課題が、授業で出たんですよ。ここは、ほれ、一応いろんな意味でおれの専門と絡むしね、黙っちゃいられないのですよ。

つーか、世の中の数値計算はLAPACK(とBLAS)に頼りすぎじゃないかね。いや、頼るのはいいんですけど、しくみを解説する資料が(ネット上に)見あたらないんですよ。あるのはLAPACKの啓蒙(使い方etc.)ばっか。

ちなみに、行列Aを
(a_11 ... a_1N)
A=(... ... ...)、
(a_N1 ... a_NN)
のようなN×N行列、同様に行列BもN×N行列としたときの、
C=A×B
(a_11*b_11+...+a_1N*b_N1 ... a_11*b_1N+...+a_1N*b_NN)
=( ... ... ... )
(a_N1*b_11+...+a_NN*b_N1 ... a_N1*b_1N+...+a_NN*b_NN)
を行列Aと行列Bの積と言います(ここまで読んだ人は書かなくても知ってるでしょうが、一応)。

さて、ここでちょっと現状を説明してみます。
まず、C言語でやってますので、行列は倍精度浮動小数点型(double)の「2次元配列」で表現することにしました。
こんな感じ:

double a[N][N];
double b[N][N];
double c[N][N];
//(中身は適当にいれるたことにして、略。)

これを使ってふつーに積をとると:

for(int i1=0;i1<N;i1++){
 for(int i2=0;i2<N;i2++){
  c[i2][i1]=0;
  for(j=0;j<N;j++){
   c[i2][i1]+=a[j][i1]*b[i2][j];
  }
 }
}

これでN=512(マシンはPPC G4 1.33GHz/memory 512MBのMac mini。以下同じ)とすると、所用時間(上記の処理だけ。以下同じ)は

26.75秒

でした。

まず、いちばん簡単そうなところから手をつけてみました。
c[i2][i1]の値をいじる、というのは、たぶん、

cの値を参照
i2の値を参照
c[i2]の値を参照
c[i2][i1]の値を参照

の様な手順になるわけで、メモリ上のあっちこっちを参照しに、非常に遠い距離を移動することになります。キャッシュミスが多発します、たぶん(知識がうすっぺらいのであんまり自信がないのです)。なので、tmpというところにとりあえず入れといて、上の作業はあとで1回だけやらせるようにしました:

double tmp;
for(int i1=0;i1<N;i1++){
 for(int i2=0;i2<N;i2++){
  tmp=0;
  for(j=0;j<N;j++){
   tmp+=a[j][i1]*b[i2][j];
  }
  c[i2][i1]=tmp;
 }
}

所要時間は24.67秒(約2秒短縮)。
だいたい8%か。ま、この程度の小細工じゃ、こんなもんでしょう。

さて、できたらtmpを参照する回数も減らしたい。
というわけで、いままでは

a[j][i1]*b[i2][j]

を1回計算するたびにtmpに足してましたが、1度に5回分計算してしまう様にしちゃいます:

for(int i1=0;i1<N;i1++){
 for(int i2=0;i2<N;i2++){
  tmp=0;
  for(j=0;j<N-5;j+=5){
   tmp+=a[j][i1]*b[i2][j]
     +a[j+1][i1]*b[i2][j+1]
     +a[j+2][i1]*b[i2][j+2]
     +a[j+3][i1]*b[i2][j+3]
     +a[j+4][i1]*b[i2][j+4];
  }
  for(;j<N;j++){
   tmp+=a[j][i1]*b[i2][j];
  }
  c[i2][i1]=tmp;
 }
}

これでも所要時間24.67秒。
ちなみに10回分まとめると22.69秒。
地味だーー。

ちなみに、LAPACKじゃなくて、LAPACKが使っているらしいBLASの方を直接使うCBLASというのがあるらしいのですが、これ使うと:

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
  N, N, N, 1.0, a[0], N, b[0], N, 0., c[0], N);

答えがあってるのかどうかチェックしてないので知りません(怪しい気もする)が、所要時間はなんと

~~0.22秒~~*1。

って、うわまじかよーーー。
なに、魔法か？え？
やだもう…。
しょうがないなあもう。BLASのコードを読みますよ。言語はFORTRANですけど。…ってFORTRANのコードなんて、見たことないんですけどー！

(つづきます: nikki-da)

mixiコメント

Ron 2007年02月02日 08:44
5回分まとめたコードだと、a[510][il]とかが計算に入ってないのでは？それとも、その分はループの後で計算する？
で、まとめないでループを回してしまう方法だと、512*512*512回のループを行っているわけで、134,217,728ステップになります。適当に考えて、かけ算で1クロック、足し算で3クロックを使っているとすると、それに4を掛けて、536,870,912クロックになります。CPUが独占できているとすると、それを1.33GHzで割って、0.40秒になります。つまり、上のアルゴリズムだと、CPUを独占できたとしてもLAPACに迫ることすら出来ないことが分かります。
それにしても遅いので、何か他に要因がないですか？仮想記憶領域に入ってしまっているとか？CPU使用率が低いとか。実行中の様子をtopコマンドなどで確認してみてください。
アルゴリズムの問題ですが、行列を列ベクトルにして計算してみたらどうでしょう？たぶん、ループ回数は減ると思います。
0.22秒のアルゴリズム、ぜひ知りたいですね。

だまん(flashingwind) 2007年02月02日 14:45
■Ronさんへ
「5回分まとめた」やつは、直後の初期化がないforループで足りない分を補ってます(たぶんOKなはず)。
あと、1命令で1クロック消費ってことは無いはずですので、そんなにかからないのではないかと思います。

だまん(flashingwind) 2007年02月02日 14:53
追伸:
http://www.netlib.org/blas/dgemm.f
をご参照下さい(FORTRANですが)。
キャッシュミスを極限まで減らす「魔法」のレシピのポイントは、配列の参照がすべて頭から1つづつである(飛ばさない)ことにあるようです。

だまん(flashingwind) 2007年02月03日 00:58
「1命令で1クロック消費ってことは無いはず」(というか、よく見ると分母分子逆のことを言いたかったんですが。)というのと、キャッシュミスの件と、どちらも何らかの資料を見たわけではなくて、1つ前のコメントにURLを載せたコードとの比較からの推測です。関連するのはおそらく
*
* Form C := alpha*A*B + beta*C.
*
(中略)
DO 80, L = 1, K
IF( B( L, J ).NE.ZERO )THEN
TEMP = ALPHA*B( L, J )
DO 70, I = 1, M
C( I, J ) = C( I, J ) + TEMP*A( I, L )
70 CONTINUE
END IF
80 CONTINUE
という部分ですが、このように、乗算して加算する部分のくり返し回数は同じで、使用する演算(インストラクション)の差もほぼ無いように見えます。