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で1秒間に約3500回実行できます。これは、内側のループ本体が1秒間に3500 * 128 * 128=5,700万回実行されることを意味します。いくつかの単純な計算しか含まれていないので、それは2.66GHzプロセッサの低い数値として私を驚かせます。
多分それはCPUパワーによってではなく、メモリ帯域幅によって制限されますか？1つの128*128floatアレイは65kBを消費するため、6つのアレイすべてがCPUのL3キャッシュ（8 MB）に簡単に収まるはずです。レジスタに何もキャッシュされていないと仮定すると、内部ループ本体で15回のメモリアクセスをカウントします。64ビットシステムでは、これは反復ごとに120バイトであるため、5,700万*120バイト=6.8GB/秒です。L3キャッシュは2.66GHzで動作するため、同じ桁数です。私の推測では、記憶は確かにボトルネックです。
これをスピードアップするために、私は次のことを試みました：
でコンパイルし
g++ -O3ます。（まあ、私は最初からこれをやっていた。）OpenMPプラグマを使用した4コア以上の並列化。同じ配列からの読み取りと同じ配列への書き込みを回避するために、Jacobiアルゴリズムに変更する必要があります。これには、2倍の反復を行う必要があり、ほぼ同じ速度の最終結果が得られます。
インデックスの代わりにポインタを使用するなど、ループ本体の実装の詳細をいじる。無効。
この男をスピードアップするための最良のアプローチは何ですか？アセンブリで内部ボディを書き直すのに役立ちますか（最初にそれを学ぶ必要があります）？代わりにこれをGPUで実行する必要があります（これは方法を知っていますが、とても面倒です）？他に明るいアイデアはありますか？
（NB私は答えに「いいえ」を取ります：「それは大幅に速くすることはできません、なぜなら...」）
更新：要求に応じて、完全なプログラムは次のとおりです。
#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
（私の「実際の」プログラムは他の多くのことも行うので、毎秒3500回の反復ではなく8000回実行します。しかしそれは代表的なものです。）
アップデート2：NaNとInfの値が遅くなる可能性があるため、単一化された値は代表的ではない可能性があると言われました。次に、サンプルコードのメモリをクリアします。ただし、実行速度には違いはありません。