c++ performance optimization floating-point g++

¿Cómo puedo mejorar el rendimiento a través de un enfoque de alto nivel al implementar ecuaciones largas en C++?



performance optimization (10)

Editar resumen

  • Mi respuesta original simplemente señaló que el código contenía muchos cálculos replicados y que muchas de las potencias involucraban factores de 1/3. Por ejemplo, pow(x, 0.1e1/0.3e1) es lo mismo que cbrt(x) .
  • Mi segunda edición fue simplemente incorrecta, y mi tercera fue extrapolada sobre esta equivocación. Esto es lo que hace que las personas tengan miedo de cambiar los resultados similares a los oráculos de los programas matemáticos simbólicos que comienzan con la letra ''M''. Eliminé (es decir, eliminé ) esas ediciones y las empujé al final de la revisión actual de esta respuesta. Sin embargo, no los eliminé. Soy humano. Es fácil para nosotros cometer un error.
  • Mi cuarta edición desarrolló una expresión muy compacta que representa correctamente la expresión enrevesada en la pregunta SI los parámetros l1 , l2 y l3 son números reales positivos y si a es un número real distinto de cero. (Todavía no hemos tenido noticias del OP sobre la naturaleza específica de estos coeficientes. Dada la naturaleza del problema, estos son supuestos razonables).
  • Esta edición intenta responder al problema genérico de cómo simplificar estas expresiones.

Lo primero es lo primero

Yo uso Maple para generar el código C ++ para evitar errores.

Maple y Mathematica a veces pierden lo obvio. Aún más importante, los usuarios de Maple y Mathematica a veces cometen errores. Sustituir "a menudo", o tal vez incluso "casi siempre", en lugar de "a veces es probablemente más cercano a la marca".

Podrías haber ayudado a Maple a simplificar esa expresión contándole los parámetros en cuestión. En el ejemplo en cuestión, sospecho que l1 , l2 y l3 son números reales positivos y que a es un número real distinto de cero. Si ese es el caso, dilo. Esos programas matemáticos simbólicos generalmente asumen que las cantidades disponibles son complejas. Restringir el dominio le permite al programa hacer suposiciones que no son válidas en los números complejos.


Cómo simplificar esos grandes problemas de los programas matemáticos simbólicos (esta edición)

Los programas matemáticos simbólicos generalmente brindan la capacidad de proporcionar información sobre los diversos parámetros. Use esa habilidad, particularmente si su problema involucra división o exponenciación. En el ejemplo en cuestión, podría haber ayudado a Maple a simplificar esa expresión diciéndole que l1 , l2 y l3 son números reales positivos y que a es un número real distinto de cero. Si ese es el caso, dilo. Esos programas matemáticos simbólicos generalmente asumen que las cantidades disponibles son complejas. Restringir el dominio permite que el programa haga suposiciones como a x b x = (ab) x . Esto es solo si b son números reales positivos y si x es real. No es válido en los números complejos.

En última instancia, esos programas matemáticos simbólicos siguen algoritmos. Ayúdalo a lo largo. Intente jugar con expandir, recopilar y simplificar antes de generar código. En este caso, podría haber recopilado aquellos términos que involucran un factor de mu y aquellos que involucran un factor de K Reducir una expresión a su "forma más simple" sigue siendo un poco un arte.

Cuando tenga un feo desorden de código generado, no lo acepte como una verdad que no debe tocar. Intenta simplificarlo tú mismo. Mire lo que tenía el programa matemático simbólico antes de generar código. Mira cómo reduje tu expresión a algo mucho más simple y mucho más rápido, y cómo la respuesta de Walter llevó la mía varios pasos más allá. No hay receta mágica. Si hubiera una receta mágica, Maple la habría aplicado y le habría dado la respuesta que dio Walter.


Sobre la pregunta específica

Estás haciendo mucha suma y resta en ese cálculo. Puede tener serios problemas si tiene términos que casi se cancelan entre sí. Está desperdiciando mucha CPU si tiene un término que domina sobre los demás.

A continuación, está desperdiciando mucha CPU al realizar cálculos repetidos. A menos que haya habilitado -ffast-math , que le permite al compilador romper algunas de las reglas del punto flotante IEEE, el compilador no (de hecho, no debe) simplificar esa expresión para usted. En su lugar, hará exactamente lo que le dijo que hiciera. Como mínimo, debe calcular l1 * l2 * l3 antes de calcular ese desorden.

Finalmente, estás haciendo muchas llamadas a pow , lo cual es extremadamente lento. Tenga en cuenta que varias de esas llamadas son de la forma (l1 * l2 * l3) (1/3) . Muchas de esas llamadas a pow podrían realizarse con una sola llamada a std::cbrt :

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

Con este,

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


Maple extrañaba lo obvio.
Por ejemplo, hay una manera mucho más fácil de escribir.

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

