diff --git a/ortools/graph/dag_constrained_shortest_path.cc b/ortools/graph/dag_constrained_shortest_path.cc
index 6b2c57e1e1..8fe8450323 100644
--- a/ortools/graph/dag_constrained_shortest_path.cc
+++ b/ortools/graph/dag_constrained_shortest_path.cc
@@ -78,8 +78,8 @@ PathWithLength ConstrainedShortestPathsOnDag(
   std::vector<NodeIndex> destinations = {destination};
   ConstrainedShortestPathsOnDagWrapper<util::StaticGraph<>>
       constrained_shortest_path_on_dag(&graph, &arc_lengths, &arc_resources,
-                                       &(*topological_order), &sources,
-                                       &destinations, &max_resources);
+                                       *topological_order, sources,
+                                       destinations, &max_resources);
 
   PathWithLength path_with_length =
       constrained_shortest_path_on_dag.RunConstrainedShortestPathOnDag();
diff --git a/ortools/graph/dag_constrained_shortest_path.h b/ortools/graph/dag_constrained_shortest_path.h
index b864d8d94f..02993acece 100644
--- a/ortools/graph/dag_constrained_shortest_path.h
+++ b/ortools/graph/dag_constrained_shortest_path.h
@@ -14,6 +14,8 @@
 #ifndef OR_TOOLS_GRAPH_DAG_CONSTRAINED_SHORTEST_PATH_H_
 #define OR_TOOLS_GRAPH_DAG_CONSTRAINED_SHORTEST_PATH_H_
 
+#include <stdbool.h>
+
 #include <cmath>
 #include <concepts>
 #include <limits>
@@ -25,13 +27,25 @@
 #include "absl/strings/str_format.h"
 #include "absl/types/span.h"
 #include "ortools/graph/dag_shortest_path.h"
