在c++中实现长方程时,如何通过高级方法提高性能?< / h1 >
How can I improve performance via a high-level approach when implementing long equations in 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
均为正实数(非零)。
我已将l1*l2*l3
替换为等效变量:J
。这确实有助于提高性能。
用cbrt(x)
代替pow(x, 0.1e1/0.3e1)
是一个很好的建议。
这将在cpu上运行,在不久的将来,这可能会在gpu上运行得更好,但现在这个选项不可用。
编辑摘要
- 我最初的回答只是指出代码包含了大量的重复计算,并且许多幂涉及1/3的因子。例如,
pow(x, 0.1e1/0.3e1)
与cbrt(x)
相同。 - 我的第二次编辑是错误的,我的第三次编辑是在这个错误的基础上推断出来的。这就是为什么人们害怕改变以字母"M"开头的符号数学程序的类似神谕的结果。我已经删除了(例如,
strike)这些编辑,并将它们推到当前版本的底部。但是,我没有删除它们。我是人类。我们很容易犯错。 - 我的第四次编辑开发了一个非常紧凑的表达式,它正确地表示了问题IF中的复杂表达式
l1
、l2
、l3
为正实数,如果a
为非零实数。(我们还没有从OP那里听到关于这些系数的具体性质。鉴于问题的性质,这些都是合理的假设。 - 这个编辑试图回答如何简化这些表达式的一般性问题。
重要的事情先做
我使用Maple来生成c++代码,以避免错误。
Maple和Mathematica有时会忽略一些显而易见的东西。更重要的是,Maple和Mathematica的用户有时会犯错误。用"often ",甚至"almost always"来代替"sometimes "可能更贴切。
你可以通过告诉Maple有问题的参数来帮助它简化表达式。在本例中,我怀疑l1
、l2
和l3
是正实数,a
是非零实数。如果是这样,就这么说。这些符号数学程序通常假设手头的数量是复杂的。限制域允许程序做出在复数中无效的假设。
如何简化符号数学程序中的大混乱(此编辑)
符号数学程序通常提供了提供各种参数信息的能力。使用这种能力,特别是当你的问题涉及除法或求幂时。在手边的示例中,您可以通过告诉Maplel1
、l2
和l3
是正实数,并且a
是非零实数来帮助Maple简化表达式。如果是这样,就这么说。这些符号数学程序通常假设手头的数量是复杂的。限制域允许程序做出假设,例如axbx=(ab)x。只有当a
和b
是正实数并且x
是实数时,才会出现这种情况。在复数中无效
最终,这些符号数学程序遵循算法。帮助它前进。在生成代码之前,尝试展开、收集和简化。在本例中,您可以收集那些涉及因子mu
和涉及因子K
的项。将一个表达简化为"最简单的形式"仍然是一门艺术。
当你得到一堆乱七八糟的生成代码时,不要把它当作你不能碰的事实。试着自己简化它。看看符号数学程序在生成代码之前有什么。看看我是如何把你的表达简化得更简单、更快的,而沃尔特的回答又使我的回答向前推进了几步。没有什么神奇的配方。如果有魔法配方的话,枫树一定会应用它,并给出沃尔特给出的答案。
关于具体问题
你在计算中做了很多加减法。如果你的项几乎相互抵消,你就会陷入大麻烦。如果你有一个项占主导地位,你会浪费大量的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
Maple确实错过了明显的。
例如,有一种更简单的方法来写
(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);
错误答案,故意保留,以示谦虚
注意这是错误的。这是错误的。
<罢工>更新罢工>
Maple确实错过了显而易见的事情。例如,有一种更简单的方法来写
(战俘(l1 * l2 * l3, -0.1 e1/0.3 e1) l1 - l2 * * l3 *战俘(l1 * l2 * l3, -0.4 e1/0.3 e1)/0.3 e1)
假设l1
、l2
、l3
为实数而非复数,提取实数立方根(而非主复方根),则上述约为零。这个零的计算重复多次。
第二更新
如果我没算错的话(有没有)我保证我的计算是正确的),问题中令人讨厌的表达就变成了
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);
我使用的是boost函数
此外,您有更多的pow
指数a
。如果a
是整数,并且在编译时已知,您还可以用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);
请注意,我还将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()
的调用,我们取消了它)。
如果碰巧a
是整数,对pow
的调用可以优化为对cbrt
的调用(加上整数幂),或者如果athird
是半整数,我们可以使用sqrt
(加上整数幂)。此外,如果碰巧l1==l2
或l1==l3
或l2==l3
可以消除对pow
的一个或两个调用。因此,如果这种可能性确实存在,那么值得将这些情况视为特殊情况。
- "many many"是多少?
- 要多久?
- 重新计算此公式时所有参数是否改变?或者可以缓存一些预先计算的值吗?
-
我已经尝试手动简化该公式,想知道它是否节省了什么?
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;
[add]我对最后三行公式做了更多的工作,并得到了这样的美丽:
T = X9 / X0 * (
(X4 * XD - XF - XG) * N1 +
(X5 * XD - XE - XG) * N2 +
(X5 * XD - XE - XF) * N3)
+ XB * (N1 + N2 + N3)
让我一步一步地展示我的作品:
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)
这可能有点简洁,但我实际上通过使用霍纳形式发现了多项式(能量函数的插值)的良好加速,这基本上将ax^3 + bx^2 + cx + d
重写为d + x(c + x(b + x(a)))
。这将避免对pow()
的大量重复调用,并阻止你做一些愚蠢的事情,比如分别调用pow(x,6)
和pow(x,7)
,而不是只调用x*pow(x,6)
。
这并不直接适用于你当前的问题,但如果你有整数次的高阶多项式,它会有所帮助。你可能不得不注意数值稳定性和溢出问题,因为操作顺序很重要(尽管通常我认为霍纳形式有助于解决这个问题,因为x^20
和x
通常相差很多数量级)。
作为一个实用提示,如果您还没有这样做,请尝试先简化maple中的表达式。您可能可以让它为您做大多数常见的子表达式消除。我不知道它对那个程序中的代码生成器有多大的影响,但我知道在Mathematica中,在生成代码之前做一个FullSimplify可以产生巨大的差异。
看起来你有很多重复的操作。
pow(l1 * l2 * l3, -0.1e1 / 0.3e1)
pow(l1 * l2 * l3, -0.4e1 / 0.3e1)
你可以预先计算这些,这样你就不会重复调用pow
函数,这可能会很昂贵。
你也可以预先计算
l1 * l2 * l3
,因为你反复使用这个词。
如果您有Nvidia CUDA显卡,您可以考虑将计算卸载到显卡上-显卡本身更适合计算复杂的计算。
https://developer.nvidia.com/how-to-cuda-c-cpp
如果不是,您可能需要考虑多个线程进行计算。
你能不能用符号表示一下?如果有向量操作,你可能真的想研究使用blas或lapack,它们在某些情况下可以并行运行操作。
可以想象(有偏离主题的风险?)您可能能够将python与numpy和/或scipy一起使用。在某种程度上,这是可能的,你的计算可能更可读。
正如您明确询问的高级优化,可能值得尝试不同的c++编译器。如今,编译器是非常复杂的优化野兽,CPU供应商可能会实现非常强大和特定的优化。但请注意,其中一些不是免费的(但可能有一个免费的学术程序)。
- GNU编译器集合是免费的,灵活的,并且可以在许多体系结构上使用
- 英特尔编译器非常快,非常昂贵,并且也可能为AMD架构产生良好的结果(我相信有一个学术计划)
- Clang编译器快速,免费,并且可能产生与GCC相似的结果(有些人说它们更快,更好,但这可能因每个应用案例而异,我建议您自己体验)
- PGI (Portland Group)不像Intel编译器那样是免费的。
- PathScale编译器可能在AMD架构上执行良好的结果
我看到过代码片段的执行速度相差2倍,这仅仅是因为改变了编译器(当然是完全优化的)。但是要注意检查输出的身份。过于激进的优化可能会导致不同的输出,这是您绝对想要避免的。
祝你好运!
- EASTL矢量<向量<int>>连续的
- 从 16UC3 到 8UC3 的高性能 OpenCV 矩阵转换
- 如何从高性能的输入迭代器返回变体?
- 编写高性能C++二传手
- 提升图形库:以高性能的方式检查vertex_descriptor的有效性
- C - 创建矢量&lt; vector&lt; double&gt;&gt;矩阵具有分配而不是inizializ
- C 字符串比较“祝您好运”&gt;“再见”
- 为什么将此对向量&lt; map&lt; int,int&gt;&gt;中的地图进行更新.失败
- 高性能程序,什么是更好的矢量数组或矢量的矢量
- 如何在不使用函数或类的情况下重复代码段,以便在C++中实现高性能循环
- C :对矢量进行排序&lt; struct&gt;(结构有2个整数)基于结构的整数之一
- C 操作员&gt;&gt;与突变器过载
- 明确的专业化“ CheckIntmap&lt;&gt;”实例化
- 是否需要使用 - &gt;运算符在C 中调用成员函数时
- 什么是模板&lt;&gt;inline bla bla
- 编辑C Qlist&lt; object*&gt; gt;QML代码和一些QML警告中的模型
- eigen :: llt&lt;eigen :: matrixxd&gt;具有不完整的类型
- 错误,包括&lt; ctype&gt;在原子上使用C 11
- 标准::矢量的高性能替代品
- 在Qt中以高性能方式将(富)文本附加到QTextEdit或QTextBrowser中