Como posso melhorar o desempenho por meio de uma abordagem de alto nível ao implementar equações longas em C ++

Estou desenvolvendo algumas simulações de engenharia. Isso envolve a implementação de algumas equações longas, como esta equação para calcular a tensão em um material semelhante a borracha:

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;

Eu uso o Maple para gerar o código C ++ para evitar erros (e economizar tempo com álgebra tediosa). Como esse código é executado milhares (senão milhões) de vezes, o desempenho é uma preocupação. Infelizmente, a matemática só simplifica até agora; as longas equações são inevitáveis.

Que abordagem posso adotar para otimizar essa implementação? Estou procurando estratégias de alto nível que devo aplicar ao implementar tais equações, não necessariamente otimizações específicas para o exemplo mostrado acima.

Estou compilando usando g ++ com --enable-optimize=-O3.

Atualizar:

Eu sei que há muitas expressões repetidas, estou supondo que o compilador iria lidar com elas; meus testes até agora sugerem que sim.

l1, l2, l3, mu, a, K são todos números reais positivos (não zero).

Tenho substituído l1*l2*l3com uma variável equivalente: J. Isso ajudou a melhorar o desempenho.

Substituir pow(x, 0.1e1/0.3e1)por cbrt(x)foi uma boa sugestão.

Ele será executado em CPUs. Em um futuro próximo, provavelmente funcionará melhor em GPUs, mas por enquanto essa opção não está disponível.

Editar resumo

• Minha resposta original apenas observou que o código continha muitos cálculos replicados e que muitos dos poderes envolviam fatores de 1/3. Por exemplo, pow(x, 0.1e1/0.3e1)é o mesmo que cbrt(x).
• Minha segunda edição estava errada, e minha terceira extrapolou sobre esse erro. Isso é o que deixa as pessoas com medo de mudar os resultados semelhantes a oráculos de programas matemáticos simbólicos que começam com a letra 'M'. Risquei (ou seja, risquei ) essas edições e empurrei-as para o final da revisão atual desta resposta. No entanto, não os eliminei. Eu sou humano. É fácil cometermos um erro.
• Minha quarta edição desenvolvido uma expressão muito compacto que representa corretamente a expressão complicada na questão IF os parâmetros l1, l2e l3são números reais positivos e se aé um número real não-zero. (Ainda não ouvimos do OP sobre a natureza específica desses coeficientes. Dada a natureza do problema, essas são suposições razoáveis.)
• Esta edição tenta responder ao problema genérico de como simplificar essas expressões.

Primeiras coisas primeiro

Eu uso o Maple para gerar o código C ++ para evitar erros.

Maple e Mathematica às vezes perdem o óbvio. Ainda mais importante, os usuários do Maple e do Mathematica às vezes cometem erros. Substituir "frequentemente", ou talvez até "quase sempre", em vez de "às vezes é provavelmente mais perto do que acertar.

Você poderia ter ajudado o Maple a simplificar essa expressão, contando sobre os parâmetros em questão. No exemplo em questão, eu suspeito que l1, l2e l3são números reais positivos e esse aé um número real diferente de zero. Se for esse o caso, diga isso. Esses programas matemáticos simbólicos normalmente assumem que as quantidades disponíveis são complexas. Restringir o domínio permite que o programa faça suposições que não são válidas nos números complexos.

Como simplificar essas grandes bagunças de programas de matemática simbólica (esta edição)

Os programas de matemática simbólica normalmente fornecem a capacidade de fornecer informações sobre os vários parâmetros. Use essa habilidade, principalmente se o seu problema envolver divisão ou exponenciação. No exemplo em mãos, você poderia ter ajudado a bordo simplificar essa expressão, dizendo-se que l1, l2e l3são números reais positivos e que aé um número real não-zero. Se for esse o caso, diga isso. Esses programas matemáticos simbólicos normalmente assumem que as quantidades disponíveis são complexas. Restringir o domínio permite que o programa faça suposições como a ^x b ^x = (ab) ^x . Isso ocorre apenas se ae bforem números reais positivos e se xfor real. Não é válido nos números complexos.

