ParallelCompLab#2

#include <iostream>
#include <memory>
#include <cstring>
#include <vector>
#include <algorithm>
#include <cmath>
#include <chrono>
#include "mpi.h"

using namespace std;

#define N 4096
#define eps 0.00001
#define lambda (1.f/(10.f*N))


class Process {
    int32_t _process_id;

public:
    Process() noexcept { MPI_Comm_rank(MPI_COMM_WORLD, (int*) &_process_id); }

    Process(const int32_t process_id) noexcept { _process_id = process_id; }

    virtual void initialize() = 0;

    virtual void run() = 0;

    virtual void finalize() = 0;

protected:
    typedef std::vector<std::vector<double>> Matrix;
    typedef std::vector<double> Vector;

    enum Command : unsigned char {
        EVALUATE_METRIC = 0,
        EVALUATE_ITERATION = 1,
        EXIT = 2
    };

    inline int32_t process_id() const noexcept { return _process_id; }

    template<typename T1>
    static void send_value(T1& a, int receiver) {
        constexpr std::uint32_t a_size = sizeof(T1);
        char buffer[a_size];
        memcpy(buffer, &a, a_size);
        MPI_Send(buffer, a_size, MPI_BYTE, receiver, 0, MPI_COMM_WORLD);
    }

    template<typename T1>
    static MPI_Status recv_value(T1& a, int sender) {
        MPI_Status result;
        constexpr std::uint32_t a_size = sizeof(T1);
        char buffer[a_size];
        MPI_Recv(buffer, a_size, MPI_BYTE, sender, 0, MPI_COMM_WORLD, &result);
        memcpy(&a, buffer, a_size);
        return result;
    }

    template<typename T1, typename T2>
    static void send_pair(const T1& a, const T2& b, int receiver) {
        constexpr std::uint32_t a_size = sizeof(T1), b_size = sizeof(T2);
        char buffer[a_size + b_size];
        memcpy(buffer, &a, a_size);
        memcpy(buffer + a_size, &b, b_size);
        MPI_Send(buffer, a_size + b_size, MPI_BYTE, receiver, 0, MPI_COMM_WORLD);
    }

    template<typename T1, typename T2>
    static MPI_Status recv_pair(T1& a, T2& b, int sender) {
        MPI_Status result;
        constexpr std::uint32_t a_size = sizeof(T1), b_size = sizeof(T2);
        char buffer[a_size + b_size];
        MPI_Recv(buffer, a_size + b_size, MPI_BYTE, sender, 0, MPI_COMM_WORLD, &result);
        memcpy(&a, buffer, a_size);
        memcpy(&b, buffer + a_size, b_size);
        return result;
    }

    template<typename ForwardIt>
    static void send_vector(ForwardIt begin, ForwardIt end, int receiver) {
        const uint32_t n = std::distance(begin, end);
        MPI_Send(&n, 1, MPI_UINT32_T, receiver, 0, MPI_COMM_WORLD);
        MPI_Send(begin.base(), n, MPI_DOUBLE, receiver, 0, MPI_COMM_WORLD);
    }

    template<typename ForwardIt>
    static MPI_Status recv_vector(ForwardIt begin, ForwardIt end, int sender) {
        MPI_Status result;
        const uint32_t n = std::distance(begin, end);
        MPI_Recv(begin.base(), n, MPI_DOUBLE, sender, 0, MPI_COMM_WORLD, &result);
        return result;
    }

    template<typename ForwardIt>
    static void send_matrix(ForwardIt begin, ForwardIt end, int receiver) {
        const uint32_t n = std::distance(begin, end), m = begin->size();
        send_pair(n, m, receiver);
        for (ForwardIt iter = begin; iter != end; iter++)
            MPI_Send(iter->data(), m, MPI_DOUBLE, receiver, 0, MPI_COMM_WORLD);
    }

    template<typename ForwardIt>
    static MPI_Status recv_matrix(ForwardIt begin, ForwardIt end, int sender) {
        MPI_Status result;
        const uint32_t n = std::distance(begin, end), m = begin->size();
        for (ForwardIt iter = begin; iter != end; iter++)
            MPI_Recv(iter->data(), m, MPI_DOUBLE, sender, 0, MPI_COMM_WORLD, &result);
        return result;
    }

