Параллельный алгоритм Гаусса средствами std::thread

#include <iostream>
#include <iomanip>
#include <cstdio>
#include <Windows.h>
#include <thread>
#include <condition_variable>
#include <mutex>
#include <atomic>
#include <cmath>
 
using namespace std;
 
condition_variable cv;
condition_variable cv2;
 
mutex cv_m;   // мьютекс для condition variable cv (ожидание i-й строки)
mutex cv2_m;  // мьютекс для condition variable cv2 (ожидание завершения i-й итерации)
mutex cout_m; // мьютекс для потока вывода
 
int i_row = 0; // индекс текущей строки, хранящейся в r
atomic_int i_counter(0); // переменная для подсчёта количества потоков, ожидающих окончания i-й итерации
 
//заполнение матрицы
void FillMatrix(double *Matrix[], int nSize)
{
    for (int i = 0; i < nSize; i++)
    {
        for (int j = 0; j < nSize; j++)
            Matrix[i][j] = nSize - fabs(i - j) - 1;
        Matrix[i][nSize] = 0;
    }
}
 
// матрица разбивается построчно, всего потоков rank_count
// в массив r сохраняется текущая строка, которую вычитаем из нижестоящих
 
void GaussThread(double *Matrix[], int nSize, double *x,
                 int my_rank, int rank_count, double *r)
{
    double tmp;
 
    // прямой проход по алгоритму Гаусса
    for (int i = 0; i < nSize; i++)
    {
        if (i % rank_count == my_rank)
        {
            // для потока, который работает с i-й строкой
            lock_guard<mutex> lk(cv_m);
 
            // вычисляем строку
            tmp = Matrix[i][i];
            for (int j = i; j <= nSize; j++)
                Matrix[i][j] /= tmp;
 
            // сохраняем её в массив r
            memcpy(r, Matrix[i], (nSize + 1) * sizeof(double));
 
            // устанавливаем номер текущей строки, остальные потоки теперь могу пользоваться массивом r
            i_row = i;
            
            cout_m.lock();
            cout << "rank = " << my_rank <<  ", setted i_row = " << i_row << endl;
            cout_m.unlock();
 
            // будим остальные потоки
            cv.notify_all();
        }
        else
        {
            // для потока, который НЕ работает с i-й строкой
            unique_lock<mutex> lk(cv_m);
            
            cout_m.lock();
            cout << "rank = " << my_rank <<  ", wait i_row..." << endl;
            cout_m.unlock();
            
            // ждём пока буфер r не будет готов, тогда же будет выполняться i_row == i
            cv.wait(lk, [&]{ return i_row == i; });
 
            // строка в буфере r готова, можем работать дальше
            cout_m.lock();
            cout << "rank = " << my_rank <<  ", getted i_row = " << i << endl;
            cout_m.unlock();
        }
 
        // этот цикл выполняется параллельно, вычисляются строки ниже текущей
        for (int k = i + my_rank; k < nSize; k += rank_count)
        {
            if (k == i)
                continue;
 
            tmp = Matrix[k][i];
            for (int j = i; j <= nSize; j++)
                Matrix[k][j] -= tmp * r[j];
        }
 
        // увеличиваем счётчик ожидающих потоков
        i_counter ++;
 
        if (i_counter < rank_count)
        {
            // если счётчик не достиг количества потоков
            unique_lock<mutex> lk2(cv2_m);
 
            cout_m.lock();
            cout << "rank = " << my_rank <<  ", wait for finish " << i_row << ", i_counter = " << i_counter << endl;
            cout_m.unlock();
 
            // засыпаем пока последний поток не сбросит счётчик
            cv2.wait(lk2, []{ return i_counter == 0; });
 
            cout_m.lock();
            cout << "rank = " << my_rank <<  ", finished, i_row = " << i_row << endl;
            cout_m.unlock();
        }
        else
        {
            // если счётчик достиг количества потоков
            lock_guard<mutex> lk(cv2_m);
            
            cout_m.lock();
            cout << "rank = " << my_rank <<  ", set i_counter = 0, i_row = " << i_row << endl;
            cout_m.unlock();
 
            // не засыпаем, сбрасываем счётчик в ноль и уведомляем остальные потоки
            i_counter = 0;
            cv2.notify_all();
        }
    }
 
    // обратный проход, только для одного потока
    if (my_rank == 0)
    {
        for (int i = nSize - 1; i >= 0; i--)
        {
            tmp = Matrix[i][nSize];
            for (int j = i + 1; j < nSize; j++)
                tmp -= Matrix[i][j] * x[j];
            x[i] = tmp / Matrix[i][i];
        }
    }
}
 
int main()
{
    setlocale(LC_ALL, "ru");
 
    int Asize;
    //cout << "Размер матрицы: ";
    //cin >> Asize;
    Asize = 10;
 
    double **A_arr = new double*[Asize];
    for (int i = 0; i < Asize; i++)
        A_arr[i] = new double[Asize + 1];
 
    FillMatrix(A_arr, Asize);
    
    double *X_arr_Thr = new double [Asize];
    double *cur_row = new double [Asize + 1];
 
    thread th1(GaussThread, A_arr, Asize, X_arr_Thr, 0, 4, cur_row);
    thread th2(GaussThread, A_arr, Asize, X_arr_Thr, 1, 4, cur_row);
    thread th3(GaussThread, A_arr, Asize, X_arr_Thr, 2, 4, cur_row);
    thread th4(GaussThread, A_arr, Asize, X_arr_Thr, 3, 4, cur_row);
    th1.join();
    th2.join();
    th3.join();
    th4.join();
    
    for (int i = 0; i < Asize; i++)
        delete A_arr[i];
    delete[] A_arr;
    delete[] X_arr_Thr;
    delete[] cur_row;
 
    return 0;
}