diff --git a/doc/bibliography.html b/doc/bibliography.html
index b541abf62..e79333f52 100644
--- a/doc/bibliography.html
+++ b/doc/bibliography.html
@@ -453,6 +453,11 @@ <H2>Bibliography</H2>
 <em>Data Structures for Weighted Matching and Nearest Common Ancestors with Linking</em><br>
 Proceedings of the First Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 434-443, 1990.
 
+<p></p><dt><a name="galil86">77</a>
+<dd>Zvi Galil<br>
+<em><a href="https://dl.acm.org/doi/pdf/10.1145/6462.6502">Efficient Algorithms for Finding Maximum Matching in Graphs</a></em><br>
+ACM Computing Surveys (CSUR), 18(1), 23-38, 1986.
+
 </dl>
 
 <br>
diff --git a/doc/maximum_weighted_matching.html b/doc/maximum_weighted_matching.html
index ab40f1f45..ca3f221a7 100644
--- a/doc/maximum_weighted_matching.html
+++ b/doc/maximum_weighted_matching.html
@@ -1,5 +1,6 @@
 <html><head><!--
      Copyright 2018 Yi Ji
+     Copyright 2025 Joris van Rantwijk
 
      Use, modification and distribution is subject to the Boost Software
      License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
@@ -20,6 +21,9 @@ <h1>
 template &lt;typename Graph, typename MateMap, typename VertexIndexMap&gt;
 void maximum_weighted_matching(const Graph&amp; g, MateMap mate, VertexIndexMap vm);
 
+template &lt;typename Graph, typename MateMap, typename VertexIndexMap, typename EdgeWeightMap&gt;
+void maximum_weighted_matching(const Graph&amp; g, MateMap mate, VertexIndexMap vm, EdgeWeightMap weights);
+
 template &lt;typename Graph, typename MateMap&gt;
 void brute_force_maximum_weighted_matching(const Graph&amp; g, MateMap mate);
 
@@ -56,77 +60,126 @@ <h1>
 <a href="../../property_map/doc/ReadWritePropertyMap.html">ReadWritePropertyMap</a>
 that maps vertices to vertices. In the mapping returned, each vertex is either mapped
 to the vertex it's matched to, or to <tt>graph_traits&lt;Graph&gt;::null_vertex()</tt> if it
-doesn't participate in the matching. If no <tt>VertexIndexMap</tt> is provided, both functions
-assume that the <tt>VertexIndexMap</tt> is provided as an internal graph property accessible
-by calling <tt>get(vertex_index, g)</tt>.
+doesn't participate in the matching.
+</p>
+
+<h3>Algorithm Description</h3>
 
 <p>
-The maximum weighted matching problem was solved by Edmonds in [<a href="bibliography.html#edmonds65:_max_weighted_match">74</a>].
-The implementation of <tt>maximum_weighted_matching</tt> followed Chapter 6, Section 10 of [<a href="bibliography.html#lawler76:_comb_opt">20</a>] and
-was written in a consistent style with <tt>edmonds_maximum_cardinality_matching</tt> because of their algorithmic similarity.
-In addition, a brute-force verifier <tt>brute_force_maximum_weighted_matching</tt> simply searches all possible matchings in any graph and selects one with the maximum weight sum.
-
-</p><h3>Algorithm Description</h3>
-
-Primal-dual method in linear programming is introduced to solve weighted matching problems. Edmonds proved that for any graph,
-the maximum number of edges in a matching is equal to the minimum capacity of an odd-set cover; this further enable us to prove a max-min duality theorem for weighted matching.
-Let <i>H<sub>k-1</sub></i> denote any graph obtained from G by contracting odd sets of three or more nodes and deleting single nodes,
-where the capacity of the family of odd sets (not necessarily a cover of G) is <i>k-1</i>. Let <i>X<sub>k</sub></i> denote any matching containing <i>k</i> edges.
-Each edge <i>(i, j)</i> has a weight <i>w<sub>ij</sub></i>. We have:
-<math>
-max<i><sub>X<sub>k</sub></sub></i> min {<i>w<sub>ij</sub>|(i, j) ϵ X<sub>k</sub></i>} = min<i><sub>H<sub>k-1</sub></sub></i> max {<i>w<sub>ij</sub>|(i, j) ϵ H<sub>k-1</sub></i>}.
-</math>
-This matching duality theorem gives an indication of how the matching problem should be formulated as a linear programming problem. That is,
-the theorem suggests a set of linear inequalities which are satisfied by any matching, and it is anticipated that these inequalities describe a convex polyhedron
-with integer vertices corresponding to feasible matchings.
+The maximum weighted matching problem was solved by Edmonds in
+[<a href="bibliography.html#edmonds65:_max_weighted_match">74</a>].
+Subsequently, Gabow [<a href="bibliography.html#gabow76">75</a>] and
+Lawler [<a href="bibliography.html#lawler76:_comb_opt">20</a>] developed methods
+to reduce the run time from <i>O(V<sup>4</sup>)</i> to <i>O(V<sup>3</sup>)</i>.
+The implementation of <tt>maximum_weighted_matching</tt> follows a description
+of the <i>O(V<sup>3</sup>)</i> algorithm by Galil [<a href="bibliography.html#galil86">77</a>].
+</p>
 
 <p>
-For <tt>maximum_weighted_matching</tt>, the management of blossoms is much more involved than in the case of <tt>max_cardinality_matching</tt>.
-It is not sufficient to record only the outermost blossoms. When an outermost blossom is expanded,
-it is necessary to know which blossom are nested immediately with it, so that these blossoms can be restored to the status of the outermost blossoms.
-When augmentation occurs, blossoms with strictly positive dual variables must be maintained for use in the next application of the labeling procedure.
+Starting from an empty matching, the algorithm repeatedly augments the matching until no further improvement is possible.
+Alternating trees and blossoms are used to find an augmenting path, in a way that is similar to maximum cardinality matching.
+A linear programming technique is used to ensure that the selected augmenting path has maximum weight.
+This involves assigning dual variables to vertices and blossoms, and computing slack values for edges.
+</p>
 
 <p>
-The outline of the algorithm is as follow:
+The outline of the algorithm is as follows:
+</p>
 
 <ol start="0">
-<li>Start with an empty matching and initialize dual variables as a half of maximum edge weight.</li>
-<li>(Labeling) Root an alternate tree at each exposed node, and proceed to construct alternate trees by labeling, using only edges with zero slack value.
+<li>Start with an empty matching and initialize dual variables to half of the maximum edge weight.</li>
+<li>(Labeling) Root a tree at each non-matched vertex, and proceed to construct alternating trees by labeling, using only edges with zero slack value.
 If an augmenting path is found, go to step 2. If a blossom is formed, go to step 3. Otherwise, go to step 4. </li>
-<li>(Augmentation) Find the augmenting path, tracing the path through shrunken blossoms. Augment the matching,
-correct labels on nodes in the augmenting path, expand blossoms with zero dual variables and remove labels from all base nodes. Go to step 1.</li>
-<li>(Blossoming) Determine the membership and base node of the new blossom and supply missing labels for all non-base nodes in the blossom.
+<li>(Augmentation) Find the augmenting path, tracing the path through shrunken blossoms. Augment the matching.
+Deconstruct and unlabel the alternating trees traversed by the augmenting path.
+Go to step 1.</li>
+<li>(Blossoming) Identify the blossoms in the alternating cycle and shrink these into a newly formed blossom.
 Return to step 1.</li>
-<li>(Revision of Dual Solution) Adjust the dual variables based on the primal-dual method. Go to step 1 or halt, accordingly.</li>
+<li>(Updating dual variables) Adjust the dual variables based on the primal-dual method.
+If this enables further growth of the alternating trees, return to step 1.
+Otherwise, halt the algorithm; the matching has maximum weight.</li>
 </ol>
 
-Note that in <tt>maximum_weighted_matching</tt>, all edge weights are multiplied by 4, so that all dual variables always remain as integers if all edge weights are integers.
-Unlike <tt>max_cardinality_matching</tt>, the initial matching and augmenting path finder are not parameterized,
-because the algorithm maintains blossoms, dual variables and node labels across all augmentations.
+<p>
+The implementation deviates from the description in [<a href="bibliography.html#galil86">77</a>] on a few points:
+</p>
+
+<ul>
+<li>
+After augmenting the matching, [<a href="bibliography.html#galil86">77</a>] expands all blossoms with zero dual.
+This is correct but potentially unnecessary; some of the expanded blossoms will be recreated during the next stage.
+The implementation holds off on the expansion of such blossoms.
+In any scenario where a blossom <i>must</i> be expanded to make progress, label T will be assigned to it.
+At that point, the blossom will be expanded through a zero-delta dual step.
+The additional cost is <i>O(V)</i> per expanded blossom, and the total number of blossom expansions is <i>O(V<sup>2</sup>)</i>,
+therefore this does not change the overall time complexity of the algorithm.
+</li>
+<li>
+The algorithm records a minimum-slack edge <i>e<sub>k,l</sub></i> for every pair of blossoms <i>B<sub>k</sub>, B<sub>l</sub></i>.
+Storing these edges in a 2-dimensional array would increase the space complexity to <i>O(V<sup>2</sup>)</i>;
+an unfortunate cost for large sparse graphs.
+The implementation of <tt>maximum_weighted_matching</tt> avoids this by maintaining, for each blossom <i>B<sub>k</sub></i>,
+a list <tt>B<sub>k</sub>.best_edge_set</tt> that contains the non-empty elements of <i>e<sub>k,l</sub></i>,
+plus edges between <i>B<sub>k</sub></i> and <i>B<sub>l</sub></i> discovered after shrinking <i>B<sub>k</sub></i>.
+The combined length of these lists is <i>O(E)</i>, since any edge is contained in at most two lists.
+The algorithm only ever examines <i>e<sub>k,l</sub></i> during the process of merging <i>B<sub>k</sub></i> into a larger blossom.
+In that case, it sequentially examines <i>e<sub>k,l</sub></i> for all <i>l</i>, taking time <i>O(V)</i>.
+The implementation handles this by examining every element of <tt>B<sub>k</sub>.best_edge_set</tt>,
+taking time <i>O(V)</i> plus a constant time per newly discovered edge.
+This does not change the overall time complexity of the algorithm.
+</li>
+<li>
+After augmenting, [<a href="bibliography.html#galil86">77</a>] removes all labels and alternating trees.
+This keeps the algorithm simple, but it causes many trees and labels to be rediscovered at the beginning of the new stage.
+To avoid such redundent work, <tt>maximum_weighted_matching</tt> removes labels only from the two alternating trees that are touched by the augmenting path.
+All other labels and trees are reused in the next stage.
+Erasing labels from a subset of the graph, requires updating the data structures that track minimum-slack edges.
+This involves examining all edges that either touch a blossom which loses its label or touch a vertex which depended on a lost label.
+These steps take <i>O(V + E)</i> per stage, adding up to <i>O(V<sup>3</sup>)</i> in total.
+Although this strategy does not improve the worst-case time complexity of the algorithm, it achieves a significant speedup in benchmarks on random graphs.
+</li>
+</ul>
 
-The algorithm's time complexity is reduced from <i>O(V<sup>4</sup>)</i> (naive implementation of [<a href="bibliography.html#edmonds65:_max_weighted_match">74</a>])
-to <i>O(V<sup>3</sup>)</i>, by a delicate labeling procedure [<a href="bibliography.html#gabow76">75</a>] to avoid re-scanning labels after revision of the dual solution.
-Special variables <i>pi, tau, gamma</i> and two arrays <i>critical_edge, tau_idx</i> are introduced for this purpose.
-Please refer to [<a href="bibliography.html#lawler76:_comb_opt">20</a>] and code comments for more implementation details.
+<p>
+If edge weights are integers, the algorithm uses only integer computations.
+This is achieved by internally multiplying edge weights by 2, which ensures that all dual variables are also integers.
+</p>
+
+<p>
+A brute-force implementation <tt>brute_force_maximum_weighted_matching</tt> is also provided.
+This algorithm simply searches all possible matchings and selects one with the maximum weight sum.
+</p>
 
 </p><h3>Where Defined</h3>
 
 <p>
 <a href="../../../boost/graph/maximum_weighted_matching.hpp"><tt>boost/graph/maximum_weighted_matching.hpp</tt></a>
+</p>
 
-</p><h3>Parameters</h3>
+<h3>Parameters</h3>
 
 IN: <tt>const Graph&amp; g</tt>
 <blockquote>
 An undirected graph. The graph type must be a model of
 <a href="VertexAndEdgeListGraph.html">Vertex and Edge List Graph</a> and
 <a href="IncidenceGraph.html">Incidence Graph</a>.
-The edge property of the graph <tt>property_map&lt;Graph, edge_weight_t&gt;</tt> must exist and have numeric value type.<br>
+The graph may not contain parallel edges.
 </blockquote>
 
 IN: <tt>VertexIndexMap vm</tt>
 <blockquote>
-Must be a model of <a href="../../property_map/doc/ReadablePropertyMap.html">ReadablePropertyMap</a>, mapping vertices to integer indices.
+Must be a model of <a href="../../property_map/doc/ReadablePropertyMap.html">ReadablePropertyMap</a>,
+mapping vertices to integer indices in the range <tt>[0, num_vertices(g))</tt>.
+If this parameter is not explicitly specified, it is obtained from the internal vertex property of the graph
+by calling <tt>get(vertex_index, g)</tt>.
+</blockquote>
+
+IN: <tt>EdgeWeightMap weights</tt>
+<blockquote>
+Must be a model of <a href="../../property_map/doc/ReadablePropertyMap.html">ReadablePropertyMap</a>, mapping edges to weights.
+Edge weights must be integers or floating point values.
+If this parameter is not explicitly specified, it is obtained from the internal edge property of the graph
+by calling <tt>get(edge_weight, g)</tt>.
 </blockquote>
 
 OUT: <tt>MateMap mate</tt>
@@ -140,24 +193,29 @@ <h3>Complexity</h3>
 
 <p>
 Let <i>m</i> and <i>n</i> be the number of edges and vertices in the input graph, respectively. Assuming the
-<tt>VertexIndexMap</tt> supplied allows constant-time lookup, the time complexity for
-<tt>maximum_weighted_matching</tt> is <i>O(n<sup>3</sup>)</i>. For <tt>brute_force_maximum_weighted_matching</tt>, the time complexity is exponential of <i>m</i>.
+<tt>VertexIndexMap</tt> and <tt>EdgeWeightMap</tt> both provide constant-time lookup, the time complexity for
+<tt>maximum_weighted_matching</tt> is <i>O(n<sup>3</sup>)</i>.
+For <tt>brute_force_maximum_weighted_matching</tt>, the time complexity is exponential in <i>m</i>.
+</p>
+
+<p>
 Note that the best known time complexity for maximum weighted matching in general graph
 is <i>O(nm+n<sup>2</sup>log(n))</i> by [<a href="bibliography.html#gabow90">76</a>], but relies on an
 efficient algorithm for solving nearest ancestor problem on trees, which is not provided in Boost C++ libraries.
-</p><p>
+</p>
 
-</p><h3>Example</h3>
+<h3>Example</h3>
 
-<p> The file <a href="../example/weighted_matching_example.cpp"><tt>example/weighted_matching_example.cpp</tt></a>
+<p>The file <a href="../example/weighted_matching_example.cpp"><tt>example/weighted_matching_example.cpp</tt></a>
 contains an example.
+</p>
 
-<br>
-</p><hr>
+<hr>
 <table>
 <tbody><tr valign="top">
-<td nowrap="nowrap">Copyright © 2018</td><td>
-Yi Ji (<a href="mailto:jiy@pku.edu.cn">jiy@pku.edu.cn</a>)<br>
+<td nowrap="nowrap">Copyright © 2018, 2025</td><td>
+Yi Ji (<a href="mailto:jiy@pku.edu.cn">jiy@pku.edu.cn</a>),
+Joris van Rantwijk
 </td></tr></tbody></table>
 
 </body></html>
diff --git a/example/Jamfile.v2 b/example/Jamfile.v2
index 30c730a2f..e316a1914 100644
--- a/example/Jamfile.v2
+++ b/example/Jamfile.v2
@@ -163,6 +163,7 @@ run vertex-name-property.cpp : $(TEST_DIR)/makefile-dependencies.dat $(TEST_DIR)
 run vf2_sub_graph_iso_example.cpp ;
 run vf2_sub_graph_iso_multi_example.cpp ;
 run visitor.cpp ;
+run weighted_matching_example.cpp ;
 run write_graphviz.cpp ;
 
 #
diff --git a/include/boost/graph/maximum_weighted_matching.hpp b/include/boost/graph/maximum_weighted_matching.hpp
index 098a9c2aa..86a133bd0 100644
--- a/include/boost/graph/maximum_weighted_matching.hpp
+++ b/include/boost/graph/maximum_weighted_matching.hpp
@@ -1,5 +1,6 @@
 //=======================================================================
 // Copyright (c) 2018 Yi Ji
+// Copyright (c) 2025 Joris van Rantwijk
 //
 // Distributed under the Boost Software License, Version 1.0.
 // (See accompanying file LICENSE_1_0.txt or copy at
