如何在 C 中改进/加速这个频率函数？答案

【问题标题】：How Can I Improve/SpeedUp This FrequentFunction in C?如何在 C 中改进/加速这个频率函数？
【发布时间】：2010-04-20 09:29:04
【问题描述】：

我怎样才能改进/加速这个频繁的功能？

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define M 10 // This is fixed
#define N 8  // This is NOT fixed

// Assumptions: 1. x, a, b and c are all arrays of 10 (M).
//              2. y and z are all matrices of 8 x 10 (N x M).
// Requirement: 1. return the value of ret;
//              2. get all elements of array c
float fnFrequentFunction(const float* x, const float* const* y, const float* const* z,
                         const float* a, const float* b, float *c, int n)
{
    register float tmp;
    register float sum;
    register float ret = 0;
    register const float* yy;
    register const float* zz;
    int i;

    for (i = 0; i < n; i++)  // M == 1, 2, 4, or 8
    {
        sum = 0;
        yy = y[i];
        zz = z[i];

        tmp = x[0] - yy[0]; sum += tmp * tmp * zz[0];
        tmp = x[1] - yy[1]; sum += tmp * tmp * zz[1];
        tmp = x[2] - yy[2]; sum += tmp * tmp * zz[2];
        tmp = x[3] - yy[3]; sum += tmp * tmp * zz[3];
        tmp = x[4] - yy[4]; sum += tmp * tmp * zz[4];
        tmp = x[5] - yy[5]; sum += tmp * tmp * zz[5];
        tmp = x[6] - yy[6]; sum += tmp * tmp * zz[6];
        tmp = x[7] - yy[7]; sum += tmp * tmp * zz[7];
        tmp = x[8] - yy[8]; sum += tmp * tmp * zz[8];
        tmp = x[9] - yy[9]; sum += tmp * tmp * zz[9];

        ret += (c[i] = log(a[i] * b[i]) + sum);
    }

    return ret;
}

// In the main function, all values are just example data.
int main()
{
    float x[M] = {0.001251f, 0.563585f, 0.193304f, 0.808741f, 0.585009f, 0.479873f, 0.350291f, 0.895962f, 0.622840f, 0.746605f};
    float* y[N];
    float* z[N];
    float a[M] = {0.870205f, 0.733879f, 0.711386f, 0.588244f, 0.484176f, 0.852962f, 0.168126f, 0.684286f, 0.072573f, 0.632160f};
    float b[M] = {0.871487f, 0.998108f, 0.798608f, 0.134831f, 0.576281f, 0.410779f, 0.402936f, 0.522935f, 0.623218f, 0.193030f};
    float c[N];

    float t1[M] = {0.864406f, 0.709006f, 0.091433f, 0.995727f, 0.227180f, 0.902585f, 0.659047f, 0.865627f, 0.846767f, 0.514359f};
    float t2[M] = {0.866817f, 0.581347f, 0.175542f, 0.620197f, 0.781823f, 0.778588f, 0.938688f, 0.721610f, 0.940214f, 0.811353f};
    int i, j;

    int n = 10000000;
    long start;

    // Initialize y, z for test example:
    for(i = 0; i < N; ++i)
    {
        y[i] = (float*)malloc(sizeof(float) * M);
        z[i] = (float*)malloc(sizeof(float) * M);

        for(j = 0; j < M; ++j)
        {
            y[i][j] = t1[j] * j;
            z[i][j] = t2[j] * j;
        }
    }


    // Speed test here:
    start = clock();
    while(--n)
        fnFrequentFunction(x, y, z, a, b, c, 8);
    printf("Time used: %ld\n", clock() - start);


    // Output the result here:
    printf("fnFrequentFunction == %f\n", fnFrequentFunction(x, y, z, a, b, c, 8));
    for(j = 0; j < N; ++j)
        printf("  c[%d] == %f\n", j, c[j]);
    printf("\n");


    // Free memory
    for(j = 0; j < N; ++j)
    {
        free(y[j]);
        free(z[j]);
    }

    return 0;
}

