寻找一种加速功能的方法

我首先重写了它，以便更容易理解：

#include <math.h> 
double Param_RE_Tterm_approx(double Tterm, double const* par) {
double value = 0.;
if (Tterm > 0.) {
// time after Che angle peak
if ( fabs(Tterm/ par[0]) >= 1.e-2) {
value += -1./(par[0])*exp(-1.*Tterm/(par[0]));
} else {
value += -1./par[0]*(1. - Tterm/par[0] + Tterm*Tterm/(par[0]*par[0]*2.) - Tterm*Tterm*Tterm/(par[0]*par[0]*par[0]*6.) );
}
if ( fabs(Tterm* par[1]) >= 1.e-2) {
value += par[2]* par[1]*pow( 1.+par[1]*Tterm, par[2]-1. );
} else {
value += par[2]*par[1]*( 1.+(par[2]-1.)*par[1]*Tterm + (par[2]-1.)*(par[2]-1.-1.)/2.*par[1]*par[1]*Tterm*Tterm + (par[2]-1.)*(par[2]-1.-1.)*(par[2]-1.-2.)/6.*par[1]*par[1]*par[1]*Tterm*Tterm*Tterm );
}
} else {
// time before Che angle peak
if ( fabs(Tterm/ par[3]) >= 1.e-2 ) {
value += -1./ par[3] *exp(-1.*Tterm/ par[3]);
} else {
value += -1./par[3]*(1. - Tterm/par[3] + Tterm*Tterm/(par[3]*par[3]*2.) - Tterm*Tterm*Tterm/(par[3]*par[3]*par[3]*6.) );
}
if ( fabs(Tterm* par[4]) >= 1.e-2 ) {
value += par[5]* par[4] *pow( 1.+ par[4]*Tterm, par[5]-1. );
} else {
value += par[5]*par[4]*( 1.+(par[5]-1.)*par[4]*Tterm + (par[5]-1.)*(par[5]-1.-1.)/2.*par[4]*par[4]*Tterm*Tterm + (par[5]-1.)*(par[5]-1.-1.)*(par[5]-1.-2.)/6.*par[4]*par[4]*par[4]*Tterm*Tterm*Tterm );
}
}
return value * 1.e9;
}

然后我们可以看看它的结构。

有两个主要分支 - Tterm负(之前)和正(之后)。 这些对应于在par数组中使用 0,1,2 或 3,4,5。

然后在每种情况下，我们都会做两件事来增加价值。 在这两种情况下，对于小情况，我们使用多项式，对于大情况，我们使用指数/幂方程。

作为猜测，这是因为多项式是小值的指数的体面近似 - 误差是可以接受的。 你应该做的是确认这个猜测——看看基于"大"幂/指数的方程的泰勒级数展开，看看它是否以某种方式与多项式一致。 或以数字方式检查。

如果是这种情况，这意味着该方程具有可接受的已知误差量。 通常，exp或pow的更快版本具有已知的最大误差量;考虑使用这些。

如果不是这种情况，仍然可能存在可接受的误差量，但泰勒级数近似可以为您提供有关可接受的误差量的"代码"信息。

我要采取的下一步是将这个等式的 8 个部分拆开。 有正/负，每个分支中的第一个和第二个value+=，然后是多项式/指数情况。

我认为 exp 占用 ~1/3 的 pow 时间是因为您的函数中有 3 次调用 pow 到 1 次调用 exp，但您可能会发现一些有趣的事情，例如"我们所有的时间实际上都在 Tterm> 0。案例"或你有什么。

现在检查呼叫站点。 您传递此函数的 Tterm 中是否有模式？ 即，您是否倾向于按大致排序顺序传递 Tterms？ 如果是这样，您可以在调用此函数之外测试要调用哪个函数，并分批进行。

简单地批量进行编译并使用优化和内联函数的主体进行编译可能会产生惊人的差异;编译器在矢量化工作方面越来越好。

如果这不起作用，您可以开始解决问题。 在现代计算机上，您可以有 4-60 个线程独立解决此问题，并且这个问题看起来您会获得近乎线性的加速。 一个基本的线程库，如TBB，将适用于这种任务。

