9

最近、ホットスポットが間違いなくこれであるプログラムをプロファイリングしていました

double d = somevalue();
double d2=d*d;
double c = 1.0/d2   // HOT SPOT

値 c のみが必要なため、値 d2 は後で使用されません。少し前に、高速逆平方根のカーマック法について読んだことがありますが、これは明らかにそうではありませんが、同様のアルゴリズムが 1/x^2 の計算に役立つかどうか疑問に思っています。

非常に正確な精度が必要です。私のプログラムが gcc -ffast-math オプションで正しい結果を出さないことを確認しました。(g++-4.5)

4

3 に答える 3

19

高速な平方根などを行うためのトリックは、精度を犠牲にすることでパフォーマンスを得ています。(まぁ、そのほとんどです。)

  1. double本当に精度が必要ですか?精度を簡単に犠牲にすることができます:

    double d = somevalue();
    float c = 1.0f / ((float) d * (float) d);
    

    この1.0f場合、 は絶対に必須です。1.0代わりに使用すると、double精度が得られます。

  2. コンパイラで「ずさんな」計算を有効にしてみましたか? を使用できる GCC-ffast-mathでは、他のコンパイラにも同様のオプションがあります。ずさんな計算は、アプリケーションには十分すぎるかもしれません。(編集:結果のアセンブリに違いは見られませんでした。)

  3. GCC を使用している場合、使用を検討しました-mrecipか? 約 12 ビットの精度しかない「逆数推定」機能がありますが、はるかに高速です。Newton-Raphson 法を使用して、結果の精度を上げることができます。この-mrecipオプションを指定すると、逆数推定とニュートン ラフソン ステップがコンパイラによって自動的に生成されますが、パフォーマンスと精度のトレードオフを微調整したい場合は、いつでも自分でアセンブリを作成できます。(Newton-Raphsonは非常に迅速に収束します。) (編集: GCC で RCPSS を生成できませんでした。以下を参照してください。)

あなたが経験している正確な問題について議論しているブログ投稿 ( source ) を見つけました。著者の結論は、Carmack メソッドのような手法は RCPSS 命令 ( -mrecipGCC のフラグが使用する) と競合しないというものです。

除算が非常に遅くなる理由は、通常、プロセッサには除算ユニットが 1 つしかなく、パイプライン化されていないことが多いためです。したがって、パイプ内でいくつかの乗算をすべて同時に実行できますが、前の除算が終了するまで除算を発行することはできません。

うまくいかないトリック

  1. カーマックの方法: 相反推定オペコードを持つ最近のプロセッサでは時代遅れです。逆数の場合、私が見た中で最高のバージョンは 1 ビットの精度しか提供しません。12 ビットのRCPSS. このトリックが平方根の逆数でうまく機能するのは偶然だと思います。二度とない偶然。

  2. 変数の再ラベル付け。コンパイラに関する限り、 と の間にほとんど違いはありませ1.0/(x*x)double x2 = x*x; 1.0/x2。最適化を最低レベルまでオンにして、2 つのバージョンで異なるコードを生成するコンパイラを見つけたら、私は驚かれることでしょう。

  3. を使用してpowいます。powライブラリ機能は総モンスターです。GCC-ffast-mathがオフになっていると、ライブラリの呼び出しはかなり高価になります。GCC を-ffast-mathオンにすると、 の場合とまったく同じアセンブリ コードが得pow(x, -2)られるため1.0/(x*x)、メリットはありません。

アップデート

倍精度浮動小数点値の逆二乗のニュートン ラフソン近似の例を次に示します。

static double invsq(double x)
{
    double y;
    int i;
    __asm__ (
        "cvtpd2ps %1, %0\n\t"
        "rcpss %0, %0\n\t"
        "cvtps2pd %0, %0"
        : "=x"(y)
        : "x"(x));
    for (i = 0; i < RECIP_ITER; ++i)
        y *= 2 - x * y;
    return y * y;
}

