如何加快这个循环（在 C 中）？答案

【问题标题】：How can I speed-up this loop (in C)?如何加快这个循环（在 C 中）？
【发布时间】：2011-02-09 09:09:52
【问题描述】：

我正在尝试在 C 中并行化一个卷积函数。这是对两个 64 位浮点数数组进行卷积的原始函数：

void convolve(const Float64 *in1,
              UInt32 in1Len,
              const Float64 *in2,
              UInt32 in2Len,
              Float64 *results)
{
    UInt32 i, j;

    for (i = 0; i < in1Len; i++) {
        for (j = 0; j < in2Len; j++) {
            results[i+j] += in1[i] * in2[j];
        }
    }
}

为了允许并发（没有信号量），我创建了一个函数来计算 results 数组中特定位置的结果：

void convolveHelper(const Float64 *in1,
                    UInt32 in1Len,
                    const Float64 *in2,
                    UInt32 in2Len,
                    Float64 *result,
                    UInt32 outPosition)
{
    UInt32 i, j;

    for (i = 0; i < in1Len; i++) {
        if (i > outPosition)
            break;
        j = outPosition - i;
        if (j >= in2Len)
            continue;
        *result += in1[i] * in2[j];
    }
}

问题是，使用convolveHelper 会使代码减慢大约 3.5 倍（在单线程上运行时）。

关于如何在保持线程安全的同时加快 convolveHelper 的任何想法？

【问题讨论】：

当然convolveHelper 比较慢。你有一定数量的乘法和加法要做，还有一定数量的数组索引要计算。您添加的所有其他内容都会使其花费更长的时间，例如第一个 if、j = 和第二个 if，更不用说函数进入/退出了。
迈克，我的问题不是“为什么这个实现慢？”；很明显，convolveHelper在单个线程上应该更慢。我重写函数的全部目的是让它并行运行，也许使用 CUDA。我在问convolveHelper 是否可以在保持线程安全的同时更有效地编写。

标签： c optimization concurrency loops performance

【解决方案1】：

您可以在循环之前计算i 的正确最小值/最大值，而不是循环中的两个if 语句。
您正在分别计算每个结果位置。相反，您可以将results 数组拆分为块，并让每个线程计算一个块。块的计算类似于convolve 函数。

【讨论】：

【解决方案2】：

最明显的帮助是预先计算循环的开始和结束索引，并删除i 和j 上的额外测试（以及它们相关的跳转）。这个：

for (i = 0; i < in1Len; i++) {
   if (i > outPosition)
     break;
   j = outPosition - i;
   if (j >= in2Len)
     continue;
   *result += in1[i] * in2[j];
}

可以改写为：

UInt32 start_i = (in2Len < outPosition) ? outPosition - in2Len + 1 : 0;
UInt32 end_i = (in1Len < outPosition) ? in1Len : outPosition + 1;

for (i = start_i; i < end_i; i++) {
   j = outPosition - i;
   *result += in1[i] * in2[j];
}

这样，条件j >= in2Len 永远不会为真，循环测试本质上是测试i < in1Len 和i < outPosition 的组合。

理论上你也可以去掉对j的赋值，把i++变成++i，但是编译器可能已经为你做了这些优化。

【讨论】：