欢迎任何建议:-)

我在我的职能中犯了一个大错误，我感到很糟糕。上面的代码是新的。我现在正在重新检查以确保这是我需要的。

【问题讨论】：

x、y、z、a 和 b 多久变化一次？
这有点像家庭作业......不是说它是@Peter，但如果是，请标记为这样。
为什么是 y 和 z 指针数组而不是实际的 2d 数组？即使你的矩阵要改变大小，它也应该是一个包含宽度、高度、元素*的结构。双重间接在这里不好。
不，这不是家庭作业，这是真正的工作:-)
二维数组不能作为参数传递（第一维大小会被编译器省略，如M in y[M][N]）。

标签： c performance optimization

【解决方案1】：

把它放在循环之外

sum = 0;

tmp = x[0] - y[0]; sum += tmp * tmp * z[0];
tmp = x[1] - y[1]; sum += tmp * tmp * z[1];
tmp = x[2] - y[2]; sum += tmp * tmp * z[2];
tmp = x[3] - y[3]; sum += tmp * tmp * z[3];
tmp = x[4] - y[4]; sum += tmp * tmp * z[4];
tmp = x[5] - y[5]; sum += tmp * tmp * z[5];
tmp = x[6] - y[6]; sum += tmp * tmp * z[6];
tmp = x[7] - y[7]; sum += tmp * tmp * z[7];
tmp = x[8] - y[8]; sum += tmp * tmp * z[8];
tmp = x[9] - y[9]; sum += tmp * tmp * z[9];

【讨论】：

++ 完全正确。这需要每次调用计算一次。
我猜我是和你同时写的，只是慢了一点。
不过，这可能是一个简化的例子。
@Peter Lee：请注意，这段代码完全独立于 M。它只是在每次迭代中不必要地对相同数据重复相同的计算。它应该只在循环之前执行一次。
@Peter Lee：那么log() 有必要吗？切换到double 不够吗？我建议尽可能避免使用像 log() 这样的复杂函数，因为它们只是 DAMN SLOW。

【解决方案2】：

此功能非常适合 SIMD 处理。查看您的编译器文档，了解与 SSE 指令相对应的内在函数。
您可以分解sum 变量的依赖链。而不是单个sum 累加器，而是交替使用两个累加器sum1 和sum2 - 一个用于偶数，一个用于奇数索引。之后把它们加起来。
这里最大的性能瓶颈是log() 函数。检查近似值是否足够。这个计算也可以向量化——我相信英特尔发布了一个高性能数学库——包括函数的向量化版本，如log()。您可能喜欢使用它。
你在这里操作floats，log() 使用double 精度。请改用logf()。它可能（或可能不会）更快。肯定不会慢。
如果您的编译器理解 C99，请在作为函数参数的指针上放置 restrict 限定符。这告诉编译器这些数组不重叠，并可能帮助它生成更高效的代码。
更改矩阵在内存中的保存方式。不要使用指向不相交内存块的指针数组，而是使用大小为 M*N 的单个数组元素。

所以，总而言之，这就是函数的外观。这是便携式C99。使用特定于编译器的 SIMD 内部函数，这可以做得更快。

更新：请注意，我更改了输入矩阵的定义方式。矩阵是一个单一的大数组。

float fnFrequentFunction(const float *restrict x, const float *restrict y,
                         const float *restrict z, const float *restrict a,
                         const float *restrict b, float *restrict c, int n)
{
    float ret = 0;
    const float *restrict yy = y; //for readability
    const float *restrict zz = z; // -||-

    for (int i = 0; i < n; i++, yy += M, zz += M)  // n == 1, 2, 4, or 8
    {
        float sum = 0;
        float sum2 = 0;

        for(int j = 0; j < 10; j += 2)
        {
            float tmp  = x[j]   - yy[j];   sum  += tmp  * tmp  * zz[j];
            float tmp2 = x[j+1] - yy[j+1]; sum2 += tmp2 * tmp2 * zz[j+1];
        }
        sum += sum2;

        ret += (c[i] = logf(a[i] * b[i]) + sum);
    }
    return ret;
}

