我正在使用 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 * 128float阵列占用 65 kB，因此所有 6 个阵列应该很容易放入 CPU 的 L3 高速缓存（即 8 MB）。假设没有缓存在寄存器中，我在内部循环体中计算了 15 次内存访问。在 64 位系统上，这是每次迭代 120 字节，因此 5700 万 * 120 字节 = 6.8 GB/s。L3 高速缓存以 2.66 GHz 运行，因此处于同一数量级。我的猜测是内存确实是瓶颈。
为了加快速度，我尝试了以下方法：
用 编译
g++ -O3。（嗯，我从一开始就这样做了。）使用 OpenMP pragma 并行化超过 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 值可能会减慢速度。现在清除示例代码中的内存。不过，执行速度对我来说并没有什么不同。