有限差分法的性能问题

Question

我写了一段 C 代码，它使用有限差分法来估计值。这是一种平均方法。我分析了代码，发现一个iterate() 函数是最慢的。

void iterate(double data[][ARRAY_SIZE], int nx, int ny, int dx, int dy)
{
    for (int i = 0; i < nx; ++i)
    {
        for (int j = 0; j < ny; ++j)
        {
            if (i % (dx + 1) == 0 && j % (dy + 1) == 0)
                continue;
            else if (i == 0 && 0 < j && j < ny)
                data[i][j] = (data[i][j - 1] + data[i][j + 1] + data[i + 1][j]) / 3;
            else if (j == 0 && 0 < i && i < nx)
                data[i][j] = (data[i - 1][j] + data[i + 1][j] + data[i][j + 1]) / 3;
            else if (i == nx - 1 && 0 < j && j < ny)
                data[i][j] = (data[i][j - 1] + data[i][j + 1] + data[i - 1][j]) / 3;
            else if (j == ny - 1 && 0 < i && i < nx)
                data[i][j] = (data[i - 1][j] + data[i + 1][j] + data[i][j - 1]) / 3;
            else
                data[i][j] = (data[i - 1][j] + data[i + 1][j] + data[i][j - 1] + data[i][j + 1]) / 4;
        }
    }
}

此循环运行缓慢，我不确定我在这里缺少什么使其变慢。有没有更好的方法来做同样的事情？

使用 400x400 double 数组进行 2000 次迭代

real    0m1.950s
user    0m1.940s
sys 0m0.004s

Answer 1

这里有一些想法：

1. 看来ny 必须等于ARRAY_SIZE。您也可以将其作为参数省略，而只使用编译时常量。
2. 除最后一个之外的所有 if/else 子句仅适用于特定的行或列。所以把它们吊起来。例如，您可以将第一行和第一列作为一维循环处理，然后再在边缘之外处理整个矩阵，然后最后处理最右边的列和底行。

最后，你的核心循环应该更像这样：

for (int i = 1; i < nx - 1; ++i)
{
    for (int j = 1; j < ARRAY_SIZE - 1; ++j)
    {
        data[i][j] = (data[i - 1][j] + data[i + 1][j] + data[i][j - 1] + data[i][j + 1]) / 4;
    }
}

Answer 2

考虑这个实现：

void iterate(double data[][ARRAY_SIZE], int nx, int ny, int dx, int dy)
{
    // because nx - 1 and ny - 1 are used
    nx--;
    ny--;
    // because dx + 1 and dy + 1 are used
    dx++;
    dy++;

    int i = 0;
    int j = 0;

    // case i == 0 && 0 < j && j < ny
    for (j = 1; j < ny; ++j)
    {
        if (j % dy)
            data[0][j] = (data[i][j - 1] + data[i][j + 1] + data[i + 1][j]) / 3.0;
    }

    j = 0;

    // case j == 0 && 0 < i && i < nx
    for (i = 1; i < nx; ++i)
    {
        if (i % dx)
            data[i][0] = (data[i - 1][j] + data[i + 1][j] + data[i][j + 1]) / 3.0;
    }

    // default case
    for (i = 1; i < nx; ++i)
    {
        for (j = 1; j < ny; ++j)
        {
            if (i % dx || j % dy)
                data[i][j] = (data[i - 1][j] + data[i + 1][j] + data[i][j - 1] + data[i][j + 1]) * 0.25;
        }
    }

    // case i == nx && 0 < j && j < ny
    for (j = 1; j < ny; ++j)
    {
        if (nx % dx || j % dy)
            data[nx][j] = (data[i][j - 1] + data[i][j + 1] + data[i - 1][j]) / 3.0;
    }

    // case j == ny && 0 < i && i < nx
    for (i = 1; i < nx; ++i)
    {
        if (ny % dy || i % dx)
            data[i][ny] = (data[i - 1][j] + data[i + 1][j] + data[i][j - 1]) / 3.0;
    }
}

主要的三点是：

1. 为双 for 循环减少内循环的操作量
2. 只需执行一次即可减少琐碎操作的数量
3. 不要混合数据类型和强制强制（使用/ 3.0 和* 0.25）

我的代码中唯一没有解释的是i % dx || j % dy 等于!(i % dx == 0 && j % dy == 0)。