diff --git a/ortools/graph/BUILD.bazel b/ortools/graph/BUILD.bazel
index 35525277bf..da456654fc 100644
--- a/ortools/graph/BUILD.bazel
+++ b/ortools/graph/BUILD.bazel
@@ -44,6 +44,18 @@ cc_library(
     ],
 )
 
+cc_library(
+    name = "flow_graph",
+    hdrs = ["flow_graph.h"],
+    deps = [
+        ":graph",
+        ":iterators",
+        "//ortools/base:stl_util",
+        "@com_google_absl//absl/log:check",
+        "@com_google_absl//absl/types:span",
+    ],
+)
+
 cc_library(
     name = "bfs",
     hdrs = ["bfs.h"],
@@ -456,7 +468,6 @@ cc_library(
         "//ortools/base",
         "//ortools/util:stats",
         "//ortools/util:zvector",
-        "@com_google_absl//absl/memory",
         "@com_google_absl//absl/strings",
     ],
 )
@@ -467,6 +478,7 @@ cc_test(
     srcs = ["generic_max_flow_test.cc"],
     deps = [
         ":ebert_graph",
+        ":flow_graph",
         ":generic_max_flow",
         ":graph",
         "//ortools/base",
@@ -533,18 +545,20 @@ cc_binary(
     name = "solve_flow_model",
     srcs = ["solve_flow_model.cc"],
     deps = [
+        ":flow_graph",
         ":flow_problem_cc_proto",
+        ":generic_max_flow",
         ":graph",
-        ":max_flow",
         ":min_cost_flow",
         "//ortools/base",
         "//ortools/base:file",
         "//ortools/base:path",
-        "//ortools/base:status_macros",
         "//ortools/base:timer",
+        "//ortools/util:file_util",
         "//ortools/util:filelineiter",
         "//ortools/util:stats",
         "@com_google_absl//absl/flags:flag",
+        "@com_google_absl//absl/log:check",
         "@com_google_absl//absl/strings",
         "@com_google_absl//absl/strings:str_format",
     ],
diff --git a/ortools/graph/flow_graph.h b/ortools/graph/flow_graph.h
new file mode 100644
index 0000000000..4bb84b7daf
--- /dev/null
+++ b/ortools/graph/flow_graph.h
@@ -0,0 +1,469 @@
+// Copyright 2010-2025 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef UTIL_GRAPH_FLOW_GRAPH_H_
+#define UTIL_GRAPH_FLOW_GRAPH_H_
+
+#include <cstddef>
+#include <cstdint>
+#include <memory>
+#include <vector>
+
+#include "absl/log/check.h"
+#include "absl/types/span.h"
+#include "ortools/base/stl_util.h"
+#include "ortools/graph/graph.h"
+#include "ortools/graph/iterators.h"
+
+namespace util {
+
+// Graph specialized for max-flow/min-cost-flow algorithms.
+// It follows the ortools/graph/graph.h interface.
+//
+// For max-flow or min-cost-flow we need a directed graph where each arc from
+// tail to head has a "reverse" arc from head to tail. In practice many input
+// graphs already have such reverse arc and it can make a big difference just to
+// reuse them.
+//
+// This is similar to ReverseArcStaticGraph but handle reverse arcs in a
+// different way. Instead of always creating a NEW reverse arc for each arc of
+// the input graph, this will detect if a reverse arc was already present in the
+// input, and do not create a new one when this is the case. In the best case,
+// this can gain a factor 2 in the final graph size, note however that the graph
+// construction is slighlty slower because of this detection.
+//
+// A side effect of reusing reverse arc is also that we cannot anymore clearly
+// separate the original arcs from the newly created one. So the algorithm needs
+// to be able to handle that.
+//
+// TODO(user): Currently only max-flow handles this graph, but not
+// min-cost-flow.
+template <typename NodeIndexType = int32_t, typename ArcIndexType = int32_t>
+class FlowGraph : public BaseGraph<NodeIndexType, ArcIndexType, false> {
+  // Note that we do NOT use negated indices for reverse arc. So we use false
+  // for the last template argument here HasNegativeReverseArcs.
+  typedef BaseGraph<NodeIndexType, ArcIndexType, false> Base;
+  using Base::arc_capacity_;
+  using Base::const_capacities_;
+  using Base::node_capacity_;
+  using Base::num_arcs_;
+  using Base::num_nodes_;
+
+ public:
+  FlowGraph() = default;
+
+  FlowGraph(NodeIndexType num_nodes, ArcIndexType arc_capacity) {
+    this->Reserve(num_nodes, arc_capacity);
+    this->FreezeCapacities();
+    this->AddNode(num_nodes - 1);
+  }
+
+  NodeIndexType Head(ArcIndexType arc) const {
+    DCHECK(is_built_);
+    DCHECK_GE(arc, 0);
+    DCHECK_LT(arc, num_arcs_);
+    return heads_[arc];
+  }
+
+  NodeIndexType Tail(ArcIndexType arc) const {
+    DCHECK(is_built_);
+    DCHECK_GE(arc, 0);
+    DCHECK_LT(arc, num_arcs_);
+
+    // Note that we could trade memory for speed if this happens to be hot.
+    // However, it is expected that most user will access arcs via
+    // for (const int arc : graph.OutgoingArcs(tail)) {}
+    // in which case all arcs already have a known tail.
+    return heads_[reverse_[arc]];
+  }
+
+  // Each arc has a reverse.
+  // If not added by the client, we have created one, see Build().
+  ArcIndexType OppositeArc(ArcIndexType arc) const {
+    DCHECK(is_built_);
+    DCHECK_GE(arc, 0);
+    DCHECK_LT(arc, num_arcs_);
+    return reverse_[arc];
+  }
+
+  // Iteration.
+  util::IntegerRange<ArcIndexType> OutgoingArcs(NodeIndexType node) const {
+    return OutgoingArcsStartingFrom(node, start_[node]);
+  }
+  util::IntegerRange<ArcIndexType> OutgoingArcsStartingFrom(
+      NodeIndexType node, ArcIndexType from) const {
+    DCHECK(is_built_);
+    DCHECK_GE(node, 0);
+    DCHECK_LT(node, num_nodes_);
+    DCHECK_GE(from, start_[node]);
+    return util::IntegerRange<ArcIndexType>(from, start_[node + 1]);
+  }
+
+  // These are needed to use with the current max-flow implementation.
+  // We don't distinguish direct from reverse arc anymore, and this is just
+  // the same as OutgoingArcs()/OutgoingArcsStartingFrom().
+  util::IntegerRange<ArcIndexType> OutgoingOrOppositeIncomingArcs(
+      NodeIndexType node) const {
+    return OutgoingArcs(node);
+  }
+  util::IntegerRange<ArcIndexType> OutgoingOrOppositeIncomingArcsStartingFrom(
+      NodeIndexType node, ArcIndexType from) const {
+    return OutgoingArcsStartingFrom(node, from);
+  }
+
+  absl::Span<const NodeIndexType> operator[](NodeIndexType node) const {
+    const int b = start_[node];
+    const size_t size = start_[node + 1] - b;
+    if (size == 0) return {};
+    return {&heads_[b], size};
+  }
+
+  void ReserveArcs(ArcIndexType bound) override {
+    Base::ReserveArcs(bound);
+    tmp_tails_.reserve(bound);
+    tmp_heads_.reserve(bound);
+  }
+
+  void AddNode(NodeIndexType node) {
+    num_nodes_ = std::max(num_nodes_, node + 1);
+  }
+
+  ArcIndexType AddArc(NodeIndexType tail, NodeIndexType head) {
+    AddNode(tail > head ? tail : head);
+    tmp_tails_.push_back(tail);
+    tmp_heads_.push_back(head);
+    return num_arcs_++;
+  }
+
+  void Build() { Build(nullptr); }
+  void Build(std::vector<ArcIndexType>* permutation);
+
+  // This influence what Build() does. If true, we will detect already existing
+  // pairs of (arc, reverse_arc) and only construct new reverse arc for the one
+  // that we couldn't match. Otherwise, we will construct a new reverse arc for
+  // each input arcs.
+  void SetDetectReverse(bool value) { option_detect_reverse_ = value; }
+
+  // This influence what Build() does. If true, the order of each outgoing arc
+  // will be sorted by their head. Otherwise, it will be in the original order
+  // with the newly created reverse arc afterwards.
+  void SetSortByHead(bool value) { option_sort_by_head_ = value; }
+
+ private:
+  // Computes the permutation that would stable_sort input, and fill start_
+  // to correspond to the beginning of each section of identical elements.
+  // This assumes input only contains element in [0, num_nodes_).
+  void FillPermutationAndStart(absl::Span<const NodeIndexType> input,
+                               absl::Span<ArcIndexType> perm);
+
+  // These are similar but tie-break identical element by "second_criteria".
+  // We have two versions. It seems filling the reverse permutation instead is
+  // slighthly faster and require less memory.
+  void FillPermutationAndStart(absl::Span<const NodeIndexType> first_criteria,
+                               absl::Span<const NodeIndexType> second_criteria,
+                               absl::Span<ArcIndexType> perm);
+  void FillReversePermutationAndStart(
+      absl::Span<const NodeIndexType> first_criteria,
+      absl::Span<const NodeIndexType> second_criteria,
+      absl::Span<ArcIndexType> reverse_perm);
+
+  // Internal helpers for the Fill*() functions above.
+  void InitializeStart(absl::Span<const NodeIndexType> input);
+  void RestoreStart();
+
+  // Different build options.
+  bool option_detect_reverse_ = true;
+  bool option_sort_by_head_ = false;
+
+  // Starts at false and set to true on Build(). This is mainly used in
+  // DCHECKs() to guarantee that the graph is just built once, and no new arcs
+  // are added afterwards.
+  bool is_built_ = false;
+
+  // This is just used before Build() to store the AddArc() data.
+  std::vector<NodeIndexType> tmp_tails_;
+  std::vector<NodeIndexType> tmp_heads_;
+
+  // First outgoing arc for a node.
+  // Indexed by NodeIndexType + a sentinel start_[num_nodes_] = num_arcs_.
+  std::unique_ptr<ArcIndexType[]> start_;
+
+  // Indexed by ArcIndexType an of size num_arcs_.
+  // Head for an arc.
+  std::unique_ptr<NodeIndexType[]> heads_;
+  // Reverse arc for an arc.
+  std::unique_ptr<ArcIndexType[]> reverse_;
+};
+
+namespace internal {
+
+// Helper to permute a given span into another one.
+template <class Key, class Value>
+void PermuteInto(absl::Span<const Key> permutation,
+                 absl::Span<const Value> input, absl::Span<Value> output) {
+  DCHECK_EQ(permutation.size(), input.size());
+  DCHECK_EQ(permutation.size(), output.size());
+  for (int i = 0; i < permutation.size(); ++i) {
+    output[permutation[i]] = input[i];
+  }
+}
+
+}  // namespace internal
+
+template <typename NodeIndexType, typename ArcIndexType>
+void FlowGraph<NodeIndexType, ArcIndexType>::InitializeStart(
+    absl::Span<const NodeIndexType> input) {
+  // Computes the number of each elements.
+  std::fill(start_.get(), start_.get() + num_nodes_, 0);
+  start_[num_nodes_] = input.size();  // sentinel.
+
+  for (const NodeIndexType node : input) start_[node]++;
+
+  // Compute the cumulative sums (shifted one element to the right).
+  int sum = 0;
+  for (int i = 0; i < num_nodes_; ++i) {
+    int temp = start_[i];
+    start_[i] = sum;
+    sum += temp;
+  }
+  DCHECK_EQ(sum, input.size());
+}
+
+template <typename NodeIndexType, typename ArcIndexType>
+void FlowGraph<NodeIndexType, ArcIndexType>::RestoreStart() {
+  // Restore in start[i] the index of the first arc with tail >= i.
+  for (int i = num_nodes_; --i > 0;) {
+    start_[i] = start_[i - 1];
+  }
+  start_[0] = 0;
+}
+
+template <typename NodeIndexType, typename ArcIndexType>
+void FlowGraph<NodeIndexType, ArcIndexType>::FillPermutationAndStart(
+    absl::Span<const NodeIndexType> input, absl::Span<ArcIndexType> perm) {
+  DCHECK_EQ(input.size(), perm.size());
+  InitializeStart(input);
+
+  // Computes the permutation.
+  // Note that this temporarily alters the start_ vector.
+  for (int i = 0; i < input.size(); ++i) {
+    perm[i] = start_[input[i]]++;
+  }
+
+  RestoreStart();
+}
+
+template <typename NodeIndexType, typename ArcIndexType>
+void FlowGraph<NodeIndexType, ArcIndexType>::FillPermutationAndStart(
+    absl::Span<const NodeIndexType> first_criteria,
+    absl::Span<const NodeIndexType> second_criteria,
+    absl::Span<ArcIndexType> perm) {
+  // First sort by second_criteria.
+  FillPermutationAndStart(second_criteria, absl::MakeSpan(perm));
+
+  // Create a temporary permuted version of first_criteria.
+  const auto tmp_storage = std::make_unique<ArcIndexType[]>(num_arcs_);
+  const auto tmp = absl::MakeSpan(tmp_storage.get(), num_arcs_);
+  internal::PermuteInto<ArcIndexType, NodeIndexType>(perm, first_criteria, tmp);
+
+  // Now sort by permuted first_criteria.
+  const auto second_perm_storage = std::make_unique<ArcIndexType[]>(num_arcs_);
+  const auto second_perm = absl::MakeSpan(second_perm_storage.get(), num_arcs_);
+  FillPermutationAndStart(tmp, second_perm);
+
+  // Update the final permutation. This guarantee that for the same
+  // first_criteria, the second_criteria will be used.
+  for (ArcIndexType& image : perm) {
+    image = second_perm[image];
+  }
+}
+
+template <typename NodeIndexType, typename ArcIndexType>
+void FlowGraph<NodeIndexType, ArcIndexType>::FillReversePermutationAndStart(
+    absl::Span<const NodeIndexType> first_criteria,
+    absl::Span<const NodeIndexType> second_criteria,
+    absl::Span<ArcIndexType> reverse_perm) {
+  // Compute the reverse perm according to the second criteria.
+  InitializeStart(second_criteria);
+  const auto r_perm_storage = std::make_unique<ArcIndexType[]>(num_arcs_);
+  const auto r_perm = absl::MakeSpan(r_perm_storage.get(), num_arcs_);
+  auto* start = start_.get();
+  for (int i = 0; i < second_criteria.size(); ++i) {
+    r_perm[start[second_criteria[i]]++] = i;
+  }
+
+  // Now radix-sort by the first criteria and combine the two permutations.
+  InitializeStart(first_criteria);
+  for (const int i : r_perm) {
+    reverse_perm[start[first_criteria[i]]++] = i;
+  }
+  RestoreStart();
+}
+
+template <typename NodeIndexType, typename ArcIndexType>
+void FlowGraph<NodeIndexType, ArcIndexType>::Build(
+    std::vector<ArcIndexType>* permutation) {
+  DCHECK(!is_built_);
+  if (is_built_) return;
+  is_built_ = true;
+
+  start_ = std::make_unique<ArcIndexType[]>(num_nodes_ + 1);
+  std::vector<ArcIndexType> perm(num_arcs_);
+
+  const int kNoExplicitReverseArc = -1;
+  std::vector<ArcIndexType> reverse(num_arcs_, kNoExplicitReverseArc);
+
+  bool fix_final_permutation = false;
+  if (option_detect_reverse_) {
+    // Step 1 we only keep a "canonical version" of each arc where tail <= head.
+    // This will allow us to detect reverse as duplicates instead.
+    std::vector<bool> was_reversed_(num_arcs_, false);
+    for (int arc = 0; arc < num_arcs_; ++arc) {
+      if (tmp_heads_[arc] < tmp_tails_[arc]) {
+        std::swap(tmp_heads_[arc], tmp_tails_[arc]);
+        was_reversed_[arc] = true;
+      }
+    }
+
+    // Step 2, compute the perm to sort by tail then head.
+    // This is something we do a few times here, we compute the permutation with
+    // a kind of radix sort by computing the degree of each node.
+    auto reverse_perm = absl::MakeSpan(perm);  // we reuse space
+    FillReversePermutationAndStart(tmp_tails_, tmp_heads_,
+                                   absl::MakeSpan(reverse_perm));
+
+    // Identify arc pairs that are reverse of each other and fill reverse for
+    // them. The others position will stay at -1.
+    int candidate_i = 0;
+    int num_filled = 0;
+    for (int i = 0; i < num_arcs_; ++i) {
+      const int arc = reverse_perm[i];
+      const int candidate_arc = reverse_perm[candidate_i];
+      if (tmp_heads_[arc] != tmp_heads_[candidate_arc] ||
+          tmp_tails_[arc] != tmp_tails_[candidate_arc]) {
+        candidate_i = i;
+        continue;
+      }
+
+      if (was_reversed_[arc] != was_reversed_[candidate_arc]) {
+        DCHECK_EQ(reverse[arc], -1);
+        DCHECK_EQ(reverse[candidate_arc], -1);
+        reverse[arc] = candidate_arc;
+        reverse[candidate_arc] = arc;
+        num_filled += 2;
+
+        // Find the next candidate without reverse if there is a gap.
+        for (; ++candidate_i < i;) {
+          if (reverse[reverse_perm[candidate_i]] == -1) break;
+        }
+        if (candidate_i == i) {
+          candidate_i = ++i;  // Advance by two.
+        }
+      }
+    }
+
+    // Create the extra reversed arcs needed.
+    {
+      const int extra_size = num_arcs_ - num_filled;
+      tmp_tails_.resize(num_arcs_ + extra_size);
+      tmp_heads_.resize(num_arcs_ + extra_size);
+      reverse.resize(num_arcs_ + extra_size);
+      int new_index = num_arcs_;
+      for (int i = 0; i < num_arcs_; ++i) {
+        // Fix the initial swap.
+        if (was_reversed_[i]) {
+          std::swap(tmp_tails_[i], tmp_heads_[i]);
+        }
+
+        if (reverse[i] != kNoExplicitReverseArc) continue;
+        reverse[i] = new_index;
+        reverse[new_index] = i;
+        tmp_tails_[new_index] = tmp_heads_[i];
+        tmp_heads_[new_index] = tmp_tails_[i];
+        ++new_index;
+      }
+      CHECK_EQ(new_index, num_arcs_ + extra_size);
+    }
+  } else {
+    // Just create a reverse for each arc.
+    tmp_tails_.resize(2 * num_arcs_);
+    tmp_heads_.resize(2 * num_arcs_);
+    reverse.resize(2 * num_arcs_);
+    for (int arc = 0; arc < num_arcs_; ++arc) {
+      const int image = num_arcs_ + arc;
+      tmp_heads_[image] = tmp_tails_[arc];
+      tmp_tails_[image] = tmp_heads_[arc];
+      reverse[image] = arc;
+      reverse[arc] = image;
+    }
+
+    // It seems better to put all the reverse first instead of last in
+    // the adjacency list, so lets do that here. Note that we need to fix
+    // the permutation returned to the user in this case.
+    //
+    // With this, we should have almost exactly the same behavior
+    // as ReverseArcStaticGraph().
+    fix_final_permutation = true;
+    for (int arc = 0; arc < num_arcs_; ++arc) {
+      const int image = num_arcs_ + arc;
+      std::swap(tmp_heads_[image], tmp_heads_[arc]);
+      std::swap(tmp_tails_[image], tmp_tails_[arc]);
+    }
+  }
+
+  num_arcs_ = tmp_heads_.size();
+  perm.resize(num_arcs_);
+
+  // Do we sort each OutgoingArcs(node) by head ?
+  // Or is it better to keep all new reverse arc before or after ?
+  //
+  // TODO(user): For now we only support sorting, or all new reverse after
+  // and keep the original arc order.
+  if (option_sort_by_head_) {
+    FillPermutationAndStart(tmp_tails_, tmp_heads_, absl::MakeSpan(perm));
+  } else {
+    FillPermutationAndStart(tmp_tails_, absl::MakeSpan(perm));
+  }
+
+  // Create the final heads_.
+  heads_ = std::make_unique<NodeIndexType[]>(num_arcs_);
+  internal::PermuteInto<ArcIndexType, NodeIndexType>(
+      perm, tmp_heads_, absl::MakeSpan(heads_.get(), num_arcs_));
+
+  // No longer needed.
+  gtl::STLClearObject(&tmp_tails_);
+  gtl::STLClearObject(&tmp_heads_);
+
+  // We now create reverse_, this needs the permutation on both sides.
+  reverse_ = std::make_unique<ArcIndexType[]>(num_arcs_);
+  for (int i = 0; i < num_arcs_; ++i) {
+    reverse_[perm[i]] = perm[reverse[i]];
+  }
+
+  if (permutation != nullptr) {
+    if (fix_final_permutation) {
+      for (int i = 0; i < num_arcs_ / 2; ++i) {
+        std::swap(perm[i], perm[num_arcs_ / 2 + i]);
+      }
+    }
+    permutation->swap(perm);
+  }
+
+  node_capacity_ = num_nodes_;
+  arc_capacity_ = num_arcs_;
+  this->FreezeCapacities();
+}
+
+}  // namespace util
+
+#endif  // UTIL_GRAPH_FLOW_GRAPH_H_
diff --git a/ortools/graph/generic_max_flow.h b/ortools/graph/generic_max_flow.h
index 5c441da43c..bd068ed29d 100644
--- a/ortools/graph/generic_max_flow.h
+++ b/ortools/graph/generic_max_flow.h
@@ -124,6 +124,7 @@
 #define OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_
 
 #include <algorithm>
+#include <cstdint>
 #include <limits>
 #include <memory>
 #include <string>
@@ -229,13 +230,12 @@ class MaxFlowStatusClass {
 //    reverse can have a positive initial capacity. This can lead to twice fewer
 //    arcs and a faster algo if the "user graph" as a lot of reverse arcs
 //    already.
-template <typename Graph>
+template <typename Graph, typename FlowType = int64_t>
 class GenericMaxFlow : public MaxFlowStatusClass {
  public:
   typedef typename Graph::NodeIndex NodeIndex;
   typedef typename Graph::ArcIndex ArcIndex;
-  // TODO(b/385094969): This should be a template parameter.
-  typedef FlowQuantity FlowQuantityT;
+  typedef FlowType FlowQuantityT;
 
   // The height of a node never excess 2 times the number of node, so we
   // use the same type as a Node index.
@@ -273,20 +273,20 @@ class GenericMaxFlow : public MaxFlowStatusClass {
   //
   // Note that this will be ignored for self-arc, so do not be surprised if you
   // get zero when reading the capacity of a self-arc back.
-  void SetArcCapacity(ArcIndex arc, FlowQuantity new_capacity);
+  void SetArcCapacity(ArcIndex arc, FlowType new_capacity);
 
   // Returns true if a maximum flow was solved.
   bool Solve();
 
   // Returns the total flow found by the algorithm.
-  FlowQuantity GetOptimalFlow() const { return node_excess_[sink_]; }
+  FlowType GetOptimalFlow() const { return node_excess_[sink_]; }
 
   // Returns the current flow on the given arc.
   //
   // Note that on a couple (arc, opposite_arc) the flow only goes in one
   // direction (where it is positive) and the other direction will have the
   // negation of that flow.
-  FlowQuantity Flow(ArcIndex arc) const {
+  FlowType Flow(ArcIndex arc) const {
     if constexpr (Graph::kHasNegativeReverseArcs) {
       if (IsArcDirect(arc)) return residual_arc_capacity_[Opposite(arc)];
       return -residual_arc_capacity_[arc];
@@ -295,7 +295,7 @@ class GenericMaxFlow : public MaxFlowStatusClass {
   }
 
   // Returns the initial capacity of an arc.
-  FlowQuantity Capacity(ArcIndex arc) const {
+  FlowType Capacity(ArcIndex arc) const {
     if constexpr (Graph::kHasNegativeReverseArcs) {
       if (!IsArcDirect(arc)) return 0;
       return residual_arc_capacity_[arc] +
@@ -401,19 +401,19 @@ class GenericMaxFlow : public MaxFlowStatusClass {
 
   // Tries to saturate all the outgoing arcs from the source that can reach the
   // sink. Most of the time, we can do that in one go, except when more flow
-  // than kMaxFlowQuantity can be pushed out of the source in which case we
-  // have to be careful. Returns true if some flow was pushed.
+  // than kMaxFlowQuantity can be pushed out of the source in which case we have
+  // to be careful. Returns true if some flow was pushed.
   bool SaturateOutgoingArcsFromSource();
 
   // Pushes flow on arc,  i.e. consumes flow on residual_arc_capacity_[arc],
   // and consumes -flow on residual_arc_capacity_[Opposite(arc)]. Updates
   // node_excess_ at the tail and head of arc accordingly.
-  void PushFlow(FlowQuantity flow, NodeIndex tail, ArcIndex arc);
+  void PushFlow(FlowType flow, NodeIndex tail, ArcIndex arc);
 
   // Same as PushFlow() but with a cached node_excess_.data(), this is faster
   // in hot loop as we remove bound checking and the pointer is constant.
-  void FastPushFlow(FlowQuantity flow, NodeIndex tail, ArcIndex arc,
-                    FlowQuantity* node_excess);
+  void FastPushFlow(FlowType flow, NodeIndex tail, ArcIndex arc,
+                    FlowType* node_excess);
 
   // Relabels a node, i.e. increases its height by the minimum necessary amount.
   // This version of Relabel is relaxed in a way such that if an admissible arc
@@ -435,14 +435,14 @@ class GenericMaxFlow : public MaxFlowStatusClass {
   void ComputeReachableNodes(NodeIndex start, std::vector<NodeIndex>* result);
 
   // Maximum manageable flow.
-  static constexpr FlowQuantity kMaxFlowQuantity =
-      std::numeric_limits<FlowQuantity>::max();
+  static constexpr FlowType kMaxFlowQuantity =
+      std::numeric_limits<FlowType>::max();
 
   // A pointer to the graph passed as argument.
   const Graph* graph_;
 
   // An array representing the excess for each node in graph_.
-  std::unique_ptr<FlowQuantity[]> node_excess_;
+  std::unique_ptr<FlowType[]> node_excess_;
 
   // An array representing the height function for each node in graph_. For a
   // given node, this is a lower bound on the shortest path length from this
@@ -473,12 +473,12 @@ class GenericMaxFlow : public MaxFlowStatusClass {
   // Using these facts enables one to only maintain residual_arc_capacity_,
   // instead of both capacity and flow, for each direct and indirect arc. This
   // reduces the amount of memory for this information by a factor 2.
-  ZVector<FlowQuantity> residual_arc_capacity_;
+  ZVector<FlowType> residual_arc_capacity_;
 
   // The initial capacity as set by SetArcCapacity(), unused if
   // `Graph::kNegativeReverseArcs`, as we can always recover the initial
   // capacity from the residual capacities of an arc and its reverse.
-  std::unique_ptr<FlowQuantity[]> initial_capacity_;
+  std::unique_ptr<FlowType[]> initial_capacity_;
 
   // An array representing the first admissible arc for each node in graph_.
   std::unique_ptr<ArcIndex[]> first_admissible_arc_;
@@ -560,9 +560,10 @@ Element PriorityQueueWithRestrictedPush<Element, IntegerPriority>::PopBack(
   return element;
 }
 
-template <typename Graph>
-GenericMaxFlow<Graph>::GenericMaxFlow(const Graph* graph, NodeIndex source,
-                                      NodeIndex sink)
+template <typename Graph, typename FlowType>
+GenericMaxFlow<Graph, FlowType>::GenericMaxFlow(const Graph* graph,
+                                                NodeIndex source,
+                                                NodeIndex sink)
     : graph_(graph),
       residual_arc_capacity_(),
       source_(source),
@@ -577,7 +578,7 @@ GenericMaxFlow<Graph>::GenericMaxFlow(const Graph* graph, NodeIndex source,
     //
     // TODO(user): Unfortunately std/absl::make_unique_for_overwrite is not
     // yet available in open-source.
-    node_excess_ = std::make_unique<FlowQuantity[]>(max_num_nodes);
+    node_excess_ = std::make_unique<FlowType[]>(max_num_nodes);
     node_potential_ = std::make_unique<NodeHeight[]>(max_num_nodes);
     first_admissible_arc_ = std::make_unique<ArcIndex[]>(max_num_nodes);
     bfs_queue_.reserve(max_num_nodes);
@@ -588,16 +589,16 @@ GenericMaxFlow<Graph>::GenericMaxFlow(const Graph* graph, NodeIndex source,
       residual_arc_capacity_.Reserve(-max_num_arcs, max_num_arcs - 1);
     } else {
       // We will need to store the initial capacity in this case.
-      initial_capacity_ = std::make_unique<FlowQuantity[]>(max_num_arcs);
+      initial_capacity_ = std::make_unique<FlowType[]>(max_num_arcs);
       residual_arc_capacity_.Reserve(0, max_num_arcs - 1);
     }
     residual_arc_capacity_.SetAll(0);
   }
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::SetArcCapacity(ArcIndex arc,
-                                           FlowQuantity new_capacity) {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::SetArcCapacity(ArcIndex arc,
+                                                     FlowType new_capacity) {
   SCOPED_TIME_STAT(&stats_);
   DCHECK_LE(0, new_capacity);
   DCHECK(IsArcDirect(arc));
@@ -618,19 +619,20 @@ void GenericMaxFlow<Graph>::SetArcCapacity(ArcIndex arc,
   }
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::GetSourceSideMinCut(
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::GetSourceSideMinCut(
     std::vector<NodeIndex>* result) {
   ComputeReachableNodes<false>(source_, result);
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::GetSinkSideMinCut(std::vector<NodeIndex>* result) {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::GetSinkSideMinCut(
+    std::vector<NodeIndex>* result) {
   ComputeReachableNodes<true>(sink_, result);
 }
 
-template <typename Graph>
-bool GenericMaxFlow<Graph>::CheckResult() const {
+template <typename Graph, typename FlowType>
+bool GenericMaxFlow<Graph, FlowType>::CheckResult() const {
   SCOPED_TIME_STAT(&stats_);
   if (node_excess_[source_] != -node_excess_[sink_]) {
     LOG(DFATAL) << "-node_excess_[source_] = " << -node_excess_[source_]
@@ -648,8 +650,8 @@ bool GenericMaxFlow<Graph>::CheckResult() const {
   }
   for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
     const ArcIndex opposite = Opposite(arc);
-    const FlowQuantity direct_capacity = residual_arc_capacity_[arc];
-    const FlowQuantity opposite_capacity = residual_arc_capacity_[opposite];
+    const FlowType direct_capacity = residual_arc_capacity_[arc];
+    const FlowType opposite_capacity = residual_arc_capacity_[opposite];
     if (direct_capacity < 0) {
       LOG(DFATAL) << "residual_arc_capacity_[" << arc
                   << "] = " << direct_capacity << " < 0";
@@ -676,8 +678,8 @@ bool GenericMaxFlow<Graph>::CheckResult() const {
   return true;
 }
 
-template <typename Graph>
-bool GenericMaxFlow<Graph>::AugmentingPathExists() const {
+template <typename Graph, typename FlowType>
+bool GenericMaxFlow<Graph, FlowType>::AugmentingPathExists() const {
   SCOPED_TIME_STAT(&stats_);
 
   // We simply compute the reachability from the source in the residual graph.
@@ -703,8 +705,9 @@ bool GenericMaxFlow<Graph>::AugmentingPathExists() const {
   return is_reached[sink_];
 }
 
-template <typename Graph>
-bool GenericMaxFlow<Graph>::CheckRelabelPrecondition(NodeIndex node) const {
+template <typename Graph, typename FlowType>
+bool GenericMaxFlow<Graph, FlowType>::CheckRelabelPrecondition(
+    NodeIndex node) const {
   DCHECK(IsActive(node));
   for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) {
     DCHECK(!IsAdmissible(node, arc, node_potential_.data()))
@@ -713,9 +716,9 @@ bool GenericMaxFlow<Graph>::CheckRelabelPrecondition(NodeIndex node) const {
   return true;
 }
 
-template <typename Graph>
-std::string GenericMaxFlow<Graph>::DebugString(absl::string_view context,
-                                               ArcIndex arc) const {
+template <typename Graph, typename FlowType>
+std::string GenericMaxFlow<Graph, FlowType>::DebugString(
+    absl::string_view context, ArcIndex arc) const {
   const NodeIndex tail = Tail(arc);
   const NodeIndex head = Head(arc);
   return absl::StrFormat(
@@ -729,8 +732,8 @@ std::string GenericMaxFlow<Graph>::DebugString(absl::string_view context,
       node_potential_[head], node_excess_[tail], node_excess_[head]);
 }
 
-template <typename Graph>
-bool GenericMaxFlow<Graph>::Solve() {
+template <typename Graph, typename FlowType>
+bool GenericMaxFlow<Graph, FlowType>::Solve() {
   status_ = NOT_SOLVED;
   InitializePreflow();
 
@@ -758,8 +761,8 @@ bool GenericMaxFlow<Graph>::Solve() {
   return true;
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::InitializePreflow() {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::InitializePreflow() {
   SCOPED_TIME_STAT(&stats_);
 
   // TODO(user): Ebert graph has an issue with nodes with no arcs, so we
@@ -807,8 +810,8 @@ void GenericMaxFlow<Graph>::InitializePreflow() {
 // potentials because of the way we cancel flow on cycle. However, we only call
 // that at the end of the algorithm, or just before a GlobalUpdate() that will
 // restore the precondition on the node potentials.
-template <typename Graph>
-void GenericMaxFlow<Graph>::PushFlowExcessBackToSource() {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::PushFlowExcessBackToSource() {
   SCOPED_TIME_STAT(&stats_);
   const NodeIndex num_nodes = graph_->num_nodes();
 
@@ -877,7 +880,7 @@ void GenericMaxFlow<Graph>::PushFlowExcessBackToSource() {
     visited[node] = true;
     index_branch.push_back(arc_stack.size() - 1);
     for (const ArcIndex arc : graph_->OutgoingArcs(node)) {
-      const FlowQuantity flow = Flow(arc);
+      const FlowType flow = Flow(arc);
       const NodeIndex head = Head(arc);
       if (flow > 0 && !stored[head]) {
         if (!visited[head]) {
@@ -895,11 +898,11 @@ void GenericMaxFlow<Graph>::PushFlowExcessBackToSource() {
 
           // Compute the maximum flow that can be canceled on the cycle and the
           // min index such that arc_stack[index_branch[i]] will be saturated.
-          FlowQuantity flow_on_cycle = flow;
+          FlowType flow_on_cycle = flow;
           int first_saturated_index = index_branch.size();
           for (int i = index_branch.size() - 1; i >= cycle_begin; --i) {
             const ArcIndex arc_on_cycle = arc_stack[index_branch[i]];
-            if (const FlowQuantity arc_flow = Flow(arc_on_cycle);
+            if (const FlowType arc_flow = Flow(arc_on_cycle);
                 arc_flow <= flow_on_cycle) {
               flow_on_cycle = arc_flow;
               first_saturated_index = i;
@@ -907,7 +910,7 @@ void GenericMaxFlow<Graph>::PushFlowExcessBackToSource() {
           }
 
           // This is just here for a DCHECK() below.
-          const FlowQuantity excess = node_excess_[head];
+          const FlowType excess = node_excess_[head];
 
           // Cancel the flow on the cycle, and set visited[node] = false for
           // the node that will be backtracked over.
@@ -949,10 +952,10 @@ void GenericMaxFlow<Graph>::PushFlowExcessBackToSource() {
     const NodeIndex node = reverse_topological_order[i];
     if (node_excess_[node] == 0) continue;
     for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) {
-      const FlowQuantity flow = Flow(arc);
+      const FlowType flow = Flow(arc);
       if (flow < 0) {
         DCHECK_GT(residual_arc_capacity_[arc], 0);
-        const FlowQuantity to_push = std::min(node_excess_[node], -flow);
+        const FlowType to_push = std::min(node_excess_[node], -flow);
         PushFlow(to_push, node, arc);
         if (node_excess_[node] == 0) break;
       }
@@ -962,8 +965,8 @@ void GenericMaxFlow<Graph>::PushFlowExcessBackToSource() {
   DCHECK_EQ(-node_excess_[source_], node_excess_[sink_]);
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::GlobalUpdate() {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::GlobalUpdate() {
   SCOPED_TIME_STAT(&stats_);
   bfs_queue_.clear();
   int queue_index = 0;
@@ -972,7 +975,7 @@ void GenericMaxFlow<Graph>::GlobalUpdate() {
   node_in_bfs_queue_[sink_] = true;
 
   // We cache this as this is a hot-loop and we don't want bound checking.
-  FlowQuantity* const node_excess = node_excess_.get();
+  FlowType* const node_excess = node_excess_.get();
 
   // We do a BFS in the reverse residual graph, starting from the sink.
   // Because all the arcs from the source are saturated (except in
@@ -1015,7 +1018,7 @@ void GenericMaxFlow<Graph>::GlobalUpdate() {
         // Note(user): I haven't seen this anywhere in the literature.
         // TODO(user): Investigate more and maybe write a publication :)
         if (node_excess[head] > 0) {
-          const FlowQuantity flow =
+          const FlowType flow =
               std::min(node_excess[head], residual_arc_capacity_[opposite_arc]);
           FastPushFlow(flow, head, opposite_arc, node_excess);
 
@@ -1062,8 +1065,8 @@ void GenericMaxFlow<Graph>::GlobalUpdate() {
   }
 }
 
-template <typename Graph>
-bool GenericMaxFlow<Graph>::SaturateOutgoingArcsFromSource() {
+template <typename Graph, typename FlowType>
+bool GenericMaxFlow<Graph, FlowType>::SaturateOutgoingArcsFromSource() {
   SCOPED_TIME_STAT(&stats_);
   const NodeIndex num_nodes = graph_->num_nodes();
 
@@ -1074,18 +1077,17 @@ bool GenericMaxFlow<Graph>::SaturateOutgoingArcsFromSource() {
 
   bool flow_pushed = false;
   for (const ArcIndex arc : graph_->OutgoingArcs(source_)) {
-    const FlowQuantity flow = residual_arc_capacity_[arc];
+    const FlowType flow = residual_arc_capacity_[arc];
 
     // This is a special IsAdmissible() condition for the source.
     if (flow == 0 || node_potential_[Head(arc)] >= num_nodes) continue;
 
     // We are careful in case the sum of the flow out of the source is greater
     // than kMaxFlowQuantity to avoid overflow.
-    const FlowQuantity current_flow_out_of_source = -node_excess_[source_];
+    const FlowType current_flow_out_of_source = -node_excess_[source_];
     DCHECK_GE(flow, 0) << flow;
     DCHECK_GE(current_flow_out_of_source, 0) << current_flow_out_of_source;
-    const FlowQuantity capped_flow =
-        kMaxFlowQuantity - current_flow_out_of_source;
+    const FlowType capped_flow = kMaxFlowQuantity - current_flow_out_of_source;
     if (capped_flow < flow) {
       // We push as much flow as we can so the current flow on the network will
       // be kMaxFlowQuantity.
@@ -1104,9 +1106,9 @@ bool GenericMaxFlow<Graph>::SaturateOutgoingArcsFromSource() {
   return flow_pushed;
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::PushFlow(FlowQuantity flow, NodeIndex tail,
-                                     ArcIndex arc) {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::PushFlow(FlowType flow, NodeIndex tail,
+                                               ArcIndex arc) {
   SCOPED_TIME_STAT(&stats_);
 
   // Update the residual capacity of the arc and its opposite arc.
@@ -1127,10 +1129,10 @@ void GenericMaxFlow<Graph>::PushFlow(FlowQuantity flow, NodeIndex tail,
   node_excess_[Head(arc)] += flow;
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::FastPushFlow(FlowQuantity flow, NodeIndex tail,
-                                         ArcIndex arc,
-                                         FlowQuantity* node_excess) {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::FastPushFlow(FlowType flow,
+                                                   NodeIndex tail, ArcIndex arc,
+                                                   FlowType* node_excess) {
   SCOPED_TIME_STAT(&stats_);
 
   DCHECK_NE(flow, 0);
@@ -1143,8 +1145,8 @@ void GenericMaxFlow<Graph>::FastPushFlow(FlowQuantity flow, NodeIndex tail,
   node_excess[Head(arc)] += flow;
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::InitializeActiveNodeContainer() {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::InitializeActiveNodeContainer() {
   SCOPED_TIME_STAT(&stats_);
   DCHECK(IsEmptyActiveNodeContainer());
   const NodeIndex num_nodes = graph_->num_nodes();
@@ -1159,8 +1161,8 @@ void GenericMaxFlow<Graph>::InitializeActiveNodeContainer() {
   }
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::RefineWithGlobalUpdate() {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::RefineWithGlobalUpdate() {
   SCOPED_TIME_STAT(&stats_);
 
   const NodeIndex num_nodes = graph_->num_nodes();
@@ -1227,13 +1229,13 @@ void GenericMaxFlow<Graph>::RefineWithGlobalUpdate() {
   }
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::Discharge(const NodeIndex node) {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::Discharge(const NodeIndex node) {
   SCOPED_TIME_STAT(&stats_);
   const NodeIndex num_nodes = graph_->num_nodes();
 
   // We cache this as this is a hot-loop and we don't want bound checking.
-  FlowQuantity* const node_excess = node_excess_.get();
+  FlowType* const node_excess = node_excess_.get();
   NodeHeight* const node_potentials = node_potential_.get();
   ArcIndex* const first_admissible_arc = first_admissible_arc_.get();
 
@@ -1250,7 +1252,7 @@ void GenericMaxFlow<Graph>::Discharge(const NodeIndex node) {
           // push the sink_, but that is handled properly in Refine().
           PushActiveNode(head);
         }
-        const FlowQuantity delta =
+        const FlowType delta =
             std::min(node_excess[node], residual_arc_capacity_[arc]);
         FastPushFlow(delta, node, arc, node_excess);
         if (node_excess[node] == 0) {
@@ -1267,8 +1269,8 @@ void GenericMaxFlow<Graph>::Discharge(const NodeIndex node) {
   }
 }
 
-template <typename Graph>
-void GenericMaxFlow<Graph>::Relabel(NodeIndex node) {
+template <typename Graph, typename FlowType>
+void GenericMaxFlow<Graph, FlowType>::Relabel(NodeIndex node) {
   SCOPED_TIME_STAT(&stats_);
   // Because we use a relaxed version, this is no longer true if the
   // first_admissible_arc_[node] was not actually the first arc!
@@ -1298,24 +1300,25 @@ void GenericMaxFlow<Graph>::Relabel(NodeIndex node) {
   first_admissible_arc_[node] = first_admissible_arc;
 }
 
-template <typename Graph>
-typename Graph::ArcIndex GenericMaxFlow<Graph>::Opposite(ArcIndex arc) const {
+template <typename Graph, typename FlowType>
+typename Graph::ArcIndex GenericMaxFlow<Graph, FlowType>::Opposite(
+    ArcIndex arc) const {
   return graph_->OppositeArc(arc);
 }
 
-template <typename Graph>
-bool GenericMaxFlow<Graph>::IsArcDirect(ArcIndex arc) const {
+template <typename Graph, typename FlowType>
+bool GenericMaxFlow<Graph, FlowType>::IsArcDirect(ArcIndex arc) const {
   return IsArcValid(arc) && arc >= 0;
 }
 
-template <typename Graph>
-bool GenericMaxFlow<Graph>::IsArcValid(ArcIndex arc) const {
+template <typename Graph, typename FlowType>
+bool GenericMaxFlow<Graph, FlowType>::IsArcValid(ArcIndex arc) const {
   return graph_->IsArcValid(arc);
 }
 
-template <typename Graph>
+template <typename Graph, typename FlowType>
 template <bool reverse>
-void GenericMaxFlow<Graph>::ComputeReachableNodes(
+void GenericMaxFlow<Graph, FlowType>::ComputeReachableNodes(
     NodeIndex start, std::vector<NodeIndex>* result) {
   // If start is not a valid node index, it can reach only itself.
   // Note(user): This is needed because source and sink are given independently
@@ -1346,8 +1349,8 @@ void GenericMaxFlow<Graph>::ComputeReachableNodes(
   *result = bfs_queue_;
 }
 
-template <typename Graph>
-FlowModelProto GenericMaxFlow<Graph>::CreateFlowModel() {
+template <typename Graph, typename FlowType>
+FlowModelProto GenericMaxFlow<Graph, FlowType>::CreateFlowModel() {
   FlowModelProto model;
   model.set_problem_type(FlowModelProto::MAX_FLOW);
   for (int n = 0; n < graph_->num_nodes(); ++n) {
diff --git a/ortools/graph/generic_max_flow_test.cc b/ortools/graph/generic_max_flow_test.cc
index 9f42eb2e69..99d22fead8 100644
--- a/ortools/graph/generic_max_flow_test.cc
+++ b/ortools/graph/generic_max_flow_test.cc
@@ -32,6 +32,7 @@
 #include "ortools/base/gmock.h"
 #include "ortools/base/logging.h"
 #include "ortools/graph/ebert_graph.h"
+#include "ortools/graph/flow_graph.h"
 #include "ortools/graph/graph.h"
 #include "ortools/linear_solver/linear_solver.h"
 
@@ -98,7 +99,7 @@ typename GenericMaxFlow<Graph>::Status MaxFlowTester(
 template <typename Graph>
 class GenericMaxFlowTest : public ::testing::Test {};
 
-typedef ::testing::Types<util::ReverseArcListGraph<>,
+typedef ::testing::Types<util::FlowGraph<>, util::ReverseArcListGraph<>,
                          util::ReverseArcStaticGraph<>,
                          util::ReverseArcMixedGraph<>>
     GraphTypes;
@@ -556,7 +557,11 @@ void FullRandomFlow(std::optional<FlowQuantity> expected_flow,
   TEST(MaxFlowListGraphTest, test_name##size) {                                \
     test_name<util::ReverseArcListGraph<>>(std::nullopt, SolveMaxFlow, size,   \
                                            size);                              \
+  }                                                                            \
+  TEST(MaxFlowNewGraphTest, test_name##size) {                                 \
+    test_name<util::FlowGraph<>>(std::nullopt, SolveMaxFlow, size, size);      \
   }
+
 // These are absl::BitGen random test, so they will always work on different
 // graphs.
 LP_AND_FLOW_TEST(FullAssignment, 300);
@@ -602,18 +607,22 @@ static void BM_FullRandomFlow(benchmark::State& state) {
 }
 
 // Note that these benchmark include the graph creation and generation...
+BENCHMARK_TEMPLATE(BM_FullRandomAssignment, util::FlowGraph<>);
 BENCHMARK_TEMPLATE(BM_FullRandomAssignment, util::ReverseArcListGraph<>);
 BENCHMARK_TEMPLATE(BM_FullRandomAssignment, util::ReverseArcStaticGraph<>);
 BENCHMARK_TEMPLATE(BM_FullRandomAssignment, util::ReverseArcMixedGraph<>);
 
+BENCHMARK_TEMPLATE(BM_PartialRandomFlow, util::FlowGraph<>);
 BENCHMARK_TEMPLATE(BM_PartialRandomFlow, util::ReverseArcListGraph<>);
 BENCHMARK_TEMPLATE(BM_PartialRandomFlow, util::ReverseArcStaticGraph<>);
 BENCHMARK_TEMPLATE(BM_PartialRandomFlow, util::ReverseArcMixedGraph<>);
 
+BENCHMARK_TEMPLATE(BM_FullRandomFlow, util::FlowGraph<>);
 BENCHMARK_TEMPLATE(BM_FullRandomFlow, util::ReverseArcListGraph<>);
 BENCHMARK_TEMPLATE(BM_FullRandomFlow, util::ReverseArcStaticGraph<>);
 BENCHMARK_TEMPLATE(BM_FullRandomFlow, util::ReverseArcMixedGraph<>);
 
+BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, util::FlowGraph<>);
 BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, util::ReverseArcListGraph<>);
 BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, util::ReverseArcStaticGraph<>);
 BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, util::ReverseArcMixedGraph<>);
diff --git a/ortools/graph/graph.h b/ortools/graph/graph.h
index 89e6dba2ac..57398b0556 100644
--- a/ortools/graph/graph.h
+++ b/ortools/graph/graph.h
@@ -79,6 +79,8 @@
 // It is however possible to do so for the outgoing arcs and the opposite
 // incoming arcs. It is why the OutgoingOrOppositeIncomingArcs() and
 // OutgoingArcs() iterations are more efficient than the IncomingArcs() one.
+// TODO(b/385094969): Once we no longer use `Next()/Ok()` for iterators, we can
+// get rid of `limit_`, which will make iteration much more efficient.
 //
 // If you know the graph size in advance, this already set the number of nodes,
 // reserve space for the arcs and check in DEBUG mode that you don't go over the
@@ -188,7 +190,7 @@ class SVector;
 // Note: The type can be unsigned, except for the graphs with reverse arcs
 // where the ArcIndexType must be signed, but not necessarily the NodeIndexType.
 template <typename NodeIndexType = int32_t, typename ArcIndexType = int32_t,
-          bool HasReverseArcs = false>
+          bool HasNegativeReverseArcs = false>
 class BaseGraph {
  public:
   // Typedef so you can use Graph::NodeIndex and Graph::ArcIndex to be generic
@@ -196,7 +198,7 @@ class BaseGraph {
   // that you define a typedef ... Graph; for readability.
   typedef NodeIndexType NodeIndex;
   typedef ArcIndexType ArcIndex;
-  static constexpr bool kHasNegativeReverseArcs = HasReverseArcs;
+  static constexpr bool kHasNegativeReverseArcs = HasNegativeReverseArcs;
 
   BaseGraph()
       : num_nodes_(0),
@@ -231,7 +233,7 @@ class BaseGraph {
   // Returns true if the given arc is a valid arc of the graph.
   // Note that the arc validity range changes for graph with reverse arcs.
   bool IsArcValid(ArcIndexType arc) const {
-    return (HasReverseArcs ? -num_arcs_ : 0) <= arc && arc < num_arcs_;
+    return (HasNegativeReverseArcs ? -num_arcs_ : 0) <= arc && arc < num_arcs_;
   }
 
   // Capacity reserved for future nodes, always >= num_nodes_.
@@ -960,46 +962,54 @@ class SVector {
 
 // BaseGraph implementation ----------------------------------------------------
 
-template <typename NodeIndexType, typename ArcIndexType, bool HasReverseArcs>
-IntegerRange<NodeIndexType>
-BaseGraph<NodeIndexType, ArcIndexType, HasReverseArcs>::AllNodes() const {
+template <typename NodeIndexType, typename ArcIndexType,
+          bool HasNegativeReverseArcs>
+IntegerRange<NodeIndexType> BaseGraph<
+    NodeIndexType, ArcIndexType, HasNegativeReverseArcs>::AllNodes() const {
   return IntegerRange<NodeIndexType>(0, num_nodes_);
 }
 
-template <typename NodeIndexType, typename ArcIndexType, bool HasReverseArcs>
+template <typename NodeIndexType, typename ArcIndexType,
+          bool HasNegativeReverseArcs>
 IntegerRange<ArcIndexType>
-BaseGraph<NodeIndexType, ArcIndexType, HasReverseArcs>::AllForwardArcs() const {
+BaseGraph<NodeIndexType, ArcIndexType, HasNegativeReverseArcs>::AllForwardArcs()
+    const {
   return IntegerRange<ArcIndexType>(0, num_arcs_);
 }
 
-template <typename NodeIndexType, typename ArcIndexType, bool HasReverseArcs>
+template <typename NodeIndexType, typename ArcIndexType,
+          bool HasNegativeReverseArcs>
 const NodeIndexType
-    BaseGraph<NodeIndexType, ArcIndexType, HasReverseArcs>::kNilNode =
+    BaseGraph<NodeIndexType, ArcIndexType, HasNegativeReverseArcs>::kNilNode =
         std::numeric_limits<NodeIndexType>::max();
 
-template <typename NodeIndexType, typename ArcIndexType, bool HasReverseArcs>
+template <typename NodeIndexType, typename ArcIndexType,
+          bool HasNegativeReverseArcs>
 const ArcIndexType
-    BaseGraph<NodeIndexType, ArcIndexType, HasReverseArcs>::kNilArc =
+    BaseGraph<NodeIndexType, ArcIndexType, HasNegativeReverseArcs>::kNilArc =
         std::numeric_limits<ArcIndexType>::max();
 
-template <typename NodeIndexType, typename ArcIndexType, bool HasReverseArcs>
-NodeIndexType
-BaseGraph<NodeIndexType, ArcIndexType, HasReverseArcs>::node_capacity() const {
+template <typename NodeIndexType, typename ArcIndexType,
+          bool HasNegativeReverseArcs>
+NodeIndexType BaseGraph<NodeIndexType, ArcIndexType,
+                        HasNegativeReverseArcs>::node_capacity() const {
   // TODO(user): Is it needed? remove completely? return the real capacities
   // at the cost of having a different implementation for each graphs?
   return node_capacity_ > num_nodes_ ? node_capacity_ : num_nodes_;
 }
 
-template <typename NodeIndexType, typename ArcIndexType, bool HasReverseArcs>
-ArcIndexType
-BaseGraph<NodeIndexType, ArcIndexType, HasReverseArcs>::arc_capacity() const {
+template <typename NodeIndexType, typename ArcIndexType,
+          bool HasNegativeReverseArcs>
+ArcIndexType BaseGraph<NodeIndexType, ArcIndexType,
+                       HasNegativeReverseArcs>::arc_capacity() const {
   // TODO(user): Same questions as the ones in node_capacity().
   return arc_capacity_ > num_arcs_ ? arc_capacity_ : num_arcs_;
 }
 
-template <typename NodeIndexType, typename ArcIndexType, bool HasReverseArcs>
+template <typename NodeIndexType, typename ArcIndexType,
+          bool HasNegativeReverseArcs>
 void BaseGraph<NodeIndexType, ArcIndexType,
-               HasReverseArcs>::FreezeCapacities() {
+               HasNegativeReverseArcs>::FreezeCapacities() {
   // TODO(user): Only define this in debug mode at the cost of having a lot
   // of ifndef NDEBUG all over the place? remove the function completely ?
   const_capacities_ = true;
@@ -1009,8 +1019,9 @@ void BaseGraph<NodeIndexType, ArcIndexType,
 
 // Computes the cumulative sum of the entry in v. We only use it with
 // in/out degree distribution, hence the Check() at the end.
-template <typename NodeIndexType, typename ArcIndexType, bool HasReverseArcs>
-void BaseGraph<NodeIndexType, ArcIndexType, HasReverseArcs>::
+template <typename NodeIndexType, typename ArcIndexType,
+          bool HasNegativeReverseArcs>
+void BaseGraph<NodeIndexType, ArcIndexType, HasNegativeReverseArcs>::
     ComputeCumulativeSum(std::vector<ArcIndexType>* v) {
   ArcIndexType sum = 0;
   for (int i = 0; i < num_nodes_; ++i) {
@@ -1026,8 +1037,9 @@ void BaseGraph<NodeIndexType, ArcIndexType, HasReverseArcs>::
 // - Put the head of the new arc #i in (*head)[i].
 // - Put in start[i] the index of the first arc with tail >= i.
 // - Update "permutation" to reflect the change, unless it is NULL.
-template <typename NodeIndexType, typename ArcIndexType, bool HasReverseArcs>
-void BaseGraph<NodeIndexType, ArcIndexType, HasReverseArcs>::
+template <typename NodeIndexType, typename ArcIndexType,
+          bool HasNegativeReverseArcs>
+void BaseGraph<NodeIndexType, ArcIndexType, HasNegativeReverseArcs>::
     BuildStartAndForwardHead(SVector<NodeIndexType>* head,
                              std::vector<ArcIndexType>* start,
                              std::vector<ArcIndexType>* permutation) {
diff --git a/ortools/graph/max_flow_test.cc b/ortools/graph/max_flow_test.cc
index ed6c7ca854..51c3f546f4 100644
--- a/ortools/graph/max_flow_test.cc
+++ b/ortools/graph/max_flow_test.cc
@@ -60,7 +60,8 @@ TEST(SimpleMaxFlowTest, EmptyWithInvalidSourceAndSink) {
 
 TEST(SimpleMaxFlowTest, TrivialGraphWithMaxCapacity) {
   SimpleMaxFlow max_flow;
-  const FlowQuantity kCapacityMax = std::numeric_limits<FlowQuantity>::max();
+  const SimpleMaxFlow::FlowQuantity kCapacityMax =
+      std::numeric_limits<SimpleMaxFlow::FlowQuantity>::max();
   EXPECT_EQ(0, max_flow.AddArcWithCapacity(0, 1, kCapacityMax));
   EXPECT_EQ(SimpleMaxFlow::OPTIMAL, max_flow.Solve(1, 0));
   EXPECT_EQ(0, max_flow.OptimalFlow());
diff --git a/ortools/graph/min_cost_flow.cc b/ortools/graph/min_cost_flow.cc
index 8e09a8e3bc..71e0e50096 100644
--- a/ortools/graph/min_cost_flow.cc
+++ b/ortools/graph/min_cost_flow.cc
@@ -53,6 +53,9 @@ GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::GenericMinCostFlow(
       alpha_(absl::GetFlag(FLAGS_min_cost_flow_alpha)),
       stats_("MinCostFlow"),
       check_feasibility_(absl::GetFlag(FLAGS_min_cost_flow_check_feasibility)) {
+  // This class assumes we have negative reverse arcs.
+  static_assert(Graph::kHasNegativeReverseArcs);
+
   const NodeIndex max_num_nodes = graph_->node_capacity();
   if (max_num_nodes > 0) {
     first_admissible_arc_.assign(max_num_nodes, Graph::kNilArc);
diff --git a/ortools/graph/python/CMakeLists.txt b/ortools/graph/python/CMakeLists.txt
index a6a394eccc..85c0474037 100644
--- a/ortools/graph/python/CMakeLists.txt
+++ b/ortools/graph/python/CMakeLists.txt
@@ -62,3 +62,9 @@ endif()
 target_link_libraries(min_cost_flow_pybind11 PRIVATE ${PROJECT_NAMESPACE}::ortools)
 add_library(${PROJECT_NAMESPACE}::min_cost_flow_pybind11 ALIAS min_cost_flow_pybind11)
 
+if(BUILD_TESTING)
+  file(GLOB PYTHON_SRCS "*_test.py")
+  foreach(FILE_NAME IN LISTS PYTHON_SRCS)
+    add_python_test(${FILE_NAME})
+  endforeach()
+endif()
diff --git a/ortools/graph/python/max_flow.cc b/ortools/graph/python/max_flow.cc
index 8cdee52f2c..b9c04cebac 100644
--- a/ortools/graph/python/max_flow.cc
+++ b/ortools/graph/python/max_flow.cc
@@ -19,7 +19,6 @@
 #include "pybind11/pybind11.h"
 #include "pybind11/stl.h"
 
-using ::operations_research::NodeIndex;
 using ::operations_research::SimpleMaxFlow;
 using ::pybind11::arg;
 
@@ -44,12 +43,12 @@ PYBIND11_MODULE(max_flow, m) {
   smf.def("flow", &SimpleMaxFlow::Flow, arg("arc"));
   smf.def("flows", pybind11::vectorize(&SimpleMaxFlow::Flow));
   smf.def("get_source_side_min_cut", [](SimpleMaxFlow* smf) {
-    std::vector<NodeIndex> result;
+    std::vector<SimpleMaxFlow::NodeIndex> result;
     smf->GetSourceSideMinCut(&result);
     return result;
   });
   smf.def("get_sink_side_min_cut", [](SimpleMaxFlow* smf) {
-    std::vector<NodeIndex> result;
+    std::vector<SimpleMaxFlow::NodeIndex> result;
     smf->GetSinkSideMinCut(&result);
     return result;
   });
diff --git a/ortools/graph/solve_flow_model.cc b/ortools/graph/solve_flow_model.cc
index 64da0cd8da..deed8253d4 100644
--- a/ortools/graph/solve_flow_model.cc
+++ b/ortools/graph/solve_flow_model.cc
@@ -21,13 +21,15 @@
 
 #include <algorithm>
 #include <cstdint>
+#include <cstdio>
 #include <cstdlib>
+#include <functional>
 #include <sstream>
 #include <string>
 #include <vector>
 
 #include "absl/flags/flag.h"
-#include "absl/status/status.h"
+#include "absl/log/check.h"
 #include "absl/strings/match.h"
 #include "absl/strings/str_format.h"
 #include "absl/strings/string_view.h"
@@ -37,17 +39,24 @@
 #include "ortools/base/helpers.h"
 #include "ortools/base/init_google.h"
 #include "ortools/base/logging.h"
+#include "ortools/base/options.h"
 #include "ortools/base/path.h"
 #include "ortools/base/timer.h"
+#include "ortools/graph/flow_graph.h"
 #include "ortools/graph/flow_problem.pb.h"
+#include "ortools/graph/generic_max_flow.h"
 #include "ortools/graph/graph.h"
-#include "ortools/graph/max_flow.h"
 #include "ortools/graph/min_cost_flow.h"
+#include "ortools/util/file_util.h"
 #include "ortools/util/filelineiter.h"
 #include "ortools/util/stats.h"
 
 ABSL_FLAG(std::string, input, "", "Input file of the problem.");
 ABSL_FLAG(std::string, output_dimacs, "", "Output problem as a dimacs file.");
+ABSL_FLAG(std::string, output_proto, "", "Output problem as a flow proto.");
+ABSL_FLAG(bool, use_flow_graph, true, "Use special kind of graph.");
+ABSL_FLAG(bool, sort_heads, false, "Sort outgoing arcs by head.");
+ABSL_FLAG(bool, detect_reverse_arcs, true, "Detect reverse arcs.");
 
 namespace operations_research {
 
@@ -194,6 +203,8 @@ void SolveMinCostFlow(const FlowModelProto& flow_model, double* loading_time,
   }
   std::vector<Graph::ArcIndex> permutation;
   graph.Build(&permutation);
+  absl::PrintF("%d,", graph.num_nodes());
+  absl::PrintF("%d,", graph.num_arcs());
 
   GenericMinCostFlow<Graph> min_cost_flow(&graph);
   for (int i = 0; i < flow_model.arcs_size(); ++i) {
@@ -220,29 +231,32 @@ void SolveMinCostFlow(const FlowModelProto& flow_model, double* loading_time,
 }
 
 // Loads a FlowModelProto proto into the MaxFlow class and solves it.
-void SolveMaxFlow(const FlowModelProto& flow_model, double* loading_time,
-                  double* solving_time) {
+template <typename GraphType>
+void SolveMaxFlow(
+    const FlowModelProto& flow_model, double* loading_time,
+    double* solving_time,
+    std::function<void(GraphType* graph)> configure_graph_options = nullptr) {
   WallTimer timer;
   timer.Start();
 
-  // Compute the number of nodes.
-  int64_t num_nodes = 0;
-  for (int i = 0; i < flow_model.arcs_size(); ++i) {
-    num_nodes = std::max(num_nodes, flow_model.arcs(i).tail() + 1);
-    num_nodes = std::max(num_nodes, flow_model.arcs(i).head() + 1);
-  }
-
   // Build the graph.
-  Graph graph(num_nodes, flow_model.arcs_size());
+  GraphType graph(flow_model.nodes_size(), flow_model.arcs_size());
   for (int i = 0; i < flow_model.arcs_size(); ++i) {
     graph.AddArc(flow_model.arcs(i).tail(), flow_model.arcs(i).head());
   }
-  std::vector<Graph::ArcIndex> permutation;
+  std::vector<typename GraphType::ArcIndex> permutation;
+
+  if (configure_graph_options != nullptr) {
+    configure_graph_options(&graph);
+  }
   graph.Build(&permutation);
 
+  absl::PrintF("%d,", graph.num_nodes());
+  absl::PrintF("%d,", graph.num_arcs());
+
   // Find source & sink.
-  Graph::NodeIndex source = -1;
-  Graph::NodeIndex sink = -1;
+  typename GraphType::NodeIndex source = -1;
+  typename GraphType::NodeIndex sink = -1;
   CHECK_EQ(2, flow_model.nodes_size());
   for (int i = 0; i < flow_model.nodes_size(); ++i) {
     if (flow_model.nodes(i).supply() > 0) {
@@ -256,9 +270,10 @@ void SolveMaxFlow(const FlowModelProto& flow_model, double* loading_time,
   CHECK_NE(sink, -1);
 
   // Create the max flow instance and set the arc capacities.
-  GenericMaxFlow<Graph> max_flow(&graph, source, sink);
+  GenericMaxFlow<GraphType> max_flow(&graph, source, sink);
   for (int i = 0; i < flow_model.arcs_size(); ++i) {
-    const Graph::ArcIndex image = i < permutation.size() ? permutation[i] : i;
+    const typename GraphType::ArcIndex image =
+        i < permutation.size() ? permutation[i] : i;
     max_flow.SetArcCapacity(image, flow_model.arcs(i).capacity());
   }
 
@@ -268,7 +283,7 @@ void SolveMaxFlow(const FlowModelProto& flow_model, double* loading_time,
 
   timer.Start();
   CHECK(max_flow.Solve());
-  CHECK_EQ(GenericMaxFlow<Graph>::OPTIMAL, max_flow.status());
+  CHECK_EQ(GenericMaxFlow<GraphType>::OPTIMAL, max_flow.status());
   *solving_time = timer.Get();
   absl::PrintF("%f,", *solving_time);
   absl::PrintF("%d", max_flow.GetOptimalFlow());
@@ -301,7 +316,8 @@ int main(int argc, char** argv) {
   TimeDistribution solving_time_distribution("Solving time summary");
 
   absl::PrintF(
-      "file_name, parsing_time, loading_time, solving_time, optimal_cost\n");
+      "file_name, parsing_time, num_nodes, num_arcs,loading_time, "
+      "solving_time, optimal_cost\n");
   for (int i = 0; i < file_list.size(); ++i) {
     const std::string file_name = file_list[i];
     absl::PrintF("%s,", file::Basename(file_name));
@@ -310,11 +326,16 @@ int main(int argc, char** argv) {
     // Parse the input as a proto.
     double parsing_time = 0;
     operations_research::FlowModelProto proto;
-    if (absl::EndsWith(file_name, ".bin")) {
+    if (absl::EndsWith(file_name, ".bin") ||
+        absl::EndsWith(file_name, ".bin.gz")) {
       ScopedWallTime timer(&parsing_time);
-      std::string raw_data;
-      CHECK_OK(file::GetContents(file_name, &raw_data, file::Defaults()));
-      proto.ParseFromString(raw_data);
+      if (absl::EndsWith(file_name, "gz")) {
+        std::string raw_data;
+        CHECK_OK(file::GetContents(file_name, &raw_data, file::Defaults()));
+        CHECK_OK(StringToProto(raw_data, &proto));
+      } else {
+        CHECK_OK(file::GetBinaryProto(file_name, &proto, file::Defaults()));
+      }
     } else {
       ScopedWallTime timer(&parsing_time);
       if (!ConvertDimacsToFlowModel(file_name, &proto)) continue;
@@ -322,6 +343,13 @@ int main(int argc, char** argv) {
     absl::PrintF("%f,", parsing_time);
     fflush(stdout);
 
+    if (!absl::GetFlag(FLAGS_output_proto).empty()) {
+      LOG(INFO) << "Dumping binary proto to '"
+                << absl::GetFlag(FLAGS_output_proto) << "'.";
+      CHECK_OK(file::SetBinaryProto(absl::GetFlag(FLAGS_output_proto), proto,
+                                    file::Defaults()));
+    }
+
     // TODO(user): improve code to convert many files.
     if (!absl::GetFlag(FLAGS_output_dimacs).empty()) {
       LOG(INFO) << "Converting first input file to dimacs format.";
@@ -344,7 +372,18 @@ int main(int argc, char** argv) {
         SolveMinCostFlow(proto, &loading_time, &solving_time);
         break;
       case FlowModelProto::MAX_FLOW:
-        SolveMaxFlow(proto, &loading_time, &solving_time);
+        if (absl::GetFlag(FLAGS_use_flow_graph)) {
+          SolveMaxFlow<util::FlowGraph<>>(
+              proto, &loading_time, &solving_time,
+              [](util::FlowGraph<>* graph) {
+                graph->SetDetectReverse(
+                    absl::GetFlag(FLAGS_detect_reverse_arcs));
+                graph->SetSortByHead(absl::GetFlag(FLAGS_sort_heads));
+              });
+        } else {
+          SolveMaxFlow<util::ReverseArcStaticGraph<>>(proto, &loading_time,
+                                                      &solving_time);
+        }
         break;
       default:
         LOG(ERROR) << "Problem type not supported: " << proto.problem_type();