プログラマーになりたい。

プログラミングや写真や本や読書会のことや、日常のこと。

トップ > 数値計算 > LAPACKに勝負を挑んでVecLibの魔法に負けた件

LAPACKに勝負を挑んでVecLibの魔法に負けた件

数値計算

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

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


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


つまり正味

int main(){
return 0;
}

ぐらいに空っぽだということ。
だとしたらこの数字も十分あり得ると思う。
ぜんぜん分かってなかった！

「LAPACKに無謀な勝負を挑むの巻」のつづき(読まなくていいですけど)。

(元ネタ: LAPACKに無謀な勝負を挑むの巻 - nikki-da)

なんなんでしょうねぇ。まったく。
勝てないことは分かってるんで、あきらめて手の内を覗いたんですよ。引き分けねらいで。で、そっくりそのまま写してみたんですが、だめでした。
Appleはいったいどんだけ最適化してるんだよ…。


Macには標準で「XCode」というプログラミング用のツール群(IDE:統合開発環境というのかな？)が付いてきますが、それにはベクトル計算用の最適化された(めちゃくちゃ速い)ライブラリ群(というかフレームワーク)である「Vector Library(vecLib)」が含まれてます。

このvecLibは

  • C-BLAS
  • C-LAPACK
  • などなど(他にもいっぱい)

を(Appleが最適化したバイナリを)含みます。


で、C-BLASで定義されている関数「cblas_dgemm」をつかって行列の積をとったら、

  • 工夫なし(before): 26秒前後
  • →C-BLAS使用(after): 0.22秒*1

と、このぐらい違いました。そこで、C-BLASや、(「C-」が付かない、Fortranで書いてある)本来のBLASのコードを参照して、入らない部分(転置とか)をばっさり切って、コピペ(コピー＆ペースト)しつつプログラムを作り直してみました。

double a[N1][M];
double b[M][N2];
double c[N1][N2];
double *A, *C, *B;
//(中略)
for(i1=0;i1<N1;i1++){
  C=c[i1];
  for(i2=0;i2<N2;i2++){
    *(C++)=0;
  }
  B=b[0];
  for(j=0;j<M;j++){
    const double temp=a[i1][j];
    if(temp!=0.0){
      //tmp=a[i1][j];
      C=c[i1];
      for(i2=0;i2<N2;i2++){
        *(C++) += temp * *(B++);
      }
    }
  }
}

結果は、

  • →C-BLASとほぼ同じ方法: 5.21秒

あれ？ポインタなんか使いたくないのに使ってみたりして、結構がんばったのになあ…。っつーことで最適化オプション(-O3)をつけてコンパイル。

やっぱだめなのです。結局、

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
 N1, N2, M, 1.0, a[0], M, b[0], N2, 0., c[0], N1);

には勝てないという結果でした。原因はよくわかりません。


vecLibは「PowerPC G4」以降のCPUに搭載された「AltiVec」(Apple的には「Velocity Engine」と呼ぶ)を利用して高速に処理するとか書いてあったりもしますが、AltiVecはどうも倍精度浮動小数点型(double型)の計算はできそうにない。並列化できるのはfloat×4個までとか。ということは1つしかないFPUをフル稼働させてるんでしょうが…。


あ、あと、Ronさんに嘘をいったことが判明しました。すいません。G4は倍精度浮動小数点型の加算・乗算など四則演算はぜんぶ3クロックサイクルで、パイプラインが3段階(1段階=1クロックサイクル)。どんなに効率が良くてもIPCは1以上ということらしいです。

mixiコメント

いわちぃくん 2007年02月04日 10:21
相変わらず奇藻男ちゃんだな。(笑)また教えてくらさいな。
次の金曜学校居ないんだっけ？遊び行くよ。

Ron 2007年02月04日 15:24
Velocicy Engineすごいな。ってか、コンピュータ遅い、とか言う前にプログラムなり、コンパイルを工夫しなさいという話ですね。プロセッサーに最適化してやれば100倍くらい早くなる可能性があると。

だまん(flashingwind) 2007年02月05日 02:03
■いわちぃくん
うるせいやい。キモイのはデフォルトですが何か。
金曜は卒研だからいるはず。何時ぐらいにいらっしゃるのかしら。

■Ronさん
そうそう。すごいんですよ。そういえば、PS3かなんか(とりあえず最近のゲーム機のどれか)も似たような(単精度までなら速いけど、doubleは駄目な)回路を積んでるという話を聞いたことがあります。

Intelのコンパイラといい、VecLibといいほんと、まったく、一体どんな魔法を使ってるのか…。ほんと、レシピを知りたいです。

いわちぃ 2007年02月05日 12:48
んじゃ、昼過ぎにでも。(笑)

なかじ 2007年02月05日 18:04
あなた学校来るのですか？
今日もいなかったですよね？
まぁもう一人の問題児も来ませんでしたけどね、笑

あと私のパソコンで行列の積を行ったところ、行列の要素数が２９４までしかできませんでした。ウィンドウズで２ＧＨｚのＲＡＭが５１２ＭＢのやつだと思うんですけど、なぜなんでしょうか？ＵＮＩＸじゃないから？

ちなみにコンパイラのせいかなと思い、gccじゃなくてborlandでやったら逆に要素数が減って２００くらいまでしかできませんでした。

このへんも素人にもわかりやすく、ご教授していただけると、有難き幸せでございます。

本当は今日学校で聞こうと思ってたんだけどね！！