【讨论】：

我将查看分析结果。是的，log() 函数与其他行相比是个问题。
我会试一试，然后告诉你结果。

【解决方案3】：

使用memoization 缓存结果。这是一种时间/空间权衡优化。

在 Perl 中使用 memoize 包很容易做到这一点，并且可能在许多其他动态语言中。在 C 中，您需要自己滚动。

使用包装函数对参数进行哈希处理，并使用它来检查值是否已经计算过。如果有，请将其退回。如果没有，则传递给原函数，并缓存返回的结果。

或者，您可以在程序启动时预先计算您的查找表，或者甚至计算一次然后将其持久化，具体取决于您的需要。

【讨论】：

@Peter Lee：扩展。希望对您有所帮助:)
如果函数被非常频繁地调用并且具有不同的值（浮点输入很可能就是这种情况），这只会浪费时间和内存。
@San Jacinto：当然。尽管 OP 要求 any 建议，但我确实根据他的需要说。即使整个算法不是，算法的某些部分也可能是可记忆的。

【解决方案4】：

上述强度降低循环外 tmp 值的建议是正确的。我什至可能会考虑将这 10 行代码放入自己的 for 循环中，因为这样可以提高代码缓存效率。

除此之外，您还想知道您的目标处理器类型。如果它具有本机 SIMD 支持、FPU、它使用哪种缓存等。还取决于通过寄存器传递的参数数量，通过组合成单个结构和按引用传递来减少参数可能会给您带来一点提升。将 vars 声明为 register 可能有帮助，也可能没有帮助。再次分析和检查汇编器输出将回答这个问题。

由于 sum 在循环之前已知，您可以在循环之后添加 M * 它的值以进行提升。这只是在里面留下了 2 个 log muls。

如果 M 始终为 8 或具有其他已知模式，您可以进行一些小循环展开，但对日志调用的收益几乎为零。

唯一需要关注的主要内容是 log()。这是如何实施的？如果您的输入范围已知，您是否可以通过表查找推出自己的更快版本。更好的是，如果有足够的可用 RAM，请列出日志产品。

只是一些想法。

【讨论】：

【解决方案5】：

您使用编译器优化吗？

在现代编译器过时的变量之前注册。如果将它们与编译器优化一起使用，甚至会损害编译器的性能。例如 gcc 简单编译提供：

Time used: 8720000

并且没有寄存器浮动：

Time used: 8710000

我知道这并不多。

我假设你做了所有这些总和是为了避免 for 循环，因为你认为这要慢得多。它不是。现代编译器也会为您进行优化。

我认为一个很大的优化是使用一个表来记录日志，如果你不介意内存，那会更快，只有当你超出范围时才使用日志。

【讨论】：

【解决方案6】：

我想知道是否将其作为缩放整数而不是浮点数可能会加快速度。我不知道数据范围，所以我不知道这是否可能

【讨论】：

【解决方案7】：

除了安德烈的回答，还可以在循环中添加一些预取：

float fnFrequentFunction(const float* x, const float* y, const float* z,
                         const float *a, const float *b, float *c, int M)
{
    register float tmp;
    register float sum;
    register float ret = 0;
    int i;
    sum = 0;

    tmp = x[0] - y[0]; sum += tmp * tmp * z[0];
    tmp = x[1] - y[1]; sum += tmp * tmp * z[1];
    tmp = x[2] - y[2]; sum += tmp * tmp * z[2];
    tmp = x[3] - y[3]; sum += tmp * tmp * z[3];
    tmp = x[4] - y[4]; sum += tmp * tmp * z[4];
    tmp = x[5] - y[5]; sum += tmp * tmp * z[5];
    tmp = x[6] - y[6]; sum += tmp * tmp * z[6];
    tmp = x[7] - y[7]; sum += tmp * tmp * z[7];
    tmp = x[8] - y[8]; sum += tmp * tmp * z[8];
    tmp = x[9] - y[9]; sum += tmp * tmp * z[9];

    for (i = 0; i < M; i++)  // M == 1, 2, 4, or 8
    {
        //----------------------------------------
        // Prefetch data into the processor's cache
        //----------------------------------------
        float a_value = a[i];
        float b_value = b[i];
        float c_value = 0.0;

        //----------------------------------------
        // Calculate using prefetched data.
        //----------------------------------------
        c_value = log(a_value * b_value) + sum;
        c[i] = c_value;
        ret += c_value;
    }

    return ret;
}

