gauss(other, finished)

// Copyright 2023 Soloninko Andrey

#include <iostream>
#include <random>
#include <vector>
#include "task_2/soloninko_a_gauss_horizontal/Gauss.h"


int get_slice(int size, int proc_count, int rank) {
    int q = size / proc_count;
    if (size % proc_count)
        q++;
    int r = proc_count * q - size;

    int slice = q;
    if (rank >= proc_count - r)
        slice = q - 1;
    return slice;
}

bool check_res(std::vector<double>& matrix, std::vector<double>& x, int size) {
    for (int i = 0; i < size; i++) {
        double sum = 0;
        for (int j = 0; j < size; j++)
            sum += matrix[i * (size + 1) + j] * x[j];

        if (std::abs(sum - matrix[i * (size + 1) + size]) > 0.0001)
            return false;
    }
    return true;
}

std::vector<double> rand_matr(int size) {
    std::random_device rd;
    std::uniform_int_distribution<int> unif(-100, 100);
    std::vector<double> matrix((size + 1) * size);
    int i, j;
    for (i = 0; i < size; i++)
        for (j = 0; j < size + 1; j++)
            matrix[i * (size + 1) + j] = unif(rd);

    return matrix;
}

std::vector<double> Gauss(std::vector<double>& matrix, int size) {
    int rank, comm_size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &comm_size);

    std::vector<double> Xi(size, 0.0);  // Инициализация вектора нулями

    int nrows = get_slice(size, comm_size, rank);
    std::vector<int> rows(nrows);
    std::vector<double> tmp(size + 1);

    for (int i = 0; i < nrows; i++)
        rows[i] = rank + comm_size * i;

    MPI_Bcast(matrix.data(), size * (size + 1), MPI_DOUBLE, 0, MPI_COMM_WORLD);

    int row = 0;
    for (int i = 0; i < size - 1; i++) {
        if (i == rows[row]) {
            MPI_Bcast(&matrix.data()[rows[row] * (size + 1)], size + 1, MPI_DOUBLE, rank, MPI_COMM_WORLD);
            for (int j = 0; j <= size; j++)
                tmp[j] = matrix[rows[row] * (size + 1) + j];
            row++;
        } else {
            MPI_Bcast(tmp.data(), size + 1, MPI_DOUBLE, i % comm_size, MPI_COMM_WORLD);
        }

        for (int j = row; j < nrows; j++) {
            double scaling = matrix[rows[j] * (size + 1) + i] / tmp[i];
            for (int k = i; k < size + 1; k++)
                matrix[rows[j] * (size + 1) + k] -= scaling * tmp[k];
        }
    }

    row = 0;
    for (int i = 0; i < size; i++) {
        if (i == rows[row]) {
            Xi[i] = matrix[i * (size + 1) + size];
            row++;
        }
    }

    row = nrows - 1;
    for (int i = size - 1; i > 0; i--) {
        if (row >= 0) {
            if (i == rows[row]) {
                Xi[i] /= matrix[i * (size + 1) + i];
                MPI_Bcast(&(Xi.data())[i], 1, MPI_DOUBLE, rank, MPI_COMM_WORLD);
                row--;
            } else {
                MPI_Bcast(&(Xi.data())[i], 1, MPI_DOUBLE, i % comm_size, MPI_COMM_WORLD);
            }
        } else {
            MPI_Bcast(&(Xi.data())[i], 1, MPI_DOUBLE, i % comm_size, MPI_COMM_WORLD);
        }

        for (int j = 0; j <= row; j++)
            Xi[rows[j]] -= matrix[rows[j] * (size + 1) + i] * Xi[i];
    }

    if (rank == 0)
        Xi[0] /= matrix[rows[row] * (size + 1)];

    return Xi;
}