c++ - C代码的半随机慢度

Question

我有一个更长的代码，在保持问题存在的同时尽可能地减少了代码。我的代码针对不同的参数值运行 MCMC 计算。对于某些值的组合，代码运行的时间要长得多，比典型情况慢大约 100 倍。然而，它不应该，因为操作的数量不依赖于参数值。

我在 AMD64 Linux 机器上运行它，glibc-2.17 在 Gentoo 上用 GCC 4.8.1 编译。编译标志看起来并不重要。我还在配备较旧的 AMD64 处理器的不同 Gentoo 机器上对其进行了测试，结果是一样的。

我做了一堆测试:

• 我已经尝试使用 Valgrind 进行调试，但没有发现内存问题或其他令人讨厌的问题。
• 其次，我尝试在修复有问题的参数值的情况下运行代码，但没有遇到速度慢的问题。
• 我尝试在迭代之间放置 sleep(4)，但没有任何改变。

当迭代命中 k = 1 时，问题显示 i = j = 0，这转化为 mu[0] = -0.05，mu[1] = -0.05 和 mu[2] = 0.05。正如我所说，对所有迭代使用此固定值可以消除我所看到的问题。

以下是消除问题的方法:

• 更改限制。
• 修复 mu[] 系数。
• 从计算中删除 dW3。
• 删除 rand()。
• 移除 q 的计算。
• 删除 s[j] 的更新。

我已经阅读了一些关于 slowpow 的内容，因此尝试通过编写我自己的版本来消除 exp。这解决了我在使用此 MWE 时遇到的问题，但在将重新实现的 exp 放置在生产代码中时却没有。

问题:半随机缓慢的原因是什么？

MWE 代码如下。非常感谢所有有关如何进行的帮助和建议。

注意:这段代码是用g++编译的，虽然它本质上是C。更改编译器不会改变任何内容。

关于分支预测:使用

删除其中一个if语句

q = exp(dW);        
q = q / (1.0 + q); 

无论dW 的值是多少都不会改变代码的行为；如果这确实是由于分支预测，则必须是由于第二个 if。

---

#include <cstdio>
#include <cstdlib>
#include <cmath>

inline int index(int const i, int const j, int const n)
{
    return (i + n) % n + ((j + n) % n) * n;
}

void get_sample(int* s, int n, double* mu)
{
    for (int i = 0; i < 10 * n * n; i++)
    {
        int j = i % (n * n); 
        int x = j % n;
        int y = (j - x) / n;

        double dW1 = mu[0] * (s[index(x - 1, y, n)] + s[index(x + 1, y, n)] + s[index(x, y - 1, n)] + s[index(x, y + 1, n)]);
        double dW2 = mu[1] * (s[index(x - 1, y - 1, n)] + s[index(x + 1, y - 1, n)] + s[index(x + 1, y + 1, n)] + s[index(x - 1, y + 1, n)]);
        double dW3 = mu[2] * (s[index(x - 1, y, n)] * s[index(x - 1, y - 1, n)] * s[index(x, y - 1, n)] + s[index(x - 1, y, n)] * s[index(x - 1, y + 1, n)] * s[index(x, y + 1, n)]
                                        + s[index(x, y + 1, n)] * s[index(x + 1, y + 1, n)] * s[index(x + 1, y, n)] + s[index(x + 1, y, n)] * s[index(x + 1, y - 1, n)] * s[index(x, y - 1, n)]);

        double dW = 2.0 * (dW1 + dW2 + dW3);

        double q;
        if (dW < 0.0)
        {   
            q = exp(dW);

            q = q / (1.0 + q); 
        }
        else
        {
            q = exp(-dW);

            q = 1.0 / (1.0 + q);
        } 

        double p = ((double) rand()) / ((double) RAND_MAX);

        if (p < q)
        {
            s[j] = 1;
        }
        else
        {
            s[j] = -1;
        }
    }
}

int main(int argc, char** argv)
{
    double mu[3];

    double limits[6] = {-0.05, 0.8, -0.05, 0.45, -0.45, 0.05};

    int s[16];

    for (int i = 0; i < 16; i++)
    {
        s[i] = -1;
    }

    for (int k = 0; k < 2; k++)
    {
        for (int j = 0; j < 2; j++)
        {
            for (int i = 0; i < 2; i++)
            {               
                mu[0] = limits[0] + ((limits[1] - limits[0]) * i);
                mu[1] = limits[2] + ((limits[3] - limits[2]) * j);
                mu[2] = limits[4] + ((limits[5] - limits[4]) * k);

                printf(" Computing (% .6lf, % .6lf, % .6lf)...\n", mu[0], mu[1], mu[2]);

                for (int sample = 0; sample < 1000; sample++)
                {
                    get_sample(s, 4, mu);
                }                           
            }
        }
    }               

    return 0;
}

Answer 1

However, it shouldn't, because the number of operations does not depend on parameter values.

不过，浮点运算的速度确实取决于参数值。如果您在计算中引入 NaN 或其他异常值(我没有检查代码)，它将大大降低您的浮点性能。

编辑:我围绕 exp() 手动分析(使用简单的 rdtsc 计数)，很容易将“好”和“坏”的情况分类。当我打印坏情况时，它都是 dW ~= 0。如果将这种情况分开，您将获得均匀的性能:

    double q;
    if (dW < -0.1e-15)
    {
        q = exp(dW);

        q = q / (1.0 + q);
    }
    else if (dW > 0.1e-15)
    {
        q = exp(-dW);

        q = 1.0 / (1.0 + q);
    }
    else
    {
        q = 0.5;
    }