OpenMP Lab3

#include <vector>
#include <iostream>
#include <omp.h>
#include <algorithm>
#include <chrono>

#define eps 0.00001

long double __denumerator_squared_st_diviation = 0;

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

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());
#pragma omp parallel for default(none) shared(lhs, rhs, result)
    for (uint64_t idx = 0; idx < lhs.size(); idx++)
        result[idx] = lhs[idx] - rhs[idx];
    return result;
}
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);
#pragma omp parallel for default(none) shared(lhs, rhs, result)
    for (std::uint64_t i = 0; i < result.size(); i++)
        for (std::uint64_t j = 0; j < lhs[i].size(); j++)
            result[i] += lhs[i][j] * rhs[j];

    return result;
}

void initialize_SLE(Matrix& mt, Vector& rhs, Vector& solution){
    const auto N = mt.size();
    for(std::uint64_t i = 0; i < mt.size(); i++)
        for(std::uint64_t j = 0; j < mt[i].size(); j++)
            mt[i][j] = (i == j) ? 2 : 1;
    for(auto& elem: rhs)
        elem = N + 1;
    __denumerator_squared_st_diviation = rhs.size()*(N+1)*(N+1);
}

long double count_metric(Matrix& mt, Vector& rhs, Vector& solution){
    long double numerator_squared_st_diviation = 0;
    const auto numerator = vector_subtraction(matrix_multiply(mt, solution), rhs);
#pragma omp parallel for reduction(+:numerator_squared_st_diviation) shared(numerator) default(none)
    for(std::uint64_t i = 0; i < numerator.size(); i++)
        numerator_squared_st_diviation += numerator[i] * numerator[i];
    return numerator_squared_st_diviation/__denumerator_squared_st_diviation;
}

void evaluate_iteration(Matrix& mt, Vector& rhs, Vector& solution){
    const double lambda = (1.f/(10.f * static_cast<float>(mt.size())));
    auto vec = vector_subtraction(matrix_multiply(mt, solution), rhs);
#pragma omp parallel for default(none) shared(lambda, vec)
    for(std::uint64_t i = 0; i < vec.size(); i++)
        vec[i] *= lambda;
#pragma omp parallel for default(none) shared(lambda, vec, solution)
    for(std::uint64_t i = 0; i < solution.size(); i++)
        solution[i] -= vec[i];
}

int main(int argc, char** argv){
    std::uint64_t N = 0;
    std::int32_t threads = 0;
    if(argc == 3){
        N = std::stoul(argv[1]);
        threads = std::stoi(argv[2]);
    }else{
        std::cerr << "Invalid number of arguments: [matrix size] [thread number]";
        exit(1);
    }
    omp_set_num_threads(threads);
    Matrix SLE_coef(N, Vector(N, 0));
    Vector SLE_rhs(N, 0), SLE_sol(N, 0);
    initialize_SLE(SLE_coef, SLE_rhs, SLE_sol);
    long double metric = count_metric(SLE_coef, SLE_rhs, SLE_sol);
    //std::cout << "Iteration #0. metric: " << metric << ", duration: 0s\n";
    uint32_t idx = 1;
    const auto start = std::chrono::system_clock::now();
    while(metric >= eps){
        evaluate_iteration(SLE_coef, SLE_rhs, SLE_sol);
        metric = count_metric(SLE_coef, SLE_rhs, SLE_sol);
        //std::cout << "Iteration #" << (idx++) << ". metric: " << metric << "ms\n";
    }
    const 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_sol.begin(), SLE_sol.end(), [](const double value) { std::cout << value << ", "; });
    std::cout << "\n}\n";
    return 0;
}