Em última análise, esses programas matemáticos simbólicos seguem algoritmos. Ajude-o junto. Tente expandir, coletar e simplificar antes de gerar o código. Nesse caso, você poderia ter coletado os termos que envolvem um fator de mue aqueles que envolvem um fator de K. Reduzir uma expressão à sua "forma mais simples" continua sendo um pouco uma arte.

Quando você receber uma bagunça feia de código gerado, não aceite isso como uma verdade que você não deve tocar. Tente simplificar você mesmo. Veja o que o programa matemático simbólico tinha antes de gerar o código. Veja como eu reduzi sua expressão a algo muito mais simples e rápido, e como a resposta de Walter levou a minha vários passos adiante. Não existe receita mágica. Se houvesse uma receita mágica, Maple a teria aplicado e dado a resposta que Walter deu.

Sobre a questão específica

Você está fazendo muitas adições e subtrações nesse cálculo. Você pode entrar em apuros se tiver termos que quase cancelam um ao outro. Você está desperdiçando muita CPU se tiver um termo que domina os outros.

Em seguida, você está desperdiçando muita CPU ao realizar cálculos repetidos. A menos que você tenha ativado -ffast-math, o que permite ao compilador quebrar algumas das regras do ponto flutuante IEEE, o compilador não irá (na verdade, não deve) simplificar essa expressão para você. Em vez disso, fará exatamente o que você disse para fazer. No mínimo, você deve calcular l1 * l2 * l3antes de computar essa bagunça.

Por fim, você está fazendo muitas chamadas para pow, o que é extremamente lento. Observe que várias dessas chamadas são no formato (l1 * l2 * l3) ^(1/3) . Muitas dessas chamadas para powpoderiam ser realizadas com uma única chamada para std::cbrt:

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

Com isso,

• X * pow(l1 * l2 * l3, 0.1e1 / 0.3e1)torna-se X * l123_pow_1_3.
• X * pow(l1 * l2 * l3, -0.1e1 / 0.3e1)torna-se X / l123_pow_1_3.
• X * pow(l1 * l2 * l3, 0.4e1 / 0.3e1)torna-se X * l123_pow_4_3.
• X * pow(l1 * l2 * l3, -0.4e1 / 0.3e1)torna-se X / l123_pow_4_3.

Maple não percebeu o óbvio.
Por exemplo, há uma maneira muito mais fácil de escrever

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

Partindo do princípio de que l1, l2e l3são reais em vez de números complexos, e que a verdadeira raiz cubo (em vez da raiz complexo princípio) estão a ser extraído, o acima reduz a

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

ou

2.0/(3.0 * l123_pow_1_3)

Usando em cbrt_l123vez de l123_pow_1_3, a expressão desagradável na pergunta se reduz a

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);

Sempre verifique, mas sempre simplifique também.

Aqui estão alguns dos meus passos para chegar ao acima:

// 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);

Resposta errada, intencionalmente mantida para humildade

Observe que isso é atingido. Está errado.

Atualizar

Maple não percebeu o óbvio. Por exemplo, há uma maneira muito mais fácil de escrever

(pow (l1 * l2 * l3, -0,1e1 / 0,3e1) - l1 * l2 * l3 * pow (l1 * l2 * l3, -0,4e1 / 0,3e1) / 0,3e1)

Partindo do princípio de que l1, l2e l3são reais em vez de números complexos, e que a verdadeira raiz cubo (em vez da raiz complexo princípio) estão a ser extraído, o acima reduz a zero. Este cálculo de zero é repetido várias vezes.

Segunda atualização

Se eu fiz as contas certas (não há garantia de que fiz as contas certas), a expressão desagradável na pergunta se reduz a

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);

O acima assume que l1, l2e l3são números reais positivos.

