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

92

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.

TylerH
fonte
32
Bem, a primeira coisa que vem à mente (a menos que o compilador o otimize sozinho) é substituir todos aqueles pow(l1 * l2 * l3, -0.1e1 / 0.3e1)por uma variável ... Você precisa avaliar seu código para ter certeza se ele executa rápido ou lento, no entanto.
SingerOfTheFall
6
Também formate o código para torná-lo mais legível - pode ajudar a identificar possibilidades de melhoria.
Ed Heal,
26
Por que todos os votos negativos e votos para fechar? Para aqueles de vocês que não gostam de programação numérica ou científica, vejam outras questões. Esta é uma boa pergunta que se adapta bem a este site. O site scicomp ainda é beta; migração não é uma boa opção. O site de revisão de código não tem olhos de sciomp suficientes. O que o OP fez acontece com bastante frequência na computação científica: construir um problema em um programa de matemática simbólica, pedir ao programa para gerar o código e não mexer no resultado porque o código gerado é uma bagunça.
David Hammen,
6
@DavidHammen, o site de revisão de código não tem olhos sciomp suficientes - soa como um problema da galinha e do ovo, e uma mentalidade que não está ajudando o CR a obter mais desses olhos. O mesmo se aplica à ideia de recusar o site beta do scicomp porque ele é beta - se todos pensassem assim, o único site a crescer seria o Stack Overflow.
Mathieu Guindon,
13
Esta questão está sendo discutida no meta aqui
NathanOliver

Respostas:

88

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 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.

David Hammen
fonte
2
Bem, a eliminação de CSE deve funcionar, independentemente da semântica relaxada (e OP indicado nos comentários que funciona). Embora, é claro, se for importante (medido), isso deve ser inspecionado (montagem gerada). Seus pontos sobre termos dominantes, simplificações de fórmulas perdidas, funções mais especializadas e os perigos de cancelamento são muito bons.
Deduplicator
3
@Deduplicator - Não com ponto flutuante. A menos que se habilite otimizações matemáticas inseguras (por exemplo, especificando -ffast-mathcom gcc ou clang), o compilador não pode confiar em pow(x,-1.0/3.0)ser igual a x*pow(x,-4.0/3.0). O último pode estar submerso, enquanto o primeiro não. Para ser compatível com o padrão de ponto flutuante, o compilador não deve otimizar esse cálculo para zero.
David Hammen,
Bem, esses são muito mais ambiciosos do que qualquer coisa que eu quis dizer.
Deduplicator
1
@Deduplicator: Como comentei em outra resposta : Você precisa -fno-math-errnodo g ++ para fazer powchamadas idênticas do CSE . (A menos que talvez possa provar que o pow não precisará definir errno?)
Peter Cordes
1
@Lefti - Aceite a resposta de Walter. É muito mais rápido. Há um problema potencial com todas essas respostas, que é o cancelamento numérico. Assumindo que o seu N1, N2e N3são não-negativo, um dos 2*N_i-(N_j+N_k)será negativo, um será positiva, eo outro será em algum lugar no meio. Isso pode facilmente resultar em problemas de cancelamento numérico.
David Hammen,
32

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.

Mariomulansky
fonte
5
Em nosso código, usar -ffast-mathleva o código a ficar instável ou dar respostas erradas. Temos um problema semelhante com os compiladores Intel e precisamos usar a -fp-model preciseopção, caso contrário, o código pode explodir ou fornecer as respostas erradas. Isso -ffast-mathpoderia ser mais rápido, mas eu recomendaria prosseguir com muito cuidado com essa opção, além dos efeitos colaterais listados em sua pergunta vinculada.
tpg2114
2
@ tpg2114: De acordo com meus testes, você só precisa-fno-math-errno que o g ++ seja capaz de içar chamadas idênticas para powfora de um loop. Essa é a parte menos "perigosa" de -ffast-math, para a maioria dos códigos.
Peter Cordes,
1
@PeterCordes Esses são resultados interessantes! Também tivemos problemas em pow ser extremamente lento e acabamos usando o dlsymhack mencionado nos comentários para obter aumentos de desempenho consideráveis, quando na verdade podíamos fazer com um pouco menos de precisão.
tpg2114
O GCC não entenderia que o pow é uma função pura? Isso provavelmente é conhecimento embutido.
usr
6
@ usr: Esse é o ponto, eu acho. nãopow é uma função pura, de acordo com o padrão, porque é para definirerrno em algumas circunstâncias. Definir sinalizadores como -fno-math-errnofazer com que ele não seja definido errno(violando o padrão), mas é uma função pura e pode ser otimizado como tal.
Nate Eldredge,
20

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.