    static Vector matrix_multiply(const Matrix& lhs, const Vector& rhs) {
        if (lhs[0].size() != rhs.size()) throw std::invalid_argument("can't multiply matrix and vector");
        Vector result(lhs.size(), 0);
        for (std::uint32_t i = 0; i < result.size(); i++)
            for (std::uint32_t j = 0; j < lhs[i].size(); j++)
                result[i] += lhs[i][j] * rhs[j];
        return result;
    }

    static Vector vector_subtraction(const Vector& lhs, const Vector& rhs) {
        if (lhs.size() != rhs.size()) throw std::invalid_argument("can't subtract vector and vector");
        Vector result(lhs.size());
        for (uint32_t idx = 0; idx < lhs.size(); idx++)
            result[idx] = lhs[idx] - rhs[idx];
        return result;
    }
};

class RootProcess : public Process {
    int _proccesses_count;
    Matrix SLE_coef;
    Vector SLE_rhs;
    Vector SLE_solution;
    long double __denumerator_squared_st_diviation = 0;
public:
    RootProcess() noexcept: Process(0), SLE_coef(N, Vector(N, 0)),
                            SLE_rhs(N, 0), SLE_solution(N, 0) {
        MPI_Comm_size(MPI_COMM_WORLD, (int*) &_proccesses_count);
    }

    void initialize() final {
        generate_SLE(SLE_coef, SLE_rhs);
        std::for_each(SLE_rhs.begin(), SLE_rhs.end(),
                      [this](const double value) { __denumerator_squared_st_diviation += value * value; });
        int slice = N / (_proccesses_count - 1);
        auto begin = SLE_coef.begin();
        auto end = min(begin + slice, SLE_coef.end());
        for (int32_t process_id = 1; process_id < _proccesses_count && begin != SLE_coef.end(); process_id++) {
            send_matrix(begin, end, process_id);
            std::uint32_t begin_idx = std::distance(SLE_coef.begin(), begin), end_idx = std::distance(SLE_coef.begin(), end);
            send_vector(SLE_rhs.begin() + begin_idx,
                        SLE_rhs.begin() + end_idx,
                        process_id);
            send_vector(SLE_solution.begin(), SLE_solution.end(), process_id);
            send_pair(begin_idx, end_idx, process_id);
            begin = end;
            end = min(end + slice, SLE_coef.end());
        }
    }

    void run() final {
        long double metric = count_metric();
        //std::cout << "Iteration #0. metric: " << metric << ", duration: 0s\n";
        uint32_t idx = 1;
        auto start = std::chrono::system_clock::now();
        while (metric > eps) {
            auto cycle_start = std::chrono::system_clock::now();
            evaluate_iteration();
            metric = count_metric();
            auto cycle_duration = std::chrono::duration_cast<std::chrono::milliseconds>(
                                                   std::chrono::system_clock::now() - cycle_start);
            //std::cout << "Iteration #" << (idx++) << ". metric: " << metric << ", duration: " << cycle_duration.count() << "ms\n";
        }
        auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now() - start);
        std::cout << "Elapsed Time: " << duration.count() << "ms. Solution: {\n";
        std::for_each(SLE_solution.begin(), SLE_solution.end(), [](const double value) { std::cout << value << ", "; });
        std::cout << "\n}\n";
    }

    void finalize() final {
        for (int32_t process_id = 1; process_id < _proccesses_count; process_id++)
            send_command(Command::EXIT, process_id);
    }

