#include<iostream>
#include <vector>
#include <cmath>
#include "mpi.h"
using namespace std;
int main(int argc, char** argv) {
MPI_Init(&argc, &argv);
    int rank, size;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
    // coefficient matrix
    vector<vector<double>> C = {
        {0.05, -0.06, -0.12, 0.14},
        {0.04, -0.12, 0.08, 0.11},
        {0.34, 0.08, -0.06, 0.14},
        {0.11, 0.12, 0, -0.03}
    };
    // constant vector
    vector<double> d = {-2.17, 1.4, -2.1, -0.8};
    // initial approximation (vector of zeros)
    vector<double> x = {0, 0, 0, 0};
    // desired accuracy
    double eps = 1e-6;
    // iteration counter
    int k = 0;
    // get the number of rows per process
    int rows_per_proc = 4 / size;
    // scatter the coefficient matrix and constant vector
    vector<double>local_C(rows_per_proc * 4);
    vector<double>local_d(rows_per_proc);
MPI_Scatter(C.data(), rows_per_proc * 4, MPI_DOUBLE, local_C.data(), rows_per_proc * 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Scatter(d.data(), rows_per_proc, MPI_DOUBLE, local_d.data(), rows_per_proc, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    // iteration process
    while (true) {
        // compute new approximation locally
        vector<double>x_new(rows_per_proc);
        for (int i = 0; i<rows_per_proc; i++) {
            double s = 0;
            for (int j = 0; j < 4; j++) {
                s += local_C[i * 4 + j] * x[j];
            }
x_new[i] = s + local_d[i];
        }
        // gather new approximation from all processes
        vector<double>global_x(4);
MPI_Allgather(x_new.data(), rows_per_proc, MPI_DOUBLE, global_x.data(), rows_per_proc, MPI_DOUBLE, MPI_COMM_WORLD);
        // check for accuracy on root process
        if (rank == 0) {
            double norm = 0;
            for (int i = 0; i< 4; i++) {
                norm += pow(global_x[i] - x[i], 2);
            }
            if (sqrt(norm) < eps) {
                break;
            }
        }
        // broadcast updated approximation to all processes
MPI_Bcast(global_x.data(), 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
        x = global_x;
        // increment counter on all processes
        k++;   }
    // print solution and number of iterations on root process
    if (rank == 0) {
cout<< "Solution: ";
        for (int i = 0; i< 4; i++) {
cout<< x[i] << " ";}
cout<<endl;
cout<< "Number of iterations: " << k <<endl;   }
MPI_Finalize();
return 0;}