diff --git a/networking_flow/README.md b/networking_flow/README.md
new file mode 100644
index 000000000000..265df176512f
--- /dev/null
+++ b/networking_flow/README.md
@@ -0,0 +1,45 @@
+# Networking Flow
+
+This directory collects algorithms for the **maximum-flow problem**: given a
+directed graph whose edges have capacities, a `source`, and a `sink`, how
+much flow can be pushed from `source` to `sink` without exceeding any edge's
+capacity?
+
+Maximum flow turns up all over the place — routing traffic through a network,
+matching people to jobs, scheduling, image segmentation, and any problem that can
+be phrased as "move as much as possible from here to there through a shared
+network." Its close relative, the **minimum cut**, finds the cheapest set of
+edges whose removal disconnects the sink from the source, and the
+[max-flow min-cut theorem](https://en.wikipedia.org/wiki/Max-flow_min-cut_theorem)
+says the two always have the same value.
+
+New to the topic? These are good starting points:
+
+* <https://en.wikipedia.org/wiki/Maximum_flow_problem>
+* <https://en.wikipedia.org/wiki/Flow_network>
+* <https://en.wikipedia.org/wiki/Max-flow_min-cut_theorem>
+
+## What's in this directory
+
+| File | Description |
+| ---- | ----------- |
+| [`ford_fulkerson.py`](ford_fulkerson.py) | The [Ford-Fulkerson](https://en.wikipedia.org/wiki/Ford%E2%80%93Fulkerson_algorithm) method, finding augmenting paths with a breadth-first search (the [Edmonds-Karp](https://en.wikipedia.org/wiki/Edmonds%E2%80%93Karp_algorithm) refinement). Uses an adjacency-matrix representation. Runs in `O(V * E^2)`. |
+| [`minimum_cut.py`](minimum_cut.py) | Finds the edges of a [minimum s-t cut](https://en.wikipedia.org/wiki/Minimum_cut) from the residual graph left behind by Ford-Fulkerson, illustrating the max-flow min-cut theorem. |
+| [`dinic.py`](dinic.py) | [Dinic's algorithm](https://en.wikipedia.org/wiki/Dinic%27s_algorithm): repeatedly build a BFS *level graph* and saturate a *blocking flow* on it. Adjacency-list based, so it handles parallel edges and sparse graphs well. Runs in `O(V^2 * E)`, or `O(E * sqrt(V))` on unit-capacity networks. |
+| [`push_relabel.py`](push_relabel.py) | The [push-relabel](https://en.wikipedia.org/wiki/Push%E2%80%93relabel_maximum_flow_algorithm) (Goldberg-Tarjan) method: instead of augmenting whole paths, it maintains a *preflow* and locally pushes excess towards the sink. With highest-label selection it runs in `O(V^2 * sqrt(E))`, and is a strong choice on dense graphs. |
+
+## Which one should I use?
+
+All four compute the same maximum-flow value; they differ in speed and in how
+the graph is represented.
+
+* **Just learning the idea?** Start with `ford_fulkerson.py` and
+  `minimum_cut.py` — the augmenting-path picture is the most intuitive.
+* **Sparse graph, or parallel edges?** Reach for `dinic.py`; the adjacency-list
+  representation and level-graph batching make it fast in practice.
+* **Dense graph?** `push_relabel.py` tends to win, because it avoids
+  re-scanning long augmenting paths.
+
+Each file is self-contained, fully type-hinted, and verified with doctests — run
+any of them directly (for example `python networking_flow/dinic.py`) to execute
+the tests.
diff --git a/networking_flow/dinic.py b/networking_flow/dinic.py
new file mode 100644
index 000000000000..f522affe304b
--- /dev/null
+++ b/networking_flow/dinic.py
@@ -0,0 +1,166 @@
+"""
+Dinic's algorithm for the maximum-flow problem.
+
+Dinic's algorithm repeatedly builds a *level graph* with a breadth-first search
+(shortest augmenting paths, measured in edges) and then, in one pass, saturates
+a *blocking flow* on that level graph using depth-first search.  Grouping the
+augmenting paths by length this way gives a much better worst case than the
+plain Ford-Fulkerson / Edmonds-Karp augmenting-path method:
+
+* Dinic's algorithm:            O(V^2 * E)
+* on unit-capacity networks:    O(E * sqrt(V))
+
+Unlike the adjacency-matrix implementations in ``ford_fulkerson.py`` and
+``minimum_cut.py`` in this directory, this version stores the graph as an
+adjacency list of residual edges, so it also handles graphs with parallel edges
+and is efficient on sparse graphs.
+
+Reference: https://en.wikipedia.org/wiki/Dinic%27s_algorithm
+"""
+
+from collections import deque
+
+
+class Dinic:
+    """
+    Maximum flow in a directed graph with non-negative integer capacities.
+
+    Add edges with :meth:`add_edge`, then call :meth:`max_flow`.
+
+    >>> g = Dinic(6)
+    >>> capacities = {
+    ...     (0, 1): 16, (0, 2): 13, (1, 2): 10, (1, 3): 12,
+    ...     (2, 1): 4, (2, 4): 14, (3, 2): 9, (3, 5): 20,
+    ...     (4, 3): 7, (4, 5): 4,
+    ... }
+    >>> for (u, v), cap in capacities.items():
+    ...     g.add_edge(u, v, cap)
+    >>> g.max_flow(0, 5)
+    23
+
+    A source with no outgoing edges (or a sink with no incoming edges) has zero
+    maximum flow:
+
+    >>> Dinic(3).max_flow(0, 2)
+    0
+
+    Parallel edges between the same pair of vertices are supported and their
+    capacities add up:
+
+    >>> h = Dinic(2)
+    >>> h.add_edge(0, 1, 3)
+    >>> h.add_edge(0, 1, 5)
+    >>> h.max_flow(0, 1)
+    8
+    """
+
+    def __init__(self, vertices: int) -> None:
+        if vertices <= 0:
+            raise ValueError("number of vertices must be positive")
+        self.size = vertices
+        # graph[vertex] holds indices into self.edges for edges leaving that vertex.
+        self.graph: list[list[int]] = [[] for _ in range(vertices)]
+        # Each edge is stored as [destination, residual_capacity].
+        # Edge i and its reverse edge i ^ 1 are always created together.
+        self.edges: list[list[int]] = []
+
+    def add_edge(self, source: int, destination: int, capacity: int) -> None:
+        """
+        Add a directed edge ``source -> destination`` with the given capacity.
+
+        >>> g = Dinic(2)
+        >>> g.add_edge(0, 1, 5)
+        >>> g.add_edge(0, 1, -1)
+        Traceback (most recent call last):
+            ...
+        ValueError: capacity must be non-negative
+        >>> g.add_edge(0, 2, 5)
+        Traceback (most recent call last):
+            ...
+        ValueError: vertex out of range
+        """
+        if capacity < 0:
+            raise ValueError("capacity must be non-negative")
+        if not (0 <= source < self.size and 0 <= destination < self.size):
+            raise ValueError("vertex out of range")
+        self.graph[source].append(len(self.edges))
+        self.edges.append([destination, capacity])
+        self.graph[destination].append(len(self.edges))
+        self.edges.append([source, 0])  # reverse edge starts saturated
+
+    def _build_level_graph(self, source: int) -> list[int]:
+        """Breadth-first search; return per-vertex levels (-1 if unreachable)."""
+        level = [-1] * self.size
+        level[source] = 0
+        queue = deque([source])
+        while queue:
+            vertex = queue.popleft()
+            for edge_index in self.graph[vertex]:
+                destination, residual = self.edges[edge_index]
+                if residual > 0 and level[destination] == -1:
+                    level[destination] = level[vertex] + 1
+                    queue.append(destination)
+        return level
+
+    def _send_flow(
+        self,
+        vertex: int,
+        pushed: int,
+        sink: int,
+        level: list[int],
+        progress: list[int],
+    ) -> int:
+        """Depth-first search that pushes a blocking flow along the level graph."""
+        if vertex == sink:
+            return pushed
+        while progress[vertex] < len(self.graph[vertex]):
+            edge_index = self.graph[vertex][progress[vertex]]
+            destination, residual = self.edges[edge_index]
+            if residual > 0 and level[destination] == level[vertex] + 1:
+                flow = self._send_flow(
+                    destination, min(pushed, residual), sink, level, progress
+                )
+                if flow > 0:
+                    self.edges[edge_index][1] -= flow
+                    self.edges[edge_index ^ 1][1] += flow
+                    return flow
+            progress[vertex] += 1
+        return 0
+
+    def max_flow(self, source: int, sink: int) -> int:
+        """
+        Return the maximum flow from ``source`` to ``sink``.
+
+        >>> g = Dinic(4)
+        >>> for (u, v), cap in {(0, 1): 3, (0, 2): 2, (1, 2): 5,
+        ...                     (1, 3): 2, (2, 3): 3}.items():
+        ...     g.add_edge(u, v, cap)
+        >>> g.max_flow(0, 3)
+        5
+        >>> g.max_flow(0, 0)
+        Traceback (most recent call last):
+            ...
+        ValueError: source and sink must be different
+        """
+        if not (0 <= source < self.size and 0 <= sink < self.size):
+            raise ValueError("vertex out of range")
+        if source == sink:
+            raise ValueError("source and sink must be different")
+        infinity = sum(capacity for _, capacity in self.edges) + 1
+        flow = 0
+        level = self._build_level_graph(source)
+        while level[sink] != -1:
+            progress = [0] * self.size
+            while True:
+                pushed = self._send_flow(source, infinity, sink, level, progress)
+                if pushed == 0:
+                    break
+                flow += pushed
+            level = self._build_level_graph(source)
+        return flow
+
+
+if __name__ == "__main__":
+    from doctest import testmod
+
+    testmod()
diff --git a/networking_flow/minimum_cut.py b/networking_flow/minimum_cut.py
index 164b45f1012a..c1f4a83b6aee 100644
--- a/networking_flow/minimum_cut.py
+++ b/networking_flow/minimum_cut.py
@@ -1,4 +1,16 @@
-# Minimum cut on Ford_Fulkerson algorithm.
+"""
+Minimum cut of a flow network via the Ford-Fulkerson algorithm.
+
+The max-flow min-cut theorem says the value of a maximum flow from the source to
+the sink equals the total capacity of the edges in a minimum s-t cut -- the
+cheapest set of edges whose removal disconnects the sink from the source.  This
+module finds those cut edges: it runs Ford-Fulkerson to build the residual
+graph, then reports every original edge that goes from a vertex still reachable
+from the source to a vertex that is not.
+
+Reference: https://en.wikipedia.org/wiki/Minimum_cut
+See also: https://en.wikipedia.org/wiki/Max-flow_min-cut_theorem
+"""
 
 test_graph = [
     [0, 16, 13, 0, 0, 0],
@@ -10,57 +22,80 @@
 ]
 
 
-def bfs(graph, s, t, parent):
-    # Return True if there is node that has not iterated.
+def bfs(graph: list[list[int]], source: int, sink: int, parent: list[int]) -> bool:
+    """
+    Return True if the ``sink`` is reachable from the ``source`` in the
+    residual ``graph``, recording the traversal tree in ``parent``.
+
+    >>> bfs(test_graph, 0, 5, [-1] * 6)
+    True
+    >>> bfs([[0, 0], [0, 0]], 0, 1, [-1, -1])
+    False
+    """
     visited = [False] * len(graph)
-    queue = [s]
-    visited[s] = True
+    queue = [source]
+    visited[source] = True
 
     while queue:
-        u = queue.pop(0)
-        for ind in range(len(graph[u])):
-            if visited[ind] is False and graph[u][ind] > 0:
-                queue.append(ind)
-                visited[ind] = True
-                parent[ind] = u
+        node = queue.pop(0)
+        for neighbor in range(len(graph[node])):
+            if visited[neighbor] is False and graph[node][neighbor] > 0:
+                queue.append(neighbor)
+                visited[neighbor] = True
+                parent[neighbor] = node
+
+    return visited[sink]
+
 
-    return visited[t]
+def mincut(graph: list[list[int]], source: int, sink: int) -> list[tuple[int, int]]:
+    """
+    Return the edges of a minimum s-t cut as ``(from, to)`` tuples.
 
+    The input ``graph`` is an adjacency matrix of capacities and is left
+    unchanged (the algorithm works on an internal copy).
 
-def mincut(graph, source, sink):
-    """This array is filled by BFS and to store path
     >>> mincut(test_graph, source=0, sink=5)
     [(1, 3), (4, 3), (4, 5)]
+
+    The capacities of the cut edges sum to the maximum flow (23 here):
+
+    >>> sum(test_graph[u][v] for u, v in mincut(test_graph, 0, 5))
+    23
+
+    A single saturated edge is its own minimum cut:
+
+    >>> mincut([[0, 7], [0, 0]], source=0, sink=1)
+    [(0, 1)]
     """
-    parent = [-1] * (len(graph))
-    max_flow = 0
+    residual = [row[:] for row in graph]  # work on a copy; keep the input intact
+    parent = [-1] * (len(residual))
     res = []
-    temp = [i[:] for i in graph]  # Record original cut, copy.
-    while bfs(graph, source, sink, parent):
-        path_flow = float("Inf")
+    while bfs(residual, source, sink, parent):
+        path_flow = float("inf")
         s = sink
 
         while s != source:
-            # Find the minimum value in select path
-            path_flow = min(path_flow, graph[parent[s]][s])
+            # Find the minimum residual capacity along the augmenting path.
+            path_flow = min(path_flow, residual[parent[s]][s])
             s = parent[s]
 
-        max_flow += path_flow
         v = sink
-
         while v != source:
             u = parent[v]
-            graph[u][v] -= path_flow
-            graph[v][u] += path_flow
+            residual[u][v] -= path_flow
+            residual[v][u] += path_flow
             v = parent[v]
 
     for i in range(len(graph)):
         for j in range(len(graph[0])):
-            if graph[i][j] == 0 and temp[i][j] > 0:
+            if graph[i][j] > 0 and residual[i][j] == 0:
                 res.append((i, j))
 
     return res
 
 
 if __name__ == "__main__":
+    from doctest import testmod
+
+    testmod()
     print(mincut(test_graph, source=0, sink=5))
diff --git a/networking_flow/push_relabel.py b/networking_flow/push_relabel.py
new file mode 100644
index 000000000000..9aa6083f2836
--- /dev/null
+++ b/networking_flow/push_relabel.py
@@ -0,0 +1,162 @@
+"""
+Push-relabel (Goldberg-Tarjan) algorithm for the maximum-flow problem.
+
+The push-relabel method takes a very different approach from the augmenting-path
+algorithms in this directory (``ford_fulkerson.py`` builds up a valid flow one
+path at a time).  Instead it works with a *preflow*, in which a vertex may
+temporarily receive more flow than it sends out.  Each active vertex either
+*pushes* its excess towards a neighbour that is one level lower, or is *relabeled*
+to a higher level so that a push becomes possible.  When no vertex other than the
+source and sink has excess, the preflow has become a maximum flow.
+
+Using the highest-label selection rule (always discharge an active vertex whose
+label is largest) this implementation runs in O(V^2 * sqrt(E)) time, which beats
+the augmenting-path methods on dense graphs.
+
+Reference: https://en.wikipedia.org/wiki/Push%E2%80%93relabel_maximum_flow_algorithm
+"""
+
+from __future__ import annotations
+
+
+class PushRelabel:
+    """
+    Maximum flow in a directed graph with non-negative integer capacities.
+
+    Add edges with :meth:`add_edge`, then call :meth:`max_flow`.
+
+    >>> g = PushRelabel(6)
+    >>> capacities = {
+    ...     (0, 1): 16, (0, 2): 13, (1, 2): 10, (1, 3): 12,
+    ...     (2, 1): 4, (2, 4): 14, (3, 2): 9, (3, 5): 20,
+    ...     (4, 3): 7, (4, 5): 4,
+    ... }
+    >>> for (u, v), cap in capacities.items():
+    ...     g.add_edge(u, v, cap)
+    >>> g.max_flow(0, 5)
+    23
+
+    It agrees with the classic four-vertex example:
+
+    >>> h = PushRelabel(4)
+    >>> for (u, v), cap in {(0, 1): 3, (0, 2): 2, (1, 2): 5,
+    ...                     (1, 3): 2, (2, 3): 3}.items():
+    ...     h.add_edge(u, v, cap)
+    >>> h.max_flow(0, 3)
+    5
+
+    Parallel edges add up, and a disconnected sink gives zero flow:
+
+    >>> p = PushRelabel(2)
+    >>> p.add_edge(0, 1, 3)
+    >>> p.add_edge(0, 1, 5)
+    >>> p.max_flow(0, 1)
+    8
+    >>> PushRelabel(3).max_flow(0, 2)
+    0
+    """
+
+    def __init__(self, vertices: int) -> None:
+        if vertices <= 0:
+            raise ValueError("number of vertices must be positive")
+        self.size = vertices
+        self.graph: list[list[int]] = [[] for _ in range(vertices)]
+        # Each edge is stored as [destination, residual_capacity].
+        self.edges: list[list[int]] = []
+
+    def add_edge(self, source: int, destination: int, capacity: int) -> None:
+        """
+        Add a directed edge ``source -> destination`` with the given capacity.
+
+        >>> g = PushRelabel(2)
+        >>> g.add_edge(0, 1, -1)
+        Traceback (most recent call last):
+            ...
+        ValueError: capacity must be non-negative
+        >>> g.add_edge(2, 0, 1)
+        Traceback (most recent call last):
+            ...
+        ValueError: vertex out of range
+        """
+        if capacity < 0:
+            raise ValueError("capacity must be non-negative")
+        if not (0 <= source < self.size and 0 <= destination < self.size):
+            raise ValueError("vertex out of range")
+        self.graph[source].append(len(self.edges))
+        self.edges.append([destination, capacity])
+        self.graph[destination].append(len(self.edges))
+        self.edges.append([source, 0])  # reverse edge starts saturated
+
+    def max_flow(self, source: int, sink: int) -> int:
+        """
+        Return the maximum flow from ``source`` to ``sink``.
+
+        >>> PushRelabel(2).max_flow(0, 0)
+        Traceback (most recent call last):
+            ...
+        ValueError: source and sink must be different
+        """
+        if not (0 <= source < self.size and 0 <= sink < self.size):
+            raise ValueError("vertex out of range")
+        if source == sink:
+            raise ValueError("source and sink must be different")
+
+        height = [0] * self.size
+        excess = [0] * self.size
+        height[source] = self.size
+
+        # Saturate every edge leaving the source to create the initial preflow.
+        for edge_index in self.graph[source]:
+            destination, residual = self.edges[edge_index]
+            if residual > 0:
+                self.edges[edge_index][1] -= residual
+                self.edges[edge_index ^ 1][1] += residual
+                excess[destination] += residual
+                excess[source] -= residual
+
+        active = [
+            v for v in range(self.size) if v not in (source, sink) and excess[v] > 0
+        ]
+
+        while active:
+            u = max(active, key=lambda v: height[v])
+            if not self._discharge(u, height):
+                # Relabel: lift u just above its lowest usable neighbour.
+                min_height = min(
+                    height[self.edges[i][0]]
+                    for i in self.graph[u]
+                    if self.edges[i][1] > 0
+                )
+                height[u] = min_height + 1
+            self._apply_pushes(u, height, excess)
+            active = [
+                v for v in range(self.size) if v not in (source, sink) and excess[v] > 0
+            ]
+
+        return excess[sink]
+
+    def _discharge(self, u: int, height: list[int]) -> bool:
+        """Return ``True`` if ``u`` has at least one admissible outgoing edge."""
+        return any(
+            self.edges[i][1] > 0 and height[self.edges[i][0]] == height[u] - 1
+            for i in self.graph[u]
+        )
+
+    def _apply_pushes(self, u: int, height: list[int], excess: list[int]) -> None:
+        """Push as much excess as possible from ``u`` along admissible edges."""
+        for edge_index in self.graph[u]:
+            if excess[u] == 0:
+                break
+            destination, residual = self.edges[edge_index]
+            if residual > 0 and height[u] == height[destination] + 1:
+                delta = min(excess[u], residual)
+                self.edges[edge_index][1] -= delta
+                self.edges[edge_index ^ 1][1] += delta
+                excess[u] -= delta
+                excess[destination] += delta
+
+
+if __name__ == "__main__":
+    from doctest import testmod
+
+    testmod()