时间序列堆叠——C算法优化

Question

我正在解决一个问题，我想堆叠在不同位置记录的时间序列并提取相干信号。繁重的工作是在 C 中完成的，使用 Python 包装器来提供更友好的界面。我已经达到了对算法的理论正确性感到满意的地步，并希望尽可能地对其进行优化。我对 C 语言的理解足以编写一些可以与 openMP 并行运行的东西，但仅此而已。

问题的优化很重要，因为我正在处理大型数据集 - 最多可堆叠 200 个时间序列，采样率高达 1000Hz，记录数月至数年。使用合理的计算设施，处理可以持续几天到几周。我在连续时间序列的块上运行此堆叠步骤，以免占用内存。

我有几个问题：

1. 有什么明显的我遗漏的东西会有所帮助（通过编译器进行优化，简化算法）？

2. 到目前为止，最显着的收获是优化标志 -Ofast - 我已经阅读并只是想了解更多为什么这更快以及它是否“安全” ' 出于我的目的？

3. 我应该在哪里（除了通过 SO 拖网）了解有关此类问题的更多信息？我的研究中有其他问题想用 C 来解决！

算法

我将每个位置的时间序列连续堆叠在一个 3-D 网格体中。完成给定单元格的完整堆栈后，我需要对结果取幂并按贡献时间序列的数量进行归一化。

#define MAX(a,b) (((a)>(b))?(a):(b))

EXPORT void migrate(double *sigPt, int32_t *indPt, double *mapPt, int32_t fsmp, int32_t lsmp, int32_t nsamp, int32_t nstation, int32_t avail, int64_t ncell, int64_t threads)
{
    double  *stnPt, *stkPt, *eStkPt;
    int32_t *ttpPt;
    int32_t ttp;
    int32_t to, tm, st;
    int64_t cell;

    #pragma omp parallel for private(cell,stkPt,eStkPt,ttpPt,st,ttp,tm) num_threads(threads)
    for (cell=0; cell<ncell; cell++)
    {
        stkPt = &mapPt[cell * (int64_t) nsamp];
        eStkPt = &mapPt[cell * (int64_t) nsamp];
        ttpPt = &indPt[cell * (int64_t) nstation];
        for(st=0; st<nstation; st++)
        {
            ttp   = MAX(0,ttpPt[st]);
            stnPt = &sigPt[st*(fsmp + lsmp + nsamp) + ttp + fsmp];
            for(tm=0; tm<nsamp; tm++)
            {
                stkPt[tm] += stnPt[tm];
            }
        }
        for(tm=0; tm<nsamp; tm++)
        {
            eStkPt[tm] = exp(stkPt[tm] / avail);
        }
    }
}

我目前正在编译：

gcc -shared -fPIC -std=gnu99 ./source.c -fopenmp -Ofast -lm -o ./output

我已阅读：

Profiling python C extensions

What GCC optimization flags and techniques are safe across CPUs?

等等。如果我重复一个问题/我的问题定义不明确，请道歉。

Answer 1

有什么明显的我遗漏的东西会有所帮助（通过编译器进行优化，简化算法）？

提议的代码看起来相当不错（GCC 应该能够对其进行矢量化，并且并行化似乎还可以）。但这里有一些可能的建议来提高性能：

• exp(stkPt[tm] / avail) 可以重写为可能更快的表达式pow(availFactor, stkPt[tm])，其中availFactor 是一个在循环外部定义并设置为exp(1.0 / avail) 的常量。正如@tim18 所建议的那样，您还应该检查此循环是否已矢量化，因为指数计算通常很慢。
• 您可以尝试使用编译器标志-march=native 以牺牲可移植性较差的二进制文件来加快代码速度（如果您不想降低旧处理器的可移植性，您可以尝试-mtune=native，它是一般不如-march=native)。关于你的处理器，这个选项可以让 GCC 使用 AVX 和 FMA 指令（在相对较新的处理器上可用），这应该可以加速你的代码。您也可以使用-mavx 和-mfma 手动启用此类功能。
• 您可以尝试调整您的算法，以便包含stkPt[tm] += stnPt[tm] 的热循环主要用于缓存中的数据。这一点非常重要。事实上，算法中最热门的部分可能是内存受限的。一个好的起点是进行平铺（例如同时处理 2 个或 4 个 nstation）。
• 如果结果精度足够好，请考虑使用单精度浮点类型而不是双精度。

不要忘记检查结果，因为这些优化可能会影响代码的准确性！

到目前为止，最重要的收获是优化标志 -Ofast - 我已经阅读并只是想更多地了解为什么它会更快以及它是否对我的目的“安全”？

使用-Ofast 是不安全的。实际上，此选项启用-ffast-math，从而启用更多选项，如-funsafe-math-optimizations 或-ffinite-math-only。因此，您的代码不应处理 NaN 值、无穷大和计算的准确性可能会低得多。关于您的期望，这可能是也可能不是问题。 请注意，此选项有助于 GCC 加快指数计算（感谢矢量化）。