@@ -10,23 +11,35 @@
 #ifndef BOOST_GRAPH_MAXIMUM_WEIGHTED_MATCHING_HPP
 #define BOOST_GRAPH_MAXIMUM_WEIGHTED_MATCHING_HPP
 
-#include <algorithm> // for std::iter_swap
-#include <boost/shared_ptr.hpp>
-#include <boost/make_shared.hpp>
-#include <boost/graph/max_cardinality_matching.hpp>
+#include <boost/assert.hpp>
+#include <boost/optional.hpp>
+#include <boost/property_map/property_map.hpp>
+#include <boost/range/iterator_range_core.hpp>
+#include <boost/graph/exception.hpp>
+#include <boost/graph/graph_concepts.hpp>
+#include <boost/graph/max_cardinality_matching.hpp>  // for empty_matching
+
+#include <algorithm>
+#include <deque>
+#include <limits>
+#include <list>
+#include <stack>
+#include <tuple>  // for std::tie
+#include <utility>  // for std::pair, std::swap
+#include <vector>
 
 namespace boost
 {
-template < typename Graph, typename MateMap, typename VertexIndexMap >
+
+template <typename Graph, typename MateMap, typename VertexIndexMap>
 typename property_traits<
-    typename property_map< Graph, edge_weight_t >::type >::value_type
+    typename property_map<Graph, edge_weight_t>::type>::value_type
 matching_weight_sum(const Graph& g, MateMap mate, VertexIndexMap vm)
 {
-    typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;
-    typedef
-        typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t;
-    typedef typename property_traits< typename property_map< Graph,
-        edge_weight_t >::type >::value_type edge_property_t;
+    using vertex_iterator_t = typename graph_traits<Graph>::vertex_iterator;
+    using vertex_descriptor_t = typename graph_traits<Graph>::vertex_descriptor;
+    using edge_property_t = typename property_traits<
+        typename property_map<Graph, edge_weight_t>::type>::value_type;
 
     edge_property_t weight_sum = 0;
     vertex_iterator_t vi, vi_end;
@@ -34,1280 +47,1445 @@ matching_weight_sum(const Graph& g, MateMap mate, VertexIndexMap vm)
     for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
     {
         vertex_descriptor_t v = *vi;
-        if (get(mate, v) != graph_traits< Graph >::null_vertex()
+        if (get(mate, v) != graph_traits<Graph>::null_vertex()
             && get(vm, v) < get(vm, get(mate, v)))
             weight_sum += get(edge_weight, g, edge(v, mate[v], g).first);
     }
     return weight_sum;
 }
 
-template < typename Graph, typename MateMap >
+template <typename Graph, typename MateMap>
 inline typename property_traits<
-    typename property_map< Graph, edge_weight_t >::type >::value_type
+    typename property_map<Graph, edge_weight_t>::type>::value_type
 matching_weight_sum(const Graph& g, MateMap mate)
 {
     return matching_weight_sum(g, mate, get(vertex_index, g));
 }
 
-template < typename Graph, typename MateMap, typename VertexIndexMap >
-class weighted_augmenting_path_finder
+// brute-force matcher searches all possible combinations of matched edges to
+// get the maximum weighted matching which can be used for testing on small
+// graphs (within dozens vertices)
+template <typename Graph, typename MateMap, typename VertexIndexMap>
+class brute_force_matching
 {
 public:
-    template < typename T > struct map_vertex_to_
+    using vertex_descriptor_t = typename graph_traits<Graph>::vertex_descriptor;
+    using vertex_iterator_t = typename graph_traits<Graph>::vertex_iterator;
+    using vertex_vec_iter_t =
+        typename std::vector<vertex_descriptor_t>::iterator;
+    using edge_iterator_t = typename graph_traits<Graph>::edge_iterator;
+    using vertex_to_vertex_map_t =
+        boost::iterator_property_map<vertex_vec_iter_t, VertexIndexMap>;
+
+    brute_force_matching(
+        const Graph& arg_g, MateMap arg_mate, VertexIndexMap arg_vm)
+    : g(&arg_g)
+    , vm(arg_vm)
+    , mate_vector(num_vertices(*g))
+    , best_mate_vector(num_vertices(*g))
+    , mate(mate_vector.begin(), vm)
+    , best_mate(best_mate_vector.begin(), vm)
     {
-        typedef boost::iterator_property_map<
-            typename std::vector< T >::iterator, VertexIndexMap >
-            type;
-    };
-    typedef typename graph::detail::VERTEX_STATE vertex_state_t;
-    typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;
-    typedef
-        typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t;
-    typedef typename std::vector< vertex_descriptor_t >::const_iterator
-        vertex_vec_iter_t;
-    typedef
-        typename graph_traits< Graph >::out_edge_iterator out_edge_iterator_t;
-    typedef typename graph_traits< Graph >::edge_descriptor edge_descriptor_t;
-    typedef typename graph_traits< Graph >::edge_iterator edge_iterator_t;
-    typedef typename property_traits< typename property_map< Graph,
-        edge_weight_t >::type >::value_type edge_property_t;
-    typedef std::deque< vertex_descriptor_t > vertex_list_t;
-    typedef std::vector< edge_descriptor_t > edge_list_t;
-    typedef typename map_vertex_to_< vertex_descriptor_t >::type
-        vertex_to_vertex_map_t;
-    typedef
-        typename map_vertex_to_< edge_property_t >::type vertex_to_weight_map_t;
-    typedef typename map_vertex_to_< bool >::type vertex_to_bool_map_t;
-    typedef typename map_vertex_to_< std::pair< vertex_descriptor_t,
-        vertex_descriptor_t > >::type vertex_to_pair_map_t;
-    typedef
-        typename map_vertex_to_< std::pair< edge_descriptor_t, bool > >::type
-            vertex_to_edge_map_t;
-    typedef typename map_vertex_to_< vertex_to_edge_map_t >::type
-        vertex_pair_to_edge_map_t;
-
-    class blossom
+        vertex_iterator_t vi, vi_end;
+        for (boost::tie(vi, vi_end) = vertices(*g); vi != vi_end; ++vi)
+            best_mate[*vi] = mate[*vi] = get(arg_mate, *vi);
+    }
+
+    template <typename PropertyMap> void find_matching(PropertyMap pm)
     {
-    public:
-        typedef boost::shared_ptr< blossom > blossom_ptr_t;
-        std::vector< blossom_ptr_t > sub_blossoms;
-        edge_property_t dual_var;
-        blossom_ptr_t father;
-
-        blossom() : dual_var(0), father(blossom_ptr_t()) {}
-
-        // get the base vertex of a blossom by recursively getting
-        // its base sub-blossom, which is always the first one in
-        // sub_blossoms because of how we create and maintain blossoms
-        virtual vertex_descriptor_t get_base() const
-        {
-            const blossom* b = this;
-            while (!b->sub_blossoms.empty())
-                b = b->sub_blossoms[0].get();
-            return b->get_base();
-        }
+        edge_iterator_t ei;
+        boost::tie(ei, ei_end) = edges(*g);
+        select_edge(ei);
 
-        // set a sub-blossom as a blossom's base by exchanging it
-        // with its first sub-blossom
-        void set_base(const blossom_ptr_t& sub)
-        {
-            for (blossom_iterator_t bi = sub_blossoms.begin();
-                 bi != sub_blossoms.end(); ++bi)
-            {
-                if (sub.get() == bi->get())
-                {
-                    std::iter_swap(sub_blossoms.begin(), bi);
-                    break;
-                }
-            }
-        }
+        vertex_iterator_t vi, vi_end;
+        for (boost::tie(vi, vi_end) = vertices(*g); vi != vi_end; ++vi)
+            put(pm, *vi, best_mate[*vi]);
+    }
+
+private:
+    const Graph* g;
+    VertexIndexMap vm;
+    std::vector<vertex_descriptor_t> mate_vector, best_mate_vector;
+    vertex_to_vertex_map_t mate, best_mate;
+    edge_iterator_t ei_end;
 
-        // get all vertices inside recursively
-        virtual std::vector< vertex_descriptor_t > vertices() const
+    void select_edge(edge_iterator_t ei)
+    {
+        if (ei == ei_end)
         {
-            std::vector< vertex_descriptor_t > all_vertices;
-            for (typename std::vector< blossom_ptr_t >::const_iterator bi
-                 = sub_blossoms.begin();
-                 bi != sub_blossoms.end(); ++bi)
+            if (matching_weight_sum(*g, mate)
+                > matching_weight_sum(*g, best_mate))
             {
-                std::vector< vertex_descriptor_t > some_vertices
-                    = (*bi)->vertices();
-                all_vertices.insert(all_vertices.end(), some_vertices.begin(),
-                    some_vertices.end());
+                vertex_iterator_t vi, vi_end;
+                for (boost::tie(vi, vi_end) = vertices(*g); vi != vi_end; ++vi)
+                    best_mate[*vi] = mate[*vi];
             }
-            return all_vertices;
+            return;
         }
-    };
 
-    // a trivial_blossom only has one vertex and no sub-blossom;
-    // for each vertex v, in_blossom[v] is the trivial_blossom that contains it
-    // directly
-    class trivial_blossom : public blossom
-    {
-    public:
-        trivial_blossom(vertex_descriptor_t v) : trivial_vertex(v) {}
-        virtual vertex_descriptor_t get_base() const { return trivial_vertex; }
+        vertex_descriptor_t v, w;
+        v = source(*ei, *g);
+        w = target(*ei, *g);
 
-        virtual std::vector< vertex_descriptor_t > vertices() const
+        select_edge(++ei);
+
+        if (mate[v] == graph_traits<Graph>::null_vertex()
+            && mate[w] == graph_traits<Graph>::null_vertex())
         {
-            std::vector< vertex_descriptor_t > all_vertices;
-            all_vertices.push_back(trivial_vertex);
-            return all_vertices;
+            mate[v] = w;
+            mate[w] = v;
+            select_edge(ei);
+            mate[v] = mate[w] = graph_traits<Graph>::null_vertex();
         }
+    }
+};
 
-    private:
-        vertex_descriptor_t trivial_vertex;
-    };
-
-    typedef boost::shared_ptr< blossom > blossom_ptr_t;
-    typedef typename std::vector< blossom_ptr_t >::iterator blossom_iterator_t;
-    typedef
-        typename map_vertex_to_< blossom_ptr_t >::type vertex_to_blossom_map_t;
-
-    weighted_augmenting_path_finder(
-        const Graph& arg_g, MateMap arg_mate, VertexIndexMap arg_vm)
-    : g(arg_g)
-    , vm(arg_vm)
-    , null_edge(std::pair< edge_descriptor_t, bool >(
-          num_edges(g) == 0 ? edge_descriptor_t() : *edges(g).first, false))
-    , mate_vector(num_vertices(g))
-    , label_S_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
-    , label_T_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
-    , outlet_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
-    , tau_idx_vector(num_vertices(g), graph_traits< Graph >::null_vertex())
-    , dual_var_vector(std::vector< edge_property_t >(
-          num_vertices(g), std::numeric_limits< edge_property_t >::min()))
-    , pi_vector(std::vector< edge_property_t >(
-          num_vertices(g), std::numeric_limits< edge_property_t >::max()))
-    , gamma_vector(std::vector< edge_property_t >(
-          num_vertices(g), std::numeric_limits< edge_property_t >::max()))
-    , tau_vector(std::vector< edge_property_t >(
-          num_vertices(g), std::numeric_limits< edge_property_t >::max()))
-    , in_blossom_vector(num_vertices(g))
-    , old_label_vector(num_vertices(g))
-    , critical_edge_vectors(num_vertices(g),
-          std::vector< std::pair< edge_descriptor_t, bool > >(
-              num_vertices(g), null_edge))
-    ,
-
-        mate(mate_vector.begin(), vm)
-    , label_S(label_S_vector.begin(), vm)
-    , label_T(label_T_vector.begin(), vm)
-    , outlet(outlet_vector.begin(), vm)
-    , tau_idx(tau_idx_vector.begin(), vm)
-    , dual_var(dual_var_vector.begin(), vm)
-    , pi(pi_vector.begin(), vm)
-    , gamma(gamma_vector.begin(), vm)
-    , tau(tau_vector.begin(), vm)
-    , in_blossom(in_blossom_vector.begin(), vm)
-    , old_label(old_label_vector.begin(), vm)
-    {
-        vertex_iterator_t vi, vi_end;
-        edge_iterator_t ei, ei_end;
-
-        edge_property_t max_weight
-            = std::numeric_limits< edge_property_t >::min();
-        for (boost::tie(ei, ei_end) = edges(g); ei != ei_end; ++ei)
-            max_weight = std::max(max_weight, get(edge_weight, g, *ei));
-
-        typename std::vector<
-            std::vector< std::pair< edge_descriptor_t, bool > > >::iterator vei;
+template <typename Graph, typename MateMap, typename VertexIndexMap>
+void brute_force_maximum_weighted_matching(
+    const Graph& g, MateMap mate, VertexIndexMap vm)
+{
+    empty_matching<Graph, MateMap>::find_matching(g, mate);
+    brute_force_matching<Graph, MateMap, VertexIndexMap> brute_force_matcher(
+        g, mate, vm);
+    brute_force_matcher.find_matching(mate);
+}
 
-        for (boost::tie(vi, vi_end) = vertices(g),
-                            vei = critical_edge_vectors.begin();
-             vi != vi_end; ++vi, ++vei)
-        {
-            vertex_descriptor_t u = *vi;
-            mate[u] = get(arg_mate, u);
-            dual_var[u] = 2 * max_weight;
-            in_blossom[u] = boost::make_shared< trivial_blossom >(u);
-            outlet[u] = u;
-            critical_edge_vector.push_back(
-                vertex_to_edge_map_t(vei->begin(), vm));
-        }
+template <typename Graph, typename MateMap>
+inline void brute_force_maximum_weighted_matching(const Graph& g, MateMap mate)
+{
+    brute_force_maximum_weighted_matching(g, mate, get(vertex_index, g));
+}
 
-        critical_edge
-            = vertex_pair_to_edge_map_t(critical_edge_vector.begin(), vm);
+namespace graph
+{
+namespace detail
+{
 
-        init();
+/** Check that vertex indices are unique and in range [0, V). */
+template <typename Graph, typename VertexIndexMap>
+void check_vertex_index_range(const Graph& g, VertexIndexMap vm)
+{
+    using index_t = typename property_traits<VertexIndexMap>::value_type;
+    using unsigned_index_t = typename std::make_unsigned<index_t>::type;
+    auto nv = num_vertices(g);
+    std::vector<bool> got_vertex(nv);
+    for (const auto& x : make_iterator_range(vertices(g)))
+    {
+        index_t i = get(vm, x);
+        if ((i < 0) || (static_cast<unsigned_index_t>(i) >= nv))
+            throw bad_graph("Invalid vertex index.");
+        if (got_vertex[i])
+            throw bad_graph("Duplicate vertex index.");
+        got_vertex[i] = true;
     }
+}
 
-    // return the top blossom where v is contained inside
-    blossom_ptr_t in_top_blossom(vertex_descriptor_t v) const
+/** Check that edge weights are valid. */
+template <typename Graph, typename EdgeWeightMap>
+void check_maximum_weighted_matching_edge_weights(
+    const Graph& g, EdgeWeightMap edge_weights)
+{
+    for (const auto& e : make_iterator_range(edges(g)))
     {
-        blossom_ptr_t b = in_blossom[v];
-        while (b->father != blossom_ptr_t())
-            b = b->father;
-        return b;
+        auto w = get(edge_weights, e);
+        auto max_weight = (std::numeric_limits<decltype(w)>::max)() / 4;
+        if (!(w <= max_weight))  // inverted logic to catch NaN
+            throw bad_graph("Edge weight exceeds maximum supported value.");
     }
+}
 
-    // check if vertex v is in blossom b
-    bool is_in_blossom(blossom_ptr_t b, vertex_descriptor_t v) const
+/** Implementation of the matching algorithm. */
+template <typename Graph, typename VertexIndexMap, typename EdgeWeightMap>
+struct maximum_weighted_matching_context
+{
+    using vertex_t = typename graph_traits<Graph>::vertex_descriptor;
+    using edge_t = typename graph_traits<Graph>::edge_descriptor;
+    using vertices_size_t = typename graph_traits<Graph>::vertices_size_type;
+    using edges_size_t = typename graph_traits<Graph>::edges_size_type;
+    using weight_t = typename property_traits<EdgeWeightMap>::value_type;
+
+    /** Ordered pair of vertices. */
+    using vertex_pair_t = std::pair<vertex_t, vertex_t>;
+
+    /**
+     * List of edges forming an alternating path or alternating cycle.
+     *
+     * The path is defined over top-level blossoms; it skips parts of the path
+     * that are internal to blossoms. Vertex pairs are oriented to match the
+     * direction of the path.
+     */
+    using alternating_path_t = std::deque<vertex_pair_t>;
+
+    /** Top-level blossoms may be labeled "S" or "T" or unlabeled. */
+    enum blossom_label_t { LABEL_NONE = 0, LABEL_S = 1, LABEL_T = 2 };
+
+    struct nontrivial_blossom_t;  // forward declaration
+
+    /**
+     * A blossom is either a single vertex, or an odd-length alternating
+     * cycle over sub-blossoms.
+     */
+    struct blossom_t
     {
-        if (v == graph_traits< Graph >::null_vertex())
-            return false;
-        blossom_ptr_t vb = in_blossom[v]->father;
-        while (vb != blossom_ptr_t())
+        /** Parent of this blossom, or "nullptr" for top-level blossom. */
+        nontrivial_blossom_t* parent;
+
+        /**
+         * Base vertex of this blossom.
+         *
+         * This is the unique vertex inside the blossom which is not matched
+         * to another vertex in the same blossom.
+         */
+        vertex_t base_vertex;
+
+        /** Label S or T or NONE, only valid for top-level blossoms. */
+        blossom_label_t label;
+
+        /** True if this is an instance of nontrivial_blossom. */
+        bool is_nontrivial_blossom;
+
+        /** Edge that attaches this blossom to the alternating tree. */
+        vertex_pair_t tree_edge;
+
+        /** Base vertex of the blossom at the root of the alternating tree. */
+        vertex_t tree_root;
+
+        /** Least-slack edge to a different S-blossom. */
+        optional<edge_t> best_edge;
+
+        /** Initialize a blossom instance. */
+        blossom_t(vertex_t base_vertex = graph_traits<Graph>::null_vertex(),
+                  bool is_nontrivial_blossom = false)
+          : parent(nullptr)
+          , base_vertex(base_vertex)
+          , label(LABEL_NONE)
+          , is_nontrivial_blossom(is_nontrivial_blossom)
+        { }
+
+        /**
+         * Cast this blossom_t instance to nontrivial_blossom_t if possible,
+         * otherwise return "nullptr".
+         */
+        nontrivial_blossom_t* nontrivial()
         {
-            if (vb.get() == b.get())
-                return true;
-            vb = vb->father;
+            // This is much faster than dynamic_cast.
+            return (is_nontrivial_blossom ?
+                static_cast<nontrivial_blossom_t*>(this) : nullptr);
         }
-        return false;
-    }
-
-    // return the base vertex of the top blossom that contains v
-    inline vertex_descriptor_t base_vertex(vertex_descriptor_t v) const
-    {
-        return in_top_blossom(v)->get_base();
-    }
 
-    // add an existed top blossom of base vertex v into new top
-    // blossom b as its sub-blossom
-    void add_sub_blossom(blossom_ptr_t b, vertex_descriptor_t v)
-    {
-        blossom_ptr_t sub = in_top_blossom(v);
-        sub->father = b;
-        b->sub_blossoms.push_back(sub);
-        if (sub->sub_blossoms.empty())
-            return;
-        for (blossom_iterator_t bi = top_blossoms.begin();
-             bi != top_blossoms.end(); ++bi)
+        const nontrivial_blossom_t* nontrivial() const
         {
-            if (bi->get() == sub.get())
-            {
-                top_blossoms.erase(bi);
-                break;
-            }
+            return (is_nontrivial_blossom ?
+                static_cast<const nontrivial_blossom_t*>(this) : nullptr);
         }
-    }
+    };
 
-    // when a top blossom is created or its base vertex getting an S-label,
-    // add all edges incident to this blossom into even_edges
-    void bloom(blossom_ptr_t b)
+    /**
+     * A non-trivial blossom is an odd-length alternating cycle over
+     * 3 or more sub-blossoms.
+     */
+    struct nontrivial_blossom_t : public blossom_t
     {
-        std::vector< vertex_descriptor_t > vertices_of_b = b->vertices();
-        vertex_vec_iter_t vi;
-        for (vi = vertices_of_b.begin(); vi != vertices_of_b.end(); ++vi)
+        struct sub_blossom_t {
+            /** Pointer to sub-blossom. */
+            blossom_t* blossom;
+
+            /** Vertex pair (x, y) where "x" is a vertex in "blossom" and
+                "y" is a vertex in the next sub-blossom. */
+            vertex_pair_t edge;
+        };
+
+        /** List of sub-blossoms, ordered along the alternating cycle. */
+        std::list<sub_blossom_t> subblossoms;
+
+        /** Dual LPP variable for this blossom. */
+        weight_t dual_var;
+
+        /** Least-slack edges to other S-blossoms. */
+        std::list<edge_t> best_edge_set;
+
+        /** Initialize a non-trivial blossom. */
+        nontrivial_blossom_t(
+            const std::vector<blossom_t*>& blossoms,
+            const std::deque<vertex_pair_t>& edges)
+          : blossom_t(blossoms.front()->base_vertex, true)
+          , dual_var(0)
         {
-            out_edge_iterator_t oei, oei_end;
-            for (boost::tie(oei, oei_end) = out_edges(*vi, g); oei != oei_end;
-                 ++oei)
-            {
-                if (target(*oei, g) != *vi && mate[*vi] != target(*oei, g))
-                    even_edges.push_back(*oei);
+            BOOST_ASSERT(blossoms.size() == edges.size());
+            BOOST_ASSERT(blossoms.size() % 2 == 1);
+            BOOST_ASSERT(blossoms.size() >= 3);
+
+            auto blossom_it = blossoms.begin();
+            auto blossom_end = blossoms.end();
+            auto edge_it = edges.begin();
+            while (blossom_it != blossom_end) {
+                subblossoms.push_back({*blossom_it, *edge_it});
+                ++blossom_it;
+                ++edge_it;
             }
         }
-    }
-
-    // assigning a T-label to a non S-vertex, along with outlet and updating pi
-    // value if updated pi[v] equals zero, augment the matching from its mate
-    // vertex
-    void put_T_label(vertex_descriptor_t v, vertex_descriptor_t T_label,
-        vertex_descriptor_t outlet_v, edge_property_t pi_v)
-    {
-        if (label_S[v] != graph_traits< Graph >::null_vertex())
-            return;
 
-        label_T[v] = T_label;
-        outlet[v] = outlet_v;
-        pi[v] = pi_v;
-
-        vertex_descriptor_t v_mate = mate[v];
-        if (pi[v] == 0)
+        /** Find the position of the specified subblossom. */
+        std::pair<vertices_size_t, typename std::list<sub_blossom_t>::iterator>
+        find_subblossom(blossom_t* child)
         {
-            label_T[v_mate] = graph_traits< Graph >::null_vertex();
-            label_S[v_mate] = v;
-            bloom(in_top_blossom(v_mate));
+            vertices_size_t pos = 0;
+            auto it = subblossoms.begin();
+            while (it->blossom != child) {
+                ++it;
+                ++pos;
+                BOOST_ASSERT(it != subblossoms.end());
+            }
+            return std::make_pair(pos, it);
         }
-    }
+    };
 
-    // get the missing T-label for a to-be-expanded base vertex
-    // the missing T-label is the last vertex of the path from outlet[v] to v
-    std::pair< vertex_descriptor_t, vertex_descriptor_t > missing_label(
-        vertex_descriptor_t b_base)
+    /** Specification of a delta step. */
+    struct delta_step_t
     {
-        vertex_descriptor_t missing_outlet = outlet[b_base];
+        /** Type of delta step: 1, 2, 3 or 4. */
+        int kind;
 
-        if (outlet[b_base] == b_base)
-            return std::make_pair(
-                graph_traits< Graph >::null_vertex(), missing_outlet);
+        /** Delta value. */
+        weight_t value;
 
-        vertex_iterator_t vi, vi_end;
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
-            old_label[*vi] = std::make_pair(label_T[*vi], outlet[*vi]);
-
-        std::pair< vertex_descriptor_t, vertex_state_t > child(
-            outlet[b_base], graph::detail::V_EVEN);
-        blossom_ptr_t b = in_blossom[child.first];
-        for (; b->father->father != blossom_ptr_t(); b = b->father)
-            ;
-        child.first = b->get_base();
+        /** Edge on which the minimum delta occurs (for delta type 2 or 3). */
+        edge_t edge;
 
-        if (child.first == b_base)
-            return std::make_pair(
-                graph_traits< Graph >::null_vertex(), missing_outlet);
+        /** Blossom on which the minimum delta occurs (for delta type 4). */
+        nontrivial_blossom_t* blossom;
+    };
 
-        while (true)
-        {
-            std::pair< vertex_descriptor_t, vertex_state_t > child_parent
-                = parent(child, true);
+    /** Similar to vector_property_map, but uses a fixed-size vector. */
+    template <typename T>
+    struct vertex_map
+    {
+        using key_type = typename property_traits<VertexIndexMap>::key_type;
+        std::vector<T> vec;
+        VertexIndexMap vm;
 
-            for (b = in_blossom[child_parent.first];
-                 b->father->father != blossom_ptr_t(); b = b->father)
-                ;
-            missing_outlet = child_parent.first;
-            child_parent.first = b->get_base();
+        vertex_map(vertices_size_t arg_size, VertexIndexMap arg_vm)
+          : vec(arg_size)
+          , vm(arg_vm)
+        { }
 
-            if (child_parent.first == b_base)
-                break;
-            else
-                child = child_parent;
+        const T& operator[](const key_type& v) const
+        {
+            return vec[get(vm, v)];
         }
-        return std::make_pair(child.first, missing_outlet);
-    }
 
-    // expand a top blossom, put all its non-trivial sub-blossoms into
-    // top_blossoms
-    blossom_iterator_t expand_blossom(
-        blossom_iterator_t bi, std::vector< blossom_ptr_t >& new_ones)
-    {
-        blossom_ptr_t b = *bi;
-        for (blossom_iterator_t i = b->sub_blossoms.begin();
-             i != b->sub_blossoms.end(); ++i)
+        T& operator[](const key_type& v)
         {
-            blossom_ptr_t sub_blossom = *i;
-            vertex_descriptor_t sub_base = sub_blossom->get_base();
-            label_S[sub_base] = label_T[sub_base]
-                = graph_traits< Graph >::null_vertex();
-            outlet[sub_base] = sub_base;
-            sub_blossom->father = blossom_ptr_t();
-            // new top blossoms cannot be pushed back into top_blossoms
-            // immediately, because push_back() may cause reallocation and then
-            // invalid iterators
-            if (!sub_blossom->sub_blossoms.empty())
-                new_ones.push_back(sub_blossom);
+            return vec[get(vm, v)];
         }
-        return top_blossoms.erase(bi);
-    }
+    };
 
-    // when expanding a T-blossom with base v, it requires more operations:
-    // supply the missing T-labels for new base vertices by picking the minimum
-    // tau from vertices of each corresponding new top-blossoms; when label_T[v]
-    // is null or we have a smaller tau from missing_label(v), replace T-label
-    // and outlet of v (but don't bloom v)
-    blossom_iterator_t expand_T_blossom(
-        blossom_iterator_t bi, std::vector< blossom_ptr_t >& new_ones)
+    /** Keep track of the least-slack edge out of a bunch of edges. */
+    struct least_slack_edge_t
     {
-        blossom_ptr_t b = *bi;
+        optional<edge_t> edge;
+        weight_t slack;
 
-        vertex_descriptor_t b_base = b->get_base();
-        std::pair< vertex_descriptor_t, vertex_descriptor_t > T_and_outlet
-            = missing_label(b_base);
+        least_slack_edge_t() : slack(0) {}
 
-        blossom_iterator_t next_bi = expand_blossom(bi, new_ones);
-
-        for (blossom_iterator_t i = b->sub_blossoms.begin();
-             i != b->sub_blossoms.end(); ++i)
+        void update(const edge_t& e, weight_t s)
         {
-            blossom_ptr_t sub_blossom = *i;
-            vertex_descriptor_t sub_base = sub_blossom->get_base();
-            vertex_descriptor_t min_tau_v
-                = graph_traits< Graph >::null_vertex();
-            edge_property_t min_tau
-                = std::numeric_limits< edge_property_t >::max();
-
-            std::vector< vertex_descriptor_t > sub_vertices
-                = sub_blossom->vertices();
-            for (vertex_vec_iter_t v = sub_vertices.begin();
-                 v != sub_vertices.end(); ++v)
+            if ((!edge.has_value()) || (s < slack))
             {
-                if (tau[*v] < min_tau)
-                {
-                    min_tau = tau[*v];
-                    min_tau_v = *v;
-                }
+                edge = e;
+                slack = s;
             }
-
-            if (min_tau < std::numeric_limits< edge_property_t >::max())
-                put_T_label(
-                    sub_base, tau_idx[min_tau_v], min_tau_v, tau[min_tau_v]);
         }
+    };
 
-        if (label_T[b_base] == graph_traits< Graph >::null_vertex()
-            || tau[old_label[b_base].second] < pi[b_base])
-            boost::tie(label_T[b_base], outlet[b_base]) = T_and_outlet;
-
-        return next_bi;
-    }
+    /** Scale integer edge weights to enable integer-only calculations. */
+    static constexpr weight_t weight_factor =
+        std::numeric_limits<weight_t>::is_integer ? 2 : 1;
 
-    // when vertices v and w are matched to each other by augmenting,
-    // we must set v/w as base vertex of any blossom who contains v/w and
-    // is a sub-blossom of their lowest (smallest) common blossom
-    void adjust_blossom(vertex_descriptor_t v, vertex_descriptor_t w)
+    /** Input graph. */
+    const Graph* g;
+    VertexIndexMap vm;
+    EdgeWeightMap edge_weights;
+
+    /** Current mate of each vertex. */
+    vertex_map<vertex_t> vertex_mate;
+
+    /**
+     * For each vertex, the trivial blossom that contains it.
+     *
+     * Pointers to blossoms must remain valid for the life time of
+     * this data structure, therefore the underlying vector must
+     * have a fixed size.
+     */
+    vertex_map<blossom_t> trivial_blossom;
+
+    /**
+     * List of non-trivial blossoms.
+     *
+     * This must be a linked list to ensure that elements can be added
+     * and removed without invalidating pointers to other elements.
+     */
+    std::list<nontrivial_blossom_t> nontrivial_blossom;
+
+    /** For each vertex, the unique top-level blossom that contains it. */
+    vertex_map<blossom_t*> vertex_top_blossom;
+
+    /** For each vertex, a variable in the dual LPP. */
+    vertex_map<weight_t> vertex_dual;
+
+    /** For T-vertex or unlabeled vertex, least-slack edge to any S-vertex. */
+    vertex_map<optional<edge_t>> vertex_best_edge;
+
+    /** Queue of S-vertices to be scanned. */
+    std::deque<vertex_t> scan_queue;
+
+    /** Initialize the matching algorithm. */
+    explicit maximum_weighted_matching_context(
+        const Graph& arg_g, VertexIndexMap arg_vm, EdgeWeightMap weights)
+      : g(&arg_g)
+      , vm(arg_vm)
+      , edge_weights(weights)
+      , vertex_mate(num_vertices(*g), arg_vm)
+      , trivial_blossom(num_vertices(*g), arg_vm)
+      , vertex_top_blossom(num_vertices(*g), arg_vm)
+      , vertex_dual(num_vertices(*g), arg_vm)
+      , vertex_best_edge(num_vertices(*g), arg_vm)
     {
-        blossom_ptr_t vb = in_blossom[v], wb = in_blossom[w],
-                      lowest_common_blossom;
-        std::vector< blossom_ptr_t > v_ancestors, w_ancestors;
+        // Vertex duals are initialized to half the maximum edge weight.
+        weight_t max_weight = 0;
+        for (const edge_t& e : make_iterator_range(edges(*g)))
+            max_weight = (std::max)(max_weight, get(weights, e));
+        weight_t init_vertex_dual = max_weight * (weight_factor / 2);
 
-        while (vb->father != blossom_ptr_t())
+        for (const vertex_t& x : make_iterator_range(vertices(*g)))
         {
-            v_ancestors.push_back(vb->father);
-            vb = vb->father;
-        }
-        while (wb->father != blossom_ptr_t())
-        {
-            w_ancestors.push_back(wb->father);
-            wb = wb->father;
+            vertex_mate[x] = null_vertex();
+            trivial_blossom[x].base_vertex = x;
+            vertex_top_blossom[x] = &trivial_blossom[x];
+            vertex_dual[x] = init_vertex_dual;
         }
+    }
 
-        typename std::vector< blossom_ptr_t >::reverse_iterator i, j;
-        i = v_ancestors.rbegin();
-        j = w_ancestors.rbegin();
-        while (i != v_ancestors.rend() && j != w_ancestors.rend()
-            && i->get() == j->get())
-        {
-            lowest_common_blossom = *i;
-            ++i;
-            ++j;
-        }
+    /** Return the null vertex. */
+    static vertex_t null_vertex()
+    {
+        return graph_traits<Graph>::null_vertex();
+    }
 
-        vb = in_blossom[v];
-        wb = in_blossom[w];
-        while (vb->father != lowest_common_blossom)
-        {
-            vb->father->set_base(vb);
-            vb = vb->father;
-        }
-        while (wb->father != lowest_common_blossom)
-        {
-            wb->father->set_base(wb);
-            wb = wb->father;
+    /** Call a function for every vertex inside the specified blossom. */
+    template <typename Func>
+    static void for_vertices_in_blossom(const blossom_t& blossom, Func func)
+    {
+        auto ntb = blossom.nontrivial();
+        if (ntb) {
+            // Visit all vertices in the non-trivial blossom.
+            // Use an explicit stack to avoid deep call chains.
+            std::vector<const nontrivial_blossom_t*> stack;
+            stack.push_back(ntb);
+            while (!stack.empty()) {
+                auto cur = stack.back();
+                stack.pop_back();
+                for (const auto& sub : cur->subblossoms) {
+                    ntb = sub.blossom->nontrivial();
+                    if (ntb) {
+                        stack.push_back(ntb);
+                    } else {
+                        func(sub.blossom->base_vertex);
+                    }
+                }
+            }
+        } else {
+            // A trivial blossom contains just one vertex.
+            func(blossom.base_vertex);
         }
     }
 
-    // every edge weight is multiplied by 4 to ensure integer weights
-    // throughout the algorithm if all input weights are integers
-    inline edge_property_t slack(const edge_descriptor_t& e) const
+    /** Mark blossom as the top-level blossom of its vertices. */
+    void update_top_level_blossom(blossom_t& blossom)
     {
-        vertex_descriptor_t v, w;
-        v = source(e, g);
-        w = target(e, g);
-        return dual_var[v] + dual_var[w] - 4 * get(edge_weight, g, e);
+        BOOST_ASSERT(!blossom.parent);
+        for_vertices_in_blossom(blossom,
+            [this,&blossom](vertex_t x)
+            {
+                vertex_top_blossom[x] = &blossom;
+            });
     }
 
-    // backtrace one step on vertex v along the augmenting path
-    // by its labels and its vertex state;
-    // boolean parameter "use_old" means whether we are updating labels,
-    // if we are, then we use old labels to backtrace and also we
-    // don't jump to its base vertex when we reach an odd vertex
-    std::pair< vertex_descriptor_t, vertex_state_t > parent(
-        std::pair< vertex_descriptor_t, vertex_state_t > v,
-        bool use_old = false) const
+    /**
+     * Calculate edge slack.
+     * The two vertices must be in different top-level blossoms.
+     */
+    weight_t edge_slack(const edge_t& e) const
     {
-        if (v.second == graph::detail::V_EVEN)
-        {
-            // a paranoid check: label_S shoule be the same as mate in
-            // backtracing
-            if (label_S[v.first] == graph_traits< Graph >::null_vertex())
-                label_S[v.first] = mate[v.first];
-            return std::make_pair(label_S[v.first], graph::detail::V_ODD);
-        }
-        else if (v.second == graph::detail::V_ODD)
-        {
-            vertex_descriptor_t w = use_old ? old_label[v.first].first
-                                            : base_vertex(label_T[v.first]);
-            return std::make_pair(w, graph::detail::V_EVEN);
-        }
-        return std::make_pair(v.first, graph::detail::V_UNREACHED);
+        vertex_t x = source(e, *g);
+        vertex_t y = target(e, *g);
+        weight_t w = get(edge_weights, e);
+        BOOST_ASSERT(vertex_top_blossom[x] != vertex_top_blossom[y]);
+        return vertex_dual[x] + vertex_dual[y] - weight_factor * w;
     }
 
-    // backtrace from vertices v and w to their free (unmatched) ancesters,
-    // return the nearest common ancestor (null_vertex if none) of v and w
-    vertex_descriptor_t nearest_common_ancestor(vertex_descriptor_t v,
-        vertex_descriptor_t w, vertex_descriptor_t& v_free_ancestor,
-        vertex_descriptor_t& w_free_ancestor) const
+    /** Clear least-slack edge tracking. */
+    void clear_best_edge()
     {
-        std::pair< vertex_descriptor_t, vertex_state_t > v_up(
-            v, graph::detail::V_EVEN);
-        std::pair< vertex_descriptor_t, vertex_state_t > w_up(
-            w, graph::detail::V_EVEN);
-        vertex_descriptor_t nca;
-        nca = w_free_ancestor = v_free_ancestor
-            = graph_traits< Graph >::null_vertex();
-
-        std::vector< bool > ancestor_of_w_vector(num_vertices(g), false);
-        std::vector< bool > ancestor_of_v_vector(num_vertices(g), false);
-        vertex_to_bool_map_t ancestor_of_w(ancestor_of_w_vector.begin(), vm);
-        vertex_to_bool_map_t ancestor_of_v(ancestor_of_v_vector.begin(), vm);
-
-        while (nca == graph_traits< Graph >::null_vertex()
-            && (v_free_ancestor == graph_traits< Graph >::null_vertex()
-                || w_free_ancestor == graph_traits< Graph >::null_vertex()))
+        for (vertex_t x : make_iterator_range(vertices(*g)))
         {
-            ancestor_of_w[w_up.first] = true;
-            ancestor_of_v[v_up.first] = true;
-
-            if (w_free_ancestor == graph_traits< Graph >::null_vertex())
-                w_up = parent(w_up);
-            if (v_free_ancestor == graph_traits< Graph >::null_vertex())
-                v_up = parent(v_up);
-
-            if (mate[v_up.first] == graph_traits< Graph >::null_vertex())
-                v_free_ancestor = v_up.first;
-            if (mate[w_up.first] == graph_traits< Graph >::null_vertex())
-                w_free_ancestor = w_up.first;
-
-            if (ancestor_of_w[v_up.first] == true || v_up.first == w_up.first)
-                nca = v_up.first;
-            else if (ancestor_of_v[w_up.first] == true)
-                nca = w_up.first;
-            else if (v_free_ancestor == w_free_ancestor
-                && v_free_ancestor != graph_traits< Graph >::null_vertex())
-                nca = v_up.first;
+            vertex_best_edge[x].reset();
+            trivial_blossom[x].best_edge.reset();
         }
 
-        return nca;
+        for (nontrivial_blossom_t& b : nontrivial_blossom)
+        {
+            b.best_edge.reset();
+            b.best_edge_set.clear();
+        }
     }
 
-    // when a new top blossom b is created by connecting (v, w), we add
-    // sub-blossoms into b along backtracing from v_prime and w_prime to
-    // stop_vertex (the base vertex); also, we set labels and outlet for each
-    // base vertex we pass by
-    void make_blossom(blossom_ptr_t b, vertex_descriptor_t w_prime,
-        vertex_descriptor_t v_prime, vertex_descriptor_t stop_vertex)
+    /** Add edge from unlabeled verter or T-vertex "y" to an S-vertex. */
+    void add_delta2_edge(vertex_t y, const edge_t& e, weight_t slack)
     {
-        std::pair< vertex_descriptor_t, vertex_state_t > u(
-            v_prime, graph::detail::V_ODD);
-        std::pair< vertex_descriptor_t, vertex_state_t > u_up(
-            w_prime, graph::detail::V_EVEN);
+        auto& best_edge = vertex_best_edge[y];
+        if ((!best_edge.has_value()) || slack < edge_slack(*best_edge))
+            best_edge = e;
+    }
 
-        for (; u_up.first != stop_vertex; u = u_up, u_up = parent(u))
+    /** Return least-slack edge between any unlabeled vertex and S-vertex. */
+    least_slack_edge_t get_least_slack_delta2_edge()
+    {
+        least_slack_edge_t best_edge;
+        for (vertex_t x : make_iterator_range(vertices(*g)))
         {
-            if (u_up.second == graph::detail::V_EVEN)
+            if (vertex_top_blossom[x]->label == LABEL_NONE)
             {
-                if (!in_top_blossom(u_up.first)->sub_blossoms.empty())
-                    outlet[u_up.first] = label_T[u.first];
-                label_T[u_up.first] = outlet[u.first];
+                const auto& vertex_edge = vertex_best_edge[x];
+                if (vertex_edge.has_value())
+                    best_edge.update(*vertex_edge, edge_slack(*vertex_edge));
             }
-            else if (u_up.second == graph::detail::V_ODD)
-                label_S[u_up.first] = u.first;
-
-            add_sub_blossom(b, u_up.first);
         }
+        return best_edge;
     }
 
-    // the design of recursively expanding augmenting path in
-    // (reversed_)retrieve_augmenting_path functions is inspired by same
-    // functions in max_cardinality_matching.hpp; except that in weighted
-    // matching, we use "outlet" vertices instead of "bridge" vertex pairs: if
-    // blossom b is the smallest non-trivial blossom that contains its base
-    // vertex v, then v and outlet[v] are where augmenting path enters and
-    // leaves b
-    void retrieve_augmenting_path(
-        vertex_descriptor_t v, vertex_descriptor_t w, vertex_state_t v_state)
+    /** Add edge between different S-blossoms. */
+    void add_delta3_edge(blossom_t& b, const edge_t& e, weight_t slack)
     {
-        if (v == w)
-            aug_path.push_back(v);
-        else if (v_state == graph::detail::V_EVEN)
-        {
-            aug_path.push_back(v);
-            retrieve_augmenting_path(label_S[v], w, graph::detail::V_ODD);
-        }
-        else if (v_state == graph::detail::V_ODD)
-        {
-            if (outlet[v] == v)
-                aug_path.push_back(v);
-            else
-                reversed_retrieve_augmenting_path(
-                    outlet[v], v, graph::detail::V_EVEN);
-            retrieve_augmenting_path(label_T[v], w, graph::detail::V_EVEN);
-        }
-    }
+        auto& best_edge = b.best_edge;
+        if ((!best_edge.has_value()) || slack < edge_slack(*best_edge))
+            best_edge = e;
 
-    void reversed_retrieve_augmenting_path(
-        vertex_descriptor_t v, vertex_descriptor_t w, vertex_state_t v_state)
-    {
-        if (v == w)
-            aug_path.push_back(v);
-        else if (v_state == graph::detail::V_EVEN)
-        {
-            reversed_retrieve_augmenting_path(
-                label_S[v], w, graph::detail::V_ODD);
-            aug_path.push_back(v);
-        }
-        else if (v_state == graph::detail::V_ODD)
-        {
-            reversed_retrieve_augmenting_path(
-                label_T[v], w, graph::detail::V_EVEN);
-            if (outlet[v] != v)
-                retrieve_augmenting_path(outlet[v], v, graph::detail::V_EVEN);
-            else
-                aug_path.push_back(v);
-        }
+        auto ntb = b.nontrivial();
+        if (ntb)
+            ntb->best_edge_set.push_back(e);
     }
 
-    // correct labels for vertices in the augmenting path
-    void relabel(vertex_descriptor_t v)
+    /** Update least-slack edge tracking after merging blossoms. */
+    void merge_delta3_blossoms(nontrivial_blossom_t& blossom)
     {
-        blossom_ptr_t b = in_blossom[v]->father;
-
-        if (!is_in_blossom(b, mate[v]))
-        { // if v is a new base vertex
-            std::pair< vertex_descriptor_t, vertex_state_t > u(
-                v, graph::detail::V_EVEN);
-            while (label_S[u.first] != u.first
-                && is_in_blossom(b, label_S[u.first]))
-                u = parent(u, true);
-
-            vertex_descriptor_t old_base = u.first;
-            if (label_S[old_base] != old_base)
-            { // if old base is not exposed
-                label_T[v] = label_S[old_base];
-                outlet[v] = old_base;
-            }
-            else
-            { // if old base is exposed then new label_T[v] is not in b,
-                // we must (i) make b2 the smallest blossom containing v but not
-                // as base vertex (ii) backtrace from b2's new base vertex to b
-                label_T[v] = graph_traits< Graph >::null_vertex();
-                for (b = b->father; b != blossom_ptr_t() && b->get_base() == v;
-                     b = b->father)
-                    ;
-                if (b != blossom_ptr_t())
+        // Build a temporary array holding the least-slack edges to
+        // other S-blossoms. The array is indexed by base vertex.
+        std::vector<least_slack_edge_t> tmp_best_edge(num_vertices(*g));
+
+        // Collect edges from sub-blossoms that were S-blossoms.
+        for (auto& sub : blossom.subblossoms)
+        {
+            blossom_t* b = sub.blossom;
+            if (b->label == LABEL_S)
+            {
+                b->best_edge.reset();
+                auto ntb = b->nontrivial();
+                if (ntb)
+                {
+                    // Use least-slack edges from subblossom.
+                    for (const edge_t& e : ntb->best_edge_set)
+                    {
+                        blossom_t* bx = vertex_top_blossom[source(e, *g)];
+                        blossom_t* by = vertex_top_blossom[target(e, *g)];
+                        BOOST_ASSERT(bx == &blossom);
+                        // Only use edges between top-level blossoms.
+                        if (bx != by)
+                        {
+                            BOOST_ASSERT(by->label == LABEL_S);
+                            tmp_best_edge[get(vm, by->base_vertex)].update(
+                                e, edge_slack(e));
+                        }
+                    }
+                    ntb->best_edge_set.clear();
+                }
+                else
                 {
-                    u = std::make_pair(b->get_base(), graph::detail::V_ODD);
-                    while (!is_in_blossom(
-                        in_blossom[v]->father, old_label[u.first].first))
-                        u = parent(u, true);
-                    label_T[v] = u.first;
-                    outlet[v] = old_label[u.first].first;
+                    // Trivial blossoms don't maintain a least-slack edge set.
+                    // Consider all incident edges.
+                    for (const edge_t& e :
+                        make_iterator_range(out_edges(b->base_vertex, *g)))
+                    {
+                        BOOST_ASSERT(source(e, *g) == b->base_vertex);
+                        vertex_t y = target(e, *g);
+                        blossom_t* by = vertex_top_blossom[y];
+                        // Only use edges to different S-blossoms.
+                        // Ignore edges with negative weight.
+                        if ((by != &blossom)
+                            && (by->label == LABEL_S)
+                            && (get(edge_weights, e) >= 0))
+                        {
+                            tmp_best_edge[get(vm, by->base_vertex)].update(
+                                e, edge_slack(e));
+                        }
+                    }
                 }
             }
         }
-        else if (label_S[v] == v || !is_in_blossom(b, label_S[v]))
-        { // if v is an old base vertex
-            // let u be the new base vertex; backtrace from u's old T-label
-            std::pair< vertex_descriptor_t, vertex_state_t > u(
-                b->get_base(), graph::detail::V_ODD);
-            while (
-                old_label[u.first].first != graph_traits< Graph >::null_vertex()
-                && old_label[u.first].first != v)
-                u = parent(u, true);
-            label_T[v] = old_label[u.first].second;
-            outlet[v] = v;
+
+        // Build a compact list of edges from the temporary array.
+        // Also find the overall least-slack edge to any other S-blossom.
+        BOOST_ASSERT(blossom.best_edge_set.empty());
+        BOOST_ASSERT(!blossom.best_edge.has_value());
+        least_slack_edge_t best_edge;
+        for (const least_slack_edge_t& item : tmp_best_edge)
+        {
+            if (item.edge.has_value())
+            {
+                blossom.best_edge_set.push_back(*item.edge);
+                best_edge.update(*item.edge, item.slack);
+            }
         }
-        else // if v is neither a new nor an old base vertex
-            label_T[v] = label_S[v];
+        blossom.best_edge = best_edge.edge;
     }
 
-    void augmenting(vertex_descriptor_t v, vertex_descriptor_t v_free_ancestor,
-        vertex_descriptor_t w, vertex_descriptor_t w_free_ancestor)
+    /** Return least-slack edge between any pair of S-blossoms. */
+    least_slack_edge_t get_least_slack_delta3_edge()
     {
-        vertex_iterator_t vi, vi_end;
-
-        // retrieve the augmenting path and put it in aug_path
-        reversed_retrieve_augmenting_path(
-            v, v_free_ancestor, graph::detail::V_EVEN);
-        retrieve_augmenting_path(w, w_free_ancestor, graph::detail::V_EVEN);
+        least_slack_edge_t best_edge;
 
-        // augment the matching along aug_path
-        vertex_descriptor_t a, b;
-        vertex_list_t reversed_aug_path;
-        while (!aug_path.empty())
+        for (vertex_t x : make_iterator_range(vertices(*g)))
         {
-            a = aug_path.front();
-            aug_path.pop_front();
-            reversed_aug_path.push_back(a);
-            b = aug_path.front();
-            aug_path.pop_front();
-            reversed_aug_path.push_back(b);
-
-            mate[a] = b;
-            mate[b] = a;
-
-            // reset base vertex for every blossom in augment path
-            adjust_blossom(a, b);
+            blossom_t* b = vertex_top_blossom[x];
+            BOOST_ASSERT(!b->parent);
+            if ((b->label == LABEL_S) && (b->best_edge.has_value()))
+                best_edge.update(*b->best_edge, edge_slack(*b->best_edge));
         }
 
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
-            old_label[*vi] = std::make_pair(label_T[*vi], outlet[*vi]);
+        return best_edge;
+    }
+
+    /** Add the vertices in a blossom to the scan queue. */
+    void add_vertices_to_scan_queue(blossom_t& blossom)
+    {
+        for_vertices_in_blossom(blossom,
+            [this](vertex_t x)
+            {
+                scan_queue.push_back(x);
+            });
+    }
 
-        // correct labels for in-blossom vertices along aug_path
-        while (!reversed_aug_path.empty())
-        {
-            a = reversed_aug_path.front();
-            reversed_aug_path.pop_front();
+    /** Trace back through the alternating trees from vertices "x" and "y". */
+    alternating_path_t trace_alternating_paths(vertex_t x, vertex_t y)
+    {
+        // Initialize a path containing only the edge (x, y).
+        alternating_path_t path;
+        path.emplace_back(x, y);
 
-            if (in_blossom[a]->father != blossom_ptr_t())
-                relabel(a);
+        // Trace alternating path from "x" to the root of the tree.
+        while (x != null_vertex())
+        {
+            blossom_t* bx = vertex_top_blossom[x];
+            x = bx->tree_edge.first;
+            if (x != null_vertex())
+                path.push_front(bx->tree_edge);
         }
 
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
+        // Trace alternating path from "y" to the root of the tree.
+        while (y != null_vertex())
         {
-            vertex_descriptor_t u = *vi;
-            if (mate[u] != graph_traits< Graph >::null_vertex())
-                label_S[u] = mate[u];
+            blossom_t* by = vertex_top_blossom[y];
+            y = by->tree_edge.first;
+            if (y != null_vertex())
+                path.emplace_back(by->tree_edge.second, y);
         }
 
-        // expand blossoms with zero dual variables
-        std::vector< blossom_ptr_t > new_top_blossoms;
-        for (blossom_iterator_t bi = top_blossoms.begin();
-             bi != top_blossoms.end();)
+        // If we found a common ancestor, trim the paths so they end there.
+        while (path.front().second == path.back().first)
         {
-            if ((*bi)->dual_var <= 0)
-                bi = expand_blossom(bi, new_top_blossoms);
-            else
-                ++bi;
+            BOOST_ASSERT(path.size() > 2);
+            path.pop_front();
+            path.pop_back();
         }
-        top_blossoms.insert(top_blossoms.end(), new_top_blossoms.begin(),
-            new_top_blossoms.end());
-        init();
+
+        // Any alternating path between S-blossoms must have odd length.
+        BOOST_ASSERT(path.size() % 2 == 1);
+
+        return path;
     }
 
-    // create a new blossom and set labels for vertices inside
-    void blossoming(vertex_descriptor_t v, vertex_descriptor_t v_prime,
-        vertex_descriptor_t w, vertex_descriptor_t w_prime,
-        vertex_descriptor_t nca)
+    /** Create a new S-blossom from an alternating cycle. */
+    void make_blossom(const alternating_path_t& path)
     {
-        vertex_iterator_t vi, vi_end;
-
-        std::vector< bool > is_old_base_vector(num_vertices(g));
-        vertex_to_bool_map_t is_old_base(is_old_base_vector.begin(), vm);
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
+        BOOST_ASSERT(path.size() % 2 == 1);
+        BOOST_ASSERT(path.size() >= 3);
+
+        // Collect pointers to sub-blossoms.
+        std::vector<blossom_t*> subblossoms;
+        subblossoms.reserve(path.size());
+        for (const vertex_pair_t& edge : path)
+            subblossoms.push_back(vertex_top_blossom[edge.first]);
+
+        // Check that the path is cyclic.
+        vertices_size_t pos = 0;
+        for (const vertex_pair_t& edge : path)
         {
-            if (*vi == base_vertex(*vi))
-                is_old_base[*vi] = true;
+            pos = (pos + 1) % subblossoms.size();
+            BOOST_ASSERT(vertex_top_blossom[edge.second] == subblossoms[pos]);
         }
 
-        blossom_ptr_t b = boost::make_shared< blossom >();
-        add_sub_blossom(b, nca);
+        // Create the new blossom.
+        nontrivial_blossom.emplace_back(subblossoms, path);
+        nontrivial_blossom_t& blossom = nontrivial_blossom.back();
 
-        label_T[w_prime] = v;
-        label_T[v_prime] = w;
-        outlet[w_prime] = w;
-        outlet[v_prime] = v;
+        // Link the subblossoms to the new parent.
+        // Insert former T-vertices into the scan queue.
+        for (blossom_t* b : subblossoms)
+        {
+            BOOST_ASSERT(!b->parent);
+            b->parent = &blossom;
+            if (b->label == LABEL_T)
+                add_vertices_to_scan_queue(*b);
+        }
 
-        make_blossom(b, w_prime, v_prime, nca);
-        make_blossom(b, v_prime, w_prime, nca);
+        // Mark vertices as belonging to the new blossom.
+        update_top_level_blossom(blossom);
 
-        label_T[nca] = graph_traits< Graph >::null_vertex();
-        outlet[nca] = nca;
+        // Assign label S to the new blossom and link to the alternating tree.
+        BOOST_ASSERT(subblossoms.front()->label == LABEL_S);
+        blossom.label = LABEL_S;
+        blossom.tree_edge = subblossoms.front()->tree_edge;
+        blossom.tree_root = subblossoms.front()->tree_root;
 
-        top_blossoms.push_back(b);
-        bloom(b);
+        // Merge least-slack edges for the S-sub-blossoms.
+        merge_delta3_blossoms(blossom);
+    }
 
-        // set gamma[b_base] = min_slack{critical_edge(b_base, other_base)}
-        // where each critical edge is updated before, by
-        // argmin{slack(old_bases_in_b, other_base)};
-        vertex_vec_iter_t i, j;
-        std::vector< vertex_descriptor_t > b_vertices = b->vertices(),
-                                           old_base_in_b, other_base;
-        vertex_descriptor_t b_base = b->get_base();
-        for (i = b_vertices.begin(); i != b_vertices.end(); ++i)
-        {
-            if (is_old_base[*i])
-                old_base_in_b.push_back(*i);
-        }
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
+    /** Expand and delete a T-blossom. */
+    void expand_t_blossom(nontrivial_blossom_t* blossom)
+    {
+        BOOST_ASSERT(!blossom->parent);
+        BOOST_ASSERT(blossom->label == LABEL_T);
+
+        // Convert sub-blossoms into top-level blossoms.
+        for (const auto& sub : blossom->subblossoms)
         {
-            if (*vi != b_base && *vi == base_vertex(*vi))
-                other_base.push_back(*vi);
+            blossom_t* b = sub.blossom;
+            BOOST_ASSERT(b->parent == blossom);
+            b->parent = nullptr;
+            b->label = LABEL_NONE;
+            update_top_level_blossom(*b);
         }
-        for (i = other_base.begin(); i != other_base.end(); ++i)
+
+        // Reconstruct the alternating tree via the sub-blossoms.
+        // Find the sub-blossom that attaches the expanding blossom to
+        // the alternating tree.
+        blossom_t* entry = vertex_top_blossom[blossom->tree_edge.second];
+
+        auto subblossom_loc = blossom->find_subblossom(entry);
+        auto entry_pos = subblossom_loc.first;
+        auto entry_it = subblossom_loc.second;
+
+        // Get the edge that attached this blossom to the alternating tree.
+        vertex_t x, y;
+        std::tie(x, y) = blossom->tree_edge;
+
+        // Reconstruct the tree in an even number of steps from entry to base.
+        auto sub_it = entry_it;
+        if (entry_pos % 2 == 0)
         {
-            edge_property_t min_slack
-                = std::numeric_limits< edge_property_t >::max();
-            std::pair< edge_descriptor_t, bool > b_vi = null_edge;
-            for (j = old_base_in_b.begin(); j != old_base_in_b.end(); ++j)
+            // Walk backward around the blossom.
+            auto sub_begin = blossom->subblossoms.begin();
+            while (sub_it != sub_begin)
             {
-                if (critical_edge[*j][*i] != null_edge
-                    && min_slack > slack(critical_edge[*j][*i].first))
-                {
-                    min_slack = slack(critical_edge[*j][*i].first);
-                    b_vi = critical_edge[*j][*i];
-                }
+                extend_tree_s_to_t(x, y);
+                --sub_it;
+                BOOST_ASSERT(sub_it != sub_begin);
+                --sub_it;
+                std::tie(y, x) = sub_it->edge;
             }
-            critical_edge[b_base][*i] = critical_edge[*i][b_base] = b_vi;
         }
-        gamma[b_base] = std::numeric_limits< edge_property_t >::max();
-        for (i = other_base.begin(); i != other_base.end(); ++i)
+        else
         {
-            if (critical_edge[b_base][*i] != null_edge)
-                gamma[b_base] = std::min(
-                    gamma[b_base], slack(critical_edge[b_base][*i].first));
+            // Walk forward around the blossom.
+            auto sub_end = blossom->subblossoms.end();
+            while (sub_it != sub_end) {
+                extend_tree_s_to_t(x, y);
+                ++sub_it;
+                BOOST_ASSERT(sub_it != sub_end);
+                std::tie(x, y) = sub_it->edge;
+                ++sub_it;
+            }
         }
+
+        // Assign label T to the base sub-blossom.
+        blossom_t* base = blossom->subblossoms.front().blossom;
+        base->label = LABEL_T;
+        base->tree_edge = std::make_pair(x, y);
+        base->tree_root = blossom->tree_root;
+
+        // Delete the expanded blossom.
+        auto blossom_it = std::find_if(
+            nontrivial_blossom.begin(),
+            nontrivial_blossom.end(),
+            [blossom](const nontrivial_blossom_t& b)
+            {
+                return (&b == blossom);
+            });
+        BOOST_ASSERT(blossom_it != nontrivial_blossom.end());
+        nontrivial_blossom.erase(blossom_it);
     }
 
-    void init()
+    void augment_blossom_rec(
+        nontrivial_blossom_t& blossom, blossom_t& entry,
+        std::stack<std::pair<nontrivial_blossom_t*, blossom_t*>>& stack)
     {
-        even_edges.clear();
+        auto subblossom_loc = blossom.find_subblossom(&entry);
+        auto entry_pos = subblossom_loc.first;
+        auto entry_it = subblossom_loc.second;
+
+        // Walk from entry to the base in an even number of steps.
+        auto sub_begin = blossom.subblossoms.begin();
+        auto sub_end = blossom.subblossoms.end();
+        auto sub_it = entry_it;
+        while ((sub_it != sub_begin) && (sub_it != sub_end)) {
+            vertex_t x, y;
+            blossom_t* bx;
+            blossom_t* by;
+
+            if (entry_pos % 2 == 0)
+            {
+                // Walk backward around the blossom.
+                --sub_it;
+                by = sub_it->blossom;
+                BOOST_ASSERT(sub_it != sub_begin);
+                --sub_it;
+                bx = sub_it->blossom;
+                std::tie(x, y) = sub_it->edge;
+            }
+            else
+            {
+                // Walk forward around the blossom.
+                ++sub_it;
+                BOOST_ASSERT(sub_it != sub_end);
+                std::tie(x, y) = sub_it->edge;
+                bx = sub_it->blossom;
+                ++sub_it;
+                by = (sub_it == sub_end) ?
+                    blossom.subblossoms.front().blossom :
+                    sub_it->blossom;
+            }
 
-        vertex_iterator_t vi, vi_end;
-        typename std::vector<
-            std::vector< std::pair< edge_descriptor_t, bool > > >::iterator vei;
+            // Pull this edge into the matching.
+            vertex_mate[x] = y;
+            vertex_mate[y] = x;
 
-        for (boost::tie(vi, vi_end) = vertices(g),
-                            vei = critical_edge_vectors.begin();
-             vi != vi_end; ++vi, ++vei)
-        {
-            vertex_descriptor_t u = *vi;
-            out_edge_iterator_t ei, ei_end;
+            // Augment through any non-trivial subblossoms touching this edge.
+            auto bx_ntb = bx->nontrivial();
+            if (bx_ntb)
+                stack.emplace(bx_ntb, &trivial_blossom[x]);
+
+            auto by_ntb = by->nontrivial();
+            if (by_ntb)
+                stack.emplace(by_ntb, &trivial_blossom[y]);
+        }
 
-            gamma[u] = tau[u] = pi[u]
-                = std::numeric_limits< edge_property_t >::max();
-            std::fill(vei->begin(), vei->end(), null_edge);
+        // Re-orient the blossom.
+        if (entry_it != sub_begin)
+            blossom.subblossoms.splice(
+                sub_begin, blossom.subblossoms, entry_it, sub_end);
+        blossom.base_vertex = entry.base_vertex;
+    }
 
-            if (base_vertex(u) != u)
-                continue;
+    void augment_blossom(nontrivial_blossom_t& blossom, blossom_t& entry)
+    {
+        // Use an explicit stack to avoid deep call chains.
+        std::stack<std::pair<nontrivial_blossom_t*, blossom_t*>> stack;
+        stack.emplace(&blossom, &entry);
 
-            label_S[u] = label_T[u] = graph_traits< Graph >::null_vertex();
-            outlet[u] = u;
+        while (!stack.empty()) {
+            nontrivial_blossom_t* outer_blossom;
+            blossom_t* inner_entry;
+            std::tie(outer_blossom, inner_entry) = stack.top();
 
-            if (mate[u] == graph_traits< Graph >::null_vertex())
-            {
-                label_S[u] = u;
-                bloom(in_top_blossom(u));
-            }
+            nontrivial_blossom_t* inner_blossom = inner_entry->parent;
+            BOOST_ASSERT(inner_blossom);
+
+            if (inner_blossom != outer_blossom)
+                stack.top() = std::make_pair(outer_blossom, inner_blossom);
+            else
+                stack.pop();
+
+            augment_blossom_rec(*inner_blossom, *inner_entry, stack);
         }
     }
 
-    bool augment_matching()
+    /** Augment the matching through the specified augmenting path. */
+    void augment_matching(const alternating_path_t& path)
     {
-        vertex_descriptor_t v, w, w_free_ancestor, v_free_ancestor;
-        v = w = w_free_ancestor = v_free_ancestor
-            = graph_traits< Graph >::null_vertex();
-        bool found_alternating_path = false;
+        BOOST_ASSERT(path.size() % 2 == 1);
 
-        // note that we only use edges of zero slack value for augmenting
-        while (!even_edges.empty() && !found_alternating_path)
+        // Process the unmatched edges on the augmenting path.
+        auto path_it = path.begin();
+        auto path_end = path.end();
+        while (path_it != path_end)
         {
-            // search for augmenting paths depth-first
-            edge_descriptor_t current_edge = even_edges.back();
-            even_edges.pop_back();
-
-            v = source(current_edge, g);
-            w = target(current_edge, g);
-
-            vertex_descriptor_t v_prime = base_vertex(v);
-            vertex_descriptor_t w_prime = base_vertex(w);
+            vertex_t x = path_it->first;
+            vertex_t y = path_it->second;
+
+            // Augment any non-trivial blossoms that touch this edge.
+            blossom_t* bx = vertex_top_blossom[x];
+            auto bx_ntb = bx->nontrivial();
+            if (bx_ntb)
+                augment_blossom(*bx_ntb, trivial_blossom[x]);
+
+            blossom_t* by = vertex_top_blossom[y];
+            auto by_ntb = by->nontrivial();
+            if (by_ntb)
+                augment_blossom(*by_ntb, trivial_blossom[y]);
+
+            // Pull this edge into the matching.
+            vertex_mate[x] = y;
+            vertex_mate[y] = x;
+
+            // Go to the next unmatched edge on the path.
+            ++path_it;
+            if (path_it == path_end)
+                break;
+            ++path_it;
+        }
+    }
 
-            // w_prime == v_prime implies that we get an edge that has been
-            // shrunk into a blossom
-            if (v_prime == w_prime)
-                continue;
+    /**
+     * Remove non-S-vertices from the scan queue.
+     * This is necessary after removing labels from S-blossoms.
+     */
+    void refresh_scan_queue()
+    {
+        std::deque<vertex_t> new_scan_queue;
+        for (const vertex_t& x : scan_queue)
+        {
+            if (vertex_top_blossom[x]->label == LABEL_S)
+                new_scan_queue.push_back(x);
+        }
+        scan_queue = std::move(new_scan_queue);
+    }
 
-            // a paranoid check
-            if (label_S[v_prime] == graph_traits< Graph >::null_vertex())
-            {
-                std::swap(v_prime, w_prime);
-                std::swap(v, w);
-            }
+    /** Remove edges to non-S-vertices from delta3 edge tracking. */
+    void refresh_delta3_blossom(blossom_t& b)
+    {
+        // Do nothing if this blossom was not tracking any delta3 edge.
+        if (!b.best_edge.has_value())
+            return;
 
-            // w_prime may be unlabeled or have a T-label; replace the existed
-            // T-label if the edge slack is smaller than current pi[w_prime] and
-            // update it. Note that a T-label is "deserved" only when pi equals
-            // zero. also update tau and tau_idx so that tau_idx becomes T-label
-            // when a T-blossom is expanded
-            if (label_S[w_prime] == graph_traits< Graph >::null_vertex())
+        auto ntb = b.nontrivial();
+        if (ntb)
+        {
+            // Remove bad edges from best_edge_set and refresh best_edge.
+            least_slack_edge_t best_edge;
+            auto it = ntb->best_edge_set.begin();
+            auto it_end = ntb->best_edge_set.end();
+            while (it != it_end)
             {
-                if (slack(current_edge) < pi[w_prime])
-                    put_T_label(w_prime, v, w, slack(current_edge));
-                if (slack(current_edge) < tau[w])
+                vertex_t y = target(*it, *g);
+                blossom_t* by = vertex_top_blossom[y];
+                BOOST_ASSERT(by != &b);
+                if (by->label == LABEL_S)
                 {
-                    if (in_blossom[w]->father == blossom_ptr_t()
-                        || label_T[w_prime] == v
-                        || label_T[w_prime]
-                            == graph_traits< Graph >::null_vertex()
-                        || nearest_common_ancestor(v_prime, label_T[w_prime],
-                               v_free_ancestor, w_free_ancestor)
-                            == graph_traits< Graph >::null_vertex())
-                    {
-                        tau[w] = slack(current_edge);
-                        tau_idx[w] = v;
-                    }
+                    best_edge.update(*it, edge_slack(*it));
+                    ++it;
                 }
-            }
-
-            else
-            {
-                if (slack(current_edge) > 0)
+                else
                 {
-                    // update gamma and critical_edges when we have a smaller
-                    // edge slack
-                    gamma[v_prime]
-                        = std::min(gamma[v_prime], slack(current_edge));
-                    gamma[w_prime]
-                        = std::min(gamma[w_prime], slack(current_edge));
-                    if (critical_edge[v_prime][w_prime] == null_edge
-                        || slack(critical_edge[v_prime][w_prime].first)
-                            > slack(current_edge))
-                    {
-                        critical_edge[v_prime][w_prime]
-                            = std::pair< edge_descriptor_t, bool >(
-                                current_edge, true);
-                        critical_edge[w_prime][v_prime]
-                            = std::pair< edge_descriptor_t, bool >(
-                                current_edge, true);
-                    }
-                    continue;
+                    // Remove edge to non-S-blossom.
+                    it = ntb->best_edge_set.erase(it);
                 }
-                else if (slack(current_edge) == 0)
+            }
+            b.best_edge = best_edge.edge;
+        }
+        else
+        {
+            // Trivial blossom does not maintain best_edge_set.
+            // If its current best_edge is invalid, recompute it.
+            vertex_t x = b.base_vertex;
+            vertex_t y = target(*b.best_edge, *g);
+            blossom_t* by = vertex_top_blossom[y];
+            BOOST_ASSERT(by != &b);
+            if (by->label != LABEL_S)
+            {
+                // Consider all incident edges to recompute best_edge.
+                least_slack_edge_t best_edge;
+                for (const edge_t& e : make_iterator_range(out_edges(x, *g)))
                 {
-                    // if nca is null_vertex then we have an augmenting path;
-                    // otherwise we have a new top blossom with nca as its base
-                    // vertex
-                    vertex_descriptor_t nca = nearest_common_ancestor(
-                        v_prime, w_prime, v_free_ancestor, w_free_ancestor);
-
-                    if (nca == graph_traits< Graph >::null_vertex())
-                        found_alternating_path
-                            = true; // to break out of the loop
-                    else
-                        blossoming(v, v_prime, w, w_prime, nca);
+                    BOOST_ASSERT(source(e, *g) == x);
+                    y = target(e, *g);
+                    by = vertex_top_blossom[y];
+                    if ((by != &b) && (by->label == LABEL_S))
+                        best_edge.update(e, edge_slack(e));
                 }
+                b.best_edge = best_edge.edge;
             }
         }
-
-        if (!found_alternating_path)
-            return false;
-
-        augmenting(v, v_free_ancestor, w, w_free_ancestor);
-        return true;
     }
 
-    // slack the vertex and blossom dual variables when there is no augmenting
-    // path found according to the primal-dual method
-    bool adjust_dual()
+    /** Recompute vertex_best_edge for the specified vertex. */
+    void recompute_vertex_best_edge(vertex_t x)
     {
-        edge_property_t delta1, delta2, delta3, delta4, delta;
-        delta1 = delta2 = delta3 = delta4
-            = std::numeric_limits< edge_property_t >::max();
-
-        vertex_iterator_t vi, vi_end;
+        least_slack_edge_t best_edge;
 
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
+        for (const edge_t& e : make_iterator_range(out_edges(x, *g)))
         {
-            delta1 = std::min(delta1, dual_var[*vi]);
-            delta4 = pi[*vi] > 0 ? std::min(delta4, pi[*vi]) : delta4;
-            if (*vi == base_vertex(*vi))
-                delta3 = std::min(delta3, gamma[*vi] / 2);
+            BOOST_ASSERT(source(e, *g) == x);
+            vertex_t y = target(e, *g);
+            blossom_t* by = vertex_top_blossom[y];
+            if (by->label == LABEL_S)
+                best_edge.update(e, edge_slack(e));
         }
+        vertex_best_edge[x] = best_edge.edge;
+    }
 
-        for (blossom_iterator_t bi = top_blossoms.begin();
-             bi != top_blossoms.end(); ++bi)
+    /** Remove the alternating trees with specified root vertices. */
+    void remove_alternating_tree(vertex_t r1, vertex_t r2)
+    {
+        // Find blossoms that are part of the specified alternating trees.
+        std::vector<blossom_t*> former_s_blossoms;
+        for (vertex_t x : make_iterator_range(vertices(*g)))
         {
-            vertex_descriptor_t b_base = (*bi)->get_base();
-            if (label_T[b_base] != graph_traits< Graph >::null_vertex()
-                && pi[b_base] == 0)
-                delta2 = std::min(delta2, (*bi)->dual_var / 2);
-        }
+            blossom_t* b = vertex_top_blossom[x];
+            if ((!b->parent)
+                && (b->label != LABEL_NONE)
+                && (b->tree_root == r1 || b->tree_root == r2))
+            {
+                // Build list of former S-blossoms.
+                if (b->label == LABEL_S)
+                    former_s_blossoms.push_back(b);
 
-        delta = std::min(std::min(delta1, delta2), std::min(delta3, delta4));
+                // Remove label.
+                b->label = LABEL_NONE;
+            }
+        }
 
-        // start updating dual variables, note that the order is important
+        vertex_map<char> blossom_recompute_best_edge(num_vertices(*g), vm);
+        vertex_map<char> vertex_recompute_best_edge(num_vertices(*g), vm);
 
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
+        // Visit former S-blossoms.
+        for (blossom_t* b : former_s_blossoms)
         {
-            vertex_descriptor_t v = *vi, v_prime = base_vertex(v);
-
-            if (label_S[v_prime] != graph_traits< Graph >::null_vertex())
-                dual_var[v] -= delta;
-            else if (label_T[v_prime] != graph_traits< Graph >::null_vertex()
-                && pi[v_prime] == 0)
-                dual_var[v] += delta;
+            // Clear best_edge tracking.
+            b->best_edge.reset();
+            auto ntb = b->nontrivial();
+            if (ntb)
+                ntb->best_edge_set.clear();
+
+            // Visit edges that touch this blossom.
+            for_vertices_in_blossom(*b,
+                [&](vertex_t x)
+                {
+                    // Mark former S-vertices.
+                    vertex_recompute_best_edge[x] = 1;
 
-            if (v == v_prime)
-                gamma[v] -= 2 * delta;
+                    for (const edge_t& e :
+                         make_iterator_range(out_edges(x, *g)))
+                    {
+                        // Mark S-blossoms adjacent to the former S-blossom.
+                        BOOST_ASSERT(source(e, *g) == x);
+                        vertex_t y = target(e, *g);
+                        blossom_t* by = vertex_top_blossom[y];
+                        if (by->label == LABEL_S)
+                            blossom_recompute_best_edge[by->base_vertex] = 1;
+
+                        // Mark non-S-vertices with least-slack edge to
+                        // former S-blossom.
+                        if (by->label != LABEL_S)
+                        {
+                            const auto& best_edge = vertex_best_edge[y];
+                            if (best_edge.has_value() && (*best_edge == e))
+                                vertex_recompute_best_edge[y] = 1;
+                        }
+                    }
+                });
         }
 
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
+        // Recompute delta3 tracking of affected S-blossoms.
+        for (vertex_t x : make_iterator_range(vertices(*g)))
         {
-            vertex_descriptor_t v_prime = base_vertex(*vi);
-            if (pi[v_prime] > 0)
-                tau[*vi] -= delta;
+            if (blossom_recompute_best_edge[x])
+                refresh_delta3_blossom(*vertex_top_blossom[x]);
         }
 
-        for (blossom_iterator_t bi = top_blossoms.begin();
-             bi != top_blossoms.end(); ++bi)
+        // Recompute vertex_best_edge of affected non-S-vertices.
+        for (vertex_t x : make_iterator_range(vertices(*g)))
         {
-            vertex_descriptor_t b_base = (*bi)->get_base();
-            if (label_T[b_base] != graph_traits< Graph >::null_vertex()
-                && pi[b_base] == 0)
-                (*bi)->dual_var -= 2 * delta;
-            if (label_S[b_base] != graph_traits< Graph >::null_vertex())
-                (*bi)->dual_var += 2 * delta;
+            if (vertex_recompute_best_edge[x])
+                recompute_vertex_best_edge(x);
         }
 
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
-        {
-            vertex_descriptor_t v = *vi;
-            if (pi[v] > 0)
-                pi[v] -= delta;
-
-            // when some T-vertices have zero pi value, bloom their mates so
-            // that matching can be further augmented
-            if (label_T[v] != graph_traits< Graph >::null_vertex()
-                && pi[v] == 0)
-                put_T_label(v, label_T[v], outlet[v], pi[v]);
-        }
+        refresh_scan_queue();
+    }
 
-        // optimal solution reached, halt
-        if (delta == delta1)
-            return false;
+    /**
+     * Extend the alternating tree via the edge from S-vertex "x"
+     * to unlabeled vertex "y".
+     *
+     * Assign label T to the blossom that contains "y", then assign
+     * label S to the blossom matched to the newly labeled T-blossom.
+     */
+    void extend_tree_s_to_t(vertex_t x, vertex_t y)
+    {
+        blossom_t* bx = vertex_top_blossom[x];
+        blossom_t* by = vertex_top_blossom[y];
+
+        BOOST_ASSERT(bx->label == LABEL_S);
+        BOOST_ASSERT(by->label == LABEL_NONE);
+        by->label = LABEL_T;
+        by->tree_edge = std::make_pair(x, y);
+        by->tree_root = bx->tree_root;
+
+        vertex_t y2 = by->base_vertex;
+        vertex_t z = vertex_mate[y2];
+        BOOST_ASSERT(z != null_vertex());
+
+        blossom_t* bz = vertex_top_blossom[z];
+        BOOST_ASSERT(bz->label == LABEL_NONE);
+        BOOST_ASSERT(!bz->best_edge.has_value());
+        bz->label = LABEL_S;
+        bz->tree_edge = std::make_pair(y2, z);
+        bz->tree_root = by->tree_root;
+        add_vertices_to_scan_queue(*bz);
+    }
 
-        // expand odd blossoms with zero dual variables and zero pi value of
-        // their base vertices
-        if (delta == delta2 && delta != delta3)
+    /**
+     * Add the edge between S-vertices "x" and "y".
+     *
+     * If the edge connects blossoms that are part of the same alternating
+     * tree, a new S-blossom is created.
+     *
+     * If the edge connects two different alternating trees, an augmenting
+     * path has been discovered. In this case the matching is augmented.
+     *
+     * @return true if the matching was augmented; otherwise false.
+     */
+    bool add_s_to_s_edge(vertex_t x, vertex_t y)
+    {
+        blossom_t* bx = vertex_top_blossom[x];
+        blossom_t* by = vertex_top_blossom[y];
+        BOOST_ASSERT(bx->label == LABEL_S);
+        BOOST_ASSERT(by->label == LABEL_S);
+        BOOST_ASSERT(bx != by);
+
+        alternating_path_t path = trace_alternating_paths(x, y);
+
+        // Check whether the path starts and ends in the same blossom.
+        blossom_t* b1 = vertex_top_blossom[path.front().first];
+        blossom_t* b2 = vertex_top_blossom[path.back().second];
+        if (b1 == b2)
         {
-            std::vector< blossom_ptr_t > new_top_blossoms;
-            for (blossom_iterator_t bi = top_blossoms.begin();
-                 bi != top_blossoms.end();)
-            {
-                const blossom_ptr_t b = *bi;
-                vertex_descriptor_t b_base = b->get_base();
-                if (b->dual_var == 0
-                    && label_T[b_base] != graph_traits< Graph >::null_vertex()
-                    && pi[b_base] == 0)
-                    bi = expand_T_blossom(bi, new_top_blossoms);
-                else
-                    ++bi;
-            }
-            top_blossoms.insert(top_blossoms.end(), new_top_blossoms.begin(),
-                new_top_blossoms.end());
+            make_blossom(path);
+            return false;
         }
+        else
+        {
+            // Remove labels the two alternating trees on the augmenting path.
+            remove_alternating_tree(bx->tree_root, by->tree_root);
+
+            // Augment matching.
+            augment_matching(path);
+            return true;
+        }
+    }
 
-        while (true)
+    /**
+     * Scan incident edges of newly labeled S-vertices.
+     *
+     * @return true if the matching was augmented; otherwise false.
+     */
+    bool scan_new_s_vertices()
+    {
+        while (!scan_queue.empty())
         {
-            // find a zero-slack critical edge (v, w) of zero gamma values
-            std::pair< edge_descriptor_t, bool > best_edge = null_edge;
-            std::vector< vertex_descriptor_t > base_nodes;
-            for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
-            {
-                if (*vi == base_vertex(*vi))
-                    base_nodes.push_back(*vi);
-            }
-            for (vertex_vec_iter_t i = base_nodes.begin();
-                 i != base_nodes.end(); ++i)
+            vertex_t x = scan_queue.front();
+            scan_queue.pop_front();
+
+            BOOST_ASSERT(vertex_top_blossom[x]->label == LABEL_S);
+
+            for (const edge_t& e : make_iterator_range(out_edges(x, *g)))
             {
-                if (gamma[*i] == 0)
+                BOOST_ASSERT(source(e, *g) == x);
+                vertex_t y = target(e, *g);
+
+                // Note: top level blossom of x may change during this loop,
+                // so we must look it up on each pass.
+                blossom_t* bx = vertex_top_blossom[x];
+                blossom_t* by = vertex_top_blossom[y];
+
+                // Ignore internal edges in blossom.
+                if (bx == by)
+                    continue;
+
+                // Ignore edges with negative slack to prevent numeric overflow.
+                if (get(edge_weights, e) < 0)
+                    continue;
+
+                weight_t slack = edge_slack(e);
+                if (slack <= 0)
                 {
-                    for (vertex_vec_iter_t j = base_nodes.begin();
-                         j != base_nodes.end(); ++j)
+                    // Tight edge.
+                    if (by->label == LABEL_NONE)
+                        extend_tree_s_to_t(x, y);
+                    else if (by->label == LABEL_S)
                     {
-                        if (critical_edge[*i][*j] != null_edge
-                            && slack(critical_edge[*i][*j].first) == 0)
-                            best_edge = critical_edge[*i][*j];
+                        bool augmented = add_s_to_s_edge(x, y);
+                        if (augmented)
+                            return true;
                     }
                 }
-            }
+                else
+                {
+                    // Track non-tight edges between S-blossoms.
+                    if (by->label == LABEL_S)
+                        add_delta3_edge(*bx, e, slack);
+                }
 
-            // if not found, continue finding other augment matching
-            if (best_edge == null_edge)
-            {
-                bool augmented = augment_matching();
-                return augmented || delta != delta1;
+                // Track all edges between S-blossom and non-S-blossom.
+                if (by->label != LABEL_S)
+                    add_delta2_edge(y, e, slack);
             }
-            // if found, determine either augmenting or blossoming
-            vertex_descriptor_t v = source(best_edge.first, g),
-                                w = target(best_edge.first, g);
-            vertex_descriptor_t v_prime = base_vertex(v),
-                                w_prime = base_vertex(w), v_free_ancestor,
-                                w_free_ancestor;
-            vertex_descriptor_t nca = nearest_common_ancestor(
-                v_prime, w_prime, v_free_ancestor, w_free_ancestor);
-            if (nca == graph_traits< Graph >::null_vertex())
-            {
-                augmenting(v, v_free_ancestor, w, w_free_ancestor);
-                return true;
-            }
-            else
-                blossoming(v, v_prime, w, w_prime, nca);
         }
 
         return false;
     }
 
-    template < typename PropertyMap > void get_current_matching(PropertyMap pm)
+    /** Calculate a delta step in the dual LPP problem. */
+    delta_step_t calc_dual_delta_step()
     {
-        vertex_iterator_t vi, vi_end;
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
-            put(pm, *vi, mate[*vi]);
-    }
+        delta_step_t delta{};
 
-private:
-    const Graph& g;
-    VertexIndexMap vm;
-    const std::pair< edge_descriptor_t, bool > null_edge;
-
-    // storage for the property maps below
-    std::vector< vertex_descriptor_t > mate_vector;
-    std::vector< vertex_descriptor_t > label_S_vector, label_T_vector;
-    std::vector< vertex_descriptor_t > outlet_vector;
-    std::vector< vertex_descriptor_t > tau_idx_vector;
-    std::vector< edge_property_t > dual_var_vector;
-    std::vector< edge_property_t > pi_vector, gamma_vector, tau_vector;
-    std::vector< blossom_ptr_t > in_blossom_vector;
-    std::vector< std::pair< vertex_descriptor_t, vertex_descriptor_t > >
-        old_label_vector;
-    std::vector< vertex_to_edge_map_t > critical_edge_vector;
-    std::vector< std::vector< std::pair< edge_descriptor_t, bool > > >
-        critical_edge_vectors;
-
-    // iterator property maps
-    vertex_to_vertex_map_t mate;
-    vertex_to_vertex_map_t label_S; // v has an S-label -> v can be an even
-                                    // vertex, label_S[v] is its mate
-    vertex_to_vertex_map_t
-        label_T; // v has a T-label -> v can be an odd vertex, label_T[v] is its
-                 // predecessor in aug_path
-    vertex_to_vertex_map_t outlet;
-    vertex_to_vertex_map_t tau_idx;
-    vertex_to_weight_map_t dual_var;
-    vertex_to_weight_map_t pi, gamma, tau;
-    vertex_to_blossom_map_t
-        in_blossom; // map any vertex v to the trivial blossom containing v
-    vertex_to_pair_map_t old_label; // <old T-label, old outlet> before
-                                    // relabeling or expanding T-blossoms
-    vertex_pair_to_edge_map_t
-        critical_edge; // an not matched edge (v, w) is critical if v and w
-                       // belongs to different S-blossoms
-
-    vertex_list_t aug_path;
-    edge_list_t even_edges;
-    std::vector< blossom_ptr_t > top_blossoms;
-};
+        // Compute delta1: minimum dual variable of any S-vertex.
+        delta.kind = 1;
+        delta.value = (std::numeric_limits<weight_t>::max)();
+        for (vertex_t x : make_iterator_range(vertices(*g)))
+        {
+            if (vertex_top_blossom[x]->label == LABEL_S)
+                delta.value = (std::min)(delta.value, vertex_dual[x]);
+        }
 
-template < typename Graph, typename MateMap, typename VertexIndexMap >
-void maximum_weighted_matching(const Graph& g, MateMap mate, VertexIndexMap vm)
-{
-    empty_matching< Graph, MateMap >::find_matching(g, mate);
-    weighted_augmenting_path_finder< Graph, MateMap, VertexIndexMap > augmentor(
-        g, mate, vm);
+        // Compute delta2: minimum slack of edge from S-vertex to unlabeled.
+        least_slack_edge_t best_edge = get_least_slack_delta2_edge();
+        if (best_edge.edge.has_value() && (best_edge.slack <= delta.value))
+        {
+            delta.kind = 2;
+            delta.edge = *best_edge.edge;
+            delta.value = best_edge.slack;
+        }
 
-    // can have |V| times augmenting at most
-    for (std::size_t t = 0; t < num_vertices(g); ++t)
-    {
-        bool augmented = false;
-        while (!augmented)
+        // Compute delta3: half minimum slack of edge between S-blossoms.
+        best_edge = get_least_slack_delta3_edge();
+        if (best_edge.edge.has_value() && (best_edge.slack / 2 <= delta.value))
+        {
+            delta.kind = 3;
+            delta.edge = *best_edge.edge;
+            delta.value = best_edge.slack / 2;
+        }
+
+        // Compute delta4: half minimum dual of a top-level T-blossom.
+        for (nontrivial_blossom_t& blossom : nontrivial_blossom)
         {
-            augmented = augmentor.augment_matching();
-            if (!augmented)
+            if ((!blossom.parent)
+                && (blossom.label == LABEL_T)
+                && (blossom.dual_var / 2 <= delta.value))
             {
-                // halt if adjusting dual variables can't bring potential
-                // augment
-                if (!augmentor.adjust_dual())
-                    break;
+                delta.kind = 4;
+                delta.blossom = &blossom;
+                delta.value = blossom.dual_var / 2;
             }
         }
-        if (!augmented)
-            break;
-    }
 
-    augmentor.get_current_matching(mate);
-}
-
-template < typename Graph, typename MateMap >
-inline void maximum_weighted_matching(const Graph& g, MateMap mate)
-{
-    maximum_weighted_matching(g, mate, get(vertex_index, g));
-}
-
-// brute-force matcher searches all possible combinations of matched edges to
-// get the maximum weighted matching which can be used for testing on small
-// graphs (within dozens vertices)
-template < typename Graph, typename MateMap, typename VertexIndexMap >
-class brute_force_matching
-{
-public:
-    typedef
-        typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t;
-    typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t;
-    typedef
-        typename std::vector< vertex_descriptor_t >::iterator vertex_vec_iter_t;
-    typedef typename graph_traits< Graph >::edge_iterator edge_iterator_t;
-    typedef boost::iterator_property_map< vertex_vec_iter_t, VertexIndexMap >
-        vertex_to_vertex_map_t;
+        return delta;
+    }
 
-    brute_force_matching(
-        const Graph& arg_g, MateMap arg_mate, VertexIndexMap arg_vm)
-    : g(arg_g)
-    , vm(arg_vm)
-    , mate_vector(num_vertices(g))
-    , best_mate_vector(num_vertices(g))
-    , mate(mate_vector.begin(), vm)
-    , best_mate(best_mate_vector.begin(), vm)
+    /** Apply a delta step to the dual LPP variables. */
+    void apply_delta_step(weight_t delta)
     {
-        vertex_iterator_t vi, vi_end;
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
-            best_mate[*vi] = mate[*vi] = get(arg_mate, *vi);
+        for (vertex_t x : make_iterator_range(vertices(*g)))
+        {
+            blossom_t* b = vertex_top_blossom[x];
+            if (b->label == LABEL_S)
+                vertex_dual[x] -= delta;
+            else if (b->label == LABEL_T)
+                vertex_dual[x] += delta;
+        }
+
+        for (nontrivial_blossom_t& blossom : nontrivial_blossom)
+        {
+            if (!blossom.parent)
+            {
+                if (blossom.label == LABEL_S)
+                    blossom.dual_var += 2 * delta;
+                else if (blossom.label == LABEL_T)
+                    blossom.dual_var -= 2 * delta;
+            }
+        }
     }
 
-    template < typename PropertyMap > void find_matching(PropertyMap pm)
+    /**
+     * Run one stage of the matching algorithm.
+     *
+     * @return True if the matching was successfully augmented;
+     *         false if no further improvement is possible.
+     */
+    bool run_stage()
     {
-        edge_iterator_t ei;
-        boost::tie(ei, ei_end) = edges(g);
-        select_edge(ei);
+        // Run substages.
+        while (true) {
+            bool augmented = scan_new_s_vertices();
+            if (augmented)
+                return true;
 
-        vertex_iterator_t vi, vi_end;
-        for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
-            put(pm, *vi, best_mate[*vi]);
-    }
+            delta_step_t delta = calc_dual_delta_step();
+            apply_delta_step(delta.value);
 
-private:
-    const Graph& g;
-    VertexIndexMap vm;
-    std::vector< vertex_descriptor_t > mate_vector, best_mate_vector;
-    vertex_to_vertex_map_t mate, best_mate;
-    edge_iterator_t ei_end;
+            if (delta.kind == 2)
+            {
+                vertex_t x = source(delta.edge, *g);
+                vertex_t y = target(delta.edge, *g);
+                if (vertex_top_blossom[x]->label != LABEL_S)
+                    std::swap(x, y);
+                extend_tree_s_to_t(x, y);
+            }
+            else if (delta.kind == 3)
+            {
+                vertex_t x = source(delta.edge, *g);
+                vertex_t y = target(delta.edge, *g);
+                augmented = add_s_to_s_edge(x, y);
+                if (augmented)
+                    return true;
+            }
+            else if (delta.kind == 4)
+            {
+                BOOST_ASSERT(delta.blossom);
+                expand_t_blossom(delta.blossom);
+            }
+            else
+            {
+                // No further improvement possible.
+                BOOST_ASSERT(delta.kind == 1);
+                return false;
+            }
+        }
+    }
 
-    void select_edge(edge_iterator_t ei)
+    /** Run the matching algorithm. */
+    void run()
     {
-        if (ei == ei_end)
+        // Assign label S to all vertices and put them in the queue.
+        for (vertex_t x : make_iterator_range(vertices(*g)))
         {
-            if (matching_weight_sum(g, mate)
-                > matching_weight_sum(g, best_mate))
-            {
-                vertex_iterator_t vi, vi_end;
-                for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
-                    best_mate[*vi] = mate[*vi];
-            }
-            return;
+            BOOST_ASSERT(vertex_mate[x] == null_vertex());
+            blossom_t* bx = vertex_top_blossom[x];
+            BOOST_ASSERT(bx->label == LABEL_NONE);
+            BOOST_ASSERT(bx->base_vertex == x);
+            bx->label = LABEL_S;
+            bx->tree_edge = std::make_pair(null_vertex(), x);
+            bx->tree_root = x;
+            scan_queue.push_back(x);
         }
 
-        vertex_descriptor_t v, w;
-        v = source(*ei, g);
-        w = target(*ei, g);
+        // Improve the solution until no further improvement is possible.
+        while (run_stage()) ;
 
-        select_edge(++ei);
+        // Clear temporary data structures.
+        scan_queue.clear();
+        clear_best_edge();
+    }
 
-        if (mate[v] == graph_traits< Graph >::null_vertex()
-            && mate[w] == graph_traits< Graph >::null_vertex())
-        {
-            mate[v] = w;
-            mate[w] = v;
-            select_edge(ei);
-            mate[v] = mate[w] = graph_traits< Graph >::null_vertex();
-        }
+    /** Copy matching to specified map. */
+    template <typename MateMap>
+    void extract_matching(MateMap mate)
+    {
+        for (vertex_t x : make_iterator_range(vertices(*g)))
+            put(mate, x, vertex_mate[x]);
     }
 };
 
-template < typename Graph, typename MateMap, typename VertexIndexMap >
-void brute_force_maximum_weighted_matching(
-    const Graph& g, MateMap mate, VertexIndexMap vm)
+} // namespace detail
+} // namespace graph
+
+/**
+ * Compute a maximum-weighted matching in an undirected weighted graph.
+ *
+ * This function takes time O(V^3).
+ * This function uses memory size O(V + E).
+ *
+ * @param g         Input graph.
+ *                  The graph type must be a model of VertexListGraph
+ *                  and EdgeListGraph and IncidenceGraph.
+ *                  The graph must be undirected.
+ *                  The graph may not contain parallel edges.
+ *
+ * @param mate      ReadWritePropertyMap, mapping vertices to vertices.
+ *                  This map returns the result of the computation.
+ *                  For each vertex "v", "mate[v]" is the vertex that is
+ *                  matched to "v", or "null_vertex()" if "v" is not matched.
+ *
+ * @param vm        ReadablePropertyMap, mapping vertices to indexes
+ *                  in range [0, num_vertices(g)).
+ *
+ * @param weights   ReadablePropertyMap, mapping edges to edge weights.
+ *                  Edge weights must be integers or floating point values.
+ *
+ * @throw boost::bad_graph  If the input graph is not valid.
+ */
+template <typename Graph, typename MateMap, typename VertexIndexMap,
+    typename EdgeWeightMap>
+void maximum_weighted_matching(
+    const Graph& g, MateMap mate, VertexIndexMap vm, EdgeWeightMap weights)
 {
-    empty_matching< Graph, MateMap >::find_matching(g, mate);
-    brute_force_matching< Graph, MateMap, VertexIndexMap > brute_force_matcher(
-        g, mate, vm);
-    brute_force_matcher.find_matching(mate);
+    BOOST_CONCEPT_ASSERT((VertexListGraphConcept<Graph>));
+    BOOST_CONCEPT_ASSERT((EdgeListGraphConcept<Graph>));
+    BOOST_CONCEPT_ASSERT((IncidenceGraphConcept<Graph>));
+    BOOST_STATIC_ASSERT(is_undirected_graph<Graph>::value);
+
+    using vertex_t = typename graph_traits<Graph>::vertex_descriptor;
+    using edge_t = typename graph_traits<Graph>::edge_descriptor;
+    BOOST_CONCEPT_ASSERT((ReadWritePropertyMapConcept<MateMap, vertex_t>));
+    BOOST_CONCEPT_ASSERT(
+        (ReadablePropertyMapConcept<VertexIndexMap, vertex_t>));
+    BOOST_CONCEPT_ASSERT((ReadablePropertyMapConcept<EdgeWeightMap, edge_t>));
+
+    graph::detail::check_vertex_index_range(g, vm);
+    graph::detail::check_maximum_weighted_matching_edge_weights(g, weights);
+
+    graph::detail::maximum_weighted_matching_context<
+        Graph, VertexIndexMap, EdgeWeightMap>
+        matching(g, vm, weights);
+    matching.run();
+    matching.extract_matching(mate);
 }
 
-template < typename Graph, typename MateMap >
-inline void brute_force_maximum_weighted_matching(const Graph& g, MateMap mate)
+/**
+ * Compute a maximum-weighted matching in an undirected weighted graph.
+ *
+ * This overloaded function obtains edge weights from "get(edge_weight, g)".
+ */
+template <typename Graph, typename MateMap, typename VertexIndexMap>
+inline void maximum_weighted_matching(
+    const Graph& g, MateMap mate, VertexIndexMap vm)
 {
-    brute_force_maximum_weighted_matching(g, mate, get(vertex_index, g));
+    maximum_weighted_matching(g, mate, vm, get(edge_weight, g));
 }
 
+/**
+ * Compute a maximum-weighted matching in an undirected weighted graph.
+ *
+ * This overloaded function obtains vertex indices from "get(vertex_index, g)"
+ * and edge weights from "get(edge_weight, g)".
+ */
+template <typename Graph, typename MateMap>
+inline void maximum_weighted_matching(const Graph& g, MateMap mate)
+{
+    maximum_weighted_matching(g, mate, get(vertex_index, g));
 }
 
-#endif
+} // namespace boost
+
+#endif // BOOST_GRAPH_MAXIMUM_WEIGHTED_MATCHING_HPP
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index b2656eb07..b2e29dd02 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -114,6 +114,7 @@ alias graph_test_regular :
     [ run king_ordering.cpp ]
     [ run matching_test.cpp ]
     [ run weighted_matching_test.cpp ]
+    [ run weighted_matching_test2.cpp ]
     [ run max_flow_test.cpp ]
     [ run boykov_kolmogorov_max_flow_test.cpp ]
     [ run cycle_ratio_tests.cpp ../build//boost_graph : $(CYCLE_RATIO_INPUT_FILE) ]
diff --git a/test/weighted_matching_test.cpp b/test/weighted_matching_test.cpp
index 2baa0adcf..30042c553 100644
--- a/test/weighted_matching_test.cpp
+++ b/test/weighted_matching_test.cpp
@@ -139,6 +139,7 @@ Graph make_graph(typename graph_traits< Graph >::vertices_size_type num_v,
 int main(int, char*[])
 {
     std::ifstream in_file("weighted_matching.dat");
+    BOOST_TEST(in_file.good());
     std::string line;
     while (std::getline(in_file, line))
     {
diff --git a/test/weighted_matching_test2.cpp b/test/weighted_matching_test2.cpp
new file mode 100644
index 000000000..3bdacdeec
--- /dev/null
+++ b/test/weighted_matching_test2.cpp
@@ -0,0 +1,448 @@
+//=======================================================================
+// Copyright (c) 2024 Joris van Rantwijk
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+//=======================================================================
+
+#include <boost/graph/maximum_weighted_matching.hpp>
+
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/adjacency_matrix.hpp>
+#include <boost/core/lightweight_test.hpp>
+#include <iostream>
+#include <limits>
+#include <sstream>
+#include <type_traits>
+#include <vector>
+
+using namespace boost;
+
+
+template <typename WeightType>
+using adj_vec_test_graph = adjacency_list<
+    vecS,  // OutEdgeList
+    vecS,  // VertexList
+    undirectedS,  // DirectedS
+    no_property,  // VertexProperty (vertex_index is implicit)
+    property<edge_weight_t, WeightType>>;  // EdgeProperty
+
+template <typename WeightType>
+using adj_list_test_graph = adjacency_list<
+    listS,  // OutEdgeList
+    listS,  // VertexList
+    undirectedS,  // DirectedS
+    property<vertex_index_t, unsigned int>,  // VertexProperty
+    property<edge_weight_t, WeightType>>;  // EdgeProperty
+
+template <typename WeightType>
+using adj_matrix_test_graph = adjacency_matrix<
+    undirectedS,  // Directed
+    no_property,  // VertexProperty (vertex_index is implicit)
+    property<edge_weight_t, WeightType>>;  // EdgeProperty
+
+
+template <typename WeightType>
+struct edge_info
+{
+    unsigned int x, y;
+    WeightType w;
+};
+
+
+// Initialize vertex_index if necessary.
+template <typename Graph>
+struct vertex_index_installer
+{
+    static void install(Graph&) {}
+};
+
+template <typename WeightType>
+struct vertex_index_installer<adj_list_test_graph<WeightType>>
+{
+    static void install(adj_list_test_graph<WeightType>& g)
+    {
+        auto vrange = vertices(g);
+        unsigned int i = 0;
+        for (auto it = vrange.first; it != vrange.second; ++it, ++i)
+            put(vertex_index, g, *it, i);
+    }
+};
+
+
+template <typename Graph, typename WeightType>
+Graph make_graph(
+    unsigned int nv,
+    const std::vector<edge_info<WeightType>>& edges_info)
+{
+    using EdgeProperty = property<edge_weight_t, WeightType>;
+    Graph g(nv);
+    vertex_index_installer<Graph>::install(g);
+    for (const auto& e : edges_info)
+        add_edge(vertex(e.x, g), vertex(e.y, g), EdgeProperty(e.w), g);
+    return g;
+}
+
+
+template <typename Graph>
+using graph_edge_weight_type = typename property_traits<
+    typename property_map<Graph, edge_weight_t>::type>::value_type;
+
+
+template <typename Graph, typename MateMap>
+std::pair<graph_edge_weight_type<Graph>, bool>
+check_matching(const Graph& g, const MateMap& mate)
+{
+    graph_edge_weight_type<Graph> matching_weight = 0;
+    unsigned int n_matched_edges = 0;
+
+    auto edge_range = edges(g);
+    for (auto it = edge_range.first; it != edge_range.second; ++it)
+    {
+        if (mate[source(*it, g)] == target(*it, g))
+        {
+            ++n_matched_edges;
+            matching_weight += get(edge_weight, g, *it);
+        }
+    }
+
+    unsigned int n_matched_vertices = 0;
+    auto vertex_range = vertices(g);
+    for (auto it = vertex_range.first; it != vertex_range.second; ++it)
+    {
+        if (mate[*it] != graph_traits<Graph>::null_vertex())
+        {
+            if (mate[mate[*it]] != *it)
+                return std::make_pair(0, false);
+            ++n_matched_vertices;
+        }
+    }
+
+    if (n_matched_vertices != 2 * n_matched_edges)
+        return std::make_pair(0, false);
+
+    return std::make_pair(matching_weight, true);
+}
+
+
+template <typename Graph>
+std::string show_graph(const Graph& g)
+{
+    std::ostringstream ss;
+    auto edge_range = edges(g);
+    for (auto it = edge_range.first; it != edge_range.second; ++it)
+    {
+        unsigned int x = get(vertex_index, g, source(*it, g));
+        unsigned int y = get(vertex_index, g, target(*it, g));
+        auto w = get(edge_weight, g, *it);
+        if (ss.tellp() > 0)
+            ss << " ";
+        ss << "{" << x << "," << y << "," << w << "}";
+    }
+    return ss.str();
+}
+
+
+template <typename Graph, typename WeightType>
+void test_matching(const Graph& g, WeightType answer)
+{
+    using vertex_index_map_t =
+        typename property_map<Graph, vertex_index_t>::type;
+    using mate_t = vector_property_map<
+        typename graph_traits<Graph>::vertex_descriptor, vertex_index_map_t>;
+
+    mate_t mate(num_vertices(g));
+    maximum_weighted_matching(g, mate);
+
+    WeightType weight;
+    bool consistent_matching;
+    std::tie(weight, consistent_matching) = check_matching(g, mate);
+    BOOST_TEST(consistent_matching);
+    if (! consistent_matching)
+    {
+        std::cout << std::endl
+                  << "Inconsistent answer for graph" << std::endl
+                  << "  " << show_graph(g) << std::endl;
+    }
+    else
+    {
+        bool same_answer;
+        if (std::numeric_limits<WeightType>::is_integer)
+        {
+            same_answer = (weight == answer);
+        }
+        else
+        {
+            WeightType max_error = std::max(1e-14, answer * 1e-14);
+            same_answer = (std::abs(weight - answer) <= max_error);
+        }
+
+        BOOST_TEST(same_answer);
+        if (! same_answer)
+        {
+            std::cout << std::endl
+                      << "Wrong answer for graph" << std::endl
+                      << "  " << show_graph(g) << std::endl
+                      << "  found weight " << weight
+                      << " while expecting " << answer << std::endl;
+        }
+    }
+}
+
+
+template <typename WeightType = long>
+void run_test_graph(
+    unsigned int nv,
+    const std::vector<edge_info<WeightType>>& edges_info,
+    WeightType answer)
+{
+    using vec_graph_t = adj_vec_test_graph<WeightType>;
+    test_matching(
+        make_graph<vec_graph_t, WeightType>(nv, edges_info), answer);
+
+    using list_graph_t = adj_list_test_graph<WeightType>;
+    test_matching(
+        make_graph<list_graph_t, WeightType>(nv, edges_info), answer);
+
+    using matrix_graph_t = adj_matrix_test_graph<WeightType>;
+    test_matching(
+        make_graph<matrix_graph_t, WeightType>(nv, edges_info), answer);
+}
+
+
+int main(int, char*[])
+{
+    // Simple test cases:
+    run_test_graph(0, {}, 0);
+    run_test_graph(1, {}, 0);
+    run_test_graph(2, {}, 0);
+    run_test_graph(2, {{0, 1, 1}}, 1);
+    run_test_graph(3, {{0, 1, 10}, {1, 2, 11}}, 11);
+    run_test_graph(4, {{0, 1, 5}, {1, 2, 11}, {2, 3, 5}}, 11);
+    run_test_graph(4, {{0, 1, 5}, {1, 2, 11}, {2, 3, 7}}, 12);
+
+    // Floating point edge weights:
+    run_test_graph<double>(4, {
+        {0, 1, 3.1415}, {1, 2, 2.7183}, {0, 2, 3.0}, {0, 3, 1.4142}},
+        4.1325);
+
+    // Negative edge weights:
+    run_test_graph(4, {
+        {0, 1, 2}, {0, 2, -2}, {1, 2, 1}, {1, 3, -1}, {2, 3, -6}},
+        2);
+
+    // Blossoms:
+    run_test_graph(4, {{0, 1, 8}, {0, 2, 9}, {1, 2, 10}, {2, 3, 7}}, 15);
+
+    run_test_graph(6, {
+        {0, 1, 8}, {0, 2, 9}, {1, 2, 10}, {2, 3, 7}, {0, 5, 5}, {3, 4, 6}},
+        21);
+
+    run_test_graph(6, {
+        {0, 1, 9}, {0, 2, 8}, {1, 2, 10}, {0, 3, 5}, {3, 4, 4}, {0, 5, 3}},
+        17);
+
+    run_test_graph(6, {
+        {0, 1, 9}, {0, 2, 8}, {1, 2, 10}, {0, 3, 5}, {3, 4, 3}, {0, 5, 4}},
+        17);
+
+    run_test_graph(6, {
+        {0, 1, 9}, {0, 2, 8}, {1, 2, 10}, {0, 3, 5}, {3, 4, 3}, {2, 5, 4}},
+        16);
+
+    run_test_graph(6, {
+        {0, 1, 9}, {0, 2, 9}, {1, 2, 10}, {1, 3, 8}, {2, 4, 8},
+        {3, 4, 10}, {4, 5, 6}},
+        23);
+
+    run_test_graph(8, {
+        {0, 1, 10}, {0, 6, 10}, {1, 2, 12}, {2, 3, 20}, {2, 4, 20},
+        {3, 4, 25}, {4, 5, 10}, {5, 6, 10}, {6, 7, 8}},
+        48);
+
+    run_test_graph(8, {
+        {0, 1, 8}, {0, 2, 8}, {1, 2, 10}, {1, 3, 12}, {2, 4, 12},
+        {3, 4, 14}, {3, 5, 12}, {4, 6, 12}, {5, 6, 14}, {6, 7, 12}},
+        44);
+
+    run_test_graph(8, {
+        {0, 1, 23}, {0, 4, 22}, {0, 5, 15}, {1, 2, 25}, {2, 3, 22},
+        {3, 4, 25}, {3, 7, 14}, {4, 6, 13}},
+        67);
+
+    run_test_graph(8, {
+        {0, 1, 19}, {0, 2, 20}, {0, 7, 8}, {1, 2, 25}, {1, 3, 18},
+        {2, 4, 18}, {3, 4, 13}, {3, 6, 7}, {4, 5, 7}},
+        47);
+
+    // Somewhat tricky test cases:
+    run_test_graph(10, {
+        {0, 1, 45}, {0, 4, 45}, {1, 2, 50}, {2, 3, 45}, {3, 4, 50},
+        {0, 5, 30}, {2, 8, 35}, {3, 7, 35}, {4, 6, 26}, {8, 9, 5}},
+        146);
+
+    run_test_graph(10, {
+        {0, 1, 45}, {0, 4, 45}, {1, 2, 50}, {2, 3, 45}, {3, 4, 50},
+        {0, 5, 30}, {2, 8, 35}, {3, 7, 26}, {4, 6, 40}, {8, 9, 5}},
+        151);
+
+    run_test_graph(10, {
+        {0, 1, 45}, {0, 4, 45}, {1, 2, 50}, {2, 3, 45}, {3, 4, 50},
+        {0, 5, 30}, {2, 8, 35}, {3, 7, 28}, {4, 6, 26}, {8, 9, 5}},
+        139);
+
+    run_test_graph(12, {
+        {0, 1, 45}, {0, 6, 45}, {1, 2, 50}, {2, 3, 45}, {3, 4, 95},
+        {3, 5, 94}, {4, 5, 94}, {5, 6, 50}, {0, 7, 30}, {2, 10, 35},
+        {4, 8, 36}, {6, 9, 26}, {10, 11, 5}},
+        241);
+
+    run_test_graph(10, {
+        {0, 1, 40}, {0, 2, 40}, {1, 2, 60}, {1, 3, 55}, {2, 4, 55},
+        {3, 4, 50}, {0, 7, 15}, {4, 6, 30}, {6, 5, 10}, {7, 9, 10},
+        {3, 8, 30}},
+        145);
+
+    run_test_graph(6, {
+        {0, 1, 2}, {0, 4, 3}, {1, 2, 7}, {1, 5, 2}, {2, 3, 9},
+        {2, 5, 4}, {3, 4, 8}, {3, 5, 4}},
+        15);
+
+    run_test_graph(6, {
+        {0, 1, 8}, {0, 2, 8}, {1, 2, 9}, {1, 3, 6}, {2, 4, 7},
+        {3, 4, 8}, {3, 5, 5}},
+        20);
+
+    run_test_graph(7, {
+        {0, 1, 7}, {0, 2, 7}, {1, 2, 9}, {0, 3, 7}, {0, 4, 7},
+        {3, 4, 9}, {5, 6, 2}},
+        20);
+
+    run_test_graph(7, {
+        {0, 1, 7}, {0, 4, 6}, {1, 2, 9}, {2, 3, 8}, {3, 4, 9},
+        {3, 5, 1}, {4, 5, 1}},
+        18);
+
+    run_test_graph(5, {
+        {0, 2, 7}, {0, 3, 3}, {1, 3, 1}, {2, 3, 5}, {2, 4, 5}},
+        8);
+
+    run_test_graph(7, {
+        {0, 1, 15}, {0, 2, 10}, {0, 3, 11}, {0, 5, 17}, {1, 4, 12},
+        {1, 5, 8}, {4, 6, 15}, {5, 6, 7}},
+        34);
+
+    run_test_graph(7, {
+        {0, 1, 12}, {0, 3, 11}, {0, 4, 11}, {1, 3, 12}, {1, 5, 11},
+        {2, 3, 14}, {2, 6, 14}, {3, 6, 11}},
+        37);
+
+    run_test_graph(7, {
+        {0, 1, 19}, {0, 3, 17}, {0, 4, 19}, {1, 2, 15}, {1, 4, 21},
+        {3, 5, 18}, {3, 6, 11}, {4, 5, 19}},
+        52);
+
+    run_test_graph(6, {
+        {0, 1, 19}, {0, 2, 19}, {0, 3, 15}, {0, 4, 17}, {1, 2, 21},
+        {3, 4, 16}, {4, 5, 10}},
+        46);
+
+    run_test_graph(6, {
+        {0, 1, 48}, {0, 3, 42}, {0, 4, 57}, {1, 3, 51}, {1, 5, 36},
+        {2, 3, 23}, {4, 5, 46}},
+        117);
+
+    run_test_graph(6, {
+        {0, 2, 7}, {0, 5, 20}, {1, 2, 50}, {1, 4, 46}, {2, 3, 35},
+        {2, 4, 8}, {2, 5, 25}, {3, 5, 47}},
+        101);
+
+    run_test_graph(5, {
+        {0, 1, 48}, {0, 2, 44}, {0, 4, 48}, {1, 4, 36}, {3, 4, 31}},
+        80);
+
+    run_test_graph(6, {
+        {3, 5, 71}, {2, 4, 36}, {3, 4, 27}, {5, 1, 29}, {2, 0, 48},
+        {1, 0, 60}, {1, 3, 11}, {4, 5, 54}},
+        167);
+
+    run_test_graph(6, {
+        {0, 2, 92}, {5, 4, 50}, {5, 0, 86}, {2, 1, 87}, {1, 3, 39},
+        {4, 0, 2}, {4, 1, 83}, {0, 1, 56}},
+        181);
+
+    // Triangles, unit weight:
+    run_test_graph(9, {
+        {0, 1, 1}, {0, 2, 1}, {1, 2, 1}, {0, 3, 1},
+        {3, 4, 1}, {3, 5, 1}, {4, 5, 1}, {4, 7, 1},
+        {6, 7, 1}, {6, 8, 1}, {7, 8, 1}},
+        4);
+
+    // Trigger known bugs in a previous version of maximum_weighted_matching:
+
+    // wrong answer
+    run_test_graph(8, {
+        {4, 7, 453}, {0, 4, 627}, {4, 6, 853}, {6, 7, 344}, {5, 6, 906},
+        {4, 5, 689}, {2, 3, 741}, {2, 7, 746}, {3, 6, 647}, {1, 5, 385},
+        {5, 7, 215}, {3, 7, 640}},
+        2405);
+
+    // wrong answer
+    run_test_graph(12, {
+        {0, 3, 448}, {0, 4, 919}, {0, 5, 918}, {0, 9, 72}, {1, 3, 830},
+        {1, 4, 687}, {1, 5, 559}, {1, 6, 580}, {1, 7, 679}, {1, 8, 627},
+        {2, 4, 585}, {2, 5, 835}, {2, 6, 822}, {2, 8, 462}, {2, 11, 380},
+        {3, 4, 328}, {3, 5, 860}, {3, 7, 297}, {3, 8, 590}, {3, 11, 235},
+        {4, 6, 692}, {4, 10, 227}, {4, 11, 354}, {5, 6, 160}, {5, 9, 400},
+        {6, 8, 410}, {6, 9, 420}, {7, 11, 924}, {9, 11, 825}, {10, 11, 790}},
+        4233);
+
+    // assertion
+    run_test_graph(9, {
+        {4, 6, 796}, {3, 4, 553}, {0, 1, 792}, {1, 7, 1000}, {4, 5, 360},
+        {5, 7, 183}, {4, 8, 694}, {0, 4, 741}, {0, 2, 483}, {5, 8, 228},
+        {7, 8, 644}, {1, 3, 236}, {6, 7, 895}, {1, 6, 913}, {1, 8, 617}},
+        2593);
+
+    // assertion
+    run_test_graph(10, {
+        {0, 6, 892}, {2, 4, 517}, {4, 7, 560}, {1, 5, 828}, {1, 7, 831},
+        {6, 7, 397}, {4, 5, 43}, {5, 6, 944}, {6, 9, 215}, {5, 7, 753},
+        {4, 6, 901}, {1, 2, 530}, {2, 8, 384}, {3, 4, 499}, {5, 8, 190}},
+        2674);
+
+    // segmentation fault
+    run_test_graph(7, {
+        {5, 6, 799}, {0, 4, 601}, {0, 1, 578}, {0, 3, 373}, {4, 6, 675},
+        {4, 5, 925}, {0, 5, 697}, {3, 4, 260}, {1, 3, 464}, {1, 5, 845},
+        {0, 2, 176}, {2, 5, 685}},
+        1938);
+
+    // segmentation fault
+    run_test_graph(8, {
+        {2, 7, 420}, {6, 7, 414}, {2, 3, 421}, {2, 6, 854}, {1, 6, 997},
+        {4, 5, 46}, {1, 2, 467}, {2, 4, 230}, {3, 7, 555}, {0, 3, 334},
+        {0, 7, 341}, {0, 6, 634}},
+        1805);
+
+    // hang
+    run_test_graph(7, {
+        {2, 3, 837}, {1, 4, 458}, {3, 4, 291}, {5, 6, 601}, {1, 2, 202},
+        {0, 4, 491}, {4, 5, 910}, {0, 1, 159}, {3, 5, 684}, {4, 6, 139},
+        {0, 2, 792}, {1, 3, 232}},
+        1934);
+
+    // hang
+    run_test_graph(12, {
+        {0, 1, 856}, {0, 4, 462}, {0, 6, 874}, {0, 7, 406}, {0, 10, 294},
+        {1, 3, 936}, {1, 5, 852}, {1, 7, 501}, {1, 8, 555}, {1, 10, 41},
+        {2, 3, 325}, {2, 6, 748}, {2, 9, 808}, {3, 5, 870}, {3, 7, 25},
+        {3, 8, 663}, {3, 9, 897}, {4, 7, 617}, {4, 9, 435}, {4, 10, 818},
+        {4, 11, 933}, {5, 6, 608}, {5, 8, 636}, {6, 9, 274}, {6, 10, 279},
+        {7, 9, 705}, {7, 10, 114}, {8, 11, 602}, {9, 10, 764}, {10, 11, 413}},
+        4599);
+
+    return boost::report_errors();
+}