+#include "ortools/graph/graph.h"
 
 namespace operations_research {
 
 // This library provides APIs to compute the constrained shortest path (CSP) on
 // a given directed acyclic graph (DAG) with resources on each arc. A CSP is a
 // shortest path on a DAG which does not exceed a set of maximum resources
-// consumption. The algorithm is exponential and has no guarantee to finish.
+// consumption. The algorithm is exponential and has no guarantee to finish. It
+// is based on bi-drectionnal search. First is a forward pass from the source to
+// nodes “somewhere in the middle” to generate forward labels, just as the
+// onedirectional labeling algorithm we discussed; then a symmetric backward
+// pass from the destination generates backward labels; and finally at each node
+// with both forward and backward labels, it joins any pair of labels to form a
+// feasible complete path. Intuitively, the number of labels grows exponentially
+// with the number of arcs in the path. The overall number of labels are then
+// expected to be smaller with shorter paths. For DAG with a topological
+// ordering, we can pick any node (usually right in the middle) as a *midpoint*
+// to stop each pass at. Then labels can be joined at only one half of the nodes
+// by considering all edges between each half.
 //
 // In the DAG, multiple arcs between the same pair of nodes is allowed. However,
 // self-loop arcs are not allowed.
@@ -114,70 +128,156 @@ class ConstrainedShortestPathsOnDagWrapper {
   //
   // Validity of arcs and topological order are DCHECKed.
   //
-  // If the number of labels in memory exceeds `max_num_created_labels` at any
-  // point in the algorithm, it returns the best path found so far, most
-  // particularly the empty path if none were found.
+  // If the number of labels in memory exceeds `max_num_created_labels / 2` at
+  // any point in each pass of the algorithm, new labels are not generated
+  // anymore and it returns the best path found so far, most particularly the
+  // empty path if none were found.
   //
-  // SUBTLE: You can modify the graph, the arc lengths and resources, the
-  // topological order, sources, destinations or the maximum resource between
-  // calls to the `RunConstrainedShortestPathOnDag()` function. That's fine.
-  // Doing so will obviously invalidate the result API of the last constrained
-  // shortest path run, which could return an upper bound, junk, or crash.
+  // IMPORTANT: You cannot modify anything except `arc_lengths` between calls to
+  // the `RunConstrainedShortestPathOnDag()` function.
+  // TODO(b/330285675): Change API to copy everything internally except the
+  // lengths.
   ConstrainedShortestPathsOnDagWrapper(
       const GraphType* graph, const std::vector<double>* arc_lengths,
       const std::vector<std::vector<double>>* arc_resources,
-      const std::vector<NodeIndex>* topological_order,
-      const std::vector<NodeIndex>* sources,
-      const std::vector<NodeIndex>* destinations,
+      absl::Span<const NodeIndex> topological_order,
+      absl::Span<const NodeIndex> sources,
+      absl::Span<const NodeIndex> destinations,
       const std::vector<double>* max_resources,
       int max_num_created_labels = 1e9);
 
   // Returns {+inf, {}, {}} if there is no constrained path of finite length
-  // from one node in `sources` to one node in `destinations`.
+  // wihtin resources constraints from one node in `sources` to one node in
+  // `destinations`.
   PathWithLength RunConstrainedShortestPathOnDag();
 
+  // For benchmarking and informational purposes, returns the number of labels
+  // generated in the call of `RunConstrainedShortestPathOnDag()`.
+  int label_count() const {
+    return lengths_from_sources_[FORWARD].size() +
+           lengths_from_sources_[BACKWARD].size();
+  }
+
  private:
-  // Returns the list of all the arcs of the shortest path from one node in
-  // `sources` ending by the arc from a given `label_index` if and only if
-  // `label_index` is between 0 and `labels_from_sources_.size() - 1`.
-  std::vector<ArcIndex> BestArcPathEndingWith(int label_index) const;
+  enum Direction {
+    FORWARD = 0,
+    BACKWARD = 1,
+  };
+
+  inline static Direction Reverse(Direction d) {
+    return d == FORWARD ? BACKWARD : FORWARD;
+  }
+
+  // A LabelPair includes the `length` of a path that can be constructed by
+  // merging the paths from two *linkable* labels corresponding to
+  // `label_index`.
+  struct LabelPair {
+    double length = 0.0;
+    int label_index[2];
+  };
+
+  void RunHalfConstrainedShortestPathOnDag(
+      const GraphType& reverse_graph, absl::Span<const double> arc_lengths,
+      absl::Span<const std::vector<double>> arc_resources,
+      absl::Span<const std::vector<double>> min_arc_resources,
+      absl::Span<const double> max_resources, int max_num_created_labels,
+      std::vector<double>& lengths_from_sources,
+      std::vector<std::vector<double>>& resources_from_sources,
+      std::vector<ArcIndex>& incoming_arc_indices_from_sources,
+      std::vector<int>& first_label, std::vector<int>& num_labels);
+
+  // Returns the arc index linking two nodes from each pass forming the best
+  // path. Returns -1 if no better path than the one found from
+  // `best_label_pair` is found.
+  ArcIndex MergeHalfRuns(
+      const GraphType& graph, absl::Span<const double> arc_lengths,
+      absl::Span<const std::vector<double>> arc_resources,
+      absl::Span<const double> max_resources,
+      const std::vector<NodeIndex> sub_node_indices[2],
+      const std::vector<double> lengths_from_sources[2],
+      const std::vector<std::vector<double>> resources_from_sources[2],
+      const std::vector<int> first_label[2],
+      const std::vector<int> num_labels[2], LabelPair& best_label_pair);
+
+  // Returns the path as list of arc indices that starts from a node in
+  // `sources` (if `direction` iS FORWARD) or `destinations` (if `direction` is
+  // BACKWARD) and ends in node represented by `best_label_index`.
+  std::vector<ArcIndex> ArcPathTo(
+      int best_label_index, const GraphType& reverse_graph,
+      absl::Span<const double> arc_lengths,
+      absl::Span<const double> lengths_from_sources,
+      absl::Span<const ArcIndex> incoming_arc_indices_from_sources,
+      absl::Span<const int> first_label,
+      absl::Span<const int> num_labels) const;
+
   // Returns the list of all the nodes implied by a given `arc_path`.
-  std::vector<NodeIndex> NodePathImpliedBy(
-      const std::vector<ArcIndex>& arc_path) const;
+  std::vector<NodeIndex> NodePathImpliedBy(absl::Span<const ArcIndex> arc_path,
+                                           const GraphType& graph) const;
+
+  static constexpr double kTolerance = 1e-6;
 
   const GraphType* const graph_;
   const std::vector<double>* const arc_lengths_;
   const std::vector<std::vector<double>>* const arc_resources_;
-  const std::vector<NodeIndex>* const topological_order_;
-  const std::vector<NodeIndex>* const sources_;
-  const std::vector<NodeIndex>* const destinations_;
   const std::vector<double>* const max_resources_;
-  int max_num_created_labels_;
+  absl::Span<const NodeIndex> sources_;
+  absl::Span<const NodeIndex> destinations_;
+  const int num_resources_;
 
-  std::vector<bool> node_is_source_;
-  std::vector<bool> node_is_destination_;
-  // Using the fact that the graph is a DAG, we can disregard any node that
-  // comes after the last destination (based on the topological order).
-  std::vector<bool> node_is_after_last_destination_;
-
-  // Data for reverse graph.
-  GraphType reverse_graph_;
-  std::vector<ArcIndex> reverse_inverse_arc_permutation_;
+  // Data about *reachable* sub-graphs split in two for bidirectional search.
+  // Reachable nodes are nodes that can be reached given the resources
+  // constraints, i.e., for each resource, the sum of the minimum resource to
+  // get to a node from a node in `sources` and to get from a node to a node in
+  // `destinations` should be less than the maximum resource. Reachable arcs are
+  // arcs linking reachable nodes.
+  //
+  // `sub_reverse_graph_[dir]` is the reachable sub-graph split in *half* with
+  // an additional linked to sources (resp. destinations) for the forward (resp.
+  // backward) direction. For the forward (resp. backward) direction, nodes are
+  // indexed using the original (resp. reverse) topological order.
+  GraphType sub_reverse_graph_[2];
+  std::vector<std::vector<double>> sub_arc_resources_[2];
+  // `sub_full_arc_indices_[dir]` has size `sub_reverse_graph_[dir].num_arcs()`
+  // such that `sub_full_arc_indices_[dir][sub_arc] = arc` where `sub_arc` is
+  // the arc in the reachable sub-graph for direction `dir` (i.e.
+  // `sub_reverse_graph[dir]`) and `arc` is the arc in the original graph (i.e.
+  // `graph`).
+  std::vector<NodeIndex> sub_full_arc_indices_[2];
+  // `sub_node_indices_[dir]` has size `graph->num_nodes()` such that
+  // `sub_node_indices[dir][node] = sub_node` where `node` is the node in the
+  // original graph (i.e. `graph`) and `sub_node` is the node in the reachable
+  // sub-graph for direction `dir` (i.e. `sub_reverse_graph[dir]`) and -1 if
+  // `node` is not present in reachable sub-graph.
+  std::vector<NodeIndex> sub_node_indices_[2];
+  // `sub_is_source_[dir][sub_dir]` has size
+  // `sub_reverse_graph_[dir].num_nodes()` such that
+  // `sub_is_source_[dir][sub_dir][sub_node]` is true if `sub_node` is a node in
+  // the reachable sub-graph for direction `dir` (i.e. `sub_reverse_graph[dir]`)
+  // which is a source (resp. destination) is `sub_dir` is FORWARD (resp.
+  // BACKWARD).
+  std::vector<bool> sub_is_source_[2][2];
+  // `sub_min_arc_resources_[dir]` has size `max_resources->size()` and
+  // `sub_min_arc_resources_[dir][r]`, `sub_reverse_graph_[dir].num_nodes()`
+  // such that `sub_min_arc_resources_[dir][r][sub_node]` is the minimum of
+  // resource r needed to get to a destination (resp. come from a source) if
+  // `dir` is FORWARD (resp. BACKWARD).
+  std::vector<std::vector<double>> sub_min_arc_resources_[2];
+  // Maximum number of labels created for each sub-graph.
+  int max_num_created_labels_[2];
 
   // Data about the last call of the RunConstrainedShortestPathOnDag()
-  // function. A Label includes the cumulative length, resources and the
-  // previous arc used in the path to get to this node.
-  struct Label {
-    double length;
-    // TODO(b/315786885): Optimize resources in Label struct.
-    std::vector<double> resources;
-    ArcIndex incoming_arc;
-  };
-  // A label is present in `labels_from_sources_` if and only if it is feasible
-  // with respect to all resources.
-  std::vector<Label> labels_from_sources_;
-  std::vector<int> node_first_label_;
-  std::vector<int> node_num_labels_;
+  // function. A path is only added to the following vectors if and only if
+  // it is feasible with respect to all resources.
+  // A Label includes the cumulative length, resources and the previous arc used
+  // in the path to get to this node.
+  // Instead of having a single vector of `Label` objects (cl/590819865), we
+  // split them into 3 vectors of more fundamental types as this improves
+  // push_back operations and memory release.
+  std::vector<double> lengths_from_sources_[2];
+  std::vector<std::vector<double>> resources_from_sources_[2];
+  std::vector<ArcIndex> incoming_arc_indices_from_sources_[2];
+  std::vector<int> node_first_label_[2];
+  std::vector<int> node_num_labels_[2];
 };
 
 std::vector<int> GetInversePermutation(absl::Span<const int> permutation);
@@ -194,30 +294,31 @@ ConstrainedShortestPathsOnDagWrapper<GraphType>::
     ConstrainedShortestPathsOnDagWrapper(
         const GraphType* graph, const std::vector<double>* arc_lengths,
         const std::vector<std::vector<double>>* arc_resources,
-        const std::vector<NodeIndex>* topological_order,
-        const std::vector<NodeIndex>* sources,
-        const std::vector<NodeIndex>* destinations,
+        absl::Span<const NodeIndex> topological_order,
+        absl::Span<const NodeIndex> sources,
+        absl::Span<const NodeIndex> destinations,
         const std::vector<double>* max_resources, int max_num_created_labels)
     : graph_(graph),
       arc_lengths_(arc_lengths),
       arc_resources_(arc_resources),
-      topological_order_(topological_order),
+      max_resources_(max_resources),
       sources_(sources),
       destinations_(destinations),
-      max_resources_(max_resources),
-      max_num_created_labels_(max_num_created_labels) {
+      num_resources_(max_resources->size()) {
   CHECK(graph_ != nullptr);
   CHECK(arc_lengths_ != nullptr);
   CHECK(arc_resources_ != nullptr);
-  CHECK(topological_order_ != nullptr);
-  CHECK(sources_ != nullptr);
-  CHECK(destinations_ != nullptr);
+  CHECK(!sources_.empty());
+  CHECK(!destinations_.empty());
   CHECK(max_resources_ != nullptr);
+  CHECK(!max_resources_->empty())
+      << "max_resources cannot be empty. Use "
+         "ortools/graph/dag_shortest_path.h instead";
   if (DEBUG_MODE) {
-    CHECK_EQ(arc_lengths_->size(), graph_->num_arcs());
-    CHECK_EQ(arc_resources_->size(), max_resources_->size());
-    for (const std::vector<double>& arcs_resource : *arc_resources_) {
-      CHECK_EQ(arcs_resource.size(), graph_->num_arcs());
+    CHECK_EQ(arc_lengths->size(), graph->num_arcs());
+    CHECK_EQ(arc_resources->size(), max_resources->size());
+    for (absl::Span<const double> arcs_resource : *arc_resources) {
+      CHECK_EQ(arcs_resource.size(), graph->num_arcs());
       for (const double arc_resource : arcs_resource) {
         CHECK(arc_resource >= 0 &&
               arc_resource != std::numeric_limits<double>::infinity() &&
@@ -225,53 +326,192 @@ ConstrainedShortestPathsOnDagWrapper<GraphType>::
             << absl::StrFormat("resource cannot be negative nor +inf nor NaN");
       }
     }
-    for (const double arc_length : *arc_lengths_) {
+    for (const double arc_length : *arc_lengths) {
       CHECK(arc_length != -std::numeric_limits<double>::infinity() &&
             !std::isnan(arc_length))
           << absl::StrFormat("length cannot be -inf nor NaN");
     }
-    CHECK_OK(TopologicalOrderIsValid(*graph_, *topological_order_))
+    CHECK_OK(TopologicalOrderIsValid(*graph, topological_order))
         << "Invalid topological order";
-    for (const double max_resource : *max_resources_) {
+    for (const double max_resource : *max_resources) {
       CHECK(max_resource >= 0 &&
             max_resource != std::numeric_limits<double>::infinity() &&
             !std::isnan(max_resource))
           << absl::StrFormat(
                  "max_resource cannot be negative not +inf nor NaN");
     }
-  }
-
-  node_is_source_.resize(graph_->num_nodes());
-  for (const NodeIndex source : *sources_) {
-    node_is_source_[source] = true;
-  }
-  node_is_destination_.resize(graph_->num_nodes());
-  for (const NodeIndex destination : *destinations_) {
-    node_is_destination_[destination] = true;
-  }
-  node_is_after_last_destination_.resize(graph_->num_nodes());
-  bool last_destination_found = false;
-  for (int i = graph_->num_nodes() - 1; i >= 0; --i) {
-    const NodeIndex node = (*topological_order_)[i];
-    if (node_is_destination_[node]) {
-      last_destination_found = true;
+    std::vector<bool> is_source(graph->num_nodes(), false);
+    for (const NodeIndex source : sources) {
+      is_source[source] = true;
+    }
+    for (const NodeIndex destination : destinations) {
+      CHECK(!is_source[destination])
+          << "A node cannot be both a source and destination";
     }
-    node_is_after_last_destination_[node] = !last_destination_found;
   }
 
-  // Reverse graph.
-  reverse_graph_ = GraphType(graph_->num_nodes(), graph_->num_arcs());
-  for (ArcIndex arc_index = 0; arc_index < graph_->num_arcs(); ++arc_index) {
-    reverse_graph_.AddArc(graph_->Head(arc_index), graph_->Tail(arc_index));
+  // Full graphs.
+  const GraphType* full_graph[2];
+  const std::vector<std::vector<double>>* full_arc_resources[2];
+  absl::Span<const NodeIndex> full_topological_order[2];
+  absl::Span<const NodeIndex> full_sources[2];
+  // Forward.
+  const int num_nodes = graph->num_nodes();
+  const int num_arcs = graph->num_arcs();
+  full_graph[FORWARD] = graph;
+  full_arc_resources[FORWARD] = arc_resources;
+  full_topological_order[FORWARD] = topological_order;
+  full_sources[FORWARD] = sources;
+  // Backward.
+  GraphType full_backward_graph(num_nodes, num_arcs);
+  for (ArcIndex arc_index = 0; arc_index < num_arcs; ++arc_index) {
+    full_backward_graph.AddArc(graph->Head(arc_index), graph->Tail(arc_index));
   }
-  std::vector<ArcIndex> permutation;
-  reverse_graph_.Build(&permutation);
-  reverse_inverse_arc_permutation_ = GetInversePermutation(permutation);
+  std::vector<ArcIndex> full_permutation;
+  full_backward_graph.Build(&full_permutation);
+  const std::vector<ArcIndex> full_inverse_arc_indices =
+      GetInversePermutation(full_permutation);
+  std::vector<std::vector<double>> backward_arc_resources(num_resources_);
+  for (int r = 0; r < num_resources_; ++r) {
+    backward_arc_resources[r] = (*arc_resources)[r];
+    util::Permute(full_permutation, &backward_arc_resources[r]);
+  }
+  std::vector<NodeIndex> full_backward_topological_order;
+  full_backward_topological_order.reserve(num_nodes);
+  for (int i = num_nodes - 1; i >= 0; --i) {
+    full_backward_topological_order.push_back(topological_order[i]);
+  }
+  full_graph[BACKWARD] = &full_backward_graph;
+  full_arc_resources[BACKWARD] = &backward_arc_resources;
+  full_topological_order[BACKWARD] = full_backward_topological_order;
+  full_sources[BACKWARD] = destinations;
 
-  // Memory allocation is done here and only once in order to avoid reallocation
-  // at each call of `RunConstrainedShortestPathOnDag()` for better performance.
-  node_first_label_.resize(graph_->num_nodes());
-  node_num_labels_.resize(graph_->num_nodes());
+  // Get the minimum resources sources -> node and node -> destination for each
+  // node.
+  std::vector<std::vector<double>> full_min_arc_resources[2];
+  for (const Direction dir : {FORWARD, BACKWARD}) {
+    full_min_arc_resources[dir].reserve(num_resources_);
+    std::vector<double> full_arc_resource = full_arc_resources[dir]->front();
+    ShortestPathsOnDagWrapper<GraphType> shortest_paths_on_dag(
+        full_graph[dir], &full_arc_resource, full_topological_order[dir]);
+    for (int r = 0; r < num_resources_; ++r) {
+      full_arc_resource = (*(full_arc_resources[dir]))[r];
+      shortest_paths_on_dag.RunShortestPathOnDag(full_sources[dir]);
+      full_min_arc_resources[dir].push_back(shortest_paths_on_dag.LengthTo());
+    }
+  }
+
+  // Get reachable subgraph.
+  std::vector<NodeIndex> sub_topological_order;
+  sub_topological_order.reserve(num_nodes);
+  for (const NodeIndex node_index : topological_order) {
+    bool is_reachable = true;
+    for (int r = 0; r < num_resources_; ++r) {
+      if (full_min_arc_resources[FORWARD][r][node_index] +
+              full_min_arc_resources[BACKWARD][r][node_index] >
+          (*max_resources)[r]) {
+        is_reachable = false;
+        break;
+      }
+    }
+    if (is_reachable) {
+      sub_topological_order.push_back(node_index);
+    }
+  }
+  const int reachable_node_count = sub_topological_order.size();
+
+  // Split sub-graphs and related information.
+  const int mid_index = reachable_node_count / 2;
+  // We split the number of labels evenly between each search (+1 for the
+  // additional source node)
+  // TODO(329206892): Could be optimized depending on the split and number of
+  // arcs in the sub-graph.
+  max_num_created_labels_[FORWARD] = max_num_created_labels / 2;
+  max_num_created_labels_[BACKWARD] =
+      max_num_created_labels - max_num_created_labels_[FORWARD];
+  for (const Direction dir : {FORWARD, BACKWARD}) {
+    absl::Span<const NodeIndex> const sub_nodes =
+        dir == FORWARD
+            ? absl::MakeSpan(sub_topological_order).subspan(0, mid_index)
+            : absl::MakeSpan(sub_topological_order)
+                  .subspan(mid_index, reachable_node_count - mid_index);
+    sub_node_indices_[dir].assign(num_nodes, -1);
+    sub_min_arc_resources_[dir].resize(num_resources_);
+    for (int r = 0; r < num_resources_; ++r) {
+      sub_min_arc_resources_[dir][r].resize(sub_nodes.size());
+    }
+    for (NodeIndex i = 0; i < sub_nodes.size(); ++i) {
+      const NodeIndex sub_node_index =
+          dir == FORWARD ? i : sub_nodes.size() - 1 - i;
+      sub_node_indices_[dir][sub_nodes[i]] = sub_node_index;
+      for (int r = 0; r < num_resources_; ++r) {
+        sub_min_arc_resources_[dir][r][sub_node_index] =
+            full_min_arc_resources[Reverse(dir)][r][sub_nodes[i]];
+      }
+    }
+    // IMPORTANT: The sub-graph has an additional node linked to sources (resp.
+    // destinations) for the forward (resp. backward) direction. This additional
+    // node is indexed with the last index. All added arcs are given to have an
+    // arc index in the original graph of -1.
+    const int sub_arcs_count = num_arcs + full_sources[dir].size();
+    sub_reverse_graph_[dir] = GraphType(sub_nodes.size() + 1, sub_arcs_count);
+    sub_arc_resources_[dir].resize(num_resources_);
+    for (int r = 0; r < num_resources_; ++r) {
+      sub_arc_resources_[dir][r].reserve(sub_arcs_count);
+    }
+    sub_full_arc_indices_[dir].reserve(sub_arcs_count);
+    const GraphType& reverse_full_graph = *(full_graph[Reverse(dir)]);
+    for (ArcIndex arc_index = 0; arc_index < num_arcs; ++arc_index) {
+      const NodeIndex from =
+          sub_node_indices_[dir][reverse_full_graph.Tail(arc_index)];
+      const NodeIndex to =
+          sub_node_indices_[dir][reverse_full_graph.Head(arc_index)];
+      if (from == -1 || to == -1) {
+        continue;
+      }
+      sub_reverse_graph_[dir].AddArc(from, to);
+      ArcIndex sub_full_arc_index;
+      if (dir == FORWARD && !full_inverse_arc_indices.empty()) {
+        sub_full_arc_index = full_inverse_arc_indices[arc_index];
+      } else {
+        sub_full_arc_index = arc_index;
+      }
+      for (int r = 0; r < num_resources_; ++r) {
+        sub_arc_resources_[dir][r].push_back(
+            (*arc_resources_)[r][sub_full_arc_index]);
+      }
+      sub_full_arc_indices_[dir].push_back(sub_full_arc_index);
+    }
+    for (const NodeIndex source : full_sources[dir]) {
+      const NodeIndex sub_source = sub_node_indices_[dir][source];
+      if (sub_source == -1) {
+        continue;
+      }
+      sub_reverse_graph_[dir].AddArc(sub_source, sub_nodes.size());
+      for (int r = 0; r < num_resources_; ++r) {
+        sub_arc_resources_[dir][r].push_back(0.0);
+      }
+      sub_full_arc_indices_[dir].push_back(-1);
+    }
+    std::vector<ArcIndex> sub_permutation;
+    sub_reverse_graph_[dir].Build(&sub_permutation);
+    for (int r = 0; r < num_resources_; ++r) {
+      util::Permute(sub_permutation, &sub_arc_resources_[dir][r]);
+    }
+    util::Permute(sub_permutation, &sub_full_arc_indices_[dir]);
+    // We add labels corresponding to the sources in the sub-graph as these will
+    // be artificially created labels due to the
+    max_num_created_labels_[dir] += full_sources[dir].size();
+  }
+
+  // Memory allocation is done here and only once in order to avoid
+  // reallocation at each call of `RunConstrainedShortestPathOnDag()` for
+  // better performance.
+  for (const Direction dir : {FORWARD, BACKWARD}) {
+    resources_from_sources_[dir].resize(num_resources_);
+    node_first_label_[dir].resize(sub_reverse_graph_[dir].size());
+    node_num_labels_[dir].resize(sub_reverse_graph_[dir].size());
+  }
 }
 
 template <class GraphType>
@@ -280,56 +520,308 @@ template <class GraphType>
 #endif
 PathWithLength ConstrainedShortestPathsOnDagWrapper<
     GraphType>::RunConstrainedShortestPathOnDag() {
-  // Caching the vector addresses allow to not fetch it on each access.
-  const absl::Span<const double> arc_lengths = *arc_lengths_;
-  const absl::Span<const double> max_resources = *max_resources_;
-
-  // Clear all labels from previous run.
-  labels_from_sources_.clear();
-
-  double best_length = std::numeric_limits<double>::infinity();
-  int best_label_index = -1;
-  for (const NodeIndex to : *topological_order_) {
-    node_first_label_[to] = labels_from_sources_.size();
-    int& num_labels = node_num_labels_[to];
-    num_labels = 0;
-    if (node_is_after_last_destination_[to]) {
-      break;
+  // Assign lengths on sub-relevant graphs.
+  std::vector<double> sub_arc_lengths[2];
+  for (const Direction dir : {FORWARD, BACKWARD}) {
+    sub_arc_lengths[dir].reserve(sub_reverse_graph_[dir].num_arcs());
+    for (ArcIndex sub_arc_index = 0;
+         sub_arc_index < sub_reverse_graph_[dir].num_arcs(); ++sub_arc_index) {
+      const ArcIndex arc_index = sub_full_arc_indices_[dir][sub_arc_index];
+      if (arc_index == -1) {
+        sub_arc_lengths[dir].push_back(0.0);
+        continue;
+      }
+      sub_arc_lengths[dir].push_back((*arc_lengths_)[arc_index]);
     }
-    if (node_is_source_[to]) {
-      labels_from_sources_.push_back(
-          {.length = 0,
-           .resources = std::vector<double>(max_resources.size()),
-           .incoming_arc = -1});
-      ++num_labels;
-      if (labels_from_sources_.size() >= max_num_created_labels_) {
-        const std::vector<ArcIndex> arc_path =
-            BestArcPathEndingWith(best_label_index);
-        return {.length = best_length,
-                .arc_path = arc_path,
-                .node_path = NodePathImpliedBy(arc_path)};
+  }
+
+  // TODO(b/329036127): Each search could be done in parallel as no data is
+  // shared except `best_label_pair_`.
+  LabelPair best_label_pair = {
+      .length = std::numeric_limits<double>::infinity(),
+      .label_index = {-1, -1}};
+  for (const Direction dir : {FORWARD, BACKWARD}) {
+    RunHalfConstrainedShortestPathOnDag(
+        /*reverse_graph=*/sub_reverse_graph_[dir],
+        /*arc_lengths=*/sub_arc_lengths[dir],
+        /*arc_resources=*/sub_arc_resources_[dir],
+        /*min_arc_resources=*/sub_min_arc_resources_[dir],
+        /*max_resources=*/*max_resources_,
+        /*max_num_created_labels=*/max_num_created_labels_[dir],
+        /*lengths_from_sources=*/lengths_from_sources_[dir],
+        /*resources_from_sources=*/resources_from_sources_[dir],
+        /*incoming_arc_indices_from_sources=*/
+        incoming_arc_indices_from_sources_[dir],
+        /*first_label=*/node_first_label_[dir],
+        /*num_labels=*/node_num_labels_[dir]);
+  }
+
+  // Check destinations within relevant half sub-graphs.
+  for (const Direction dir : {FORWARD, BACKWARD}) {
+    absl::Span<const NodeIndex> destinations =
+        dir == FORWARD ? destinations_ : sources_;
+    for (const NodeIndex dst : destinations) {
+      const NodeIndex sub_dst = sub_node_indices_[dir][dst];
+      if (sub_dst == -1) {
+        continue;
+      }
+      const int num_labels_dst = node_num_labels_[dir][sub_dst];
+      if (num_labels_dst == 0) {
+        continue;
+      }
+      const int first_label_dst = node_first_label_[dir][sub_dst];
+      for (int label_index = first_label_dst;
+           label_index < first_label_dst + num_labels_dst; ++label_index) {
+        const double length_dst = lengths_from_sources_[dir][label_index];
+        if (length_dst < best_label_pair.length) {
+          best_label_pair.length = length_dst;
+          best_label_pair.label_index[dir] = label_index;
+        }
       }
     }
-    for (const ArcIndex reverse_arc_index : reverse_graph_.OutgoingArcs(to)) {
-      const NodeIndex from = reverse_graph_.Head(reverse_arc_index);
-      const ArcIndex arc_index =
-          reverse_inverse_arc_permutation_.empty()
-              ? reverse_arc_index
-              : reverse_inverse_arc_permutation_[reverse_arc_index];
-      const double arc_length = arc_lengths[arc_index];
+  }
+
+  const ArcIndex merging_arc_index = MergeHalfRuns(
+      /*graph=*/*graph_, /*arc_lengths=*/*arc_lengths_,
+      /*arc_resources=*/*arc_resources_,
+      /*max_resources=*/*max_resources_,
+      /*sub_node_indices=*/sub_node_indices_,
+      /*lengths_from_sources=*/lengths_from_sources_,
+      /*resources_from_sources=*/resources_from_sources_,
+      /*first_label=*/node_first_label_,
+      /*num_labels=*/node_num_labels_, /*best_label_pair=*/best_label_pair);
+
+  std::vector<ArcIndex> arc_path;
+  for (const Direction dir : {FORWARD, BACKWARD}) {
+    for (const ArcIndex sub_arc_index : ArcPathTo(
+             /*best_label_index=*/best_label_pair.label_index[dir],
+             /*reverse_graph=*/sub_reverse_graph_[dir],
+             /*arc_lengths=*/sub_arc_lengths[dir],
+             /*lengths_from_sources=*/lengths_from_sources_[dir],
+             /*incoming_arc_indices_from_sources=*/
+             incoming_arc_indices_from_sources_[dir],
+             /*first_label=*/node_first_label_[dir],
+             /*num_labels=*/node_num_labels_[dir])) {
+      const ArcIndex arc_index = sub_full_arc_indices_[dir][sub_arc_index];
+      if (arc_index == -1) {
+        break;
+      }
+      arc_path.push_back(arc_index);
+    }
+    if (dir == FORWARD && merging_arc_index != -1) {
+      absl::c_reverse(arc_path);
+      arc_path.push_back(merging_arc_index);
+    }
+  }
+
+  // Clear all labels from the next run.
+  for (const Direction dir : {FORWARD, BACKWARD}) {
+    lengths_from_sources_[dir].clear();
+    for (int r = 0; r < num_resources_; ++r) {
+      resources_from_sources_[dir][r].clear();
+    }
+    incoming_arc_indices_from_sources_[dir].clear();
+  }
+  return {.length = best_label_pair.length,
+          .arc_path = arc_path,
+          .node_path = NodePathImpliedBy(arc_path, *graph_)};
+}
+
+template <class GraphType>
+#if __cplusplus >= 202002L
+  requires DagGraphType<GraphType>
+#endif
+void ConstrainedShortestPathsOnDagWrapper<GraphType>::
+    RunHalfConstrainedShortestPathOnDag(
+        const GraphType& reverse_graph, absl::Span<const double> arc_lengths,
+        absl::Span<const std::vector<double>> arc_resources,
+        absl::Span<const std::vector<double>> min_arc_resources,
+        absl::Span<const double> max_resources,
+        const int max_num_created_labels,
+        std::vector<double>& lengths_from_sources,
+        std::vector<std::vector<double>>& resources_from_sources,
+        std::vector<ArcIndex>& incoming_arc_indices_from_sources,
+        std::vector<int>& first_label, std::vector<int>& num_labels) {
+  // Initialize source node.
+  const NodeIndex source_node = reverse_graph.num_nodes() - 1;
+  first_label[source_node] = 0;
+  num_labels[source_node] = 1;
+  lengths_from_sources.push_back(0);
+  for (int r = 0; r < num_resources_; ++r) {
+    resources_from_sources[r].push_back(0);
+  }
+  incoming_arc_indices_from_sources.push_back(-1);
+
+  std::vector<double> lengths_to;
+  std::vector<std::vector<double>> resources_to(num_resources_);
+  std::vector<ArcIndex> incoming_arc_indices_to;
+  std::vector<int> label_indices_to;
+  std::vector<double> resources(num_resources_);
+  for (NodeIndex to = 0; to < source_node; ++to) {
+    lengths_to.clear();
+    for (int r = 0; r < num_resources_; ++r) {
+      resources_to[r].clear();
+    }
+    incoming_arc_indices_to.clear();
+    for (const ArcIndex reverse_arc_index : reverse_graph.OutgoingArcs(to)) {
+      const NodeIndex from = reverse_graph.Head(reverse_arc_index);
+      const double arc_length = arc_lengths[reverse_arc_index];
       DCHECK(arc_length != -std::numeric_limits<double>::infinity());
       if (arc_length == std::numeric_limits<double>::infinity()) {
         continue;
       }
-      for (int i = node_first_label_[from];
-           i < node_first_label_[from] + node_num_labels_[from]; ++i) {
-        const Label& label_from = labels_from_sources_[i];
+      for (int label_index = first_label[from];
+           label_index < first_label[from] + num_labels[from]; ++label_index) {
         bool path_is_feasible = true;
-        for (int r = 0; r < max_resources.size(); ++r) {
-          DCHECK_GE((*arc_resources_)[r][arc_index], 0.0);
-          // TODO(b/314756645): Include lower bound on the resources to
-          // increase pruning.
-          if (label_from.resources[r] + (*arc_resources_)[r][arc_index] >
+        for (int r = 0; r < num_resources_; ++r) {
+          DCHECK_GE(arc_resources[r][reverse_arc_index], 0.0);
+          resources[r] = resources_from_sources[r][label_index] +
+                         arc_resources[r][reverse_arc_index];
+          if (resources[r] + min_arc_resources[r][to] > max_resources[r]) {
+            path_is_feasible = false;
+            break;
+          }
+        }
+        if (!path_is_feasible) {
+          continue;
+        }
+        lengths_to.push_back(lengths_from_sources[label_index] + arc_length);
+        for (int r = 0; r < num_resources_; ++r) {
+          resources_to[r].push_back(resources[r]);
+        }
+        incoming_arc_indices_to.push_back(reverse_arc_index);
+      }
+    }
+    // Sort labels lexicographically with lengths then resources.
+    label_indices_to.clear();
+    label_indices_to.reserve(lengths_to.size());
+    for (int i = 0; i < lengths_to.size(); ++i) {
+      label_indices_to.push_back(i);
+    }
+    absl::c_sort(label_indices_to, [&](const int i, const int j) {
+      if (lengths_to[i] < lengths_to[j]) return true;
+      if (lengths_to[i] > lengths_to[j]) return false;
+      for (int r = 0; r < num_resources_; ++r) {
+        if (resources_to[r][i] < resources_to[r][j]) return true;
+        if (resources_to[r][i] > resources_to[r][j]) return false;
+      }
+      return i < j;
+    });
+
+    first_label[to] = lengths_from_sources.size();
+    int& num_labels_to = num_labels[to];
+    // Reset the number of labels to zero otherwise it holds the previous run
+    // result.
+    num_labels_to = 0;
+    for (int i = 0; i < label_indices_to.size(); ++i) {
+      // Check if label "i" on node `to` is dominated by any other label.
+      const int label_i_index = label_indices_to[i];
+      bool label_i_is_dominated = false;
+      for (int j = 0; j < i - 1; ++j) {
+        const int label_j_index = label_indices_to[j];
+        if (lengths_to[label_i_index] <= lengths_to[label_j_index]) continue;
+        bool label_j_dominates_label_i = true;
+        for (int r = 0; r < num_resources_; ++r) {
+          if (resources_to[r][label_i_index] <=
+              resources_to[r][label_j_index]) {
+            label_j_dominates_label_i = false;
+            break;
+          }
+        }
+        if (label_j_dominates_label_i) {
+          label_i_is_dominated = true;
+          break;
+        }
+      }
+      if (label_i_is_dominated) continue;
+      lengths_from_sources.push_back(lengths_to[label_i_index]);
+      for (int r = 0; r < num_resources_; ++r) {
+        resources_from_sources[r].push_back(resources_to[r][label_i_index]);
+      }
+      incoming_arc_indices_from_sources.push_back(
+          incoming_arc_indices_to[label_i_index]);
+      ++num_labels_to;
+      if (lengths_from_sources.size() >= max_num_created_labels) {
+        return;
+      }
+    }
+  }
+}
+
+template <class GraphType>
+#if __cplusplus >= 202002L
+  requires DagGraphType<GraphType>
+#endif
+typename GraphType::ArcIndex
+ConstrainedShortestPathsOnDagWrapper<GraphType>::MergeHalfRuns(
+    const GraphType& graph, absl::Span<const double> arc_lengths,
+    absl::Span<const std::vector<double>> arc_resources,
+    absl::Span<const double> max_resources,
+    const std::vector<NodeIndex> sub_node_indices[2],
+    const std::vector<double> lengths_from_sources[2],
+    const std::vector<std::vector<double>> resources_from_sources[2],
+    const std::vector<int> first_label[2], const std::vector<int> num_labels[2],
+    LabelPair& best_label_pair) {
+  const std::vector<NodeIndex>& forward_sub_node_indices =
+      sub_node_indices[FORWARD];
+  absl::Span<const double> forward_lengths = lengths_from_sources[FORWARD];
+  const std::vector<std::vector<double>>& forward_resources =
+      resources_from_sources[FORWARD];
+  absl::Span<const int> forward_first_label = first_label[FORWARD];
+  absl::Span<const int> forward_num_labels = num_labels[FORWARD];
+  const std::vector<NodeIndex>& backward_sub_node_indices =
+      sub_node_indices[BACKWARD];
+  absl::Span<const double> backward_lengths = lengths_from_sources[BACKWARD];
+  const std::vector<std::vector<double>>& backward_resources =
+      resources_from_sources[BACKWARD];
+  absl::Span<const int> backward_first_label = first_label[BACKWARD];
+  absl::Span<const int> backward_num_labels = num_labels[BACKWARD];
+  ArcIndex merging_arc_index = -1;
+  // TODO(b/328745703): These operations could be run in parallel with a mutex
+  // over `best_label_pair`.
+  for (ArcIndex arc_index = 0; arc_index < graph.num_arcs(); ++arc_index) {
+    const NodeIndex sub_from = forward_sub_node_indices[graph.Tail(arc_index)];
+    if (sub_from == -1) {
+      continue;
+    }
+    const NodeIndex sub_to = backward_sub_node_indices[graph.Head(arc_index)];
+    if (sub_to == -1) {
+      continue;
+    }
+    const int num_labels_from = forward_num_labels[sub_from];
+    if (num_labels_from == 0) {
+      continue;
+    }
+    const int num_labels_to = backward_num_labels[sub_to];
+    if (num_labels_to == 0) {
+      continue;
+    }
+    const double arc_length = arc_lengths[arc_index];
+    DCHECK(arc_length != -std::numeric_limits<double>::infinity());
+    if (arc_length == std::numeric_limits<double>::infinity()) {
+      continue;
+    }
+    const int first_label_from = forward_first_label[sub_from];
+    const int first_label_to = backward_first_label[sub_to];
+    for (int label_to_index = first_label_to;
+         label_to_index < first_label_to + num_labels_to; ++label_to_index) {
+      const double length_to = backward_lengths[label_to_index];
+      if (arc_length + length_to >= best_label_pair.length) {
+        continue;
+      }
+      for (int label_from_index = first_label_from;
+           label_from_index < first_label_from + num_labels_from;
+           ++label_from_index) {
+        const double length_from = forward_lengths[label_from_index];
+        if (length_from + arc_length + length_to >= best_label_pair.length) {
+          continue;
+        }
+        bool path_is_feasible = true;
+        for (int r = 0; r < num_resources_; ++r) {
+          DCHECK_GE(arc_resources[r][arc_index], 0.0);
+          if (forward_resources[r][label_from_index] +
+                  arc_resources[r][arc_index] +
+                  backward_resources[r][label_to_index] >
               max_resources[r]) {
             path_is_feasible = false;
             break;
@@ -338,67 +830,53 @@ PathWithLength ConstrainedShortestPathsOnDagWrapper<
         if (!path_is_feasible) {
           continue;
         }
-        Label label_to = label_from;
-        label_to.length += arc_length;
-        for (int r = 0; r < max_resources.size(); ++r) {
-          label_to.resources[r] += (*arc_resources_)[r][arc_index];
-        }
-        label_to.incoming_arc = arc_index;
-        labels_from_sources_.push_back(label_to);
-        ++num_labels;
-        if (node_is_destination_[to]) {
-          if (best_length > label_to.length) {
-            best_length = label_to.length;
-            best_label_index = labels_from_sources_.size() - 1;
-          }
-        }
-        if (labels_from_sources_.size() >= max_num_created_labels_) {
-          const std::vector<ArcIndex> arc_path =
-              BestArcPathEndingWith(best_label_index);
-          return {.length = best_length,
-                  .arc_path = arc_path,
-                  .node_path = NodePathImpliedBy(arc_path)};
-        }
+        best_label_pair.length = length_from + arc_length + length_to;
+        best_label_pair.label_index[FORWARD] = label_from_index;
+        best_label_pair.label_index[BACKWARD] = label_to_index;
+        merging_arc_index = arc_index;
       }
     }
   }
-  const std::vector<ArcIndex> arc_path =
-      BestArcPathEndingWith(best_label_index);
-  return {.length = best_length,
-          .arc_path = arc_path,
-          .node_path = NodePathImpliedBy(arc_path)};
+  return merging_arc_index;
 }
 
-template <typename GraphType>
+template <class GraphType>
 #if __cplusplus >= 202002L
   requires DagGraphType<GraphType>
 #endif
 std::vector<typename GraphType::ArcIndex>
-ConstrainedShortestPathsOnDagWrapper<GraphType>::BestArcPathEndingWith(
-    const int label_index) const {
-  if (label_index < 0 || label_index >= labels_from_sources_.size()) {
+ConstrainedShortestPathsOnDagWrapper<GraphType>::ArcPathTo(
+    const int best_label_index, const GraphType& reverse_graph,
+    absl::Span<const double> arc_lengths,
+    absl::Span<const double> lengths_from_sources,
+    absl::Span<const ArcIndex> incoming_arc_indices_from_sources,
+    absl::Span<const int> first_label, absl::Span<const int> num_labels) const {
+  if (best_label_index == -1) {
     return {};
   }
-  int current_label_index = label_index;
+  int current_label_index = best_label_index;
   std::vector<ArcIndex> arc_path;
-  while (labels_from_sources_[current_label_index].incoming_arc != -1) {
-    const Label& label = labels_from_sources_[current_label_index];
-    const NodeIndex node = graph_->Tail(label.incoming_arc);
-    arc_path.push_back(label.incoming_arc);
-    for (int i = node_first_label_[node];
-         i < node_first_label_[node] + node_num_labels_[node]; ++i) {
-      // Since all labels of `labels_from_sources_` are feasible, we can pick
-      // any previous label which satisfies the length equality (irrespective of
-      // the resources it consumes).
-      if (labels_from_sources_[i].length +
-              (*arc_lengths_)[label.incoming_arc] ==
-          label.length) {
-        current_label_index = i;
+  for (int i = 0; i < reverse_graph.num_nodes(); ++i) {
+    const ArcIndex current_arc_index =
+        incoming_arc_indices_from_sources[current_label_index];
+    if (current_arc_index == -1) {
+      break;
+    }
+    arc_path.push_back(current_arc_index);
+    const NodeIndex sub_node = reverse_graph.Head(current_arc_index);
+    const double current_length = lengths_from_sources[current_label_index];
+    for (int label_index = first_label[sub_node];
+         label_index < first_label[sub_node] + num_labels[sub_node];
+         ++label_index) {
+      if (std::abs(lengths_from_sources[label_index] +
+                   arc_lengths[current_arc_index] - current_length) <=
+          kTolerance) {
+        current_label_index = label_index;
         break;
       }
     }
   }
-  absl::c_reverse(arc_path);
+  CHECK_EQ(incoming_arc_indices_from_sources[current_label_index], -1);
   return arc_path;
 }
 
@@ -408,16 +886,16 @@ template <typename GraphType>
 #endif
 std::vector<typename GraphType::NodeIndex>
 ConstrainedShortestPathsOnDagWrapper<GraphType>::NodePathImpliedBy(
-    const std::vector<ArcIndex>& arc_path) const {
+    absl::Span<const ArcIndex> arc_path, const GraphType& graph) const {
   if (arc_path.empty()) {
     return {};
   }
   std::vector<NodeIndex> node_path;
   node_path.reserve(arc_path.size() + 1);
   for (const ArcIndex arc_index : arc_path) {
-    node_path.push_back(graph_->Tail(arc_index));
+    node_path.push_back(graph.Tail(arc_index));
   }
-  node_path.push_back(graph_->Head(arc_path.back()));
+  node_path.push_back(graph.Head(arc_path.back()));
   return node_path;
 }
 
diff --git a/ortools/graph/dag_constrained_shortest_path_test.cc b/ortools/graph/dag_constrained_shortest_path_test.cc
index 8d7df9fba7..89a47b87ce 100644
--- a/ortools/graph/dag_constrained_shortest_path_test.cc
+++ b/ortools/graph/dag_constrained_shortest_path_test.cc
@@ -22,6 +22,7 @@
 #include "absl/algorithm/container.h"
 #include "absl/log/check.h"
 #include "absl/random/random.h"
+#include "absl/strings/str_cat.h"
 #include "absl/types/span.h"
 #include "benchmark/benchmark.h"
 #include "gtest/gtest.h"
@@ -30,6 +31,7 @@
 #include "ortools/graph/dag_shortest_path.h"
 #include "ortools/graph/graph.h"
 #include "ortools/graph/io.h"
+#include "ortools/math_opt/cpp/math_opt.h"
 
 namespace operations_research {
 namespace {
@@ -154,6 +156,21 @@ TEST(ConstrainedShortestPathOnDagTest, GraphWithInefficientEdge) {
                 /*node_path=*/ElementsAre(source, a, destination)));
 }
 
+TEST(ConstrainedShortestPathOnDagTest, NoResources) {
+  const int source = 0;
+  const int destination = 1;
+  const int a = 2;
+  const int num_nodes = 3;
+  const std::vector<ArcWithLengthAndResources> arcs_with_length_and_resources =
+      {{source, a, 5.0, {}}, {a, destination, 3.0, {}}};
+
+  EXPECT_DEATH(
+      ConstrainedShortestPathsOnDag(num_nodes, arcs_with_length_and_resources,
+                                    source, destination,
+                                    /*max_resources=*/{}),
+      "ortools/graph/dag_shortest_path.h");
+}
+
 TEST(ConstrainedShortestPathOnDagTest, Cycle) {
   const int source = 0;
   const int destination = 1;
@@ -278,19 +295,21 @@ TEST(ConstrainedShortestPathsOnDagWrapperTest,
   const int num_nodes = 3;
   util::ListGraph<> graph(num_nodes, /*arc_capacity=*/2);
   std::vector<double> arc_lengths;
-  std::vector<std::vector<double>> arc_resources;
+  std::vector<std::vector<double>> arc_resources(1);
   graph.AddArc(source_1, source_2);
   arc_lengths.push_back(-7.0);
+  arc_resources[0].push_back(2.0);
   graph.AddArc(source_2, destination);
   arc_lengths.push_back(3.0);
+  arc_resources[0].push_back(3.0);
   const std::vector<int> topological_order = {source_1, source_2, destination};
   const std::vector<int> sources = {source_1, source_2};
   const std::vector<int> destinations = {destination};
-  const std::vector<double> max_resources = {};
+  const std::vector<double> max_resources = {6.0};
   ConstrainedShortestPathsOnDagWrapper<util::ListGraph<>>
       constrained_shortest_path_on_dag(&graph, &arc_lengths, &arc_resources,
-                                       &topological_order, &sources,
-                                       &destinations, &max_resources);
+                                       topological_order, sources, destinations,
+                                       &max_resources);
 
   EXPECT_THAT(
       constrained_shortest_path_on_dag.RunConstrainedShortestPathOnDag(),
@@ -324,8 +343,8 @@ TEST(ConstrainedShortestPathsOnDagWrapperTest, MultipleDestinations) {
   const std::vector<double> max_resources = {6.0};
   ConstrainedShortestPathsOnDagWrapper<util::ListGraph<>>
       constrained_shortest_path_on_dag(&graph, &arc_lengths, &arc_resources,
-                                       &topological_order, &sources,
-                                       &destinations, &max_resources);
+                                       topological_order, sources, destinations,
+                                       &max_resources);
 
   EXPECT_THAT(
       constrained_shortest_path_on_dag.RunConstrainedShortestPathOnDag(),
@@ -364,8 +383,8 @@ TEST(ConstrainedShortestPathsOnDagWrapperTest, UpdateArcsLength) {
   const std::vector<double> max_resources = {6.0, 12.0};
   ConstrainedShortestPathsOnDagWrapper<util::ListGraph<>>
       constrained_shortest_path_on_dag(&graph, &arc_lengths, &arc_resources,
-                                       &topological_order, &sources,
-                                       &destinations, &max_resources);
+                                       topological_order, sources, destinations,
+                                       &max_resources);
 
   EXPECT_THAT(
       constrained_shortest_path_on_dag.RunConstrainedShortestPathOnDag(),
@@ -397,30 +416,22 @@ TEST(ConstrainedShortestPathsOnDagWrapperTest, LimitMaximumNumberOfLabels) {
   arc_resources[0].push_back(1.0);
   const std::vector<int> topological_order = {source, a, destination};
   const std::vector<int> sources = {source};
-  const std::vector<int> destinations = {a, destination};
+  const std::vector<int> destinations = {destination};
   const std::vector<double> max_resources = {6.0};
   ConstrainedShortestPathsOnDagWrapper<util::ListGraph<>>
       constrained_shortest_path_on_dag_with_one_label(
-          &graph, &arc_lengths, &arc_resources, &topological_order, &sources,
-          &destinations, &max_resources, /*max_num_created_labels=*/1);
-  ConstrainedShortestPathsOnDagWrapper<util::ListGraph<>>
-      constrained_shortest_path_on_dag_with_two_labels(
-          &graph, &arc_lengths, &arc_resources, &topological_order, &sources,
-          &destinations, &max_resources, /*max_num_created_labels=*/2);
+          &graph, &arc_lengths, &arc_resources, topological_order, sources,
+          destinations, &max_resources, /*max_num_created_labels=*/1);
   ConstrainedShortestPathsOnDagWrapper<util::ListGraph<>>
       constrained_shortest_path_on_dag_with_three_labels(
-          &graph, &arc_lengths, &arc_resources, &topological_order, &sources,
-          &destinations, &max_resources, /*max_num_created_labels=*/3);
+          &graph, &arc_lengths, &arc_resources, topological_order, sources,
+          destinations, &max_resources, /*max_num_created_labels=*/3);
 
   EXPECT_THAT(constrained_shortest_path_on_dag_with_one_label
                   .RunConstrainedShortestPathOnDag(),
               FieldsAre(/*length=*/kInf,
                         /*arc_path=*/IsEmpty(),
                         /*node_path=*/IsEmpty()));
-  EXPECT_THAT(constrained_shortest_path_on_dag_with_two_labels
-                  .RunConstrainedShortestPathOnDag(),
-              FieldsAre(/*length=*/5.0, /*arc_path=*/ElementsAre(0),
-                        /*node_path=*/ElementsAre(source, a)));
   EXPECT_THAT(constrained_shortest_path_on_dag_with_three_labels
                   .RunConstrainedShortestPathOnDag(),
               FieldsAre(/*length=*/3.0, /*arc_path=*/ElementsAre(0, 1),
@@ -478,11 +489,67 @@ std::vector<double> GenerateRandomIntegerValues(
   return arc_values;
 }
 
-// TODO(b/316203403): Remplace bounds with correct value computed using MIP.
+double SolveConstrainedShortestPathUsingIntegerProgramming(
+    const util::StaticGraph<>& graph, absl::Span<const double> arc_lengths,
+    absl::Span<const std::vector<double>> arc_resources,
+    absl::Span<const double> max_resources,
+    absl::Span<const util::StaticGraph<>::NodeIndex> sources,
+    absl::Span<const util::StaticGraph<>::NodeIndex> destinations) {
+  using NodeIndex = util::StaticGraph<>::NodeIndex;
+  using ArcIndex = util::StaticGraph<>::ArcIndex;
+
+  math_opt::Model model;
+  std::vector<math_opt::Variable> arc_variables;
+  std::vector<math_opt::LinearExpression> flow_conservation(graph.num_nodes(),
+                                                            0.0);
+  for (ArcIndex arc_index = 0; arc_index < graph.num_arcs(); ++arc_index) {
+    arc_variables.push_back(model.AddBinaryVariable(absl::StrCat(
+        arc_index, "_", graph.Tail(arc_index), "->", graph.Head(arc_index))));
+    model.set_objective_coefficient(arc_variables[arc_index],
+                                    arc_lengths[arc_index]);
+    flow_conservation[graph.Head(arc_index)] -= arc_variables[arc_index];
+    flow_conservation[graph.Tail(arc_index)] += arc_variables[arc_index];
+  }
+
+  math_opt::LinearExpression all_sources;
+  math_opt::LinearExpression all_destinations;
+  for (NodeIndex node_index = 0; node_index < graph.num_nodes(); ++node_index) {
+    math_opt::LinearExpression net_flow = 0;
+    if (absl::c_linear_search(sources, node_index)) {
+      const math_opt::Variable s = model.AddBinaryVariable();
+      all_sources += s;
+      net_flow += s;
+    }
+    if (absl::c_linear_search(destinations, node_index)) {
+      const math_opt::Variable t = model.AddBinaryVariable();
+      all_destinations += t;
+      net_flow -= t;
+    }
+    model.AddLinearConstraint(flow_conservation[node_index] == net_flow);
+  }
+  model.AddLinearConstraint(all_sources == 1);
+  model.AddLinearConstraint(all_destinations == 1);
+  for (int r = 0; r < max_resources.size(); ++r) {
+    math_opt::LinearExpression variable_resources;
+    for (ArcIndex arc_index = 0; arc_index < graph.num_arcs(); ++arc_index) {
+      variable_resources +=
+          arc_resources[r][arc_index] * arc_variables[arc_index];
+    }
+    model.AddLinearConstraint(variable_resources <= max_resources[r]);
+  }
+  const absl::StatusOr<math_opt::SolveResult> result =
+      math_opt::Solve(model, math_opt::SolverType::kCpSat, {});
+  CHECK_OK(result.status())
+      << util::GraphToString(graph, util::PRINT_GRAPH_ARCS);
+  CHECK_OK(result->termination.EnsureIsOptimal())
+      << util::GraphToString(graph, util::PRINT_GRAPH_ARCS);
+  return result->objective_value();
+}
+
 TEST(ConstrainedShortestPathsOnDagWrapperTest,
      RandomizedStressTestSingleResource) {
   absl::BitGen bit_gen;
-  const int kNumTests = 10000;
+  const int kNumTests = 50;
   for (int test = 0; test < kNumTests; ++test) {
     const int num_nodes = absl::Uniform(bit_gen, 2, 12);
     const int num_arcs = absl::Uniform(
@@ -495,55 +562,23 @@ TEST(ConstrainedShortestPathsOnDagWrapperTest,
                                                    /*max_value=*/10.0,
                                                    /*start_to_end_value=*/1.0);
     const std::vector<int> sources = {0};
-    std::vector<int> destinations(num_nodes);
-    absl::c_iota(destinations, 0);
+    const std::vector<int> destinations = {num_nodes - 1};
     const std::vector<double> max_resources = {15.0};
     ConstrainedShortestPathsOnDagWrapper<util::StaticGraph<>>
         constrained_shortest_path_on_dag(&graph, &arc_lengths, &arc_resources,
-                                         &topological_order, &sources,
-                                         &destinations, &max_resources);
-
-    // Run DAG shortest path on `arc_length` as a LB.
-    ShortestPathsOnDagWrapper<util::StaticGraph<>> shortest_path_length(
-        &graph, &arc_lengths, topological_order);
-
-    // Run DAG shortest path on `arc_resources` as a UB.
-    ShortestPathsOnDagWrapper<util::StaticGraph<>> shortest_path_resources(
-        &graph, &(arc_resources[0]), topological_order);
-    shortest_path_resources.RunShortestPathOnDag(sources);
-
-    // Precompute the expected reached nodes: any node whose resource is <=
-    // max_resources[0].
-    std::vector<int> expected_reached_nodes;
-    for (int node = 0; node < num_nodes; ++node) {
-      if (shortest_path_resources.LengthTo(node) <= max_resources[0]) {
-        expected_reached_nodes.push_back(node);
-      }
-    }
-
-    const int kNumSamples = 20;
+                                         topological_order, sources,
+                                         destinations, &max_resources);
+    const int kNumSamples = 5;
     for (int _ = 0; _ < kNumSamples; ++_) {
-      // Draw random lengths and recompute *un*constrained shortest paths.
       arc_lengths = GenerateRandomIntegerValues(graph);
-      shortest_path_length.RunShortestPathOnDag(sources);
-
       const PathWithLength path_with_length =
           constrained_shortest_path_on_dag.RunConstrainedShortestPathOnDag();
 
-      EXPECT_FALSE(path_with_length.node_path.empty());
-
-      const int best_destination = path_with_length.node_path.back();
-      const std::vector<int> ub_arc_path =
-          shortest_path_resources.ArcPathTo(best_destination);
-      double ub_length = 0.0;
-      for (const int arc : ub_arc_path) {
-        ub_length += arc_lengths[arc];
-      }
-
-      EXPECT_GE(path_with_length.length,
-                shortest_path_length.LengthTo(best_destination))
-          << best_destination;
-      EXPECT_LE(path_with_length.length, ub_length) << best_destination;
+      EXPECT_NEAR(path_with_length.length,
+                  SolveConstrainedShortestPathUsingIntegerProgramming(
+                      graph, arc_lengths, arc_resources, max_resources, sources,
+                      destinations),
+                  1e-5);
 
       ASSERT_FALSE(HasFailure())
           << DUMP_VARS(num_nodes, num_arcs, arc_lengths) << "\n With graph :\n "
@@ -578,25 +613,33 @@ void BM_RandomDag(benchmark::State& state) {
   const std::vector<double> max_resources = {num_nodes * 0.2};
   ConstrainedShortestPathsOnDagWrapper<util::StaticGraph<>>
       constrained_shortest_path_on_dag(&graph, &arc_lengths, &arc_resources,
-                                       &topological_order, &sources,
-                                       &destinations, &max_resources);
+                                       topological_order, sources, destinations,
+                                       &max_resources);
 
+  int total_label_count = 0;
   for (auto _ : state) {
     // Pick a arc lengths scenario at random.
     arc_lengths =
         arc_lengths_scenarios[absl::Uniform(bit_gen, 0, num_scenarios)];
     const PathWithLength path_with_length =
         constrained_shortest_path_on_dag.RunConstrainedShortestPathOnDag();
+    total_label_count += constrained_shortest_path_on_dag.label_count();
     CHECK_GE(path_with_length.length, 0.0);
     CHECK_LE(path_with_length.length, 10000.0);
   }
-  state.SetItemsProcessed(state.iterations() * (num_nodes + num_arcs));
+  state.SetItemsProcessed(std::max(1, total_label_count));
 }
 
 BENCHMARK(BM_RandomDag)
     ->ArgPair(1 << 6, 4)
     ->ArgPair(1 << 6, 16)
-    ->ArgPair(1 << 9, 4);
+    ->ArgPair(1 << 9, 4)
+    ->ArgPair(1 << 9, 16)
+    ->ArgPair(1 << 12, 4)
+    ->ArgPair(1 << 12, 16)
+    ->ArgPair(1 << 15, 4)
+    ->ArgPair(1 << 15, 16)
+    ->ArgPair(1 << 18, 4);
 
 // -----------------------------------------------------------------------------
 // Debug tests.
@@ -654,6 +697,17 @@ TEST(ConstrainedShortestPathOnDagTest, NegativeMaxResource) {
                "negative");
 }
 
+TEST(ConstrainedShortestPathOnDagTest, SourceIsDestination) {
+  const int source = 0;
+  const int num_nodes = 1;
+
+  EXPECT_DEATH(
+      ConstrainedShortestPathsOnDag(
+          num_nodes, /*arcs_with_length_and_resources=*/{}, source, source,
+          /*max_resources=*/{0.0}),
+      "source and destination");
+}
+
 TEST(ConstrainedShortestPathsOnDagWrapperTest, ValidateTopologicalOrder) {
   const int source = 0;
   const int destination = 1;
@@ -670,8 +724,8 @@ TEST(ConstrainedShortestPathsOnDagWrapperTest, ValidateTopologicalOrder) {
   const std::vector<double> max_resources = {0.0};
 
   EXPECT_DEATH(ConstrainedShortestPathsOnDagWrapper<util::ListGraph<>>(
-                   &graph, &arc_lengths, &arc_resources, &topological_order,
-                   &sources, &destinations, &max_resources),
+                   &graph, &arc_lengths, &arc_resources, topological_order,
+                   sources, destinations, &max_resources),
                "Invalid topological order");
 }
 #endif  // NDEBUG