Suponiendo que l1 , l2 y l3 son números reales en lugar de números complejos, y que se debe extraer la raíz del cubo real (en lugar de la raíz compleja principal), lo anterior se reduce a

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

o

2.0/(3.0 * l123_pow_1_3)

Usando cbrt_l123 lugar de l123_pow_1_3 , la expresión desagradable en la pregunta se reduce 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);

Siempre verifique dos veces, pero siempre simplifique también.


Estos son algunos de mis pasos para llegar a lo anterior:

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


Respuesta incorrecta, intencionalmente guardada para humildad

Tenga en cuenta que esto está afectado. Está incorrecto.

Actualizar

Maple extrañaba lo obvio. Por ejemplo, hay una manera mucho más fácil de escribir.

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

Suponiendo que l1 , l2 e l3 son números reales en lugar de números complejos, y que se debe extraer la raíz del cubo real (en lugar de la raíz compleja principal), lo anterior se reduce a cero. Este cálculo de cero se repite muchas veces.

Segunda actualización

Si he hecho bien las matemáticas (no hay garantía de que haya hecho bien las matemáticas), la expresión desagradable en la pregunta se reduce 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);

Lo anterior supone que l1 , l2 y l3 son números reales positivos.

Estoy desarrollando algunas simulaciones de ingeniería. Esto implica implementar algunas ecuaciones largas como esta ecuación para calcular la tensión en un material similar al caucho:

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;

Utilizo Maple para generar el código C ++ para evitar errores (y ahorrar tiempo con álgebra tediosa). Como este código se ejecuta miles (si no millones) de veces, el rendimiento es preocupante. Lamentablemente, las matemáticas solo se simplifican hasta ahora; Las ecuaciones largas son inevitables.

¿Qué enfoque puedo tomar para optimizar esta implementación? Estoy buscando estrategias de alto nivel que debería aplicar al implementar tales ecuaciones, no necesariamente optimizaciones específicas para el ejemplo que se muestra arriba.

Estoy compilando usando g ++ con --enable-optimize=-O3 .

Actualizar:

Sé que hay muchas expresiones repetidas, estoy usando el supuesto de que el compilador las manejaría; mis pruebas hasta ahora sugieren que sí.

l1, l2, l3, mu, a, K son todos números reales positivos (no cero).

He reemplazado l1*l2*l3 con una variable equivalente: J Esto ayudó a mejorar el rendimiento.

Reemplazar pow(x, 0.1e1/0.3e1) con cbrt(x) fue una buena sugerencia.

Esto se ejecutará en las CPU. En el futuro cercano, esto probablemente funcionaría mejor en las GPU, pero por ahora esa opción no está disponible.

  1. ¿Cuántos son "muchos muchos"?
  2. ¿Cuánto tiempo se tarda?
  3. ¿Cambian TODOS los parámetros entre el recálculo de esta fórmula? ¿O puede almacenar en caché algunos valores calculados previamente?

  4. He intentado una simplificación manual de esa fórmula, ¿me gustaría saber si guarda algo?

    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;

[AGREGADO] He trabajado un poco más en la última fórmula de tres líneas y me puse a esta belleza:

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

Déjame mostrarte mi trabajo, paso a paso:

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)

Como ha preguntado explícitamente sobre optimizaciones de alto nivel, podría valer la pena probar diferentes compiladores de C ++. Hoy en día, los compiladores son bestias de optimización muy complejas y los proveedores de CPU pueden implementar optimizaciones muy potentes y específicas. Pero tenga en cuenta que algunos de ellos no son gratuitos (pero puede haber un programa académico gratuito).

  • La colección del compilador GNU es gratuita, flexible y está disponible en muchas arquitecturas.
  • Los compiladores Intel son muy rápidos, muy caros y también pueden producir buenos resultados para las arquitecturas AMD (creo que hay un programa académico)
  • Los compiladores de Clang son rápidos, gratuitos y pueden producir resultados similares a los de GCC (algunas personas dicen que son más rápidos, mejor, pero esto puede diferir para cada caso de aplicación, sugiero que haga sus propias experiencias)
  • PGI (Portland Group) no es gratuito como los compiladores de Intel.
  • Los compiladores PathScale pueden obtener buenos resultados en arquitecturas AMD

He visto que los fragmentos de código difieren en la velocidad de ejecución en un factor de 2, solo cambiando el compilador (con optimizaciones completas, por supuesto). Pero tenga cuidado de verificar la identidad de la salida. Una optimización agresiva puede conducir a resultados diferentes, que es algo que definitivamente desea evitar.

¡Buena suerte!

Esto puede ser un poco breve, pero en realidad he encontrado una buena aceleración para los polinomios (interpolación de funciones de energía) usando la Forma de Horner, que básicamente reescribe ax^3 + bx^2 + cx + d como d + x(c + x(b + x(a))) . Esto evitará muchas llamadas repetidas a pow() y le impedirá hacer cosas tontas como llamar por separado a pow(x,6) y pow(x,7) lugar de simplemente hacer x*pow(x,6) .

