// solve in parallel Ax = b, there A is a coefficientsMatrix and b is a freeMembersVector
std::vector<double> solveParallelSeidel(matrix &a,
                                        std::vector<double> &b,
                                        double eps,
                                        int approximationsCount) {
    int size;
    int rank;
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Status status;

    int    root             = 0;
    int    tag              = 0;
    int    procNumber       = 0;
    int    indexInProc      = 0;
    int elementsCounter      = 0;
    int    n                = static_cast<int>(b.size());
    int    elementsPerBlock = n / size;
    int    remainElements   = n % size + elementsPerBlock;
    int    elementsCount    = (rank ==  size - 1 ? remainElements : elementsPerBlock);
    double currSum          = 0.;
    double globalSum        = 0.;
    double oldNorm          = 0.;
    double newNorm          = 0.;
    double currEps          = 0.;
    char endOfProgram       = 0;

    std::vector<double> resultX;
    std::vector<double> x(elementsCount);
    matrix              currA;


    // Init data on all Proc
    if (rank == root) {
        resultX.resize(n);
        for (int i = 1; i < size; i++) {
            elementsCounter = (i ==  size - 1 ? remainElements : elementsPerBlock);
            MPI_Send(a[i * elementsPerBlock], elementsCounter * n, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
        }

    } else {
        currA.init(n, elementsCount);
        MPI_Recv (currA[0], elementsCount * n, MPI_DOUBLE, root, tag, MPI_COMM_WORLD, &status);
    }

    //main process
    do {
        newNorm = 0.;

        for (int i = 0; i < n; i++) {
            currSum = 0;
            procNumber = i / elementsPerBlock;
            if (procNumber == size) procNumber--;
            indexInProc = i - procNumber * elementsPerBlock;

            for (int j = 0; j < elementsCount; j++) {
                currSum -= currA[i][j] * x[j];
            }
            //MPI_Reduce(&currSum, &globalSum, 1, MPI_DOUBLE, MPI_SUM, root, MPI_COMM_WORLD);

            if (rank == root) {
                newNorm += globalSum * globalSum;
                resultX[i] = (globalSum + b[i]) / a[i][i];
                globalSum = 0.;
                if (procNumber != root) {
                    MPI_Send (&resultX[i], 1, MPI_DOUBLE, procNumber, tag, MPI_COMM_WORLD);
                } else {
                    x[indexInProc] = resultX[i];
                }
            } else if (rank == procNumber) {
                MPI_Recv (&x[indexInProc], 1, MPI_DOUBLE, root, tag, MPI_COMM_WORLD, &status);
            }
        }
        if (rank == root) {
            std::cerr << std::endl;
            currEps = abs(sqrt(newNorm) - sqrt(oldNorm));
            approximationsCount--;
            oldNorm = newNorm;
            if (currEps < eps || !approximationsCount) endOfProgram = 1;
            MPI_Bcast (&endOfProgram, 1, MPI_BYTE, root, MPI_COMM_WORLD);
        }
    } while (!endOfProgram);

    return resultX;