private:
    long double count_metric() const noexcept {
        for (int32_t process_id = 1; process_id < _proccesses_count; process_id++){
            send_command(Command::EVALUATE_METRIC, process_id);
            send_vector(SLE_solution.begin(), SLE_solution.end(), process_id);
        }
        long double numerator = 0;
        for (int32_t process_id = 1; process_id < _proccesses_count; process_id++) {
            long double received_numerator;
            recv_value(received_numerator, process_id);
            numerator += received_numerator;
        }
        return std::sqrt(numerator) / std::sqrt(__denumerator_squared_st_diviation);
    }

    void evaluate_iteration() {
        for (int32_t process_id = 1; process_id < _proccesses_count; process_id++)
            send_command(Command::EVALUATE_ITERATION, process_id);
        std::uint32_t solution_idx = 0;
        for (int32_t process_id = 1; process_id < _proccesses_count; process_id++) {
            std::uint32_t n;
            recv_value(n, process_id);
            Vector received_solution(n, 0);
            recv_vector(received_solution.begin(), received_solution.end(), process_id);
            auto slice = N / (_proccesses_count - 1);
            for (std::uint32_t idx = slice * (process_id-1); idx < slice*process_id; idx++, solution_idx++)
                SLE_solution[solution_idx] = received_solution[idx];
        }
    }

    static void generate_SLE(Matrix& mt, Vector& rhs) noexcept {
        for (uint64_t i = 0; i < mt.size(); i++)
            for (uint64_t j = 0; j < mt[i].size(); j++)
                mt[i][j] = (i == j ? 2 : 1);
        for (auto& elem: rhs)
            elem = N + 1;
    }

    static void send_command(const Command command, const int receiver) noexcept {
        MPI_Send((uint8_t*) &command, 1, MPI_UINT8_T, receiver, 0, MPI_COMM_WORLD);
    }
};

class ChildProcess : public Process {
    Matrix SLE_coef;
    Vector SLE_rhs;
    Vector SLE_solution;
    std::pair<uint32_t, uint32_t> solution_range;
public:
    void initialize() final {
        uint32_t n, m;
        recv_pair(n, m, 0);
        SLE_coef.resize(n, Vector(m, 0));
        recv_matrix(SLE_coef.begin(), SLE_coef.end(), 0);
        recv_value(n, 0);
        SLE_rhs.resize(n, 0);
        recv_vector(SLE_rhs.begin(), SLE_rhs.end(), 0);
        recv_value(n, 0);
        SLE_solution.resize(n, 0);
        recv_vector(SLE_solution.begin(), SLE_solution.end(), 0);
        recv_pair(solution_range.first, solution_range.second, 0);
    }

    void run() final {
        bool is_cancellation_requested{};
        while (!is_cancellation_requested) {
            Command current_command;
            recv_command(current_command);
            switch (current_command) {
                case EVALUATE_METRIC: {
                    uint32_t n;
                    recv_value(n, 0);
                    recv_vector(SLE_solution.begin(), SLE_solution.end(), 0);
                    const auto numerator = vector_subtraction(matrix_multiply(SLE_coef, SLE_solution), SLE_rhs);
                    long double numerator_squared_st_diviation = 0, denumerator_squared_st_diviation = 0;
                    std::for_each(numerator.begin(), numerator.end(),
                                  [&numerator_squared_st_diviation](const double value) {
                                      numerator_squared_st_diviation += value * value;
                                  });
                    send_value(numerator_squared_st_diviation, 0);
                    break;
                }
                case EVALUATE_ITERATION: {
                    auto vec = vector_subtraction(matrix_multiply(SLE_coef, SLE_solution), SLE_rhs);
                    std::for_each(vec.begin(), vec.end(), [](double& value) { value *= lambda; });
                    for (uint32_t i = solution_range.first; i < solution_range.second; i++)
                        SLE_solution[i] -= vec[i-solution_range.first];
                    send_vector(SLE_solution.begin(), SLE_solution.end(), 0);
                    break;
                }
                case EXIT: {
                    is_cancellation_requested = true;
                    break;
                }
            }
        }
    }

    void finalize() final {

    }

private:
    static MPI_Status recv_command(Command& command) noexcept {
        MPI_Status result;
        MPI_Recv((uint8_t*) &command, 1, MPI_UINT8_T, 0, 0, MPI_COMM_WORLD, &result);
        return result;
    }
};

void resolve_process_type(std::unique_ptr<Process>& ptr) {
    int process_id;
    MPI_Comm_rank(MPI_COMM_WORLD, (int*) &process_id);
    if (process_id == 0)
        ptr = std::move(std::make_unique<RootProcess>());
    else
        ptr = std::move(std::make_unique<ChildProcess>());
}

int main(int argc, char** argv) {
    MPI_Init(&argc, &argv);
    std::unique_ptr<Process> process_ptr;
    resolve_process_type(process_ptr);
    process_ptr->initialize();
    process_ptr->run();
    process_ptr->finalize();
    MPI_Finalize();
    return 0;
}