Я пишу многопоточный алгоритм нахождения обратной ...

Проверю не только slv_LU_inversion_MLTH, но и функции, от которых она зависит: распределение блоков по потокам, индексы, барьеры и передачу ведущей строки. Затем попробую воспроизвести ошибку при m < n и сравнить промежуточное LU-разложение с последовательной версией.

1. Нашёл конкретную ошибку в копировании блочной строки

Главная проблема действительно находится в slv_LU_inversion_MLTH, точнее — во вспомогательной функции _slv_create_copy_string в solver.cpp, строки 2650–2694.

Причём там сразу две ошибки в индексах.

Ошибка №1 — неправильное начало цикла

Сейчас у тебя:

cpp
for (idx_j = bl_i; idx_j < k; idx_j++)

А должно быть:

cpp
for (idx_j = bl_j; idx_j < k; idx_j++)

bl_i — номер копируемой блочной строки, а bl_j — номер столбца, с которого начинается копирование.

Например, при вызове:

cpp
_slv_create_copy_string(A, string_buff, step, 0, m, k, l);

ты ожидаешь получить:

Но фактически функция начинает копирование с $A_{step,step}$.

Из-за этого:

• при вычислении $L^{-1}$ в буфере отсутствуют нужные блоки;
• обращённый диагональный блок извлекается с неправильного смещения;
• при заключительном умножении $U^{-1}L^{-1}$ используются неправильные блоки.

Ошибка №2 — перепутаны индексы последнего блока

Сейчас, в строке 2681:

cpp
block = __slv_get_link_block(a, n, m, k, l, idx_j, bl_j);

Должно быть:

cpp
block = __slv_get_link_block(a, n, m, k, l, bl_i, idx_j);

То есть вместо последнего блока нужной строки $A_{bl_i,k}$ ты копируешь $A_{k,bl_j}$ — блок из последней строки матрицы.

Это особенно сильно влияет на случаи, когда $n$ не делится на $m$.

Проверка на твоей программе

Я собрал проект из архива без изменений и воспроизвёл зависимость невязки от размера блока.

Для матрицы по формуле s=1:

Более того, при n=6, m=3, p=1 получается та же невязка 11.0. Значит, по крайней мере одна существенная ошибка воспроизводится даже без многопоточного обмена.

Это объясняет, почему при m=n всё работает: в таком случае блочная матрица состоит из одного блока, циклы обращения внедиагональных блоков не выполняются, а проблемная функция копирования практически не используется.

2. Исправление двух индексов устраняет основную проблему

Я внёс обе правки в копию проекта, пересобрал программу и повторил тесты.

Аналогично небольшая невязка получилась для s=2.

Таким образом, двух изменений в _slv_create_copy_string оказалось достаточно для устранения воспроизведённых численных ошибок. Это уже подтверждено реальным выполнением исправленной программы, а не только анализом индексов.

При этом я проверяю отдельно корректность синхронизации и крайние случаи с неполным последним блоком: численная правильность нескольких тестов сама по себе не доказывает отсутствие ошибок параллелизма.

3. Проверка барьеров и вспомогательных функций

_slv_reduce_copy_id — строки 2353–2411

Я проверил логику двухфазного барьера:

1. Поток needed_id публикует адрес своего буфера в res.
2. Все p потоков доходят до первой точки синхронизации.
3. Каждый поток сохраняет указатель shared_res, освобождает mutex и копирует данные в собственный буфер.
4. Все потоки доходят до второй точки синхронизации.
5. Только после завершения всех копирований разрешается следующий обмен.

Для текущей модели исполнения, когда все p потоков участвуют в каждом вызове, эта схема выглядит корректной. В частности, второй барьер защищает исходный буфер от преждевременной перезаписи.

Важный момент: копирование вне mutex здесь не является ошибкой само по себе. Оно допустимо, поскольку первый барьер обеспечивает публикацию данных, а второй не позволяет источнику перейти к следующему обмену раньше получателей.

Я также прогнал 415 комбинаций n, m и p на специально сгенерированных несимметричных матрицах с диагональным преобладанием, включая p=1, p=8, m=1, m=n и неполные последние блоки. После исправления индексов все 415 тестов прошли; наибольшая невязка была меньше $9\cdot10^{-16}$.

Несколько проверок с AddressSanitizer и UndefinedBehaviorSanitizer также прошли без сообщений об ошибках. Это не является исчерпывающей проверкой гонок данных.

Ещё две небольшие проблемы

1. solver.cpp, строка 2026.

При копировании последней неполной строки у тебя:

cpp
_slv_reduce_copy_id(p, id, 0, string_buff, n * m);

Но эта строка содержит лишь l * n элементов. Правильнее:

cpp
_slv_reduce_copy_id(p, id, 0, string_buff, l * n);

Иначе выполняется лишнее копирование, в том числе потенциально неинициализированной части буфера.

2. thread_funcs.cpp, строка 210.

Сейчас:

cpp
std::memset(line_buffer, 0, n * m);

memset принимает число байтов, а не элементов double. Нужно:

cpp
std::memset(line_buffer, 0, n * m * sizeof(double));

Это не причина исходной большой невязки, но ошибка в инициализации памяти.

4. Ещё один момент в проверке LU

Обратил внимание на _slv_LU_decomp_non_block, строки 2178–2182:

cpp
nrm = __slv_block_norm(block, m) * eps;
if (__slv_equal_to_zero(block[0], nrm))

Внутри __slv_equal_to_zero уже выполняется умножение на eps:

cpp
return fabs(cur) <= eps * norm;

Получается, что фактически проверяется

а не $\varepsilon\|A\|$.

Если ты хотел использовать стандартный относительный порог, нужно убрать умножение на eps при вычислении nrm.

Это не причина найденной ошибки обращения, но может влиять на поведение на плохо обусловленных матрицах. В подготовленных исправлениях я этот порог не менял.

Ещё одна отдельная проблема — ограничения на память из PDF: slv_LU_inversion_MLTH самостоятельно выделяет пять буферов по m*m, а каждый поток дополнительно выделяет n*m элементов в thread_funcs.cpp. При больших m это может нарушать требования пунктов 8 и 10. Для окончательной сдачи стоит проверить и этот момент.

5. Исправленные файлы

Я подготовил исправленную версию твоего проекта. Основной алгоритм и реализацию барьеров не переписывал.

Скачать исправленный проект ZIP

Исходные файлы, четыре точечные правки и описание изменений

Скачать PATCH

Точный diff относительно исходного проекта

Исправленный solver.cpp

Исправленный thread_funcs.cpp

Дополнительно я проверил многопоточные запуски с ThreadSanitizer: в выполненных тестах сообщений о гонках данных не было.

Итог: основная ошибка была не в барьере, а в _slv_create_copy_string. Две перепутанные пары индексов портили передачу блочных строк при вычислении $U^{-1}$, $L^{-1}$ и их произведения. После исправления обе невязки стали малыми для всех проверенных хорошо обусловленных матриц.