Index: lemon/bipartite_matching.h
===================================================================
--- lemon/bipartite_matching.h	(revision 2463)
+++ lemon/bipartite_matching.h	(revision 2466)
@@ -70,6 +70,6 @@
     ///
     /// Constructor of the algorithm. 
-    MaxBipartiteMatching(const BpUGraph& _graph) 
-      : anode_matching(_graph), bnode_matching(_graph), graph(&_graph) {}
+    MaxBipartiteMatching(const BpUGraph& graph) 
+      : _matching(graph), _rmatching(graph), _reached(graph), _graph(&graph) {}
 
     /// \name Execution control
@@ -88,11 +88,13 @@
     /// It initalizes the data structures and creates an empty matching.
     void init() {
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        anode_matching[it] = INVALID;
-      }
-      for (BNodeIt it(*graph); it != INVALID; ++it) {
-        bnode_matching[it] = INVALID;
-      }
-      matching_size = 0;
+      for (ANodeIt it(*_graph); it != INVALID; ++it) {
+        _matching.set(it, INVALID);
+      }
+      for (BNodeIt it(*_graph); it != INVALID; ++it) {
+        _rmatching.set(it, INVALID);
+	_reached.set(it, -1);
+      }
+      _size = 0;
+      _phase = -1;
     }
 
@@ -103,19 +105,22 @@
     /// the matching than from the initial empty matching.
     void greedyInit() {
-      matching_size = 0;
-      for (BNodeIt it(*graph); it != INVALID; ++it) {
-        bnode_matching[it] = INVALID;
-      }
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        anode_matching[it] = INVALID;
-        for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
-          if (bnode_matching[graph->bNode(jt)] == INVALID) {
-            anode_matching[it] = jt;
-            bnode_matching[graph->bNode(jt)] = jt;
-            ++matching_size;
+      _size = 0;
+      for (BNodeIt it(*_graph); it != INVALID; ++it) {
+        _rmatching.set(it, INVALID);
+	_reached.set(it, 0);
+      }
+      for (ANodeIt it(*_graph); it != INVALID; ++it) {
+        _matching[it] = INVALID;
+        for (IncEdgeIt jt(*_graph, it); jt != INVALID; ++jt) {
+          if (_rmatching[_graph->bNode(jt)] == INVALID) {
+            _matching.set(it, jt);
+	    _rmatching.set(_graph->bNode(jt), jt);
+	    _reached.set(it, -1);
+            ++_size;
             break;
           }
         }
       }
+      _phase = 0;
     }
 
@@ -125,18 +130,21 @@
     template <typename MatchingMap>
     void matchingInit(const MatchingMap& mm) {
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        anode_matching[it] = INVALID;
-      }
-      for (BNodeIt it(*graph); it != INVALID; ++it) {
-        bnode_matching[it] = INVALID;
-      }
-      matching_size = 0;
-      for (UEdgeIt it(*graph); it != INVALID; ++it) {
+      for (ANodeIt it(*_graph); it != INVALID; ++it) {
+        _matching.set(it, INVALID);
+      }
+      for (BNodeIt it(*_graph); it != INVALID; ++it) {
+        _rmatching.set(it, INVALID);
+	_reached.set(it, 0);
+      }
+      _size = 0;
+      for (UEdgeIt it(*_graph); it != INVALID; ++it) {
         if (mm[it]) {
-          ++matching_size;
-          anode_matching[graph->aNode(it)] = it;
-          bnode_matching[graph->bNode(it)] = it;
-        }
-      }
+          ++_size;
+          _matching.set(_graph->aNode(it), it);
+          _rmatching.set(_graph->bNode(it), it);
+	  _reached.set(it, 0);
+        }
+      }
+      _phase = 0;
     }
 
@@ -146,47 +154,78 @@
     /// \return %True when the given map contains really a matching.
     template <typename MatchingMap>
