diff --git a/src/edmonds.cpp b/src/edmonds.cpp
index 6918b43..7849e88 100644
--- a/src/edmonds.cpp
+++ b/src/edmonds.cpp
@@ -103,10 +103,10 @@ void maximum_matching_from_initial_matching(Graph & graph)
       for(NodeId neighbor_id : graph.node(id).neighbors())
       {
          check_integrity(graph);
-         std::cout << "Check passed" << std::endl;
+         //std::cout << "Check passed" << std::endl;
          if (graph.is_out_of_forest(neighbor_id))
          {
-            std::cout << "Grow forest" << std::endl;
+            //std::cout << "Grow forest" << std::endl;
             // Grow Forest
             graph.node(neighbor_id).phi = id;
          }
@@ -117,7 +117,7 @@ void maximum_matching_from_initial_matching(Graph & graph)
 
             if (x_path.back() != y_path.back())
             {
-               std::cout << "Augment" << std::endl;
+               //std::cout << "Augment" << std::endl;
                // Paths are disjoint -> augment
                graph.node(x_path.front()).matched_neighbor = y_path.front();
                graph.node(y_path.front()).matched_neighbor = x_path.front();
@@ -140,7 +140,7 @@ void maximum_matching_from_initial_matching(Graph & graph)
             }
             else
             {
-               std::cout << "Contract blossom" << std::endl;
+               //std::cout << "Contract blossom" << std::endl;
                // Paths are not disjoint -> shrink blossom
                size_t distance_from_x = x_path.size() - 1;
                size_t distance_from_y = y_path.size() - 1;
@@ -181,7 +181,10 @@ void maximum_matching_from_initial_matching(Graph & graph)
                   graph.node(y_path.front()).phi = x_path.front();
                }
 
-               // Update root indices
+               // Update root indices. We have to do this for all vertices v with
+               // ρ(v) in the paths from x or y to r
+               // We update ρ(v) first for the paths themselves, and then 'contract' ρ
+               // by updating ρ(v) to r for all vertices where ρ(ρ(v)) = r
                for(size_t i = 0; i <= distance_from_x; ++i)
                {
                   graph.node(x_path[i]).rho = blossom_root_id;
@@ -190,6 +193,14 @@ void maximum_matching_from_initial_matching(Graph & graph)
                {
                   graph.node(y_path[i]).rho = blossom_root_id;
                }
+               for(NodeId node_id = 0; node_id < graph.num_nodes(); ++node_id)
+               {
+                  if (graph.rho(graph.rho(node_id)) == blossom_root_id)
+                  {
+                     graph.node(node_id).rho = blossom_root_id;
+                  }
+               }
+               check_integrity(graph);
             }
          }
       }