diff --git a/bnb/push_relabel.cpp b/bnb/push_relabel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a5692f8cdfdaf5d954977df467cb387010684593
--- /dev/null
+++ b/bnb/push_relabel.cpp
@@ -0,0 +1,120 @@
+#include <limits>
+#include <gp-bnb/push_relabel.hpp>
+
+
+push_relabel::push_relabel(const graph &g, std::vector<node_id> sources, std::vector<node_id> sinks) 
+    : g_(g), sources_(sources), sinks_(sinks) {
+
+    int num_nodes = g_.num_nodes();
+
+    // Initialize all vectors
+    labels_.assign(num_nodes + 1, 0);
+    flow_.assign(num_nodes + 1, std::vector<int>(num_nodes + 1, 0));
+    excess_.assign(num_nodes + 1, 0);
+    sources_and_sinks_.assign(num_nodes + 1, 0);
+
+    // Sources marked as 1, sinks marked as -1. All the other nodes are 0.
+    for (auto s : sources_)
+        sources_and_sinks_[s] = 1;
+    for (auto t : sinks_) 
+        sources_and_sinks_[t] = -1;
+
+};
+
+
+void push_relabel::push(node_id u, node_id v) {
+
+    int d = std::min(excess_[u], 1 - flow_[u][v]);
+    flow_[u][v] += d;
+    flow_[v][u] -= d;
+    excess_[u] -= d;
+    excess_[v] += d;
+
+    if (d && excess_[v] == d) {
+        excess_nodes_.push(v);
+    }
+       
+}
+
+
+void push_relabel::relabel(node_id u) {
+
+    int d = std::numeric_limits<int>::max();
+
+    std::vector<node_id> neighbors = g_.get_adjacency(u);
+
+    for (auto v : neighbors) {
+        if (1 - flow_[u][v] > 0)
+            d = std::min(d, labels_[v]);
+    }
+
+    if (d < std::numeric_limits<int>::max())
+        labels_[u] = d + 1;
+
+}
+
+
+void push_relabel::discharge(node_id u) {
+
+    while (excess_[u] > 0) {
+
+        std::vector<node_id> neighbors = g_.get_adjacency(u);
+
+        bool push_possible = false;
+        for (auto v : neighbors) {
+            if (1 - flow_[u][v] > 0 && labels_[u] > labels_[v]) {
+                push(u, v);
+                push_possible = true;
+            }
+        }
+
+        if (!push_possible)
+             relabel(u);
+    }
+
+}
+
+/* Preflow: Labels for source nodes - |V|.
+Excess for source nodes - inf.
+Push all outgoing edges from source nodes. */
+void push_relabel::preflow() {
+
+    int num_nodes = g_.num_nodes();
+
+    for (node_id s : sources_)  {
+        labels_[s] = num_nodes;
+        excess_[s] = std::numeric_limits<int>::max();
+    }
+
+    for (node_id s : sources_) {
+        std::vector<node_id> neighbors = g_.get_adjacency(s);
+        for (auto n : neighbors) {
+            if (sources_and_sinks_[n] != 1) // if n is not a source node
+                push(s, n);
+        }
+    }
+}
+
+void push_relabel::run() {
+
+    preflow();
+
+    while (!excess_nodes_.empty()) {
+        int u = excess_nodes_.front();
+        excess_nodes_.pop();
+
+        // if u is not a source and not a sink
+        if (sources_and_sinks_[u] == 0)
+            discharge(u);
+    }
+
+    // Flow value is the sum of excesses in sinks
+    flow_value_ = 0;
+    for (node_id t : sinks_)
+        flow_value_ += excess_[t];
+        
+}
+
+int push_relabel::get_max_flow() const {
+    return flow_value_;
+}
\ No newline at end of file
diff --git a/include/gp-bnb/push_relabel.hpp b/include/gp-bnb/push_relabel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f9e88f8fbd5155a15435ac3474b198e601652bf5
--- /dev/null
+++ b/include/gp-bnb/push_relabel.hpp
@@ -0,0 +1,54 @@
+#ifndef PUSHRELABEL_H_
+#define PUSHRELABEL_H_
+
+#include <vector>
+#include <map>
+#include <queue>
+#include <gp-bnb/graph.hpp>
+
+class push_relabel {
+private:
+	const graph &g_;
+
+	std::vector<node_id> sources_;
+	std::vector<node_id> sinks_;
+	std::vector<int> sources_and_sinks_;
+
+    std::vector<int> labels_, excess_;
+    std::vector<std::vector<int>> flow_;
+    std::queue<node_id> excess_nodes_;
+
+    int flow_value_;
+
+    void push(node_id u, node_id v);
+
+	void preflow();
+
+    void relabel(node_id u);
+
+    void discharge(node_id u);
+
+public:
+	/**
+	 * Constructs an instance of the EdmondsKarp algorithm for the given graph, source and sink
+	 * @param graph The graph.
+	 * @param source The source node.
+	 * @param sink The sink node.
+	 */
+	push_relabel(const graph &g, std::vector<node_id> sources, std::vector<node_id> sinks);
+
+	/**
+	 * Computes the maximum flow, executes the EdmondsKarp algorithm.
+	 * For unweighted, undirected Graphs.
+	 */
+	void run();
+
+	/**
+	 * Returns the value of the maximum flow from source to sink.
+	 *
+	 * @return The maximum flow value
+	 */
+	int get_max_flow() const;
+};
+
+#endif /* PUSHRELABEL_H_ */
\ No newline at end of file