[解決済み] C++で長方形方程式を実装する際に、高レベルのアプローチでパフォーマンスを向上させるには？

質問

私はあるエンジニアリング・シミュレーションを開発しています。これは、ゴムのような材料の応力を計算するために、この方程式のようないくつかの長い方程式を実装することを含んでいます。

T = (
    mu * (
            pow(l1 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a
            * (
                pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
                - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1
            ) * pow(l1 * l2 * l3, 0.1e1 / 0.3e1) / l1
            - pow(l2 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l1 / 0.3e1
            - pow(l3 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l1 / 0.3e1
        ) / a
    + K * (l1 * l2 * l3 - 0.1e1) * l2 * l3
) * N1 / l2 / l3

+ (
    mu * (
        - pow(l1 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l2 / 0.3e1
        + pow(l2 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a
        * (
            pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
            - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1
        ) * pow(l1 * l2 * l3, 0.1e1 / 0.3e1) / l2
        - pow(l3 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l2 / 0.3e1
    ) / a
    + K * (l1 * l2 * l3 - 0.1e1) * l1 * l3
) * N2 / l1 / l3

+ (
    mu * (
        - pow(l1 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l3 / 0.3e1
        - pow(l2 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a / l3 / 0.3e1
        + pow(l3 * pow(l1 * l2 * l3, -0.1e1 / 0.3e1), a) * a
        * (
            pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
            - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1
        ) * pow(l1 * l2 * l3, 0.1e1 / 0.3e1) / l3
    ) / a
+ K * (l1 * l2 * l3 - 0.1e1) * l1 * l2
) * N3 / l1 / l2;

私は間違いを避けるために（そして面倒な代数の時間を節約するために）C++のコードを生成するためにMapleを使用しています。このコードは何千回も (何百万回も) 実行されるため、パフォーマンスが懸念されます。残念ながら、数学はここまでしか単純化できず、長い方程式は避けられません。

この実装を最適化するために、どのようなアプローチを取ることができますか? 私は、このような方程式を実装する際に適用すべき高レベルの戦略を探しているのであって、必ずしも上記の例に対する特定の最適化ではありません。

私は g++ でコンパイルしています。 --enable-optimize=-O3 .

更新してください。

私は多くの繰り返される式があることを知っています。私はコンパイラがこれらを処理するという仮定を使用しています。

l1, l2, l3, mu, a, K はすべて正の実数です (ゼロではありません)。

を置き換えた l1*l2*l3 を同等の変数に置き換えます。 J . これは、パフォーマンスの向上に役立ちました。

置き換えの pow(x, 0.1e1/0.3e1) を cbrt(x) は良い提案でした。

これはCPU上で実行されます。近い将来、これはGPU上でより良く実行されると思われますが、今のところそのオプションはありません。

どのように解決するのですか？

編集の概要

• 私の最初の回答は、コードに多くの複製された計算が含まれており、多くの累乗が 1/3 の係数を含んでいることを指摘したにすぎません。たとえば pow(x, 0.1e1/0.3e1) と同じです。 cbrt(x) .
• 私の2回目の編集はただ間違っており、3回目の編集はこの間違っていることを外挿したものでした。これは、人々が 'M' で始まる記号的な数学プログラムからの神託のような結果を変更することを恐れているものです。私は打ち消した（つまり。 <ストライク ストライク ) それらの編集を取り消し、この回答の現在のリビジョンの一番下に押し込みました。しかし、私はそれらを削除しませんでした。私は人間です。間違いを犯すのは簡単なことです。
• 私の4回目の編集は、質問の複雑な式を正しく表す非常にコンパクトな式を開発しました。 IF パラメータ l1 , l2 そして l3 は正の実数であり、もし a が0でない実数である場合。(これらの係数の具体的な性質については、まだOPから聞いていない。問題の性質を考えると、これらは妥当な仮定です)。
• この編集は、これらの式をどのように簡略化するかという一般的な問題に答えようとするものです。

はじめに

C++のコードを生成するのにMapleを使うのは、ミスを防ぐためです。

MapleとMathematicaは時々明らかな間違いをします。さらに重要なことは，MapleとMathematicaのユーザが間違いを犯すことがあることです．時々」の代わりに「時々」、あるいは「ほとんどいつも」を代用することが、おそらく的を得ていると思われます。

あなたは、問題のパラメータについて伝えることによって、Mapleがその式を簡略化するのを助けることができたでしょう。手元の例では、私は l1 , l2 そして l3 は正の実数であり a はゼロでない実数である。もしそうであれば、そう教えてあげてください。これらの記号計算プログラムは通常、手元の数量が複素数であると仮定します。領域を制限することで、プログラムは複素数では有効でない仮定をすることができます。

記号計算プログラムからの大きな混乱を簡略化する方法 (この記事です)

記号数学プログラムは通常、様々なパラメータに関する情報を提供する機能を備えています。特に、問題が除算や指数を含む場合は、その機能を使用してください。手元の例では、Maple が式を単純化するのを助けるために、次のように伝えることができたでしょう。 l1 , l2 と l3 は正の実数であり a はゼロでない実数である。もしそうであれば、そう教えてあげてください。これらの記号演算プログラムは通常、手元の数量が複素数であることを前提としています。領域を制限することで、プログラムは次のような仮定をすることができます。 ^x b ^x =(ab) ^x . これは、以下の場合にのみ a と b は正の実数であり、もし x は実数です。複素数では成立しません。

結局のところ、それらの記号的な数学プログラムはアルゴリズムに従うのです。それを手助けしてください。コードを生成する前に、展開、収集、および簡略化で遊んでみてください。この例では、以下の因子を含むこれらの項を収集することができました。 mu の因子を含むものと K . 式を最も単純な形にすることは、ちょっとした芸術です。

生成されたコードが醜いことになっても、それを「触れてはいけない真実」として受け止めてはいけません。自分で単純化することを試みましょう。コードを生成する前に、記号数学のプログラムが持っていたものを見てください。私があなたの式をもっと単純でもっと速いものにしたかを見てください。 ウォルターの答え は私の式をさらに数歩進めたのです。魔法のレシピはありません。もし魔法のレシピがあれば、メイプルはそれを適用してウォルターが出したような答えを出したでしょう。

具体的な質問について

あなたはその計算で多くの足し算と引き算をしています。互いにほぼ相殺される項があると、深い問題に陥る可能性があります。1つの項が他を圧倒している場合、多くのCPUを浪費していることになります。

次に、繰り返し計算を行うことによって、多くの CPU を浪費しています。有効にしていない限り -ffast-math を有効にしていない限り、コンパイラはその式を単純化することはありません（実際、してはいけません）。その代わり、あなたが指示したとおりの処理をします。最低でも l1 * l2 * l3 を計算する必要があります。

最後に、あなたはたくさんの呼び出しを pow を大量に呼び出しており、これは非常に遅いです。これらの呼び出しのいくつかは (l1*l2*l3) という形式であることに注意してください。 ^(1/3) . これらの呼び出しの多くは pow への呼び出しを一回で済ませることができます。 std::cbrt :

l123 = l1 * l2 * l3;
l123_pow_1_3 = std::cbrt(l123);
l123_pow_4_3 = l123 * l123_pow_1_3;

これを使って

• X * pow(l1 * l2 * l3, 0.1e1 / 0.3e1) となる X * l123_pow_1_3 .
• X * pow(l1 * l2 * l3, -0.1e1 / 0.3e1) となる X / l123_pow_1_3 .
• X * pow(l1 * l2 * l3, 0.4e1 / 0.3e1) となる X * l123_pow_4_3 .
• X * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) となる X / l123_pow_4_3 .

楓は明らかなことを見逃していました。

例えば、もっと簡単な書き方があります。

(pow(l1 * l2 * l3, -0.1e1 / 0.3e1) - l1 * l2 * l3 * pow(l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1)

と仮定すると l1 , l2 そして l3 が複素数ではなく実数であり、（原理的な複素数根ではなく）実数の立方根が抽出されるとすると、上記は以下のようになります。

2.0/(3.0 * pow(l1 * l2 * l3, 1.0/3.0))

または

2.0/(3.0 * l123_pow_1_3)

使用方法 cbrt_l123 の代わりに l123_pow_1_3 に変更すると、質問中の厄介な式は次のようになります。

l123 = l1 * l2 * l3; 
cbrt_l123 = cbrt(l123);
T = 
  mu/(3.0*l123)*(  pow(l1/cbrt_l123,a)*(2.0*N1-N2-N3)
                 + pow(l2/cbrt_l123,a)*(2.0*N2-N3-N1)
                 + pow(l3/cbrt_l123,a)*(2.0*N3-N1-N2))
 +K*(l123-1.0)*(N1+N2+N3);

常にダブルチェック、しかし常にシンプルに。

以下に、私が上記に至るまでに行ったことの一部を紹介します。

// Step 0: Trim all whitespace.
T=(mu*(pow(l1*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a*(pow(l1*l2*l3,-0.1e1/0.3e1)-l1*l2*l3*pow(l1*l2*l3,-0.4e1/0.3e1)/0.3e1)*pow(l1*l2*l3,0.1e1/0.3e1)/l1-pow(l2*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l1/0.3e1-pow(l3*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l1/0.3e1)/a+K*(l1*l2*l3-0.1e1)*l2*l3)*N1/l2/l3+(mu*(-pow(l1*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l2/0.3e1+pow(l2*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a*(pow(l1*l2*l3,-0.1e1/0.3e1)-l1*l2*l3*pow(l1*l2*l3,-0.4e1/0.3e1)/0.3e1)*pow(l1*l2*l3,0.1e1/0.3e1)/l2-pow(l3*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l2/0.3e1)/a+K*(l1*l2*l3-0.1e1)*l1*l3)*N2/l1/l3+(mu*(-pow(l1*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l3/0.3e1-pow(l2*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a/l3/0.3e1+pow(l3*pow(l1*l2*l3,-0.1e1/0.3e1),a)*a*(pow(l1*l2*l3,-0.1e1/0.3e1)-l1*l2*l3*pow(l1*l2*l3,-0.4e1/0.3e1)/0.3e1)*pow(l1*l2*l3,0.1e1/0.3e1)/l3)/a+K*(l1*l2*l3-0.1e1)*l1*l2)*N3/l1/l2;

// Step 1:
//   l1*l2*l3 -> l123
//   0.1e1 -> 1.0
//   0.4e1 -> 4.0
//   0.3e1 -> 3
l123 = l1 * l2 * l3;
T=(mu*(pow(l1*pow(l123,-1.0/3),a)*a*(pow(l123,-1.0/3)-l123*pow(l123,-4.0/3)/3)*pow(l123,1.0/3)/l1-pow(l2*pow(l123,-1.0/3),a)*a/l1/3-pow(l3*pow(l123,-1.0/3),a)*a/l1/3)/a+K*(l123-1.0)*l2*l3)*N1/l2/l3+(mu*(-pow(l1*pow(l123,-1.0/3),a)*a/l2/3+pow(l2*pow(l123,-1.0/3),a)*a*(pow(l123,-1.0/3)-l123*pow(l123,-4.0/3)/3)*pow(l123,1.0/3)/l2-pow(l3*pow(l123,-1.0/3),a)*a/l2/3)/a+K*(l123-1.0)*l1*l3)*N2/l1/l3+(mu*(-pow(l1*pow(l123,-1.0/3),a)*a/l3/3-pow(l2*pow(l123,-1.0/3),a)*a/l3/3+pow(l3*pow(l123,-1.0/3),a)*a*(pow(l123,-1.0/3)-l123*pow(l123,-4.0/3)/3)*pow(l123,1.0/3)/l3)/a+K*(l123-1.0)*l1*l2)*N3/l1/l2;

// Step 2:
//   pow(l123,1.0/3) -> cbrt_l123
//   l123*pow(l123,-4.0/3) -> pow(l123,-1.0/3)
//   (pow(l123,-1.0/3)-pow(l123,-1.0/3)/3) -> 2.0/(3.0*cbrt_l123)
//   *pow(l123,-1.0/3) -> /cbrt_l123
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T=(mu*(pow(l1/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l1-pow(l2/cbrt_l123,a)*a/l1/3-pow(l3/cbrt_l123,a)*a/l1/3)/a+K*(l123-1.0)*l2*l3)*N1/l2/l3+(mu*(-pow(l1/cbrt_l123,a)*a/l2/3+pow(l2/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l2-pow(l3/cbrt_l123,a)*a/l2/3)/a+K*(l123-1.0)*l1*l3)*N2/l1/l3+(mu*(-pow(l1/cbrt_l123,a)*a/l3/3-pow(l2/cbrt_l123,a)*a/l3/3+pow(l3/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l3)/a+K*(l123-1.0)*l1*l2)*N3/l1/l2;

// Step 3:
//   Whitespace is nice.
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  (mu*( pow(l1/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l1
       -pow(l2/cbrt_l123,a)*a/l1/3
       -pow(l3/cbrt_l123,a)*a/l1/3)/a
   +K*(l123-1.0)*l2*l3)*N1/l2/l3
 +(mu*(-pow(l1/cbrt_l123,a)*a/l2/3
       +pow(l2/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l2
       -pow(l3/cbrt_l123,a)*a/l2/3)/a
   +K*(l123-1.0)*l1*l3)*N2/l1/l3
 +(mu*(-pow(l1/cbrt_l123,a)*a/l3/3
       -pow(l2/cbrt_l123,a)*a/l3/3
       +pow(l3/cbrt_l123,a)*a*2.0/(3.0*cbrt_l123)*cbrt_l123/l3)/a
   +K*(l123-1.0)*l1*l2)*N3/l1/l2;

// Step 4:
//   Eliminate the 'a' in (term1*a + term2*a + term3*a)/a
//   Expand (mu_term + K_term)*something to mu_term*something + K_term*something
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  (mu*( pow(l1/cbrt_l123,a)*2.0/(3.0*cbrt_l123)*cbrt_l123/l1
       -pow(l2/cbrt_l123,a)/l1/3
       -pow(l3/cbrt_l123,a)/l1/3))*N1/l2/l3
 +K*(l123-1.0)*l2*l3*N1/l2/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l2/3
       +pow(l2/cbrt_l123,a)*2.0/(3.0*cbrt_l123)*cbrt_l123/l2
       -pow(l3/cbrt_l123,a)/l2/3))*N2/l1/l3
 +K*(l123-1.0)*l1*l3*N2/l1/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l3/3
       -pow(l2/cbrt_l123,a)/l3/3
       +pow(l3/cbrt_l123,a)*2.0/(3.0*cbrt_l123)*cbrt_l123/l3))*N3/l1/l2
 +K*(l123-1.0)*l1*l2*N3/l1/l2;

// Step 5:
//   Rearrange
//   Reduce l2*l3*N1/l2/l3 to N1 (and similar)
//   Reduce 2.0/(3.0*cbrt_l123)*cbrt_l123/l1 to 2.0/3.0/l1 (and similar)
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  (mu*( pow(l1/cbrt_l123,a)*2.0/3.0/l1
       -pow(l2/cbrt_l123,a)/l1/3
       -pow(l3/cbrt_l123,a)/l1/3))*N1/l2/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l2/3
       +pow(l2/cbrt_l123,a)*2.0/3.0/l2
       -pow(l3/cbrt_l123,a)/l2/3))*N2/l1/l3
 +(mu*(-pow(l1/cbrt_l123,a)/l3/3
       -pow(l2/cbrt_l123,a)/l3/3
       +pow(l3/cbrt_l123,a)*2.0/3.0/l3))*N3/l1/l2
 +K*(l123-1.0)*N1
 +K*(l123-1.0)*N2
 +K*(l123-1.0)*N3;

// Step 6:
//   Factor out mu and K*(l123-1.0)
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  mu*(  ( pow(l1/cbrt_l123,a)*2.0/3.0/l1
         -pow(l2/cbrt_l123,a)/l1/3
         -pow(l3/cbrt_l123,a)/l1/3)*N1/l2/l3
      + (-pow(l1/cbrt_l123,a)/l2/3
         +pow(l2/cbrt_l123,a)*2.0/3.0/l2
         -pow(l3/cbrt_l123,a)/l2/3)*N2/l1/l3
      + (-pow(l1/cbrt_l123,a)/l3/3
         -pow(l2/cbrt_l123,a)/l3/3
         +pow(l3/cbrt_l123,a)*2.0/3.0/l3)*N3/l1/l2)
 +K*(l123-1.0)*(N1+N2+N3);

// Step 7:
//   Expand
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  mu*( pow(l1/cbrt_l123,a)*2.0/3.0/l1*N1/l2/l3
      -pow(l2/cbrt_l123,a)/l1/3*N1/l2/l3
      -pow(l3/cbrt_l123,a)/l1/3*N1/l2/l3
      -pow(l1/cbrt_l123,a)/l2/3*N2/l1/l3
      +pow(l2/cbrt_l123,a)*2.0/3.0/l2*N2/l1/l3
      -pow(l3/cbrt_l123,a)/l2/3*N2/l1/l3
      -pow(l1/cbrt_l123,a)/l3/3*N3/l1/l2
      -pow(l2/cbrt_l123,a)/l3/3*N3/l1/l2
      +pow(l3/cbrt_l123,a)*2.0/3.0/l3*N3/l1/l2)
 +K*(l123-1.0)*(N1+N2+N3);

// Step 8:
//   Simplify.
l123 = l1 * l2 * l3;
cbrt_l123 = cbrt(l123);
T =
  mu/(3.0*l123)*(  pow(l1/cbrt_l123,a)*(2.0*N1-N2-N3)
                 + pow(l2/cbrt_l123,a)*(2.0*N2-N3-N1)
                 + pow(l3/cbrt_l123,a)*(2.0*N3-N1-N2))
 +K*(l123-1.0)*(N1+N2+N3);

誤答、謙遜のため意図的に残しています。

これは削除されていることに注意してください。間違っているのです。

更新

l123 = l1 * l2 * l3; 
cbrt_l123_inv = 1.0 / cbrt(l123);
nasty_expression =
    K * (l123 - 1.0) * (N1 + N2 + N3) 
    - (  pow(l1 * cbrt_l123_inv, a) * (N2 + N3) 
       + pow(l2 * cbrt_l123_inv, a) * (N1 + N3) 
       + pow(l3 * cbrt_l123_inv, a) * (N1 + N2)) * mu / (3.0*l123);