diff --git a/benches/ford_fulkerson.rs b/benches/ford_fulkerson.rs
new file mode 100644
index 0000000..ca02711
--- /dev/null
+++ b/benches/ford_fulkerson.rs
@@ -0,0 +1,24 @@
+#![feature(test)]
+extern crate petgraph;
+extern crate test;
+
+use petgraph::algo::ford_fulkerson;
+use petgraph::prelude::{Graph, NodeIndex};
+use test::Bencher;
+
+#[bench]
+fn ford_fulkerson_bench(bench: &mut Bencher) {
+    static NODE_COUNT: usize = 1_000;
+    let mut g: Graph<usize, usize> = Graph::new();
+    let nodes: Vec<NodeIndex<_>> = (0..NODE_COUNT).map(|i| g.add_node(i)).collect();
+    for i in 0..NODE_COUNT - 1 {
+        g.add_edge(nodes[i], nodes[i + 1], 1);
+    }
+    bench.iter(|| {
+        let _flow = ford_fulkerson(
+            &g,
+            NodeIndex::from(0),
+            NodeIndex::from(g.node_count() as u32 - 1),
+        );
+    });
+}
diff --git a/src/algo/ford_fulkerson.rs b/src/algo/ford_fulkerson.rs
new file mode 100644
index 0000000..3b2127f
--- /dev/null
+++ b/src/algo/ford_fulkerson.rs
@@ -0,0 +1,196 @@
+use std::{collections::VecDeque, ops::Sub};
+
+use crate::{
+    data::DataMap,
+    visit::{
+        EdgeCount, EdgeIndexable, IntoEdges, IntoEdgesDirected, NodeCount, NodeIndexable, VisitMap,
+        Visitable,
+    },
+};
+
+use super::{EdgeRef, PositiveMeasure};
+use crate::prelude::Direction;
+
+fn residual_capacity<N>(
+    network: N,
+    edge: N::EdgeRef,
+    vertex: N::NodeId,
+    flow: N::EdgeWeight,
+) -> N::EdgeWeight
+where
+    N: NodeIndexable + IntoEdges,
+    N::EdgeWeight: Sub<Output = N::EdgeWeight> + PositiveMeasure,
+{
+    if vertex == edge.source() {
+        // backward edge
+        flow
+    } else if vertex == edge.target() {
+        // forward edge
+        return *edge.weight() - flow;
+    } else {
+        let end_point = NodeIndexable::to_index(&network, vertex);
+        panic!("Illegal endpoint {}", end_point);
+    }
+}
+
+/// Gets the other endpoint of graph edge, if any, otherwise panics.
+fn other_endpoint<N>(network: N, edge: N::EdgeRef, vertex: N::NodeId) -> N::NodeId
+where
+    N: NodeIndexable + IntoEdges,
+{
+    if vertex == edge.source() {
+        edge.target()
+    } else if vertex == edge.target() {
+        edge.source()
+    } else {
+        let end_point = NodeIndexable::to_index(&network, vertex);
+        panic!("Illegal endpoint {}", end_point);
+    }
+}
+
+/// Tells whether there is an augmented path in the graph
+fn has_augmented_path<N>(
+    network: N,
+    source: N::NodeId,
+    destination: N::NodeId,
+    edge_to: &mut [Option<N::EdgeRef>],
+    flows: &[N::EdgeWeight],
+) -> bool
+where
+    N: NodeCount + IntoEdgesDirected + NodeIndexable + EdgeIndexable + Visitable,
+    N::EdgeWeight: Sub<Output = N::EdgeWeight> + PositiveMeasure,
+{
+    let mut visited = network.visit_map();
+    let mut queue = VecDeque::new();
+    visited.visit(source);
+    queue.push_back(source);
+
+    while let Some(vertex) = queue.pop_front() {
+        let out_edges = network.edges_directed(vertex, Direction::Outgoing);
+        let in_edges = network.edges_directed(vertex, Direction::Incoming);
+        for edge in out_edges.chain(in_edges) {
+            let next = other_endpoint(&network, edge, vertex);
+            let edge_index: usize = EdgeIndexable::to_index(&network, edge.id());
+            let residual_cap = residual_capacity(&network, edge, next, flows[edge_index]);
+            if !visited.is_visited(&next) && (residual_cap > N::EdgeWeight::zero()) {
+                visited.visit(next);
+                edge_to[NodeIndexable::to_index(&network, next)] = Some(edge);
+                if destination == next {
+                    return true;
+                }
+                queue.push_back(next);
+            }
+        }
+    }
+    false
+}
+
+fn adjust_residual_flow<N>(
+    network: N,
+    edge: N::EdgeRef,
+    vertex: N::NodeId,
+    flow: N::EdgeWeight,
+    delta: N::EdgeWeight,
+) -> N::EdgeWeight
+where
+    N: NodeIndexable + IntoEdges,
+    N::EdgeWeight: Sub<Output = N::EdgeWeight> + PositiveMeasure,
+{
+    if vertex == edge.source() {
+        // backward edge
+        flow - delta
+    } else if vertex == edge.target() {
+        // forward edge
+        flow + delta
+    } else {
+        let end_point = NodeIndexable::to_index(&network, vertex);
+        panic!("Illegal endpoint {}", end_point);
+    }
+}
+
+/// \[Generic\] Ford-Fulkerson algorithm.
+///
+/// Computes the [maximum flow][ff] of a weighted directed graph.
+///
+/// If it terminates, it returns the maximum flow and also the computed edge flows.
+///
+/// [ff]: https://en.wikipedia.org/wiki/Ford%E2%80%93Fulkerson_algorithm
+///
+/// # Example
+/// ```rust
+/// use petgraph::Graph;
+/// use petgraph::algo::ford_fulkerson;
+/// // Example from CLRS book
+/// let mut graph = Graph::<u8, u8>::new();
+/// let source = graph.add_node(0);
+/// let _ = graph.add_node(1);
+/// let _ = graph.add_node(2);
+/// let _ = graph.add_node(3);
+/// let _ = graph.add_node(4);
+/// let destination = graph.add_node(5);
+/// graph.extend_with_edges(&[
+///    (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),
+/// ]);
+/// let (max_flow, _) = ford_fulkerson(&graph, source, destination);
+/// assert_eq!(23, max_flow);
+/// ```
+pub fn ford_fulkerson<N>(
+    network: N,
+    source: N::NodeId,
+    destination: N::NodeId,
+) -> (N::EdgeWeight, Vec<N::EdgeWeight>)
+where
+    N: NodeCount
+        + EdgeCount
+        + IntoEdgesDirected
+        + EdgeIndexable
+        + NodeIndexable
+        + DataMap
+        + Visitable,
+    N::EdgeWeight: Sub<Output = N::EdgeWeight> + PositiveMeasure,
+{
+    let mut edge_to = vec![None; network.node_count()];
+    let mut flows = vec![N::EdgeWeight::zero(); network.edge_count()];
+    let mut max_flow = N::EdgeWeight::zero();
+    while has_augmented_path(&network, source, destination, &mut edge_to, &flows) {
+        let mut path_flow = N::EdgeWeight::max();
+
+        // Find the bottleneck capacity of the path
+        let mut vertex = destination;
+        let mut vertex_index = NodeIndexable::to_index(&network, vertex);
+        while let Some(edge) = edge_to[vertex_index] {
+            let edge_index = EdgeIndexable::to_index(&network, edge.id());
+            let residual_capacity = residual_capacity(&network, edge, vertex, flows[edge_index]);
+            // Minimum between the current path flow and the residual capacity.
+            path_flow = if path_flow > residual_capacity {
+                residual_capacity
+            } else {
+                path_flow
+            };
+            vertex = other_endpoint(&network, edge, vertex);
+            vertex_index = NodeIndexable::to_index(&network, vertex);
+        }
+
+        // Update the flow of each edge along the path
+        let mut vertex = destination;
+        let mut vertex_index = NodeIndexable::to_index(&network, vertex);
+        while let Some(edge) = edge_to[vertex_index] {
+            let edge_index = EdgeIndexable::to_index(&network, edge.id());
+            flows[edge_index] =
+                adjust_residual_flow(&network, edge, vertex, flows[edge_index], path_flow);
+            vertex = other_endpoint(&network, edge, vertex);
+            vertex_index = NodeIndexable::to_index(&network, vertex);
+        }
+        max_flow = max_flow + path_flow;
+    }
+    (max_flow, flows)
+}
diff --git a/src/algo/mod.rs b/src/algo/mod.rs
index d49ad9c..6c0a408 100644
--- a/src/algo/mod.rs
+++ b/src/algo/mod.rs
@@ -10,6 +10,7 @@
 pub mod dominators;
 pub mod feedback_arc_set;
 pub mod floyd_warshall;