Esto no es directamente aplicable a su problema actual, pero si tiene polinomios de alto orden con potencias enteras, puede ayudar. Es posible que tenga que tener cuidado con los problemas de estabilidad numérica y desbordamiento, ya que el orden de las operaciones es importante para eso (aunque en general creo que Horner Form ayuda para esto, ya que x^20 x suelen estar separados por muchos órdenes de magnitud).

También como un consejo práctico, si aún no lo ha hecho, primero intente simplificar la expresión en arce. Probablemente pueda lograr que haga la mayor parte de la eliminación de subexpresión común por usted. No sé cuánto afecta al generador de código en ese programa en particular, pero sé que en Mathematica hacer un FullSimplify antes de generar el código puede resultar en una gran diferencia.

Lo primero que debe tener en cuenta es que el pow es realmente costoso, por lo que debe deshacerse de esto tanto como sea posible. Al escanear la expresión, veo muchas repeticiones de pow(l1 * l2 * l3, -0.1e1 / 0.3e1) y pow(l1 * l2 * l3, -0.4e1 / 0.3e1) . Por lo tanto, esperaría una gran ganancia al precalcular esos:

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

donde estoy usando la función boost pow .

Además, tienes más pow con el exponente a . Si a es Integer y se conoce en el momento del compilador, también puede reemplazarlos con boost::math::pow<a>(...) para obtener un mayor rendimiento. También sugeriría reemplazar términos como a / l1 / 0.3e1 con a / (l1 * 0.3e1) ya que la multiplicación es más rápida que la división.

Finalmente, si usa g ++ puede usar el -ffast-math que permite que el optimizador sea más agresivo en la transformación de ecuaciones. Lea acerca de lo que realmente hace esta bandera , ya que tiene efectos secundarios.

Parece que tienes muchas operaciones repetidas en curso.

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

Puede precalcularlos para no llamar repetidamente a la función pow , lo que puede ser costoso.

También podrías precalutar

l1 * l2 * l3

mientras usa ese término repetidamente.

Por casualidad, ¿podría proporcionar el cálculo simbólicamente? Si hay operaciones vectoriales, es posible que desee investigar utilizando blas o lapack, que en algunos casos puede ejecutar operaciones en paralelo.

Es concebible (¿con el riesgo de estar fuera de tema?) Que pueda usar Python con numpy y / o scipy. En la medida de lo posible, sus cálculos podrían ser más legibles.

Si tiene una tarjeta gráfica Nvidia CUDA, podría considerar descargar los cálculos en la tarjeta gráfica, que es más adecuada para cálculos computacionalmente complicados.

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

Si no, es posible que desee considerar múltiples hilos para los cálculos.

Woah, qué expresión tan infernal. Crear la expresión con Maple en realidad fue una elección subóptima aquí. El resultado es simplemente ilegible.

  1. eligió hablar nombres de variables (no l1, l2, l3, sino, por ejemplo, altura, anchura, profundidad, si eso es lo que significan). Entonces es más fácil para usted comprender su propio código.
  2. calcule subterms, que usa varias veces, por adelantado y almacene los resultados en variables con nombres de voz.
  3. Usted menciona que la expresión se evalúa muchas veces. Supongo que solo unos pocos parámetros varían en el bucle más interno. Calcule todos los subterráneos invariantes antes de ese ciclo. Repita para el segundo ciclo interno y así sucesivamente hasta que todos los invariantes estén fuera del ciclo.

Teóricamente, el compilador debería poder hacer todo eso por usted, pero a veces no puede hacerlo, por ejemplo, cuando la anidación del bucle se extiende sobre múltiples funciones en diferentes unidades de compilación. De todos modos, eso le dará un código mucho mejor legible, comprensible y mantenible.

La respuesta de David Hammen es buena, pero aún lejos de ser óptima. Continuemos con su última expresión (al momento de escribir esto)

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 se puede optimizar aún más. En particular, podemos evitar la llamada a cbrt() y una de las llamadas a pow() si explotamos algunas identidades matemáticas. Hagamos esto nuevamente paso a paso.

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

Tenga en cuenta que también he optimizado 2.0*N1 a N1+N1 etc. A continuación, solo podemos hacer dos llamadas a 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);

Dado que las llamadas a pow() son, con mucho, la operación más costosa aquí, vale la pena reducirlas lo más posible (la siguiente operación costosa fue la llamada a cbrt() , que eliminamos).

Si por casualidad a es entero, las llamadas a pow podrían optimizarse para llamadas a cbrt (más potencias enteras), o si un athird es medio entero, podemos usar sqrt (más potencias enteras). Además, si por casualidad l1==l2 o l1==l3 o l2==l3 se pueden eliminar una o ambas llamadas a pow . Por lo tanto, vale la pena considerar estos casos especiales si tales posibilidades existen de manera realista.