正是我想要的！谢谢泰勒:)
start_i 需要登录 outPosition - inLen 经常会被否定
在这种情况下，真正应该做的是测试outPosition > inLen是否设置start_i = 0，因为让i为负数是没有意义的。
嗯.. 我刚试过你的代码，但它不起作用:(虽然我跑得很快......
另外，不要这么快用“它不起作用”来驳回答案。这是你的程序，不是我的，所以我不会花很多时间坐在这里进行心理调试。该方法是有效的：预先计算循环边界。如果我用来计算循环边界公式的代数略有偏差，那么您应该能够弄清楚并修复它。

【解决方案3】：

让 convolve 帮助器在更大的集合上工作，计算多个结果，使用一个短的外循环。

并行化的关键是在线程之间的工作分配之间找到一个很好的平衡。使用的线程数不要超过 CPU 内核数。

在所有线程之间平均分配工作。遇到这种问题，各个线程工作的复杂度应该是一样的。

【讨论】：

【解决方案4】：

时域中的卷积变成傅里叶域中的乘法。我建议你获取一个快速的 FFT 库（如 FFTW）并使用它。您将从 O(n^2) 变为 O(n log n)。

算法优化几乎总是优于微优化。

【讨论】：

我知道，但我要解决的问题特别需要时域卷积。
是的，但是 Thomas 的论点是从输入数据的时域到频域的转换，使用 O(n) 中的乘法进行卷积，然后再转换回时域比优化更有效O(n^2) 卷积算法。您不必更改表示数据的方式，只需将其转换为执行卷积的一部分即可。
是的，运行时间是个问题。我正在看看人类是否可以区分不同音频卷积的结果。在两次测试之间等待半小时很烦人；）
@splicer - 相关地，您将“泄露”的光谱相乘，然后再转换回来。泄漏是另一个域中时域波形的忠实表示。该算法将信号零填充到 newLen = 2*max(in1Len, in2Len)，计算 FFT，将频谱相乘（复数乘法），逆 FFT，取第一个 in1Len+in2Len-1 个样本作为答案。这就是 Matlab/Octave 在内部实现它们的相关函数的方式。
@splicer - 它可以用数学方法显示，但我不能:) 这是一个推理，看看你是否觉得这令人信服：首先，假设你正在卷积的两个波形实际上是部分两个较长的波形是周期性脉冲，脉冲之间有死区。其次，看到线性卷积脉冲与循环卷积其周期性“父母”的一个周期相同，因为存在死区。第三，每个信号的一个周期的DFT是精确的（纹波和泄漏等），因此它们的乘积和逆变换也是精确的。希望有帮助！

【解决方案5】：

除非您的数组非常大，否则使用线程实际上不会有太大帮助，因为启动线程的开销将大于循环的成本。但是，让我们假设您的数组很大，并且线程是一个净赢。在这种情况下，我会执行以下操作：

忘记你当前的convolveHelper，这太复杂了，没有多大帮助。

将循环内部拆分为线程函数。 IE。只做

for (j = 0; j < in2Len; j++) {
    results[i+j] += in1[i] * in2[j];
}

进入它自己的函数，该函数将i 与其他所有参数一起作为参数。

让convolve 的主体简单地启动线程。为获得最大效率，请使用信号量来确保创建的线程数永远不会超过内核数。

【讨论】：

是的，我的数组相当大（每个 256KB）。您的方法的问题在于 += 操作不是线程安全的：i 和 j 的许多组合可以产生给定的 outPosition。

【解决方案6】：

答案在于简单数学而非多线程（已更新）

这就是为什么...

考虑 ab + ac

U 可以将其优化为 a*(b+c)（一少乘法）

在你的例子中，in2Len 在 inner-loop 中有不必要的乘法。可以消除的。

因此，如下修改代码应该会给我们 reqd 卷积：

(注意：以下代码返回circle-convolution，必须将其展开以获得linear-convolution结果。

void convolve(const Float64 *in1,
              UInt32 in1Len,
              const Float64 *in2,
              UInt32 in2Len,
              Float64 *results)
{
    UInt32 i, j;

    for (i = 0; i < in1Len; i++) {

        for (j = 0; j < in2Len; j++) {
            results[i+j] += in2[j];
        }

        results[i] = results[i] * in1[i];

    }
}

这应该给 U 带来最大的性能提升。试试我们看看！！

祝你好运！！

CVS @ 2600 赫兹

【讨论】：

这将给出完全不同的结果。它不执行卷积。
当然，但是当 O(n log n) 是一个解决方案时，您仍在优化 O(n^2) 算法。当 n 增长时，O(n log n) 总是击败 O(n^2)。也许您的数据不够大，无法遇到这种转变，O(n^2) 的微优化仍然可以改善结果。

【解决方案7】：

我终于弄清楚了如何正确地预先计算开始/结束索引（Tyler McHenry 和 interjay 给出的建议）：

if (in1Len > in2Len) {
    if (outPosition < in2Len - 1) {
        start = 0;
        end = outPosition + 1;
    } else if (outPosition >= in1Len) {
        start = 1 + outPosition - in2Len;
        end = in1Len;
    } else {
        start = 1 + outPosition - in2Len;
        end = outPosition + 1;
    }
} else {
    if (outPosition < in1Len - 1) {
        start = 0;
        end = outPosition + 1;
    } else if (outPosition >= in2Len) {
        start = 1 + outPosition - in2Len;
        end = in1Len;
    } else {
        start = 0;
        end = in1Len;
    }
}

for (i = start; i < end; i++) {
    *result = in1[i] * in2[outPosition - i];
}

不幸的是，预先计算索引会导致执行时间没有明显减少 :(

【讨论】：