Я написал этот SOR-код. Не беспокойтесь о том, что делает этот алгоритм, это не проблема. Но только ради полноты: он может решить линейную систему уравнений, в зависимости от того, насколько хорошо кондиционирована система.
Я запускаю его с плохой обработанной матрицей строк 2097152 (которая никогда не сходится) с не более 7 ненулевыми столбцами на строку.
Перевод: внешний цикл do-while будет выполнять 10000 итераций (значение, которое я передаю как max_iters), средний for будет выполнять итерации 2097152, разделенные на куски work_line, разделенные между потоками OpenMP. Самый внутренний цикл for будет содержать 7 итераций, за исключением очень немногих случаев (менее 1%), где он может быть меньше.
Существует зависимость данных между потоками в значениях массива sol. Каждая итерация середины for обновляет один элемент, но считывает до 6 других элементов массива. Поскольку SOR не является точным алгоритмом, при чтении он может иметь любое из предыдущего или текущего значения в этой позиции (если вы знакомы с решателями, это Гаусс-Сидель, который терпит поведение Якоби в некоторых местах ради parallelism).
typedef struct{
size_t size;
unsigned int *col_buffer;
unsigned int *row_jumper;
real *elements;
} Mat;
int work_line;
// Assumes there are no null elements on main diagonal
unsigned int solve(const Mat* matrix, const real *rhs, real *sol, real sor_omega, unsigned int max_iters, real tolerance)
{
real *coefs = matrix->elements;
unsigned int *cols = matrix->col_buffer;
unsigned int *rows = matrix->row_jumper;
int size = matrix->size;
real compl_omega = 1.0 - sor_omega;
unsigned int count = 0;
bool done;
do {
done = true;
#pragma omp parallel shared(done)
{
bool tdone = true;
#pragma omp for nowait schedule(dynamic, work_line)
for(int i = 0; i < size; ++i) {
real new_val = rhs[i];
real diagonal;
real residual;
unsigned int end = rows[i+1];
for(int j = rows[i]; j < end; ++j) {
unsigned int col = cols[j];
if(col != i) {
real tmp;
#pragma omp atomic read
tmp = sol[col];
new_val -= coefs[j] * tmp;
} else {
diagonal = coefs[j];
}
}
residual = fabs(new_val - diagonal * sol[i]);
if(residual > tolerance) {
tdone = false;
}
new_val = sor_omega * new_val / diagonal + compl_omega * sol[i];
#pragma omp atomic write
sol[i] = new_val;
}
#pragma omp atomic update
done &= tdone;
}
} while(++count < max_iters && !done);
return count;
}
Как вы можете видеть, внутри параллельной области нет блокировки, поэтому, чему они нас всегда учат, это проблема 100% параллельной. Это не то, что я вижу на практике.
Все мои тесты выполнялись на процессоре Intel (R) Xeon (R) E5-2670 v2 @2.50GHz, 2 процессора, по 10 ядер каждый, с поддержкой гиперпотока, суммируя до 40 логических ядер.
В моих первых запусках, work_line был зафиксирован на 2048, а количество потоков варьировалось от 1 до 40 (всего 40 прогонов). Это график с временем выполнения каждого запуска (в секундах x число потоков):
Удивление было логарифмической кривой, поэтому я подумал, что, поскольку рабочая строка была настолько большой, общие кэши были не очень хорошо использованы, поэтому я выкопал этот виртуальный файл /sys/devices/system/cpu/cpu0/cache/index0/coherency_line_size, который сказал мне, что этот процессор L1-кеш синхронизирует обновления в группах по 64 байта (8 удваивается в массиве sol). Поэтому я установил work_line в 8:
Затем я подумал, что 8 было слишком низким, чтобы избежать киосков NUMA и установить work_line в 16:
Во время выполнения вышеизложенного я подумал: "Кто я такой, чтобы предсказать, что work_line хорошо? Позволяет просто увидеть...", и планировать запуск каждого work_line с 8 до 2048, шаг 8 (т.е. каждый несколько строк кеша, от 1 до 256). Результаты для 20 и 40 потоков (в секундах x размер раскола среднего цикла for, разделенный между потоками):
Я считаю, что случаи с низким work_line плохо страдают от синхронизации кеширования, в то время как более крупный work_line не дает преимуществ за пределами определенного количества потоков (я предполагаю, что путь памяти является узким местом). Очень грустно, что проблема, которая кажется на 100% параллельной, представляет такое плохое поведение на реальной машине. Поэтому, прежде чем я убежден, что многоядерные системы - очень хорошо проданная ложь, я сначала прошу вас:
Как я могу сделать этот масштаб шкалы линейно равным числу ядер? Что мне не хватает? Есть ли что-то в этой проблеме, которое делает это не так хорошо, как кажется сначала?
Обновление
Следуя рекомендациям, я тестировал как с static, так и dynamic планированием, но удалял чтение/запись атома массивом sol. Для справки, синяя и оранжевая линии совпадают с предыдущим графиком (только до work_line = 248;). Желтые и зеленые линии - новые. Для чего я мог видеть: static делает существенную разницу для низкого work_line, но после 96 преимущества dynamic перевешивают его накладные расходы, делая это быстрее. Атомные операции не имеют никакого значения.