だまん(flashingwind) 2007年02月06日 02:15
main()の中で、配列を宣言してるなら、スタック領域にとられるのですぐに限界が来ます。この原因であれば、「static」をつけて

    int main(){
      static double a[N][N];
      (略)
    }

とかやるか、あるいは

    double a[N][N];
    int main(){
      (略)
    }

とすると解決するかも。

kazusan 2007年02月11日 15:18
通りすがりの者です。

VecLib は ATLAS というフリーのパッケージを元にして作られています。
http://math-atlas.sourceforge.net/
にあります。ソースも公開されていますがVecLib の方が少しだけ速いかもしれません。
それから、double a[N][N] という形式は FORTRAN 向けのライブラリでは動作しない可能性が高いです。double a[N*N] にしましょう。

だまん(flashingwind) 2007年02月11日 15:42
コメントありがとうございます。

MacのVecLiBのBLASはATLAS系のものなんですか。知りませんでした。

というか、この日の日記を書いたあとで、BLASはインターフェイスみたいなものらしいと知ったので、本文にはあやしい記述がありますが、お許し下さい。

http://www.okada.jp.org/RWiki/?ATLAS+Rのビルド
(またはhttp://www.okada.jp.org/RWiki/?ATLAS%2BR%A4%CE%A5%D3%A5%EB%A5%C9
)
の最後の方なんかを見ると、ATLASのが速い場合もあるみたいですね。速いのはSVDなんでBLASじゃなくてLAPAC的な部分かもしれないですけど。

ちなみに、今回はBLASに関しては「ATLASよりも速い」とかいう話しを聞いて(とりあえず信じることにして)、GotoBLASをちらっと見たんですが肝心な部分はCPUごとにアセンブラで書かれていて、結局、真面目に読む気になりませんでした…。

kazusan 2007年02月11日 15:56
ほとんどのケースで生 LAPACK を使っているはずです。ブロッキングサイズをちょいといじってあるだけです。
GotoBLAS の方は読みやすさよりも速さ優先です。ただし、すべてのルーチンでやっていることは同じです。アセンブラ上での表現が違うだけ。作った本人が言うんだから間違いなし。

だまん(flashingwind) 2007年02月11日 19:48
> 作った本人が…
たしかに！(つまりいま気づいたんですが…。おそくてすいません)

そもそもが「ブロッキング」の言葉も知らずにコードを見に行ったんで(というか今もろくにわかってませんが)、読みやすさとか以前にアセンブラが読めても理解できたのかは疑問なんですが…。

しかし、速さ優先だとはいいつつも、アセンブラであの行数はすごいですね(本筋と関係ないけど…)。

kazusan 2007年02月12日 00:53
> アセンブラであの行数はすごいですね

実は本当に重要な部分は 100 行くらい。その他の部分は主に端数処理に関するもの。

> 「行列の積(行列のかけ算)」だけに絞って、どうやったら速くできるかという課題

出題者の意図は不明だけど、実際問題としてこの課題設定はかなり無謀に近い。おそらく意図としては、行った操作に対してどのように性能が変化するかを観察する、というところにあると思うんだけど、この関数は極限まで性能が出る特徴があるので、考察しづらいし、なにしろ間違ったことを覚える可能性が高い。

本当ならば、
for(j=0;j<N;j++){
tmp+=a[j]*b[j];
}
をどうやったら速くすることができるか、の方が勉強になると思う。

参考文献を挙げておきます。
・寒川光、「RISC 超高速化プログラミング技法」、共立出版
日本語で書かれたおそらく唯一のまともな最適化に関する本です。ただし、内容はちょっと古くて現在のＣＰＵの特性には合わない部分が多いのと、絶版なので入手は難しいということかな？

だまん(flashingwind) 2007年02月12日 02:42
すいません。なんかいろいろとアドバイスをいただいてしまい恐縮です。

> 出題者の意図
は、あいにくはっきりとは憶えてないですが、ざっくり言うと、数値計算でよく使われる計算だから、とかそんな感じだったような。普段はアーキテクチャの基礎、用語説明とかを中心にやっていますので、深い意図はないかもしれません。前回の(日記の直後の)授業であっさり「答え合わせ」らしきことをして、そのあとはIntelのコンパイラで最適化の効果を見せて、終わりでしたし。

とりあえず教えていただいた本を読んでみます(Amazon経由で中古がありました)。せっかく時間さいていただいておいてって感じですが、じつは個人的には、行列の積よりも(ちょうど例のような)内積の方をよく使いますので、ちょうどいい機会ですからそのへんを中心に追加でみてみます。どうもありがとうございました。

あと、ついでっぽくて申し訳ないですが、「マイミク申請」してもいいでしょうか。

kazusan 2007年02月12日 06:15
> Intelのコンパイラで最適化の効果を見せて

実はコンパイラって本気は出さないんです。Intel には MKL という商品があるので、コンパイラで完全に最適化ができてしまうとこちらが売れなくなる。だからコンパイラは「程ほどの」性能が出せるように調節してあります。

内積が理解できたら、行列ベクトル積、そして行列行列積へと段階を踏んだら良いと思います。ちなみに、G4 ではピークが 1Flop/cycle よりちょっと下の性能になるはずです。要するに理論値性能の半分です。これが、行列行列積になると理論値性能に近づくようになります。

> 「マイミク申請」してもいいでしょうか。

どうぞ。

*1:この数字はよく考えるとあり得ない気がする

Creative Commons License ©2007-2021 IIDA Munenori.