+pub mod ford_fulkerson;
 pub mod isomorphism;
 pub mod k_shortest_path;
 pub mod matching;
@@ -36,6 +37,7 @@
 pub use dijkstra::dijkstra;
 pub use feedback_arc_set::greedy_feedback_arc_set;
 pub use floyd_warshall::floyd_warshall;
+pub use ford_fulkerson::ford_fulkerson;
 pub use isomorphism::{
     is_isomorphic, is_isomorphic_matching, is_isomorphic_subgraph, is_isomorphic_subgraph_matching,
     subgraph_isomorphisms_iter,
@@ -836,3 +838,28 @@
     }
 );
 impl_unit_measure!(f32, f64);
+
+/// Some measure of positive numbers, assuming positive
+/// float-pointing numbers
+pub trait PositiveMeasure: Measure + Copy {
+    fn zero() -> Self;
+    fn max() -> Self;
+}
+
+macro_rules! impl_positive_measure(
+    ( $( $t:ident ),* )=> {
+        $(
+            impl PositiveMeasure for $t {
+                fn zero() -> Self {
+                    0 as $t
+                }
+                fn max() -> Self {
+                    std::$t::MAX
+                }
+            }
+
+        )*
+    }
+);
+
+impl_positive_measure!(u8, u16, u32, u64, u128, usize, f32, f64);
diff --git a/tests/ford_fulkerson.rs b/tests/ford_fulkerson.rs
new file mode 100644
index 0000000..0b7881d
--- /dev/null
+++ b/tests/ford_fulkerson.rs
@@ -0,0 +1,121 @@
+use petgraph::algo::ford_fulkerson;
+use petgraph::prelude::Graph;
+
+#[test]
+fn test_ford_fulkerson() {
+    // Example from https://downey.io/blog/max-flow-ford-fulkerson-algorithm-explanation/
+    let mut graph = Graph::<usize, u16>::new();
+    let source = graph.add_node(0);
+    let _ = graph.add_node(1);
+    let _ = graph.add_node(2);
+    let destination = graph.add_node(3);
+    graph.extend_with_edges(&[(0, 1, 3), (0, 2, 2), (1, 2, 5), (1, 3, 2), (2, 3, 3)]);
+    let (max_flow, _) = ford_fulkerson(&graph, source, destination);
+    assert_eq!(5, max_flow);
+
+    // Example from https://brilliant.org/wiki/ford-fulkerson-algorithm/
+    let mut graph = Graph::<usize, f32>::new();
+    let source = graph.add_node(0);
+    let _ = graph.add_node(1);
+    let _ = graph.add_node(2);
+    let _ = graph.add_node(3);
+    let _ = graph.add_node(4);
+    let destination = graph.add_node(5);
+    graph.extend_with_edges(&[
+        (0, 1, 4.),
+        (0, 2, 3.),
+        (1, 3, 4.),
+        (2, 4, 6.),
+        (3, 2, 3.),
+        (3, 5, 2.),
+        (4, 5, 6.),
+    ]);
+    let (max_flow, _) = ford_fulkerson(&graph, source, destination);
+    assert_eq!(7.0, max_flow);
+
+    // Example from https://cp-algorithms.com/graph/edmonds_karp.html
+    let mut graph = Graph::<usize, f32>::new();
+    let source = graph.add_node(0);
+    let _ = graph.add_node(1);
+    let _ = graph.add_node(2);
+    let _ = graph.add_node(3);
+    let _ = graph.add_node(4);
+    let destination = graph.add_node(5);
+    graph.extend_with_edges(&[
+        (0, 1, 7.),
+        (0, 2, 4.),
+        (1, 3, 5.),
+        (1, 4, 3.),
+        (2, 1, 3.),
+        (2, 4, 2.),
+        (3, 5, 8.),
+        (4, 3, 3.),
+        (4, 5, 5.),
+    ]);
+    let (max_flow, _) = ford_fulkerson(&graph, source, destination);
+    assert_eq!(10.0, max_flow);
+
+    // Example from https://www.programiz.com/dsa/ford-fulkerson-algorithm (corrected: result not 6 but 5)
+    let mut graph = Graph::<u8, f32>::new();
+    let source = graph.add_node(0);
+    let _ = graph.add_node(1);
+    let _ = graph.add_node(2);
+    let _ = graph.add_node(3);
+    let _ = graph.add_node(4);
+    let destination = graph.add_node(5);
+    graph.extend_with_edges(&[
+        (0, 1, 8.),
+        (0, 2, 3.),
+        (1, 3, 9.),
+        (2, 3, 7.),
+        (2, 4, 4.),
+        (3, 5, 2.),
+        (4, 5, 5.),
+    ]);
+    let (max_flow, _) = ford_fulkerson(&graph, source, destination);
+    assert_eq!(5.0, max_flow);
+
+    let mut graph = Graph::<u8, u8>::new();
+    let source = graph.add_node(0);
+    let _ = graph.add_node(1);
+    let _ = graph.add_node(2);
+    let _ = graph.add_node(3);
+    let _ = graph.add_node(4);
+    let destination = graph.add_node(5);
+    graph.extend_with_edges(&[
+        (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),
+    ]);
+    let (max_flow, _) = ford_fulkerson(&graph, source, destination);
+    assert_eq!(23, max_flow);
+
+    // Example taken from https://medium.com/@jithmisha/solving-the-maximum-flow-problem-with-ford-fulkerson-method-3fccc2883dc7
+    let mut graph = Graph::<u8, u8>::new();
+    let source = graph.add_node(0);
+    let _ = graph.add_node(1);
+    let _ = graph.add_node(2);
+    let _ = graph.add_node(3);
+    let _ = graph.add_node(4);
+    let destination = graph.add_node(5);
+    graph.extend_with_edges(&[
+        (0, 1, 10),
+        (0, 2, 10),
+        (1, 2, 2),
+        (1, 3, 4),
+        (1, 4, 8),
+        (2, 4, 9),
+        (3, 5, 10),
+        (4, 3, 6),
+        (4, 5, 10),
+    ]);
+    let (max_flow, _) = ford_fulkerson(&graph, source, destination);
+    assert_eq!(19, max_flow);
+}
diff --git a/tests/quickcheck.rs b/tests/quickcheck.rs
index 1f7a1af..3b9ae69 100644
--- a/tests/quickcheck.rs
+++ b/tests/quickcheck.rs
@@ -23,7 +23,7 @@
 use rand::Rng;
 
 use petgraph::algo::{
-    bellman_ford, condensation, dijkstra, find_negative_cycle, floyd_warshall,
+    bellman_ford, condensation, dijkstra, find_negative_cycle, floyd_warshall, ford_fulkerson,
     greedy_feedback_arc_set, greedy_matching, is_cyclic_directed, is_cyclic_undirected,
     is_isomorphic, is_isomorphic_matching, k_shortest_path, kosaraju_scc, maximum_matching,
     min_spanning_tree, page_rank, tarjan_scc, toposort, Matching,
@@ -35,7 +35,7 @@
 use petgraph::operator::complement;
 use petgraph::prelude::*;
 use petgraph::visit::{
-    EdgeFiltered, EdgeRef, IntoEdgeReferences, IntoEdges, IntoNeighbors, IntoNodeIdentifiers,
+    EdgeFiltered, EdgeIndexable, IntoEdgeReferences, IntoEdges, IntoNeighbors, IntoNodeIdentifiers,
     IntoNodeReferences, NodeCount, NodeIndexable, Reversed, Topo, VisitMap, Visitable,
 };
 use petgraph::EdgeType;
@@ -1330,3 +1330,49 @@
         true
     }
 }