-    void checkedMatchingInit(const MatchingMap& mm) {
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        anode_matching[it] = INVALID;
-      }
-      for (BNodeIt it(*graph); it != INVALID; ++it) {
-        bnode_matching[it] = INVALID;
-      }
-      matching_size = 0;
-      for (UEdgeIt it(*graph); it != INVALID; ++it) {
+    bool checkedMatchingInit(const MatchingMap& mm) {
+      for (ANodeIt it(*_graph); it != INVALID; ++it) {
+        _matching.set(it, INVALID);
+      }
+      for (BNodeIt it(*_graph); it != INVALID; ++it) {
+        _rmatching.set(it, INVALID);
+	_reached.set(it, 0);
+      }
+      _size = 0;
+      for (UEdgeIt it(*_graph); it != INVALID; ++it) {
         if (mm[it]) {
-          ++matching_size;
-          if (anode_matching[graph->aNode(it)] != INVALID) {
+          ++_size;
+          if (_matching[_graph->aNode(it)] != INVALID) {
             return false;
           }
-          anode_matching[graph->aNode(it)] = it;
-          if (bnode_matching[graph->aNode(it)] != INVALID) {
+          _matching.set(_graph->aNode(it), it);
+          if (_matching[_graph->bNode(it)] != INVALID) {
             return false;
           }
-          bnode_matching[graph->bNode(it)] = it;
-        }
+          _matching.set(_graph->bNode(it), it);
+	  _reached.set(_graph->bNode(it), -1);
+        }
+      }
+      _phase = 0;
+      return true;
+    }
+
+  private:
+    
+    bool _find_path(Node anode, int maxlevel,
+		    typename Graph::template BNodeMap<int>& level) {
+      for (IncEdgeIt it(*_graph, anode); it != INVALID; ++it) {
+	Node bnode = _graph->bNode(it); 
+	if (level[bnode] == maxlevel) {
+	  level.set(bnode, -1);
+	  if (maxlevel == 0) {
+	    _matching.set(anode, it);
+	    _rmatching.set(bnode, it);
+	    return true;
+	  } else {
+	    Node nnode = _graph->aNode(_rmatching[bnode]);
+	    if (_find_path(nnode, maxlevel - 1, level)) {
+	      _matching.set(anode, it);
+	      _rmatching.set(bnode, it);
+	      return true;
+	    }
+	  }
+	}
       }
       return false;
     }
 
+  public:
+
     /// \brief An augmenting phase of the Hopcroft-Karp algorithm
     ///
     /// It runs an augmenting phase of the Hopcroft-Karp
-    /// algorithm. This phase finds maximum count of edge disjoint
-    /// augmenting paths and augments on these paths. The algorithm
-    /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is
-    /// \f$ O(e) \f$ long.
+    /// algorithm. This phase finds maximal edge disjoint augmenting
+    /// paths and augments on these paths. The algorithm consists at
+    /// most of \f$ O(\sqrt{n}) \f$ phase and one phase is \f$ O(e)
+    /// \f$ long.
     bool augment() {
 
-      typename Graph::template ANodeMap<bool> areached(*graph, false);
-      typename Graph::template BNodeMap<bool> breached(*graph, false);
-
-      typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
-
-      std::vector<Node> queue, bqueue;
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        if (anode_matching[it] == INVALID) {
+      ++_phase;
+      
+      typename Graph::template BNodeMap<int> _level(*_graph, -1);
+      typename Graph::template ANodeMap<bool> _found(*_graph, false);
+
+      std::vector<Node> queue, aqueue;
+      for (BNodeIt it(*_graph); it != INVALID; ++it) {
+        if (_rmatching[it] == INVALID) {
           queue.push_back(it);
-          areached[it] = true;
+          _reached.set(it, _phase);
+	  _level.set(it, 0);
         }
       }
@@ -194,134 +233,106 @@
       bool success = false;
 
+      int level = 0;
       while (!success && !queue.empty()) {
-        std::vector<Node> newqueue;
+        std::vector<Node> nqueue;
         for (int i = 0; i < int(queue.size()); ++i) {
-          Node anode = queue[i];
-          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
-            Node bnode = graph->bNode(jt);
-            if (breached[bnode]) continue;
-            breached[bnode] = true;
-            bpred[bnode] = jt;
-            if (bnode_matching[bnode] == INVALID) {
-              bqueue.push_back(bnode);
+          Node bnode = queue[i];
+          for (IncEdgeIt jt(*_graph, bnode); jt != INVALID; ++jt) {
+            Node anode = _graph->aNode(jt);
+            if (_matching[anode] == INVALID) {
+
+	      if (!_found[anode]) {
+		if (_find_path(anode, level, _level)) {
+		  ++_size;
+		}
+		_found.set(anode, true);
+	      }
               success = true;
             } else {           
-              Node newanode = graph->aNode(bnode_matching[bnode]);
-              if (!areached[newanode]) {
-                areached[newanode] = true;
-                newqueue.push_back(newanode);
+              Node nnode = _graph->bNode(_matching[anode]);
+              if (_reached[nnode] != _phase) {
+                _reached.set(nnode, _phase);
+                nqueue.push_back(nnode);
+		_level.set(nnode, level + 1);
               }
             }
           }
         }
-        queue.swap(newqueue);
-      }
-
-      if (success) {
-
-        typename Graph::template ANodeMap<bool> aused(*graph, false);
-        
-        for (int i = 0; i < int(bqueue.size()); ++i) {
-          Node bnode = bqueue[i];
-
-          bool used = false;
-
-          while (bnode != INVALID) {
-            UEdge uedge = bpred[bnode];
-            Node anode = graph->aNode(uedge);
-            
-            if (aused[anode]) {
-              used = true;
-              break;
-            }
-            
-            bnode = anode_matching[anode] != INVALID ? 
-              graph->bNode(anode_matching[anode]) : INVALID;
-            
-          }
-          
-          if (used) continue;
-
-          bnode = bqueue[i];
-          while (bnode != INVALID) {
-            UEdge uedge = bpred[bnode];
-            Node anode = graph->aNode(uedge);
-            
-            bnode_matching[bnode] = uedge;
-
-            bnode = anode_matching[anode] != INVALID ? 
-              graph->bNode(anode_matching[anode]) : INVALID;
-            
-            anode_matching[anode] = uedge;
-
-            aused[anode] = true;
-          }
-          ++matching_size;
-
-        }
-      }
+	++level;
+        queue.swap(nqueue);
+      }
+      
       return success;
     }
-
-    /// \brief An augmenting phase of the Ford-Fulkerson algorithm
-    ///
-    /// It runs an augmenting phase of the Ford-Fulkerson
-    /// algorithm. This phase finds only one augmenting path and 
-    /// augments only on this paths. The algorithm consists at most 
-    /// of \f$ O(n) \f$ simple phase and one phase is at most 
-    /// \f$ O(e) \f$ long.
-    bool simpleAugment() {
-
-      typename Graph::template ANodeMap<bool> areached(*graph, false);
-      typename Graph::template BNodeMap<bool> breached(*graph, false);
-
-      typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
-
-      std::vector<Node> queue;
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        if (anode_matching[it] == INVALID) {
+  private:
+    
+    void _find_path_bfs(Node anode,
+			typename Graph::template ANodeMap<UEdge>& pred) {
+      while (true) {
+	UEdge uedge = pred[anode];
+	Node bnode = _graph->bNode(uedge);
+
+	UEdge nedge = _rmatching[bnode];
+
+	_matching.set(anode, uedge);
+	_rmatching.set(bnode, uedge);
+
+	if (nedge == INVALID) break;
+	anode = _graph->aNode(nedge);
+      }
+    }
+
+  public:
+
+    /// \brief An augmenting phase with single path augementing
+    ///
+    /// This phase finds only one augmenting paths and augments on
+    /// these paths. The algorithm consists at most of \f$ O(n) \f$
+    /// phase and one phase is \f$ O(e) \f$ long.
+    bool simpleAugment() { 
+      ++_phase;
+      
+      typename Graph::template ANodeMap<UEdge> _pred(*_graph);
+
+      std::vector<Node> queue, aqueue;
+      for (BNodeIt it(*_graph); it != INVALID; ++it) {
+        if (_rmatching[it] == INVALID) {
           queue.push_back(it);
-          areached[it] = true;
-        }
-      }
-
-      while (!queue.empty()) {
-        std::vector<Node> newqueue;
+          _reached.set(it, _phase);
+        }
+      }
+
+      bool success = false;
+
+      int level = 0;
+      while (!success && !queue.empty()) {
+        std::vector<Node> nqueue;
         for (int i = 0; i < int(queue.size()); ++i) {
-          Node anode = queue[i];
-          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
-            Node bnode = graph->bNode(jt);
-            if (breached[bnode]) continue;
-            breached[bnode] = true;
-            bpred[bnode] = jt;
-            if (bnode_matching[bnode] == INVALID) {
-              while (bnode != INVALID) {
-                UEdge uedge = bpred[bnode];
-                anode = graph->aNode(uedge);
-                
-                bnode_matching[bnode] = uedge;
-                
-                bnode = anode_matching[anode] != INVALID ? 
-                  graph->bNode(anode_matching[anode]) : INVALID;
-                
-                anode_matching[anode] = uedge;
-                
-              }
-              ++matching_size;
-              return true;
+          Node bnode = queue[i];
+          for (IncEdgeIt jt(*_graph, bnode); jt != INVALID; ++jt) {
+            Node anode = _graph->aNode(jt);
+            if (_matching[anode] == INVALID) {
+	      _pred.set(anode, jt);
+	      _find_path_bfs(anode, _pred);
+	      ++_size;
+	      return true;
             } else {           
-              Node newanode = graph->aNode(bnode_matching[bnode]);
-              if (!areached[newanode]) {
-                areached[newanode] = true;
-                newqueue.push_back(newanode);
+              Node nnode = _graph->bNode(_matching[anode]);
+              if (_reached[nnode] != _phase) {
+		_pred.set(anode, jt);
+		_reached.set(nnode, _phase);
+                nqueue.push_back(nnode);
               }
             }
           }
         }
-        queue.swap(newqueue);
+	++level;
+        queue.swap(nqueue);
       }
       
-      return false;
-    }
+      return success;
+    }
+
+
 
     /// \brief Starts the algorithm.
@@ -351,4 +362,25 @@
     ///@{
 
+    /// \brief Return true if the given uedge is in the matching.
+    /// 
+    /// It returns true if the given uedge is in the matching.
+    bool matchingEdge(const UEdge& edge) const {
+      return _matching[_graph->aNode(edge)] == edge;
+    }
+
+    /// \brief Returns the matching edge from the node.
+    /// 
+    /// Returns the matching edge from the node. If there is not such
+    /// edge it gives back \c INVALID.
+    /// \note If the parameter node is a B-node then the running time is
+    /// propotional to the degree of the node.
+    UEdge matchingEdge(const Node& node) const {
+      if (_graph->aNode(node)) {
+        return _matching[node];
+      } else {
+        return _rmatching[node];
+      }
+    }
+
     /// \brief Set true all matching uedge in the map.
     /// 
@@ -358,10 +390,10 @@
     template <typename MatchingMap>
     int quickMatching(MatchingMap& mm) const {
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        if (anode_matching[it] != INVALID) {
-          mm.set(anode_matching[it], true);
-        }
-      }
-      return matching_size;
+      for (ANodeIt it(*_graph); it != INVALID; ++it) {
+        if (_matching[it] != INVALID) {
+          mm.set(_matching[it], true);
+        }
+      }
+      return _size;
     }
 
@@ -372,8 +404,8 @@
     template <typename MatchingMap>
     int matching(MatchingMap& mm) const {
-      for (UEdgeIt it(*graph); it != INVALID; ++it) {
-        mm.set(it, it == anode_matching[graph->aNode(it)]);
-      }
-      return matching_size;
+      for (UEdgeIt it(*_graph); it != INVALID; ++it) {
+        mm.set(it, it == _matching[_graph->aNode(it)]);
+      }
+      return _size;
     }
 
@@ -385,8 +417,8 @@
     template<class MatchingMap>
     int aMatching(MatchingMap& mm) const {
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        mm.set(it, anode_matching[it]);
-      }
-      return matching_size;
+      for (ANodeIt it(*_graph); it != INVALID; ++it) {
+        mm.set(it, _matching[it]);
+      }
+      return _size;
     }
 
@@ -398,34 +430,8 @@
     template<class MatchingMap>
     int bMatching(MatchingMap& mm) const {
-      for (BNodeIt it(*graph); it != INVALID; ++it) {
-        mm.set(it, bnode_matching[it]);
-      }
-      return matching_size;
-    }
-
-    /// \brief Return true if the given uedge is in the matching.
-    /// 
-    /// It returns true if the given uedge is in the matching.
-    bool matchingEdge(const UEdge& edge) const {
-      return anode_matching[graph->aNode(edge)] == edge;
-    }
-
-    /// \brief Returns the matching edge from the node.
-    /// 
-    /// Returns the matching edge from the node. If there is not such
-    /// edge it gives back \c INVALID.
-    UEdge matchingEdge(const Node& node) const {
-      if (graph->aNode(node)) {
-        return anode_matching[node];
-      } else {
-        return bnode_matching[node];
-      }
-    }
-
-    /// \brief Gives back the number of the matching edges.
-    ///
-    /// Gives back the number of the matching edges.
-    int matchingSize() const {
-      return matching_size;
+      for (BNodeIt it(*_graph); it != INVALID; ++it) {
+        mm.set(it, _rmatching[it]);
+      }
+      return _size;
     }
 
@@ -439,46 +445,15 @@
     int coverSet(CoverMap& covering) const {
 
-      typename Graph::template ANodeMap<bool> areached(*graph, false);
-      typename Graph::template BNodeMap<bool> breached(*graph, false);
-      
-      std::vector<Node> queue;
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        if (anode_matching[it] == INVALID) {
-          queue.push_back(it);
-        }
-      }
-
-      while (!queue.empty()) {
-        std::vector<Node> newqueue;
-        for (int i = 0; i < int(queue.size()); ++i) {
-          Node anode = queue[i];
-          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
-            Node bnode = graph->bNode(jt);
-            if (breached[bnode]) continue;
-            breached[bnode] = true;
-            if (bnode_matching[bnode] != INVALID) {
-              Node newanode = graph->aNode(bnode_matching[bnode]);
-              if (!areached[newanode]) {
-                areached[newanode] = true;
-                newqueue.push_back(newanode);
-              }
-            }
-          }
-        }
-        queue.swap(newqueue);
-      }
-
       int size = 0;
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        covering.set(it, !areached[it] && anode_matching[it] != INVALID);
-        if (!areached[it] && anode_matching[it] != INVALID) {
-          ++size;
-        }
-      }
-      for (BNodeIt it(*graph); it != INVALID; ++it) {
-        covering.set(it, breached[it]);
-        if (breached[it]) {
-          ++size;
-        }
+      for (ANodeIt it(*_graph); it != INVALID; ++it) {
+	bool cn = _matching[it] != INVALID && 
+	  _reached[_graph->bNode(_matching[it])] == _phase;
+        covering.set(it, cn);
+        if (cn) ++size;
+      }
+      for (BNodeIt it(*_graph); it != INVALID; ++it) {
+	bool cn = _reached[it] != _phase;
+        covering.set(it, cn);
+        if (cn) ++size;
       }
       return size;
@@ -486,5 +461,5 @@
 
     /// \brief Gives back a barrier on the A-nodes
-    
+    ///    
     /// The barrier is s subset of the nodes on the same side of the
     /// graph, which size minus its neighbours is exactly the
@@ -494,41 +469,12 @@
     void aBarrier(BarrierMap& barrier) const {
 
-      typename Graph::template ANodeMap<bool> areached(*graph, false);
-      typename Graph::template BNodeMap<bool> breached(*graph, false);
-      
-      std::vector<Node> queue;
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        if (anode_matching[it] == INVALID) {
-          queue.push_back(it);
-        }
-      }
-
-      while (!queue.empty()) {
-        std::vector<Node> newqueue;
-        for (int i = 0; i < int(queue.size()); ++i) {
-          Node anode = queue[i];
-          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
-            Node bnode = graph->bNode(jt);
-            if (breached[bnode]) continue;
-            breached[bnode] = true;
-            if (bnode_matching[bnode] != INVALID) {
-              Node newanode = graph->aNode(bnode_matching[bnode]);
-              if (!areached[newanode]) {
-                areached[newanode] = true;
-                newqueue.push_back(newanode);
-              }
-            }
-          }
-        }
-        queue.swap(newqueue);
-      }
-
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        barrier.set(it, areached[it] || anode_matching[it] == INVALID);
+      for (ANodeIt it(*_graph); it != INVALID; ++it) {
+        barrier.set(it, _matching[it] == INVALID || 
+		    _reached[_graph->bNode(_matching[it])] != _phase);
       }
     }
 
     /// \brief Gives back a barrier on the B-nodes
-    
+    ///    
     /// The barrier is s subset of the nodes on the same side of the
     /// graph, which size minus its neighbours is exactly the
@@ -538,37 +484,14 @@
     void bBarrier(BarrierMap& barrier) const {
 
-      typename Graph::template ANodeMap<bool> areached(*graph, false);
-      typename Graph::template BNodeMap<bool> breached(*graph, false);
-      
-      std::vector<Node> queue;
-      for (ANodeIt it(*graph); it != INVALID; ++it) {
-        if (anode_matching[it] == INVALID) {
-          queue.push_back(it);
-        }
-      }
-
-      while (!queue.empty()) {
-        std::vector<Node> newqueue;
-        for (int i = 0; i < int(queue.size()); ++i) {
-          Node anode = queue[i];
-          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
-            Node bnode = graph->bNode(jt);
-            if (breached[bnode]) continue;
-            breached[bnode] = true;
-            if (bnode_matching[bnode] != INVALID) {
-              Node newanode = graph->aNode(bnode_matching[bnode]);
-              if (!areached[newanode]) {
-                areached[newanode] = true;
-                newqueue.push_back(newanode);
-              }
-            }
-          }
-        }
-        queue.swap(newqueue);
-      }
-
-      for (BNodeIt it(*graph); it != INVALID; ++it) {
-        barrier.set(it, !breached[it]);
-      }
+      for (BNodeIt it(*_graph); it != INVALID; ++it) {
+        barrier.set(it, _reached[it] == _phase);
+      }
+    }
+
+    /// \brief Gives back the number of the matching edges.
+    ///
+    /// Gives back the number of the matching edges.
+    int matchingSize() const {
+      return _size;
     }
 
