The main difference between Augmenting path algorithms and Push-Relabel is that:

• Augmenting path algorithms maintain a flow(capacity constraints and conservation constraints are valid) and augment it some s-t path in each iteration to increase its value. We iterate untill there exists no s-t path in the current residual network.
• Push-Relabel algorithm takes a different approach, it works with pre-flow(conservation constraint is violated — the amount of flow into a vertex can exceed the amount of flow exiting it, the difference(flow_in-flow_out) in the amount present at that vertex is called excess, excess >= 0). This algorithm starts with a pre flow and terminates with an actual flow, eliminates all the excesses. We start with s-t disconnected and work towards a feasible flow.

To eliminate excess at a vertex it has the following function:

Push(v):
    choose an outgoing edge (v, w) of v in Residual graph G_f (if any)
    // push as much flow as possible
    let fl = min{excess(v), capacilty of edge(v, w)}
    push fl units of flow along (v, w)

Choose an edge(what type of edge will be clear later) and push as much as we can.

We maintain heights to route flow from source to sink, otherwise, if we have a cycle in the network we might end up pushing around it indefinitely. To route the excess flow systematically we maintain these heights, heights are always non-negative, and thus we can visualize edges going uphill, downhill or at the same level. We maintain these invariants, let's denote height by h then h(src)=|V|, h(sink)=0 and for every edge (v, w) of the current residual network G_f (with positive residual capacity), h(v) <= h(w) + 1, the third one says that if we have to go downhill through the edges we only go down 1 step at a time(height doesn't go down quickly).

The motivation for invariants: We maintain that s-t is disconnected as h(s)=|V| and h(t)=0 if there is a path then it is of length atmost |V|-1.

We initialize the network with the following:

h(s) = |V|
h(v) = 0 for V-{s}
flow = capacity for (s,*) // edges going out of source 
flow = 0 for E-{(s,*)}

The Pre-flow for edges going out of s is at full capacity so as to maintain the third invariant so that we don't have edges going out of s with positive residual capacity.

We restrict the above Push() operation — flow is only allowed to be pushed downhill in the residual network. This restriction allows us to maintain invariants, if we push uphill then we get reverse arcs in the residual network which violate h(v) <= h(w) + 1.

The main iteration in this algorithm is to choose a vertex v in V-{s,t} with excess[v]>0, if there are many such vertices, choose any vertex with maximum height. Now if we have an edge going out of v in G_f with h(v)=h(w)+1, then Push() over this edge (v,w).

While there exists v in V-{s,t} with excess[v]>0:
    Choose any vertex v with a maximum height
    if there is an outgoing edge (v, w) of v in Gf with h(v) = h(w) + 1
        Push(v)
    else
        Relabel(v) 

The idea behind choosing maximum height is to select one furthest away from the sink t. The hope is that working on those nodes will move the flow towards the sink, and aggregate them along the way.

Here is the Relabel(v) algorithm

Relabel(v):
    h[v]++;

One can observe that the same vertex will be relabelled until we find a vertex w such that h(v)=h(w)+1 so we can rewrite the Relabel(v) function as follows:

Relabel(v):
    d=INF
    for all edges (v,w) with positive residual capacity
        d=min(d,h[w])
    h[v]=d+1

About the running time $$$O(n^3)$$$:

• If the vertex v has positive excess in the pre-flow f, then there is a path v->s in the residual network G_f. — Clearly, the excess must have come from the source only so there should be a path to s.

• In the push-relabel algorithm, every vertex always has height at most 2|V|. — As the height of the source is |V|, and there is a path from v to s and height invariant. Therefore, the push-relabel algorithm performs $$$O(n^2)$$$ relabels

• Now if we think when the Push will flow will be bounded by Capacity(Saturating pushes) and when it will be bounded by Excess(Non-saturating pushes), One can prove that Between two saturating pushes on the same edge (v, w) in the same direction, each of v, w is relabeled at least twice. and Between any two relabel operations, there are at most n non-saturating pushes.

Here is implementation from CP-algorithms

Here is an implementation from Stanford Notebook with Gap Heuristic.

[UPD -- Gap Heuristic]

Let g be an integer and 0 < g < n. Suppose at a certain stage of the algorithm there are no nodes with height/distance label g but there are nodes v with g < h(v) < n. Then the sink is not reachable from any of these nodes. Therefore, the labels of such nodes may be increased to n. (Note that these nodes will never be active.) If for every i we maintain linked lists of nodes with the distance label i, the overhead of detecting the gap is very small. Most work done by the gap relabeling heuristic is useful: it involves processing the nodes determined to be disconnected from the sink.

You can take a look at my implementation tested with SPOJ/MTOTALF and Kattis/maxflow for better understanding.

 #include <bits/stdc++.h>
 
using namespace std;
 
typedef long long ll;
 
struct edge{
    int from,to; ll cap,f;
};
 
vector<edge> ed;
vector<vector<int> > adj,hs;
vector<int> h;vector<ll> ex;
int N,S,T;
 
void init(int n,int s,int t){
    N=n;S=s;T=t; // S - source, T - sink
    h=vector<int>(N); // height
    ex=vector<ll>(N); // excess
    adj=vector<vector<int>>(N);
    hs=vector<vector<int>>(2*N+1); // for finding max height with excess
}
 
void add_edge(int from,int to,int cap){
    if(from==to) return;
    adj[from].push_back(ed.size());
    ed.push_back({from,to,cap,0});
    adj[to].push_back(ed.size());
    ed.push_back({to,from,0,0});
}
 
void push(int id){
    int v=ed[id].from,to=ed[id].to;
    ll fl=min(ex[v],ed[id].cap-ed[id].f);
    if(to!=S && to!=T && !ex[to] && fl) hs[h[to]].push_back(to);
    ed[id].f+=fl;
    ed[id^1].f-=fl;
    ex[v]-=fl;
    ex[to]+=fl;
}
 
int relabel(int v){
    h[v]=2*N; // Cannot exceed this value
    for(int id:adj[v])
        if(ed[id].cap>ed[id].f)
            h[v]=min(h[v],h[ed[id].to]+1);
    hs[h[v]].push_back(v);
    return h[v];
}
 
ll max_flow(){
    // initialization and preflow
    h[S]=N;
    for(int id:adj[S]){
        ex[S]+=ed[id].cap;
        push(id);
    }
 
    if(hs[0].size())
    for(int hi=0;hi>=0;){
        // find a vertex v with ex[v]>0 && max height
        int v=hs[hi].back();hs[hi].pop_back();
        while(ex[v] >0){
            for(int i=0;i<(int)adj[v].size() && ex[v];i++){
                int id=adj[v][i];edge e=ed[id];
                if(e.cap-e.f>0 && h[v]==h[e.to]+1)
                    push(id);
            }
            if(ex[v]) hi=relabel(v);
        }
        while(hi>=0 && hs[hi].empty()) --hi;
    }
    // Calculate flow
    ll res=0;
    for(int id:adj[S]) res+=ed[id].f;
    return res;
}

int main(){
    // blah blah
    init(n,s,t);
    ll flow=max_flow();
    return 0;
}

Here is a submission from zscoder

If you have better implementations(with some heuristics), you can share in comments. Learn and Let Learn.