#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <mpi.h>

using namespace std;

int main(int argc, char** argv)
{
    int rank, commSize;

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &commSize);


    double** A, * x, * b, * xnew, eps, min, max, tau = 0.01, norm = 0, norm_b = 0, sum = 0;
    int i, j, n, t = 0;

    cin >> n;
    cin >> eps; // 10^-5

    A = (double**)malloc(n * sizeof(double*)); // матрица указателей на указатели-заполняем матрицу по строчкам, у нас двумерная,поэтому **
    x = (double*)calloc(n, sizeof(double));
    b = (double*)malloc(n * sizeof(double));
    xnew = (double*)malloc(n * sizeof(double));

    for (i = 0; i < n; i++)
        A[i] = (double*)malloc(n * sizeof(double));//выделение памяти для указателей

    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
            cin >> A[i][j];
    for(i = 0; i < n; i++)
        cin >> b[i];

    for (i = 0; i < n; i++)
        sum += b[i] * b[i];
    norm_b = sqrt(sum); // норма b
    double norm_0;

    //

    int *index = new int[commSize];
    int *count = new int[commSize];

    int rest_rows = n;
    int row_number = n / commSize;

    count[0] = row_number * n;
    index[0] = 0;

    for (int i = 1; i < commSize; i++)
    {
        rest_rows -= row_number;
        row_number = rest_rows / (commSize - i);
        count[i] = row_number * n;
        index[i] = index[i - 1] + count[i - 1];
    }


    double *cur_row = new double[count[rank]];

    MPI_Scatterv(A, count, index, MPI_DOUBLE, cur_row, count[rank], MPI_DOUBLE, 0, MPI_COMM_WORLD);


    double *res_part = new double[count[rank]];

    for (int i = 0; i < count[rank] / n; i++)
    {
        res_part[i] = 0;

        for (int j = 0; j < n; j++)
        {
            res_part[i] += cur_row[j + i * n] * x[j];
        }
    }

    for (int i = 0; i < commSize; i++)
    {
        index[i] /= n;
        count[i] /= n;
    }

    MPI_Allgatherv(res_part, count[rank], MPI_DOUBLE, xnew, count, index, MPI_DOUBLE, MPI_COMM_WORLD);


    delete[] index;
    delete[] cur_row;
    delete[] res_part;
    delete[] count;

    //

    for (i = 0; i < n; i++)
        xnew[i] = xnew[i] - b[i];
    for (i = 0; i < n; i++)
        sum += xnew[i] * xnew[i];
    norm_0 = sqrt(sum); // ||Ax_0 - b||
    for (i = 0; i < n; i++)
        xnew[i] = -tau * xnew[i];
    for (i = 0; i < n; i++)
        xnew[i] = x[i] + xnew[i];
    for (i = 0; i < n; i++)
        x[i] = xnew[i];

    sum = 0;
    for(i = 0; i < n; i++)
    {
        xnew[i] = 0;// обнуляем
        for (j = 0; j < n; j++)
            xnew[i] += A[i][j] * x[j];
    }
    for (i = 0; i < n; i++)
        xnew[i] = xnew[i] - b[i];
    for (i = 0; i < n; i++)
        sum += xnew[i] * xnew[i];
    norm = sqrt(sum); // ||Ax_1 - b||

    if (norm_0 <= norm)
        tau = -tau;

    norm /= norm_b;
    while (norm >= eps)
    {
        sum = 0;
        for (i = 0; i < n; i++)
        {
            xnew[i] = 0;// обнуляем
            for (j = 0; j < n; j++)
                xnew[i] += A[i][j] * x[j];
        }
        for (i = 0; i < n; i++)
            xnew[i] = xnew[i] - b[i];
        for (i = 0; i < n; i++)
            sum += xnew[i] * xnew[i];
        norm = sqrt(sum)/norm_b;
        for (i = 0; i < n; i++)
            xnew[i] = -tau * xnew[i];
        for (i = 0; i < n; i++)
            xnew[i] = x[i] + xnew[i];
        for (i = 0; i < n; i++)
            x[i] = xnew[i];
    }

    cout << "Решения уравнений: " << "\n";
    for (i = 0; i < n; i++)
        cout << "x[" << i << "] = " << xnew[i] << "\n";

    MPI_Finalize();
    return 0;
}
//
//3
//0.001
//1 0 0 1
//0 1 0 2
//0 0 1 3