A primeira coisa a notar é que powé muito caro, então você deve se livrar disso o máximo possível. Examinando a expressão, vejo muitas repetições de pow(l1 * l2 * l3, -0.1e1 / 0.3e1)e pow(l1 * l2 * l3, -0.4e1 / 0.3e1). Portanto, eu esperaria um grande ganho da pré-computação desses:

 const double c1 = pow(l1 * l2 * l3, -0.1e1 / 0.3e1);
const double c2 = boost::math::pow<4>(c1);

onde estou usando a função boost pow .

Além disso, você tem mais alguns powcom expoente a. Se afor inteiro e conhecido no momento do compilador, você também pode substituí-los boost::math::pow<a>(...)para obter mais desempenho. Também sugiro substituir termos como a / l1 / 0.3e1por, a / (l1 * 0.3e1)pois a multiplicação é mais rápida do que a divisão.

Finalmente, se você usar g ++, poderá usar o -ffast-mathsinalizador que permite que o otimizador seja mais agressivo na transformação de equações. Leia sobre o que este sinalizador realmente faz , pois tem efeitos colaterais.

Uau, que expressão infernal. Criar a expressão com Maple na verdade foi uma escolha abaixo do ideal aqui. O resultado é simplesmente ilegível.

1. escolha nomes de variáveis ​​falantes (não l1, l2, l3, mas, por exemplo, altura, largura, profundidade, se é isso que significam). Assim, será mais fácil entender seu próprio código.
2. calcule subtermos, que você usa várias vezes, antecipadamente e armazene os resultados em variáveis ​​com nomes falados.
3. Você menciona que a expressão é avaliada muitas vezes. Eu acho, apenas alguns parâmetros variam no loop mais interno. Calcule todos os subtermos invariantes antes desse loop. Repita para o segundo loop interno e assim por diante até que todos os invariantes estejam fora do loop.

Teoricamente, o compilador deveria ser capaz de fazer tudo isso para você, mas às vezes não pode - por exemplo, quando o aninhamento de loop se espalha por várias funções em unidades de compilação diferentes. De qualquer forma, isso lhe dará um código muito melhor legível, compreensível e sustentável.

A resposta de David Hammen é boa, mas ainda está longe de ser ideal. Vamos continuar com sua última expressão (no momento em que escrevo isso)

auto l123 = l1 * l2 * l3;
auto 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);

que pode ser otimizado ainda mais. Em particular, podemos evitar a chamada para cbrt()e uma das chamadas para pow()se explorar algumas identidades matemáticas. Vamos fazer isso novamente, passo a passo.

// step 1 eliminate cbrt() by taking the exponent into pow()
auto l123 = l1 * l2 * l3;
auto athird = 0.33333333333333333 * a; // avoid division
T = mu/(3.0*l123)*(  (N1+N1-N2-N3)*pow(l1*l1/(l2*l3),athird)
                   + (N2+N2-N3-N1)*pow(l2*l2/(l1*l3),athird)
                   + (N3+N3-N1-N2)*pow(l3*l3/(l1*l2),athird))
  + K*(l123-1.0)*(N1+N2+N3);

Observe que eu também otimizei 2.0*N1para N1+N1etc. Em seguida, podemos fazer com apenas duas chamadas para pow().

// step 2  eliminate one call to pow
auto l123 = l1 * l2 * l3;
auto athird = 0.33333333333333333 * a;
auto pow_l1l2_athird = pow(l1/l2,athird);
auto pow_l1l3_athird = pow(l1/l3,athird);
auto pow_l2l3_athird = pow_l1l3_athird/pow_l1l2_athird;
T = mu/(3.0*l123)*(  (N1+N1-N2-N3)* pow_l1l2_athird*pow_l1l3_athird
                   + (N2+N2-N3-N1)* pow_l2l3_athird/pow_l1l2_athird
                   + (N3+N3-N1-N2)/(pow_l1l3_athird*pow_l2l3_athird))
  + K*(l123-1.0)*(N1+N2+N3);