残念ながら、RECIP_ITER=1私のコンピューターのベンチマークでは、単純なバージョンよりもわずかに遅くなりました (~5%) 1.0/(x*x)。反復がゼロの場合は高速 (2 倍) ですが、12 ビットの精度しか得られません。12ビットで十分かどうかはわかりません。

ここでの問題の 1 つは、マイクロ最適化が小さすぎることだと思います。この規模では、コンパイラの作成者はアセンブリ ハッカーとほぼ同等の立場にあります。より大きな全体像があれば、それを高速化する方法がわかるかもしれません。

たとえば、それ-ffast-mathが望ましくない精度の低下を引き起こしたとあなたは言いました。これは、使用しているアルゴリズムの数値安定性の問題を示している可能性があります。floatアルゴリズムを正しく選択すれば、多くの問題はの代わりに解決できますdouble。(もちろん、24 ビット以上が必要な場合もあります。わかりません。)

RCPSSこれらのいくつかを並行して計算したい場合、この方法が優れていると思います。

于 2012-03-16T11:57:17.507 に答える
5

はい、確かに何かを試してみることができます。一般的なアイデアをいくつか紹介します。詳細を記入してください。

まず、Carmack のルートが機能する理由を見てみましょう。

通常の方法でx  =  M  × 2 Eと書きます。ここで、IEEE float がバイアスによって指数オフセットを格納することを思い出してください。e が指数フィールドを表す場合 e = Bias +  E  ≥ 0 となります。並べ替えると、E  = e − Bias が得られます。

次に逆平方根: x −1/2  =  M -1/2  × 2 E /2 . 新しい指数フィールドは次のとおりです。

       e'  = バイアス −  E /2 = 3/2 バイアス − e/2

ビットをいじると、シフトによってeから値e /2 を取得できます。3/2バイアスは単なる定数です。

さらに、仮数Mは 1.0 +  x ( x < 1) として格納され、 M -1/2を 1 + x/2 として 近似できます。繰り返しますが、xのみが 2 進数で格納されているという事実は、単純なビット シフトによって 2 で除算することを意味します。


ここでx −2を見ます: これはM −2  × 2 −2 Eに等しく、指数体を探しています:

       e'  = バイアス − 2  E  = 3 バイアス − 2  e

繰り返しますが、3 Bias は単なる定数であり、ビットシフトによって e から 2 e を得ること できます仮数に関しては、(1 + x) −2を 1 − 2  xで近似できるため、問題は x から 2 x を求めること なります


Carmack の魔法の浮動小数点フィドリングは、実際にはすぐに結果を計算するわけではないことに注意してください。むしろ、従来の反復計算の開始点として使用される、非常に正確な推定値を生成します。しかし、推定値が非常に優れているため、許容できる結果を得るために必要な反復回数はごくわずかです。

于 2012-03-16T11:57:09.930 に答える
1

現在のプログラムでは、ホットスポットを特定しました - 良いです。1/d^2 を高速化する代わりに、1/d^2 を頻繁に計算しないようにプログラムを変更するオプションがあります。内側のループから巻き上げることはできますか? 1/d^2 を計算する d の値はいくつありますか? 必要なすべての値を事前に計算してから、結果を調べていただけますか? これは 1/d^2 の場合は少し面倒ですが、1/d^2 がより大きなコードの一部である場合は、このトリックをそれに適用する価値があるかもしれません。精度を下げると、十分な答えが得られないとおっしゃっています。より良い動作を提供する可能性のあるコードを言い換えることができる方法はありますか? 数値解析は非常に微妙なので、いくつかのことを試して何が起こったのかを確認する価値があるかもしれません。

もちろん、理想的には、何年にもわたる研究に基づいて最適化されたルーチンを見つけることです.lapackまたはlinpackにリンクできるものはありますか?

于 2012-03-16T18:21:34.200 に答える