您也可以尝试展开循环：

float a_value = 0.0;
float b_value = 0.0;
float c_value = 0.0;
--M;
switch (M)
{
    case 7:
        a_value = a[M];
        b_value = b[M];
        c_value = log(a_value * b_value) + sum;
        c[M] = c_value;
    ret += c_value;
    --M;
    case 6:
        a_value = a[M];
        b_value = b[M];
        c_value = log(a_value * b_value) + sum;
        c[M] = c_value;
    ret += c_value;
    --M;
    case 5:
        a_value = a[M];
        b_value = b[M];
        c_value = log(a_value * b_value) + sum;
        c[M] = c_value;
    ret += c_value;
    --M;
    case 4:
        a_value = a[M];
        b_value = b[M];
        c_value = log(a_value * b_value) + sum;
        c[M] = c_value;
    ret += c_value;
    --M;
    case 3:
        a_value = a[M];
        b_value = b[M];
        c_value = log(a_value * b_value) + sum;
        c[M] = c_value;
    ret += c_value;
    --M;
    case 2:
        a_value = a[M];
        b_value = b[M];
        c_value = log(a_value * b_value) + sum;
        c[M] = c_value;
    ret += c_value;
    --M;
    case 1:
        a_value = a[M];
        b_value = b[M];
        c_value = log(a_value * b_value) + sum;
        c[M] = c_value;
    ret += c_value;
    --M;
    case 0:
        a_value = a[M];
        b_value = b[M];
        c_value = log(a_value * b_value) + sum;
        c[M] = c_value;
    ret += c_value;
    break;
}

查看展开的版本，您可以将“+ sum”从“循环”中取出并在末尾添加为： ret += (M + 1) * sum; 因为sum 不会改变。

最后，另一种选择是一次执行所有乘法运算，然后进行所有log 计算，然后总结所有内容：

float product[8];
for (i = 0; i < M; ++i)
{
  product[i] = a[i] * b[i];
}
for (i = 0; i < M; ++i)
{
  c[i] = log(product);
  ret += c[i];
}
ret += M * sum;

【讨论】：

sum 被添加到存储在c 中的值中，因此无法从循环中删除。
我感觉很糟糕，我在之前的函数中出错了。请看更新的。对不起，伙计们。是的。 (sum) 不能从循环中删除，因为数组 c 需要它。
不幸的是，这都是错误的。数据足够小以适合 L1 缓存，因此无需预取。展开在这里是没有用的，因为这是 FP 繁重的代码，整数 ALU 一直在松弛。循环开销在这里基本上是免费的。展开只会浪费 L1 代码缓存空间，从而损害性能。

【解决方案8】：

如果你在 a 和 b 没有改变的情况下多次调用它，那么将 a 和 b 合并到 logab 中 logab[i] = log(a[i] * b[i]) 因为 a 和 b 没有被使用其他任何地方。

【讨论】：

【解决方案9】：

这似乎是一个高斯混合模型计算。几年前，我致力于优化同样的算法，该算法被用作语音处理程序的一部分。我调查了许多优化，就像您尝试做的那样，但从未发现任何使用直接 C 来获得超过百分之几的东西。我最大的收获来自使用 SIMD 指令重新编码基本的 GMM 内核。由于这仍然不能提供我想要的性能，所以下一步（也是最后一步）是使用 Nvidia GPU。这种方法很有效，但对它进行编程本身就是一件令人头疼的事情。

很抱歉，我无法提供更多帮助，但如果您坚持使用常规 CPU，我认为您只会获得名义上的速度。

【讨论】：