Uma vez que as chamadas para pow()são de longe a operação mais cara aqui, vale a pena reduzi-las ao máximo (a próxima operação cara foi a chamada para cbrt(), que eliminamos).

Se por acaso afor inteiro, as chamadas para powpoderiam ser otimizadas para chamadas para cbrt(mais potências inteiras), ou se athirdfor meio-inteiro, podemos usar sqrt(mais potências inteiras). Além disso, se por acaso l1==l2ou l1==l3ou l2==l3uma ou ambas as chamadas para powpode ser eliminado. Portanto, vale a pena considerá-los como casos especiais se tais chances existirem de forma realista.

1. Quantos são "muitos muitos"?
2. Quanto tempo leva?
3. Do TODAS parâmetros mudar entre recálculo desta fórmula? Ou você pode armazenar em cache alguns valores pré-calculados?

4. Já tentei uma simplificação manual dessa fórmula, gostaria de saber se salva alguma coisa?

   C1 = -0.1e1 / 0.3e1;
   C2 =  0.1e1 / 0.3e1;
   C3 = -0.4e1 / 0.3e1;
   
   X0 = l1 * l2 * l3;
   X1 = pow(X0, C1);
   X2 = pow(X0, C2);
   X3 = pow(X0, C3);
   X4 = pow(l1 * X1, a);
   X5 = pow(l2 * X1, a);
   X6 = pow(l3 * X1, a);
   X7 = a / 0.3e1;
   X8 = X3 / 0.3e1;
   X9 = mu / a;
   XA = X0 - 0.1e1;
   XB = K * XA;
   XC = X1 - X0 * X8;
   XD = a * XC * X2;
   
   XE = X4 * X7;
   XF = X5 * X7;
   XG = X6 * X7;
   
   T = (X9 * ( X4 * XD - XF - XG) / l1 + XB * l2 * l3) * N1 / l2 / l3 
     + (X9 * (-XE + X5 * XD - XG) / l2 + XB * l1 * l3) * N2 / l1 / l3 
     + (X9 * (-XE - XF + X6 * XD) / l3 + XB * l1 * l2) * N3 / l1 / l2;

[ADICIONADO] Trabalhei um pouco mais na fórmula das três últimas linhas e consegui essa beleza:

T = X9 / X0 * (
      (X4 * XD - XF - XG) * N1 + 
      (X5 * XD - XE - XG) * N2 + 
      (X5 * XD - XE - XF) * N3)
  + XB * (N1 + N2 + N3)

Deixe-me mostrar meu trabalho, passo a passo:

T = (X9 * (X4 * XD - XF - XG) / l1 + XB * l2 * l3) * N1 / l2 / l3 
  + (X9 * (X5 * XD - XE - XG) / l2 + XB * l1 * l3) * N2 / l1 / l3 
  + (X9 * (X5 * XD - XE - XF) / l3 + XB * l1 * l2) * N3 / l1 / l2;


T = (X9 * (X4 * XD - XF - XG) / l1 + XB * l2 * l3) * N1 / (l2 * l3) 
  + (X9 * (X5 * XD - XE - XG) / l2 + XB * l1 * l3) * N2 / (l1 * l3) 
  + (X9 * (X5 * XD - XE - XF) / l3 + XB * l1 * l2) * N3 / (l1 * l2);

T = (X9 * (X4 * XD - XF - XG) + XB * l1 * l2 * l3) * N1 / (l1 * l2 * l3) 
  + (X9 * (X5 * XD - XE - XG) + XB * l1 * l2 * l3) * N2 / (l1 * l2 * l3) 
  + (X9 * (X5 * XD - XE - XF) + XB * l1 * l2 * l3) * N3 / (l1 * l2 * l3);

T = (X9 * (X4 * XD - XF - XG) + XB * X0) * N1 / X0 
  + (X9 * (X5 * XD - XE - XG) + XB * X0) * N2 / X0 
  + (X9 * (X5 * XD - XE - XF) + XB * X0) * N3 / X0;