+
+fn sum_flows<N, F: std::iter::Sum + Copy>(
+    gr: &Graph<N, F>,
+    flows: &[F],
+    node: NodeIndex,
+    dir: Direction,
+) -> F {
+    gr.edges_directed(node, dir)
+        .map(|edge| flows[EdgeIndexable::to_index(&gr, edge.id())])
+        .sum::<F>()
+}
+
+quickcheck! {
+    // 1. (Capacity)
+    //    The flows should be <= capacities
+    // 2. (Flow conservation)
+    //    For every internal node (i.e a node different from the
+    //    source node and the destination (or sink) node), the sum
+    //    of incoming flows (i.e flows of incoming edges) is equal
+    //    to the sum of the outgoing flows (i.e flows of outgoing edges).
+    // 3. (Maximum flow)
+    //    It is equal to the sum of the destination node incoming flows and
+    //    also the sum of the outgoing flows of the source node.
+    fn test_ford_fulkerson_flows(gr: Graph<usize, u32>) -> bool {
+        if gr.node_count() <= 1 || gr.edge_count() == 0 {
+            return true;
+        }
+        let source = NodeIndex::from(0);
+        let destination = NodeIndex::from(gr.node_count() as u32 / 2);
+        let (max_flow, flows) = ford_fulkerson(&gr, source, destination);
+        let capacity_constraint = flows
+            .iter()
+            .enumerate()
+            .all(|(ix, flow)| flow <= gr.edge_weight(EdgeIndexable::from_index(&gr, ix)).unwrap());
+        let flow_conservation_constraint = (0..gr.node_count()).all(|ix| {
+            let node = NodeIndexable::from_index(&gr, ix);
+            if (node != source) && (node != destination){
+            sum_flows(&gr, &flows, node, Direction::Outgoing)
+                == sum_flows(&gr, &flows, node, Direction::Incoming)
+            } else {true}
+        });
+        let max_flow_constaint = (sum_flows(&gr, &flows, source, Direction::Outgoing) == max_flow)
+            && (sum_flows(&gr, &flows, destination, Direction::Incoming) == max_flow);
+        return capacity_constraint && flow_conservation_constraint && max_flow_constaint;
+    }
+}