Я пишу разреженный матричный решатель, используя метод Гаусса-Зейделя.По профилированию я определил, что около половины моей программы тратится внутри решателя.Критически важная часть производительности выглядит следующим образом:
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.Я должен перейти на алгоритм Якоби, чтобы избежать чтения и записи в тот же массив.Это требует, чтобы я выполнял вдвое больше итераций, что приводит к чистому результату с примерно одинаковой скоростью.
Работа с деталями реализации тела цикла, например использование указателей вместо индексов.Без эффекта.
Какой лучший способ ускорить этого парня?Поможет ли это переписать внутреннее тело в сборке (я должен сначала изучить это)?Должен ли я запустить это на GPU вместо этого (что я знаю, как это сделать, но это такая проблема)?Какие-нибудь другие яркие идеи?
(NB I do принимает «нет» за ответ, например: «это нельзя сделать значительно быстрее, потому что ...»)
Обновление : по запросу, вот полная программа:
#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 могут замедлить процесс.Теперь очищаем память в примере кода.Тем не менее, для меня нет никакой разницы в скорости выполнения.