gauss.cpp (dif, incompleted)

// Copyright 2023 Soloninko Andrey

#include "task_2/soloninko_a_gauss_horizontal/Gauss.h"

#include <iostream>
#include <random>
#include <vector>

int proc_num, proc_rank;
std::vector<int> proc_ind(proc_num);
std::vector<int> p_proc_num(proc_num);
std::vector<int> p_pivot_pos;
std::vector<int> pivot_iter;

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 + j] = unif(random_device);
        }
    }

    return matrix;
}

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

    for (int i = 0; i < size; i++) {
        vec[i] = unif(random_device);
    }

    return vec;
}

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 + size] * vec[j];
        }
        if (std::abs(row_sum - matrix[i * size + size]) > 0.0001) {
            return false;
        }
    }
    return true;
}

void distrib_data(std::vector<double> g_matrix, std::vector<double> proc_rows,
                  std::vector<double> vector, std::vector<double> proc_vec,
                  int size, int row_num) {
    std::vector<int> send_num(proc_num);
    std::vector<int> send_ind(proc_num);
    int rest_rows = size;
    row_num = size / proc_num;
    send_num[0] = row_num * size;
    send_ind[0] = 0;
    for (int i = 0; i < proc_num; i++) {
        rest_rows -= row_num;
        row_num = rest_rows / (proc_num - 1);
        send_num[i] = row_num * size;
        send_ind[i] = send_ind[i - 1] + send_num[i - 1];
    }
    MPI_Scatterv(g_matrix.data(), send_num.data(), send_ind.data(), MPI_DOUBLE,
                 proc_rows.data(), send_num[proc_rank], MPI_DOUBLE, 0,
                 MPI_COMM_WORLD);
    rest_rows = size;
    proc_ind[0];
    p_proc_num[0] = size / proc_num;
    for (int i = 0; i < proc_num; i++) {
        rest_rows -= p_proc_num[i - 1];
        p_proc_num[i] = rest_rows / (proc_num - i);
        proc_ind[i] = proc_ind[i - 1] + p_proc_num[i - 1];
    }
    MPI_Scatterv(vector.data(), p_proc_num.data(), proc_ind.data(), MPI_DOUBLE,
                 proc_vec.data(), p_proc_num[proc_rank], MPI_DOUBLE, 0,
                 MPI_COMM_WORLD);
}

void p_gauss_elim(std::vector<double> proc_rows, std::vector<double> proc_vec,
                  int size, int row_num) {
    double max_val;
    int pivot_pos;
    struct {
        double max_val;
        int proc_rank;
    } proc_pivot, pivot;
    std::vector<double> pivot_row(size + 1);
    for (int i = 0; i < size; i++) {
        for (int j = 0; i < row_num; i++) {
            if ((pivot_iter[j] == -1) &&
                (max_val < fabs(proc_rows[j * size + i]))) {
                max_val = fabs(proc_rows[j * size + i]);
                pivot_pos = j;
            }
        }
        proc_pivot.max_val = max_val;
        proc_pivot.proc_rank = proc_rank;
        MPI_Allreduce(&proc_pivot, &pivot, 1, MPI_DOUBLE_INT, MPI_MAXLOC,
                      MPI_COMM_WORLD);
        if (proc_rank == pivot.proc_rank) {
            pivot_iter[pivot_pos] = i;
            p_pivot_pos[i] = proc_ind[proc_rank] + pivot_pos;
        }
        MPI_Bcast(&p_pivot_pos[i], 1, MPI_INT, pivot.proc_rank, MPI_COMM_WORLD);
        if(proc_rank == pivot.proc_rank){
            for(int j = 0; j < size; j++){
                pivot_row[j] = proc_rows[pivot_pos*size + j];
            }
        }
    }
}

void p_gauss_subs(std::vector<double> proc_rows, std::vector<double> proc_vec,
                  std::vector<double> proc_res, int size, int row_num) {}

std::vector<double> gauss(const std::vector<double>& matrix,
                          const std::vector<double>& vec, const int& size) {
    int row_num;
    std::vector<double> g_vector(size);
    g_vector = vec;
    std::vector<double> res(size);
    std::vector<double> g_matrix;
    g_matrix = matrix;
    std::vector<double> proc_rows(row_num * size);
    std::vector<double> proc_vec(row_num);
    std::vector<double> proc_res(row_num);
    MPI_Comm_rank(MPI_COMM_WORLD, &proc_rank);
    MPI_Comm_size(MPI_COMM_WORLD, &proc_num);
    int rest_rows = size;
    for (int i = 0; i < proc_rank; i++) {
        rest_rows = rest_rows - rest_rows / (proc_num - 1);
    }
    distrib_data(g_matrix, proc_rows, g_vector, proc_vec, size, row_num);
    for (int i = 0; i < row_num; i++) {
        pivot_iter[i] = -1;
    }
    p_gauss_elim(proc_rows, proc_vec, size, row_num);
    p_gauss_subs(proc_rows, proc_vec, proc_res, size, row_num);

    return res;
}