cdonat
fonte
8
"o compilador deve fazer isso, mas às vezes não", é a chave aqui. além da legibilidade, é claro.
Javier,
3
Se o compilador não é obrigado a fazer algo, então assumir isso quase sempre está errado.
edmz
4
Escolha novamente os nomes das variáveis ​​faladas - muitas vezes essa regra legal não se aplica quando você está fazendo matemática. Ao examinar o código que supostamente implementa um algoritmo em um periódico científico, prefiro ver os símbolos no código exatamente aqueles usados ​​no periódico. Normalmente, isso significa nomes extremamente curtos, possivelmente com um subscrito.
David Hammen,
8
"O resultado é simplesmente ilegível" - por que isso é um problema? Você não se importaria se a saída de linguagem de alto nível de um gerador de lexer ou analisador fosse "ilegível" (por humanos). O que importa aqui é que a entrada para o gerador de código (Maple) seja legível e verificável. A única coisa a não fazer é editar o código gerado manualmente, se quiser ter certeza de que não há erros.
alephzero,
3
@DavidHammen: Bem, nesse caso, os de uma letra são os "nomes falados". Por exemplo, ao fazer geometria em um sistema de coordenadas cartesianas bidimensionais, xe nãoy são variáveis ​​de uma única letra sem sentido, são palavras inteiras com uma definição precisa e um significado bem e amplamente compreendido.
Jörg W Mittag
17

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.

Walter
fonte
@gnat Agradeço sua edição (pensei em fazer isso sozinho), mas teria achado mais justo, se a resposta de David também ligasse a esta. Por que você também não edita a resposta de David de maneira semelhante?
Walter
1
Só editei porque vi você mencionando isso explicitamente; Eu reli a resposta de David e não consegui encontrar nenhuma referência para sua resposta ali. Tento evitar edições onde não esteja 100% claro que as coisas que adicionei correspondem às intenções do autor
mosquito
1
@Walter - Minha resposta agora está vinculada à sua.
David Hammen,
1
Certamente não fui eu. Votei sua resposta há alguns dias. Eu também recebi um downvote flyby aleatório em minha resposta. Coisas simplesmente acontecem às vezes.
David Hammen
1
Você e eu recebemos votos negativos insignificantes cada um. Veja todos os votos negativos sobre a questão! Até agora, a questão recebeu 16 votos negativos. Ele também recebeu 80 votos positivos que mais do que compensaram todos os votos negativos.
David Hammen,
12
  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)
Vlad Feinstein
fonte
2
Isso é perceptível, hein? :) FORTRAN, IIRC, foi projetado para cálculos de fórmulas eficientes ("FOR" é para fórmulas).
Vlad Feinstein,
A maioria dos códigos F77 que vi se parecia com isso (por exemplo, BLAS e NR). Que bom que Fortran 90-> 2008 existe :)
Kyle Kanos
Sim. Se você está traduzindo uma fórmula, que maneira melhor do que FORmulaTRANslation?
Brian Drummond,
1
Sua 'otimização' ataca o lugar errado. Os bits caros são as chamadas para std::pow(), das quais você ainda tem 6, 3 vezes mais do que o necessário. Em outras palavras, seu código é 3 vezes mais lento do que possível.
Walter,
7

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.

neocpp
fonte
A forma de Horner é bastante padrão para codificar polinômios e não tem nenhuma relevância para a questão.
Walter
1
Isso pode ser verdade dado seu exemplo, mas você notará que ele disse "equações deste tipo". Achei que a resposta seria útil se o pôster tivesse polinômios em seu sistema. Eu particularmente notei que geradores de código para programas CAS, como Mathematica e Maple, tendem a NÃO dar a forma de Horner a menos que você especificamente peça por ela; eles são padronizados para a maneira como você normalmente escreveria um polinômio como um humano.
neocpp,
3

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.

NathanOliver
fonte
6
Aposto que o otimizador já faz isso para você ... embora torne pelo menos o código mais legível.
Karoly Horvath
Eu fiz isso, mas não acelerou nada. Achei que isso acontecia porque a otimização do compilador já estava cuidando disso.
armazenar l1 * l2 * l3 acelera as coisas, não sei por que com a otimização do compilador
porque o compilador às vezes simplesmente não consegue fazer algumas otimizações ou encontrá-las em conflito com outras opções.
Javier,
1
Na verdade, o compilador não deve fazer essas otimizações a menos que -ffast-mathesteja habilitado e, conforme observado em um comentário de @ tpg2114, essa otimização pode criar resultados extremamente instáveis.
David Hammen,
0

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.

user3791372
fonte
10
Essa resposta é ortogonal à pergunta em questão. Embora as GPUs tenham muitos e muitos processadores, elas são bem lentas em comparação com a FPU incorporada à CPU. Realizar um único cálculo serial com uma GPU é uma grande perda. A CPU tem que preencher o pipeline até a GPU, esperar que a GPU lenta execute aquela única tarefa e, em seguida, descarregar o resultado. Embora as GPUs sejam absolutamente fantásticas quando o problema em questão é massivamente paralelizável, elas são absolutamente atrozes quando se trata de executar tarefas seriais.
David Hammen,
1
Na pergunta original: "Como este código é executado muitas vezes, o desempenho é uma preocupação.". Isso é mais um do que "muitos". O op pode enviar os cálculos de uma maneira encadeada.
user3791372
0

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.

Fred Mitchell
fonte
0

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!

matemática
fonte