gauss.cpp

// Copyright 2023 Soloninko Andrey

#include "task_2/soloninko_a_gauss_horizontal/Gauss.h"

#include <random>
#include <vector>

/*
int GetRow(const int &rank, const int &comm_size, const int &size) {
    int segment = size / comm_size;
    if (size % comm_size) {
        segment++;
    }
    int m = comm_size * segment - size;
    int row = segment;
    if (rank >= comm_size - m) {
        row = segment - 1;
    }
    return row;
}
*/

std::vector<double> set_matrix_rand(const int size) {
    std::vector<double> matrix(size * (size + 1));
    std::random_device random_device;
    std::uniform_int_distribution<int> unif(-100, 100);

    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++) {
            matrix[i * (size + 1) + j] = unif(random_device);
        }
    }

    return matrix;
}

bool check_res(const std::vector<double>& matrix, const std::vector<double>& vec,
               const int& size) {
    for (int i = 0; i < size; i++) {
        double row_sum = 0;
        for (int j = 0; j < size; j++) {
            row_sum += matrix[i * (size + 1) + size] * vec[j];
        }
        if (std::abs(row_sum - matrix[i * (size + 1) + size]) > 0.0001) {
            return false;
        }
    }
    return true;
}

std::vector<double> gauss(const std::vector<double>& matrix, int matrix_size) {
    int rank;
    int proc_c;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &proc_c);
    std::vector<double> matrix_g = matrix;
    matrix_g = set_matrix_rand(matrix_size);

    std::vector<double> res(matrix_size);

    int segment = matrix_size / proc_c;
    if (matrix_size % proc_c != 0) {
        segment++;
    }

    std::vector<int> rows(segment);
    std::vector<double> tmp(matrix_size + 1);

    for (int i = 0; i < segment; i++) rows[i] = rank + proc_c * i;

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

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

    row = 0;
    for (int i = 0; i < matrix_size; i++) {
        res[i] = 0;
        if (i == rows[row]) {
            res[i] = matrix_g[i * (matrix_size + 1) + matrix_size];
            row++;
        }
    }

    row = proc_c - 1;
    for (int i = matrix_size - 1; i > 0; i--) {
        if (row >= 0) {
            if (i == rows[row]) {
                res[i] /= matrix_g[i * (matrix_size + 1) + i];
                MPI_Bcast(&(res.data())[i], 1, MPI_DOUBLE, rank,
                          MPI_COMM_WORLD);
                row--;
            } else {
                MPI_Bcast(&(res.data())[i], 1, MPI_DOUBLE, i % proc_c,
                          MPI_COMM_WORLD);
            }
        } else {
            MPI_Bcast(&(res.data())[i], 1, MPI_DOUBLE, i % proc_c,
                      MPI_COMM_WORLD);
        }

        for (int j = 0; j <= row; j++)
            res[rows[j]] -= matrix_g[rows[j] * (matrix_size + 1) + i] * res[i];
    }

    if (rank == 0) res[0] /= matrix_g[rows[row] * (matrix_size + 1)];

    return res;
}