T = X9 * (X4 * XD - XF - XG) * N1 / X0 + XB * N1 
  + X9 * (X5 * XD - XE - XG) * N2 / X0 + XB * N2
  + X9 * (X5 * XD - XE - XF) * N3 / X0 + XB * N3;


T = X9 * (X4 * XD - XF - XG) * N1 / X0 
  + X9 * (X5 * XD - XE - XG) * N2 / X0
  + X9 * (X5 * XD - XE - XF) * N3 / X0
  + XB * (N1 + N2 + N3)

Isso pode ser um pouco conciso, mas eu realmente encontrei uma boa aceleração para polinômios (interpolação de funções de energia) usando Horner Form, que basicamente reescreve ax^3 + bx^2 + cx + dcomo d + x(c + x(b + x(a))). Isso evitará muitas chamadas repetidas para pow()e o impedirá de fazer coisas bobas, como ligar separadamente pow(x,6)e em pow(x,7)vez de apenas fazer x*pow(x,6).

Isso não se aplica diretamente ao seu problema atual, mas se você tiver polinômios de alta ordem com potências inteiras, isso pode ajudar. Você deve ter cuidado com a estabilidade numérica e os problemas de estouro, já que a ordem das operações é importante para isso (embora, em geral, eu realmente ache que a Forma de Horner ajuda nisso, já que x^20e xgeralmente estão em muitas ordens de magnitude).

Também como uma dica prática, se ainda não o fez, tente primeiro simplificar a expressão em maple. Você provavelmente pode fazer com que ele faça a maior parte da eliminação de subexpressão comum para você. Não sei o quanto isso afeta o gerador de código naquele programa em particular, mas sei que no Mathematica fazer um FullSimplify antes de gerar o código pode resultar em uma grande diferença.

Parece que há muitas operações repetidas em andamento.

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

Você pode pré-calculá-los para não chamar repetidamente a powfunção, o que pode ser caro.

Você também pode pré-calcular

l1 * l2 * l3

conforme você usa esse termo repetidamente.

Se você tiver uma placa de vídeo Nvidia CUDA, pode considerar transferir os cálculos para a placa de vídeo - que por si só é mais adequada para cálculos complicados computacionalmente.

https://developer.nvidia.com/how-to-cuda-c-cpp

Caso contrário, você pode querer considerar vários threads para cálculos.

Por acaso, você poderia fornecer o cálculo simbolicamente. Se houver operações de vetor, você pode realmente querer investigar o uso de blas ou lapack, que em alguns casos pode executar operações em paralelo.

É concebível (correndo o risco de sair do assunto?) Que você possa usar o python com numpy e / ou scipy. Na medida do possível, seus cálculos podem ser mais legíveis.

Como você perguntou explicitamente sobre otimizações de alto nível, pode valer a pena tentar diferentes compiladores C ++. Hoje em dia, os compiladores são bestas de otimização muito complexas e os fornecedores de CPU podem implementar otimizações muito poderosas e específicas. Mas observe que alguns deles não são gratuitos (mas pode haver um programa acadêmico gratuito).

• A coleção de compiladores GNU é gratuita, flexível e está disponível em muitas arquiteturas
• Os compiladores Intel são muito rápidos, muito caros e também podem produzir bons resultados para as arquiteturas AMD (acredito que haja um programa acadêmico)
• Os compiladores do Clang são rápidos, gratuitos e podem produzir resultados semelhantes ao GCC (algumas pessoas dizem que eles são mais rápidos, melhores, mas isso pode ser diferente para cada caso de aplicação, sugiro fazer suas próprias experiências)
• PGI (Portland Group) não é gratuito como os compiladores Intel.
• Os compiladores PathScale podem apresentar bons resultados em arquiteturas AMD

Já vi trechos de código diferirem na velocidade de execução por um fator de 2, apenas mudando o compilador (com otimizações completas, é claro). Mas esteja ciente de verificar a identidade da saída. A otimização agressiva pode levar a resultados diferentes, algo que você definitivamente deseja evitar.

Boa sorte!