vDSPについてのメモ - 3

vDSPについてのメモの続き。ベンチマークをしてみよう。そして、最後に「なんじゃそりゃあ」をひとつ。

vDSPでのベクタユニットの単純ベンチマーク

とりあえず一番単純なベンチマークをやってみる。ふたつの1次元配列を用意してそのかけ算を計算してみよう。vDSPには

void vDSP_vmul(const float input1[],
   vDSP_Stride stride1,
   const float input2[],
   vDSP_Stride stride2,
   float result[],
   vDSP_Stride strideResult,
   vDSP_Length size);

というふたつの配列の要素ごとのかけ算をする関数がある。倍精度浮動小数点数用にvDSP_vmulD()がある。こっちは引数がdoubleの配列になっただけで同じことをする。

ところで、ふつうにCのforループで書けばスカラユニットの逐次計算に展開されるはず。つまり

void    byScalarUnit(float *f1, float *f2, float *f3, int num)
{
    int i;
    for (i = 0 ; i < num ; i ++)
        f3[i] = f1[i] * f2[i];
}

はさっきのvDSP_Strideが1の場合のvDSP_vmulと同じことをスカラユニットを使ってやるはず。これを比較してみよう。

まず、単精度不動小数点数の配列に対して

void    floatCompare()
{
    int     num = 1024 * 1024;                     // (1)
    int     ntimes = 50;                           // (2)
    struct  rtms  tstart, tscalar, tvector;// (*)
    float   *f1 = (float *)malloc(num * sizeof(float));
    float   *f2 = (float *)malloc(num * sizeof(float));
    float   *f3 = (float *)malloc(num * sizeof(float));
    int     i;
    
    for (i = 0 ; i < num ; i ++) {                 // (4)
        f1[i] = (float)i;
        f2[i] = (float)(num - i);
    }
    tStart(&tstart);                      // (*)
    for (i = 0 ; i < ntimes ; i ++)
        byScalarUnit(f1, f2, f3,  num);            // (5)
    tscalar = tMeasure(&tstart);          // (*)
    tStart(&tstart);                      // (*)
    for (i = 0 ; i < ntimes ; i ++)
        vDSP_vmul(f1, 1, f2, 1, f3, 1, num);       // (6)
    tvector = tMeasure(&tstart);          // (*)
    fprintf(stderr, "Scalar : %d\tvector : %d\tRatio = %f\n",
            tscalar.t.tms_utime, tvector.t.tms_utime,
            ((double)tscalar.t.tms_utime) / tvector.t.tms_utime);
    free(f1);
    free(f2);
    free(f3);
}

まず要素の数がnum個の配列をつくる。あまり小さいとCPUのキャッシュの中に収まってしまって実際の場合とかけ離れてしまうのでここでは(1)のように1M要素にして(4)のところで適当な値を入れている。あまり計算時間が短いと時計の精度も悪いので同じかけ算をntimes回繰り返すことにした。(5)のループはスカラユニットのかけ算、(6)のループはvDSPを直接呼んでいる。時間を計るのはunixの古いシステムコールを使っている。(*)のある行が時間を計るしかけ。これについては別途メモを作る。とりあえず、プロセスのユーザ時間で計測ができるようにした。

全く同じことを倍精度でやるdoubleCompare()と言う関数（ntimesは20回に減らした）も書いてmain()で呼ぶ。

int main (int argc, const char * argv[])
{
    floatCompare();
    doubleCompare();
    return 0;
}

こんだけ。

結果

これをPowerBook G4 1.33GHzでビルドして実行してみた。Xcode3.0のgcc4.0.1でビルド構成はppcのreleaseのすっぴん。で、これがそのコンソールアウト。

[Session started at 2008-08-14 17:06:34 +0900.]
Scalar : 302	vector : 152	Ratio = 1.986842
Scalar : 154	vector : 156	Ratio = 0.987179
The Debugger has exited with status 0.

時間の単位はよくわからないけど、昔のままなら10ミリ秒。だから単精度ではスカラユニットを使うと約3秒で、vDSPだと約1.5秒という意味。他に重いプロセスはなくて、実質的にこれがCPUを占有して、単精度不動小数点の場合でだいたい2倍速いという結果。何回か繰り返してみたけど、それほど変わらなかった。

PowerPCのベクタユニットは倍精度浮動小数点演算の機能はないので、vDSPを呼んでもパフォーマンスは変わらない（オーバーヘッドと思われる分だけvDSPの方が遅い）。純粋に単純計算の比較ならま、こんなもんだろ。

Intel iMacでの結果

さてこれを女房が新しく買ったiMac2.4GHzでやってみた。Xcode3.0のgcc4.0.1でビルド構成はi386のrelease。外側ループのntimesの値は10倍にした。

[Session started at 2008-08-15 07:58:27 +0900.]
Scalar : 152	vector : 141	Ratio = 1.078014
Scalar : 116	vector : 120	Ratio = 0.966667
The Debugger has exited with status 0.

