void GetBorders(const Matrix& m,
    V& border_x0,
    V& border_x1,
    V& border_y0,
    V& border_y1
    ) {
        for (int j = 0; j < m.M; ++j) {
            border_x0[j] = m.cat(0, j);
        }

        for (int j = 0; j < m.M; ++j) {
            border_x1[j] = m.cat(m.N - 1, j);
        }

        for (int i = 0; i < m.N; ++i) {
            border_y0[i] = m.cat(i, 0);
        }

        for (int i = 0; i < m.N; ++i) {
            border_y1[i] = m.cat(i, m.M - 1);
        }
}


void SyncBorder(V& s, V& r, int other_i, int other_j) {

    if (other_i < 0 || other_i >= lattice_n) {
        r.resize(r.size(), 0.0);
        return;
    }

    if (other_j < 0 || other_j >= lattice_m) {
        r.resize(r.size(), 0.0);
        return;
    }



    int me = world_rank;
    int dst = ComputeLatticeId(other_i, other_j);


    MPI_Send(s.data(), s.size(), MPI_DOUBLE, dst, 0, MPI_COMM_WORLD);
    MPI_Status status;

        MPI_Recv(r.data(), r.size(), MPI_DOUBLE, dst, 0, MPI_COMM_WORLD, &status);

}

double SyncDouble(double value) {
    double result;
    MPI_Allreduce(&value, &result, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
    return result;
}


Matrix Solve(LocalOperator &op, int max_iter, double eps=1e-6) {
    Matrix w(op.block_h, op.block_w, 0.0);

    V border_x0_s(op.block_w, 0);
    V border_x1_s(op.block_w, 0);

    V border_y0_s(op.block_h, 0);
    V border_y1_s(op.block_h, 0);

    V border_x0_r(op.block_w, 0);
    V border_x1_r(op.block_w, 0);

    V border_y0_r(op.block_h, 0);
    V border_y1_r(op.block_h, 0);

    for (int iter = 0; iter < max_iter; ++iter) {
        GetBorders(w, border_x0_s, border_x1_s, border_y0_s, border_y1_s);

        SyncBorder(border_x0_s, border_x0_r, my_i - 1, my_j);
        SyncBorder(border_x1_s, border_x1_r, my_i + 1, my_j);
        SyncBorder(border_y0_s, border_y0_r, my_i, my_j - 1);
        SyncBorder(border_y1_s, border_y1_r, my_i, my_j + 1);


        Matrix Aw = op.Call(w, border_x0_r, border_x1_r, border_y0_r, border_y1_r);



        Matrix r = Aw - op.rhs;
        GetBorders(r, border_x0_s, border_x1_s, border_y0_s, border_y1_s);

        SyncBorder(border_x0_s, border_x0_r, my_i - 1, my_j);
        SyncBorder(border_x1_s, border_x1_r, my_i + 1, my_j);
        SyncBorder(border_y0_s, border_y0_r, my_i, my_j - 1);
        SyncBorder(border_y1_s, border_y1_r, my_i, my_j + 1);
        Matrix Ar = op.Call(r, border_x0_r, border_x1_r, border_y0_r, border_y1_r);




        double rAr = SyncDouble(op.Dot(r, Ar));
        double ArAr = SyncDouble(op.Dot(Ar, Ar));
        double rr = SyncDouble(op.Dot(r, r));

        double tau = rAr/ArAr;


        w = w - r * tau;
        double d = sqrt(rr) * tau;
        double delta = SyncDouble(op.Norm2(w - op.phi));


        if (d  < eps) {
            break;
        }



    }
    MPI_Barrier(MPI_COMM_WORLD);

    return w;
}