Я пишу разреженный матричный решатель, используя метод Гаусса-Зейделя. По профилированию я решил, что около половины моего времени программы проводится внутри решателя. Критическая часть производительности выглядит следующим образом:
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, это может быть запущено около 3500 раз в секунду на Intel i7 920. Это означает, что тело внутреннего цикла работает 3500 * 128 * 128 = 57 миллионов раз в секунду. Поскольку задействована только какая-то простая арифметика, это поражает меня как низкое число для процессора с частотой 2,66 ГГц.
Возможно, это не ограничено мощностью процессора, а пропускной способностью памяти? Ну, один массив 128 * 128 float питается 65 кБ, поэтому все 6 массивов должны легко вписываться в кеш процессора L3 (что составляет 8 МБ). Предполагая, что ничего не кэшируется в регистрах, я подсчитываю 15 обращений к памяти во внутреннем цикле. В 64-битной системе это составляет 120 байт на итерацию, поэтому 57 миллионов * 120 байт = 6,8 ГБ/с. Кэш L3 работает на частоте 2,66 ГГц, поэтому он имеет тот же порядок величины. Я предполагаю, что память действительно является узким местом.
Чтобы ускорить это, я попытался выполнить следующее:
-
Скомпилируйте с помощью
g++ -O3. (Ну, я делал это с самого начала.) -
Параллелизация более 4 ядер с использованием OpenMP-прагм. Я должен перейти к алгоритму Якоби, чтобы избежать чтения и записи в тот же массив. Это требует, чтобы я делал в два раза больше итераций, что приводило к чистому результату примерно той же скорости.
-
Реализация деталей реализации тела цикла, например, с использованием указателей вместо индексов. Нет эффекта.
Какой лучший способ ускорить этот парень? Помогло бы оно переписать внутреннее тело в собрании (я должен был бы это узнать первым)? Должен ли я запускать это на графическом процессоре (что я знаю, как это сделать, но это такая проблема)? Любые другие яркие идеи?
(N.B. Я принимаю "нет" для ответа, как в: "это невозможно сделать значительно быстрее, потому что..." )
Обновить: по запросу, здесь полная программа:
#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 могут замедлять работу. Теперь очистка памяти в примере кода. Тем не менее, это не имеет никакого значения для меня в скорости выполнения.