下一步，如果您要获得大量数据并且需要进行大量处理，则可以将其填充到GPU上并在那里解决。 可悲的是，GPU<->RAM的通信很小，因此简单地在GPU上使用此函数进行数学运算并使用RAM来回读取/写入不会给您带来太多性能。 但是，如果GPU上可以做更多的工作，那可能是值得的。

总时间中只有 10% 被此函数中非 pow 或 exp 的部分使用。

如果您的函数性能瓶颈是 exp()、pow() 执行，请考虑在计算中使用向量指令。所有现代处理器都至少支持 SSE2 指令集，因此这种方法肯定会提供至少 ~2 倍的速度，因为您的计算可以很容易地矢量化。

我建议您使用这个 c++ 矢量化库，它包含所有标准数学函数(例如 exp 和 pow)，并允许在不使用汇编语言的情况下以 OOP 风格编写代码。我用过几次，它必须完美地解决你的问题。

如果你有GPU，你还应该考虑尝试cuda框架，因为同样，你的问题可以被完美地矢量化。此外，如果这个函数被调用477+百万次，GPU将真正消除你的问题......

(部分优化：)

最长的表达式有

• 常用子表达式
• 多项式评估了昂贵的方法。

预定义这些(也许将它们添加到 par[])：

a = par[5]*par[4];
b =   (par[5]-1.);
c = b*(par[5]-2.)/2.;
d = c*(par[5]-3.)/3.;

然后，例如，最长的表达式变为：

e = par[4]*Tterm;
value += a*(((d*e + c)*e + b)*e + 1.);

并简化其余部分。

如果表达式是曲线拟合近似，为什么不也使用

value += -1./(*par)*exp(-1.*Tterm/(*par));

您还应该询问是否需要所有 477M 迭代。

如果您想探索批处理/更多优化机会，以便在依赖于这些值的计算中进行融合，请尝试使用 Halide

我在这里用 Halide 重写了您的程序：

#include <Halide.h>
using namespace Halide;
class ParamReTtermApproxOpt : public Generator<ParamReTtermApproxOpt>
{
public:
Input<Buffer<float>> tterm{"tterm", 1};
Input<Buffer<float>> par{"par", 1};
Input<int> ncpu{"ncpu"};
Output<Buffer<float>> output{"output", 1};
Var x;
Func par_inv;
void generate() {
// precompute 1 / par[x]
par_inv(x) = fast_inverse(par(x));
// after che peak
Expr after_che_peak = tterm(x) > 0;
Expr first_term = -par_inv(0) * fast_exp(-tterm(x) * par_inv(0));
Expr second_term = par(2) * par(1) * fast_pow(1 + par(1) * tterm(x), par(2) - 1);
// before che peak
Expr third_term = -par_inv(3) * fast_exp(-tterm(x) * par_inv(3));
Expr fourth_term = par(5) * par(4) * fast_pow(1 + par(4) * tterm(x), par(5) - 1);
// final value
output(x) = 1.e9f * select(after_che_peak, first_term + second_term,
third_term + fourth_term);
}
void schedule() {
par_inv.bound(x, 0, 6);
par_inv.compute_root();
Var xo, xi;
// break x into two loops, one for ncpu tasks
output.split(x, xo, xi, output.extent() / ncpu)
// mark the task loop parallel
.parallel(xo)
// vectorize each thread's computation for 8-wide vector lanes
.vectorize(xi, 8);
output.print_loop_nest();
}
};
HALIDE_REGISTER_GENERATOR(ParamReTtermApproxOpt, param_re_tterm_approx_opt)

我可以在我的Surface Book上在略多于一秒的时间内运行477,000,000次迭代(ncpu=4)。批处理在这里非常重要，因为它可以实现矢量化。

请注意，使用双精度算术编写的等效程序比浮点数算法慢得多(20x)。虽然 Halide 不提供双打fast_版本，所以这可能不是苹果对苹果。无论如何，我会检查您是否需要额外的精度。