あれ？ほとんどかわんない（どうでもいいけど、PowerBook G4は単純計算では1/10のパフォーマンスしかないということ？ちょっと寂しい）。どういうこと？

しょうがない、生まれてはじめてXcodeの「アセンブラコードの表示」メニューで、アセンブラを表示させてみると、PowerPCの方は

	fmuls	f0,f1,f13

などとスカラユニットの浮動小数点レジスタを使う命令を生成しているのに対して、Intelの方は

	mulss	(%ebx,%eax,4), %xmm0

などのSSE命令が生成されている。

最適化レベル（GCC_OPTIMIZATION_LEVEL）を「最適化なし（-O0）」にしてやりなおすと

[Session started at 2008-08-15 08:02:18 +0900.]
Scalar : 336	vector : 139	Ratio = 2.417266
Scalar : 151	vector : 120	Ratio = 1.258333
The Debugger has exited with status 0.

となって、単精度の方はG4と同じ程度の差がついた。倍精度の方は中途半端な差しかない。これはメモリアクセスが律速してるのか。よくわからない。

これ以上調べるのも面倒なのでもうやめるけど、ようするにgccはIntelでは可能な場合にベクタユニットを使うコードを出力するということのようだ。

gcc最適化とvDSP

gccはgnuの中心的な存在で、この品質の高さがgnuのあらゆるソフトウェアを支えている。gccは長年に渡って改善が続けられて、叩き上げられてきた。gccのIntel CPU向けのコードジェネレータは自動的にSSEのコードを出力するようになっているようである。きっとこんなことはLinuxなどのMacOS X以外のプラットフォームでgccを使ってる人たちにはごく当たり前のことなんだろうな。gccのリリースノートをちゃんと見ればきっと書いてあるんだろう。

vDSPのほうは実はかなり古いライブラリでAccelerateフレームワークに取り込まれる前に、PowerPC G4のAltiVecベクトルユニットを使いこなすために独立したライブラリとしてOS9.1/10.0の時代から存在していた。あきらかに当時のコンパイラの最適化はベクタユニットを使うところまで考慮できず、ライブラリを用意してプログラマに任せた、というのがvDSPの位置づけだろう。

ということで、少なくともIntel Macで今回のベンチマークのような単純計算にvDSPを使う意味はまったくない。なんじゃそりゃ。

gccでのベクトル化

もちろん、gccがなんでもかんでもベクトルユニットを使うようになっているわけではない。Auto-vectorization in GCCにはベクトル化できるループとできないループの例が書いてある。例えば

  for (i=0; i<256; i++){
    a[i] = b[i] + c[i];
  }

はできまっせ、これは今回のと同じで、実際にちゃんとベクトル化できてる。もっと複雑な

for (i = 0; i < N/2; i++){
  a[i] = b[2*i+1] * c[2*i+1] - b[2*i] * c[2*i];
  d[i] = b[2*i] * c[2*i+1] + b[2*i+1] * c[2*i];
}

もOKと書いてある。また、

   while (n--){
      *p++ = *q++;
   }

もOKで、ループに入る前にnの値が確定しているからできるということらしい。

一方でできないループとして

while (*p != NULL) {
  *q++ = *p++;
}

が書いてある。これはループに入るときに回数がわからないかららしい。そらそうだわな。また、

for (i = 0; i < N; i++) {
  a[i] = i;
}

もダメとなってるけど、これはなんでだろう。

この最後のふたつなんかはごく普通にCで書くので、gccのベクトル化も限られたことに対応できるだけだと考えられる。このgccのベクトル化はSSA（Static Single Assignment）という解析手法がgccに導入されたことで実現できたらしい。変数の名前をつけ直すだけで最適化が実現しやすくなるなんて、なかなかすごいアイデアじゃん。昔のFORTRANでの積み重ねがこういうところに生かされてるのね。ドキュメントを見てると面白い。

ここでの結論

これはgcc4.0かららしいのでXcodeユーザにとってはMacOS X10.5から、ということになる。でもこういった最適化は、最終的なコードジェネレーションとは無関係の、ずっと上流の話のはず。なんでIntelのSSEでできて、PowerPCのAltiVecでできないの？コードジェネレータがそういう最適化情報を利用していないということやろな。同じOSでこういう差がでるのは美しくないけど、もうPowerPCに工数をかけたくないということの現れか。

ということで、もともとvDSPはFFTや畳み込みを計算するためのライブラリで、他の単純計算の関数などはそのユーティリティという位置づけであろうから、vDSP全部がいらないわけではないけど、Intel Macでは存在意義のなくなった関数が含まれているということになる。しばらくは互換性のために残るが、MacがPowerPCのサポートをやめてgccの最適化レベルが上がると、いずれは他のライブラリと統合してAccelerateフレームワークをアンブレラでなくしていくことになるのだろう。今は過渡的な時代ということか。いつの時代も過渡的だったような気はするけど。