@@ -577,9 +500,13 @@
   private:
 
-    ANodeMatchingMap anode_matching;
-    BNodeMatchingMap bnode_matching;
-    const Graph *graph;
-
-    int matching_size;
+    typename BpUGraph::template ANodeMap<UEdge> _matching;
+    typename BpUGraph::template BNodeMap<UEdge> _rmatching;
+
+    typename BpUGraph::template BNodeMap<int> _reached;
+
+    int _phase;
+    const Graph *_graph;
+
+    int _size;
   
   };
Index: lemon/pr_bipartite_matching.h
===================================================================
--- lemon/pr_bipartite_matching.h	(revision 2463)
+++ lemon/pr_bipartite_matching.h	(revision 2466)
@@ -60,4 +60,8 @@
   public:
 
+    /// Constructor
+
+    /// Constructor
+    ///
     PrBipartiteMatching(const Graph &g) :
       _g(g),
@@ -196,12 +200,13 @@
       }
       
-      _matching_size = _node_num;
-      for(ANodeIt n(_g);n!=INVALID;++n)
-	if(_matching[n]==INVALID) _matching_size--;
-	else if (_cov[_g.bNode(_matching[n])]>1) {
+      for(ANodeIt n(_g);n!=INVALID;++n) {
+	if (_matching[n]==INVALID)continue;
+	if (_cov[_g.bNode(_matching[n])]>1) {
 	  _cov[_g.bNode(_matching[n])]--;
-	  _matching_size--;
 	  _matching[n]=INVALID;
-	}
+	} else {
+	  ++_matching_size;
+	}
+      }
     }
 
@@ -262,4 +267,5 @@
 	  return false;
       }
+      _matching_size = _node_num;
       return true;
     }
