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;
저는 Maple을 사용하여 실수를 방지하고 지루한 대수로 시간을 절약하기 위해 C ++ 코드를 생성합니다. 이 코드는 수천 번 (수백만은 아니더라도) 실행되므로 성능이 문제가됩니다. 불행히도 수학은 지금까지만 단순화되었습니다. 긴 방정식은 피할 수 없습니다.
이 구현을 최적화하기 위해 어떤 접근 방식을 취할 수 있습니까? 위의 예에 대한 특정 최적화가 아니라 이러한 방정식을 구현할 때 적용해야하는 높은 수준의 전략을 찾고 있습니다.
나는 g ++를 사용하여 --enable-optimize=-O3.
최신 정보:
반복되는 표현이 많다는 것을 알고 있으며 컴파일러가이를 처리 할 것이라는 가정을 사용하고 있습니다. 지금까지 내 테스트에 따르면 그렇습니다.
l1, l2, l3, mu, a, K 모두 양의 실수 (0이 아님)입니다.
l1*l2*l3동등한 변수 로 대체 했습니다 : J. 이것은 성능 향상에 도움이되었습니다.
교체 pow(x, 0.1e1/0.3e1)로하는 것은 cbrt(x)좋은 제안했다.
이것은 CPU에서 실행될 것입니다. 가까운 장래에 GPU에서 더 잘 실행될 것 같지만 현재로서는 해당 옵션을 사용할 수 없습니다.
요약 수정
- 내 원래 대답은 코드에 많은 복제 계산이 포함되어 있고 많은 힘이 1/3의 요소와 관련이 있다는 점에 주목했습니다. 예를 들어
pow(x, 0.1e1/0.3e1)는cbrt(x). - 두 번째 편집은 틀렸고 세 번째 편집은이 잘못을 외삽했습니다. 이것이 사람들이 문자 'M'으로 시작하는 상징적 수학 프로그램에서 오라클과 같은 결과를 변경하는 것을 두려워하게 만드는 것입니다. 나는 그 편집을
스트라이크(즉,스트라이크) 하고이 답변의 현재 개정판의 맨 아래로 밀어 넣었습니다. 그러나 나는 그들을 삭제하지 않았다. 나는 인간이다. 우리는 실수하기 쉽습니다. - 나의 네 번째 편집 올바르게 문제의 복잡한 표현을 나타내는 매우 컴팩트 한 표정으로 개발 IF 매개 변수가
l1,l2그리고l3긍정적 인 실제 숫자와 경우에a비 제로 실수입니다. (이러한 계수의 특정 특성에 대해서는 아직 OP로부터 듣지 못했습니다. 문제의 특성을 고려할 때 이는 합리적인 가정입니다.) - 이 편집은 이러한 표현을 단순화하는 방법에 대한 일반적인 문제에 대한 답변을 시도합니다.
먼저 첫 번째 것들
실수를 피하기 위해 Maple을 사용하여 C ++ 코드를 생성합니다.
Maple과 Mathematica는 때때로 명백한 것을 놓칩니다. 더욱 중요한 것은 Maple 및 Mathematica 사용자가 때때로 실수를한다는 것입니다. "가끔"또는 "거의 항상"을 "때때로"로 대체하는 것이 아마도 마크에 더 가깝습니다.
Maple이 문제의 매개 변수에 대해 말함으로써 해당 표현을 단순화하는 데 도움이 될 수 있습니다. 손의 예에서, 그 의심 l1, l2그리고 l3있습니다 긍정적 인 실수를하고는 a비 제로 실수입니다. 그럴 경우 그렇게 말하십시오. 이러한 상징적 수학 프로그램은 일반적으로 당면한 양이 복잡하다고 가정합니다. 도메인을 제한하면 프로그램이 복소수에서 유효하지 않은 가정을 할 수 있습니다.
기호 수학 프로그램에서 이러한 큰 혼란을 단순화하는 방법 (이 편집)
기호 수학 프로그램은 일반적으로 다양한 매개 변수에 대한 정보를 제공하는 기능을 제공합니다. 특히 문제가 나눗셈이나 지수화와 관련된 경우 그 능력을 사용하십시오. 손의 예에서, 당신은 단풍 나무가 그것을 말함으로써 그 표현을 단순화 도움이 수 l1, l2및 l3포지티브 실수를하고는 a비 제로 실수입니다. 그럴 경우 그렇게 말하십시오. 이러한 상징적 수학 프로그램은 일반적으로 당면한 양이 복잡하다고 가정합니다. 도메인을 제한하면 프로그램이 a x b x = (ab) x 와 같은 가정을 할 수 있습니다 . 이 경우에만입니다 a및 b긍정적 인 실수하고있는 경우는 x진짜입니다. 복소수에서는 유효하지 않습니다.
궁극적으로 이러한 상징적 수학 프로그램은 알고리즘을 따릅니다. 함께 도와주세요. 코드를 생성하기 전에 확장, 수집 및 단순화를 시도해보십시오. 이 경우, 당신은 그의 요소를 포함하는 용어 수집 한 수 mu와의 요소를 포함하는 사람들을 K. 표현을 "가장 단순한 형태"로 줄이는 것은 약간의 예술로 남아 있습니다.
생성 된 코드가 엉망이 될 때 만져서는 안된다는 사실로 받아들이지 마십시오. 직접 단순화하십시오. 기호 수학 프로그램이 코드를 생성하기 전에 무엇을 가지고 있는지 살펴보십시오. 내가 어떻게 당신의 표현을 훨씬 더 간단하고 훨씬 빠르게 줄 였는지, 그리고 Walter의 대답 이 어떻게 내 여러 단계를 더 발전 시켰 는지 보세요. 마법의 레시피는 없습니다. 마법의 레시피가 있다면 메이플은 그것을 적용하고 월터가 준 대답을했을 것입니다.
구체적인 질문에 대해
당신은 그 계산에서 많은 덧셈과 뺄셈을하고 있습니다. 서로 거의 취소되는 조건이 있으면 심각한 문제에 빠질 수 있습니다. 하나의 용어가 다른 용어보다 우세한 경우 많은 CPU를 낭비하고 있습니다.
다음으로 반복 계산을 수행하여 많은 CPU를 낭비하고 있습니다. -ffast-math컴파일러가 IEEE 부동 소수점의 일부 규칙을 위반할 수 있도록 하는 을 활성화하지 않는 한 컴파일러는 해당 표현식을 단순화하지 않습니다 (사실 로선 안됩니다). 대신 당신이 말한대로 정확하게 할 것입니다. 최소한 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);
오답, 의도적으로 겸손하게 유지
이것은 위험하다는 점에 유의하십시오. 틀렸어.
최신 정보
메이플은 명백한 것을 놓쳤습니다. 예를 들어 훨씬 더 쉽게 작성할 수있는 방법이 있습니다.
(pow (l1 * l2 * l3, -0.1e1 / 0.3e1)-l1 * l2 * l3 * pow (l1 * l2 * l3, -0.4e1 / 0.3e1) / 0.3e1)
그 가정 l1, l2그리고 l3복잡한 숫자보다 실제가 아니라이며, 실제 큐브 루트 (보다는 원칙 복잡한 루트)입니다 추출 할 것을, 위에서 제로로 줄일 수 있습니다. 이 0의 계산은 여러 번 반복됩니다.
두 번째 업데이트
내가 수학을 제대로했다면 (내가 수학을 제대로 했다는 보장 은 없다 ), 질문의 불쾌한 표현은 다음과 같이 줄어든다.
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);
위는 가정
l1,
l2그리고
l3긍정적 인 실제 숫자입니다.
가장 먼저 주목할 점 pow은 정말 비싸기 때문에 가능한 한 많이 제거해야한다는 것입니다. 표현을 훑어 보면 pow(l1 * l2 * l3, -0.1e1 / 0.3e1)과의 반복이 많이 보입니다 pow(l1 * l2 * l3, -0.4e1 / 0.3e1). 그래서 나는 그것들을 미리 계산함으로써 큰 이득을 기대할 것입니다.
const double c1 = pow(l1 * l2 * l3, -0.1e1 / 0.3e1);
const double c2 = boost::math::pow<4>(c1);
부스트 파워 기능을 사용하고 있습니다.
또한 powexponent를 사용 하면 더 많은 것을 얻을 수 있습니다 a. aInteger이고 컴파일러 시간에 알려진 경우 boost::math::pow<a>(...)추가 성능을 얻기 위해 이를로 바꿀 수도 있습니다 . 곱셈이 나눗셈보다 빠르기 때문에 a / l1 / 0.3e1와 a / (l1 * 0.3e1)같은 용어를 바꾸는 것이 좋습니다 .
마지막으로, g ++를 사용하는 -ffast-math경우 옵티마이 저가 방정식을보다 적극적으로 변환 할 수 있도록하는 플래그를 사용할 수 있습니다 . 이 플래그 는 부작용이 있기 때문에 실제로 무엇을하는지 읽어보십시오 .
와, 정말 대단한 표정 이군요. 실제로 Maple로 표현을 만드는 것은 여기서 차선책이었습니다. 결과는 단순히 읽을 수 없습니다.
- 말하는 변수 이름을 선택했습니다 (l1, l2, l3이 아니라 높이, 너비, 깊이, 그것이 의미하는 경우). 그러면 자신의 코드를 더 쉽게 이해할 수 있습니다.
- 여러 번 사용하는 하위 용어를 계산하고 결과를 말하는 이름이있는 변수에 저장합니다.
- 당신은 표현이 매우 여러 번 평가된다는 것을 언급합니다. 대부분의 내부 루프에서는 매개 변수가 거의 변하지 않는 것 같습니다. 해당 루프 이전의 모든 불변 하위 용어를 계산하십시오. 모든 불변이 루프 외부에있을 때까지 두 번째 내부 루프에 대해 반복합니다.
이론적으로 컴파일러는 모든 작업을 수행 할 수 있어야하지만 때로는 불가능합니다. 예를 들어 루프 중첩이 다른 컴파일 단위의 여러 함수에 걸쳐 퍼져있는 경우입니다. 어쨌든, 그것은 당신에게 훨씬 더 읽기 쉽고, 이해하기 쉬우 며 유지 가능한 코드를 줄 것입니다.
David Hammen의 대답 은 좋지만 여전히 최적과는 거리가 멀습니다. 그의 마지막 표현을 계속합시다 (이 글을 쓰는 시점에서)
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);
더 최적화 할 수 있습니다. 특히, 일부 수학적 정체성을 악용 cbrt()하는 pow()경우 호출 및 호출 중 하나 를 피할 수 있습니다 . 이 작업을 단계별로 다시 해보겠습니다.
// 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);
또한 etc로 최적화 2.0*N1했습니다 N1+N1. 다음으로 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);
에 대한 호출 pow()은 여기에서 가장 비용이 많이 드는 작업이므로 가능한 한 많이 줄이는 것이 cbrt()좋습니다 (다음으로 비용이 많이 드는 작업은에 대한 호출 이며 제거했습니다).
If by any chance a is integer, the calls to pow could be optimized to calls to cbrt (plus integer powers), or if athird is half-integer, we can use sqrt (plus integer powers). Furthermore, if by any chance l1==l2 or l1==l3 or l2==l3 one or both calls to pow can be eliminated. So, it's worth to consider these as special cases if such chances realistically exist.
- How many is "many many"?
- How long does it take?
- Do ALL parameters change between recalculation of this formula? Or can you cache some pre-calculated values?
I've attempted a manual simplification of that formula, would like to know if it saves anything?
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;
[ADDED] I've worked some more on the last three-lines formula and got it down to this beauty:
T = X9 / X0 * (
(X4 * XD - XF - XG) * N1 +
(X5 * XD - XE - XG) * N2 +
(X5 * XD - XE - XF) * N3)
+ XB * (N1 + N2 + N3)
Let me show my work, step by step:
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)
This may be a little terse, but I've actually found good speedup for polynomials (interpolation of energy functions) by using Horner Form, which basically rewrites ax^3 + bx^2 + cx + d as d + x(c + x(b + x(a))). This will avoid a lot of repeated calls to pow() and stops you from doing silly things like separately calling pow(x,6) and pow(x,7) instead of just doing x*pow(x,6).
This is not directly applicable to your current problem, but if you have high order polynomials with integer powers it can help. You might have to watch out for numerical stability and overflow issues since the order of operations is important for that (although in general I actually think Horner Form helps for this, since x^20 and x are usually many orders of magnitude apart).
Also as a practical tip, if you haven't done so already, try to simplify the expression in maple first. You can probably get it to do most of the common subexpression elimination for you. I don't know how much it affects the code generator in that program in particular, but I know in Mathematica doing a FullSimplify before generating the code can result in a huge difference.
It looks like you have a lot of repeated operations going on.
pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
pow(l1 * l2 * l3, -0.4e1 / 0.3e1)
You could pre-calculate those so you are not repeatedly calling the pow function which can be expensive.
You could also pre-calutate
l1 * l2 * l3
as you use that term repeatedly.
If you have a Nvidia CUDA graphics card, you could consider offloading the calculations to the graphics card - which itself is more suitable for computationally complicated calculations.
https://developer.nvidia.com/how-to-cuda-c-cpp
If not, you may want to consider multiple threads for calculations.
By any chance, could you supply the calculation symbolically. If there are vector operations, you might really want to investigate using blas or lapack which in some cases can run operations in parallel.
It is conceivable (at the risk of being off-topic?) that you might be able to use python with numpy and/or scipy. To the extent that it was possible, your calculations might be more readable.
As you explicitly asked about high level optimizations, it might be worth trying different C++ compilers. Nowadays, compilers are very complex optimization beasts and CPU vendors might implement very powerful and specific optimizations. But please note, some of them are not free (but there might be a free academic program).
- GNU compiler collection is free, flexible, and available on many architectures
- Intel compilers are very fast, very expensive, and may also produce good results for AMD architectures (I believe there is an academic program)
- Clang compilers are fast, free, and might produce similar results to GCC (some people say they are faster, better, but this might differ for each application case, I suggest to make your own experiences)
- PGI (Portland Group) is not free as the Intel compilers.
- PathScale compilers might perform good results on AMD architectures
I've seen code snippets differ in execution speed by the factor of 2, only by changing the compiler (with full optimizations of course). But be aware of checking the identity of the output. To aggressive optimization might lead to different output, which is something you definitely want to avoid.
Good Luck!
'programing' 카테고리의 다른 글
| Scala를 사용하는 Joda Time의 클래스 중단 오류 (0) | 2020.09.03 |
|---|---|
| JSON 유형 내부의 배열 요소 쿼리 (0) | 2020.09.03 |
| C에서 stderr로 어떻게 인쇄 할 수 있습니까? (0) | 2020.09.03 |
| 두 개의 슬래시로 시작하는 URI… 어떻게 동작합니까? (0) | 2020.09.03 |
| Android xml 오류 : RelativeLayout (@ id / LinearLayout_acc, @ id / ProgressBar_statusScreen)과 함께 "주어진 이름과 일치하는 리소스를 찾을 수 없음" (0) | 2020.09.03 |