我正在使用 Gauss-Seidel 方法编写稀疏矩阵求解器。通过分析，我确定我的程序大约一半的时间都花在了求解器中。性能关键部分如下:

size_t ic = d_ny + 1, iw = d_ny, ie = d_ny + 2, is = 1, in = 2 * d_ny + 1;
for (size_t y = 1; y < d_ny - 1; ++y) {
    for (size_t x = 1; x < d_nx - 1; ++x) {
        d_x[ic] = d_b[ic]
                - d_w[ic] * d_x[iw] - d_e[ic] * d_x[ie]
                - d_s[ic] * d_x[is] - d_n[ic] * d_x[in];
        ++ic; ++iw; ++ie; ++is; ++in;
    }
    ic += 2; iw += 2; ie += 2; is += 2; in += 2;
}

所有涉及的数组都是 float 类型。实际上，它们不是数组，而是具有重载 [] 运算符的对象，(我认为)应该对其进行优化，但定义如下:

inline float &operator[](size_t i) { return d_cells[i]; }
inline float const &operator[](size_t i) const { return d_cells[i]; }

对于 d_nx = d_ny = 128，这可以在 Intel i7 920 上每秒运行大约 3500 次。这意味着内部循环体每次运行 3500 * 128 * 128 = 5700 万次第二。由于只涉及一些简单的算术运算，这让我觉得对于 2.66 GHz 处理器来说这是一个很小的数字。

也许它不受 CPU 能力的限制，而是受内存带宽的限制？好吧，一个 128 * 128 float 数组占用 65 kB，因此所有 6 个数组都应该很容易放入 CPU 的 L3 缓存(即 8 MB)。假设没有缓存在寄存器中，我在内部循环体中计算了 15 次内存访问。在 64 位系统上，这是每次迭代 120 字节，因此 5700 万 * 120 字节 = 6.8 GB/s。 L3 高速缓存以 2.66 GHz 运行，因此处于同一数量级。我的猜测是内存确实是瓶颈。

为了加快速度，我尝试了以下方法:

• 使用 g++ -O3 编译。 (嗯，我从一开始就是这样做的。)

• 使用 OpenMP 编译指示并行化超过 4 个内核。我必须更改为 Jacobi 算法以避免读取和写入同一个数组。这需要我进行两倍的迭代，从而得到大致相同速度的最终结果。

• 摆弄循环体的实现细节，例如使用指针而不是索引。没有效果。

加速这个家伙的最佳方法是什么？在汇编中重写内部主体是否有帮助(我必须先学习)？我应该在 GPU 上运行它吗(我知道该怎么做，但这太麻烦了)？还有其他好主意吗？

(注意，我确实回答“不”，例如:“它不能做得更快，因为......”)

更新:根据要求，这是一个完整的程序:

#include <iostream>
#include <cstdlib>
#include <cstring>

using namespace std;

size_t d_nx = 128, d_ny = 128;
float *d_x, *d_b, *d_w, *d_e, *d_s, *d_n;

void step() {
    size_t ic = d_ny + 1, iw = d_ny, ie = d_ny + 2, is = 1, in = 2 * d_ny + 1;
    for (size_t y = 1; y < d_ny - 1; ++y) {
        for (size_t x = 1; x < d_nx - 1; ++x) {
            d_x[ic] = d_b[ic]
                - d_w[ic] * d_x[iw] - d_e[ic] * d_x[ie]
                - d_s[ic] * d_x[is] - d_n[ic] * d_x[in];
            ++ic; ++iw; ++ie; ++is; ++in;
        }
        ic += 2; iw += 2; ie += 2; is += 2; in += 2;
    }
}

void solve(size_t iters) {
    for (size_t i = 0; i < iters; ++i) {
        step();
    }
}

void clear(float *a) {
    memset(a, 0, d_nx * d_ny * sizeof(float));
}

int main(int argc, char **argv) {
    size_t n = d_nx * d_ny;
    d_x = new float[n]; clear(d_x);
    d_b = new float[n]; clear(d_b);
    d_w = new float[n]; clear(d_w);
    d_e = new float[n]; clear(d_e);
    d_s = new float[n]; clear(d_s);
    d_n = new float[n]; clear(d_n);
    solve(atoi(argv[1]));
    cout << d_x[0] << endl; // prevent the thing from being optimized away
}

我编译运行如下:

$ g++ -o gstest -O3 gstest.cpp
$ time ./gstest 8000
0

real    0m1.052s
user    0m1.050s
sys     0m0.010s

(它每秒执行 8000 次而不是 3500 次迭代，因为我的“真实”程序也执行许多其他操作。但它具有代表性。)

更新 2:有人告诉我，未初始化的值可能不具有代表性，因为 NaN 和 Inf 值可能会减慢速度。现在清除示例代码中的内存。不过，执行速度对我来说并没有什么不同。

最佳答案

几个想法:

1. 使用 SIMD。您可以一次将每个数组中的 4 个 float 加载到 SIMD 寄存器中(例如 Intel 上的 SSE，PowerPC 上的 VMX)。这样做的缺点是某些 d_x 值将是“陈旧的”，因此您的收敛速度会受到影响(但不如 jacobi 迭代那么糟糕)；很难说加速是否抵消了它。

2. 使用SOR .它很简单，不会增加太多计算量，并且可以很好地提高您的收敛速度，即使对于相对保守的松弛值(比如 1.5)也是如此。

3. 使用共轭渐变。如果这是用于流体模拟的投影步骤(即强制不可压缩性)，您应该能够应用 CG 并获得更好的收敛速度。一个好的预处理器会带来更多帮助。

4. 使用专门的求解器。如果线性系统来自Poisson equation ，您可以比使用基于 FFT 的方法的共轭梯度做得更好。

如果您能详细说明您要解决的系统是什么样的，我可能会就 #3 和 #4 提供更多建议。

关于c++ - 如何加快我的稀疏矩阵求解器？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/2388196/