lemon/bipartite_matching.h
changeset 2466 feb7974cf4ec
parent 2463 19651a04d056
child 2488 da94e3b332f3
     1.1 --- a/lemon/bipartite_matching.h	Sat Aug 25 10:12:03 2007 +0000
     1.2 +++ b/lemon/bipartite_matching.h	Tue Aug 28 13:58:54 2007 +0000
     1.3 @@ -69,8 +69,8 @@
     1.4      /// \brief Constructor.
     1.5      ///
     1.6      /// Constructor of the algorithm. 
     1.7 -    MaxBipartiteMatching(const BpUGraph& _graph) 
     1.8 -      : anode_matching(_graph), bnode_matching(_graph), graph(&_graph) {}
     1.9 +    MaxBipartiteMatching(const BpUGraph& graph) 
    1.10 +      : _matching(graph), _rmatching(graph), _reached(graph), _graph(&graph) {}
    1.11  
    1.12      /// \name Execution control
    1.13      /// The simplest way to execute the algorithm is to use
    1.14 @@ -87,13 +87,15 @@
    1.15      ///
    1.16      /// It initalizes the data structures and creates an empty matching.
    1.17      void init() {
    1.18 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
    1.19 -        anode_matching[it] = INVALID;
    1.20 +      for (ANodeIt it(*_graph); it != INVALID; ++it) {
    1.21 +        _matching.set(it, INVALID);
    1.22        }
    1.23 -      for (BNodeIt it(*graph); it != INVALID; ++it) {
    1.24 -        bnode_matching[it] = INVALID;
    1.25 +      for (BNodeIt it(*_graph); it != INVALID; ++it) {
    1.26 +        _rmatching.set(it, INVALID);
    1.27 +	_reached.set(it, -1);
    1.28        }
    1.29 -      matching_size = 0;
    1.30 +      _size = 0;
    1.31 +      _phase = -1;
    1.32      }
    1.33  
    1.34      /// \brief Initalize the data structures.
    1.35 @@ -102,21 +104,24 @@
    1.36      /// matching.  From this matching sometimes it is faster to get
    1.37      /// the matching than from the initial empty matching.
    1.38      void greedyInit() {
    1.39 -      matching_size = 0;
    1.40 -      for (BNodeIt it(*graph); it != INVALID; ++it) {
    1.41 -        bnode_matching[it] = INVALID;
    1.42 +      _size = 0;
    1.43 +      for (BNodeIt it(*_graph); it != INVALID; ++it) {
    1.44 +        _rmatching.set(it, INVALID);
    1.45 +	_reached.set(it, 0);
    1.46        }
    1.47 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
    1.48 -        anode_matching[it] = INVALID;
    1.49 -        for (IncEdgeIt jt(*graph, it); jt != INVALID; ++jt) {
    1.50 -          if (bnode_matching[graph->bNode(jt)] == INVALID) {
    1.51 -            anode_matching[it] = jt;
    1.52 -            bnode_matching[graph->bNode(jt)] = jt;
    1.53 -            ++matching_size;
    1.54 +      for (ANodeIt it(*_graph); it != INVALID; ++it) {
    1.55 +        _matching[it] = INVALID;
    1.56 +        for (IncEdgeIt jt(*_graph, it); jt != INVALID; ++jt) {
    1.57 +          if (_rmatching[_graph->bNode(jt)] == INVALID) {
    1.58 +            _matching.set(it, jt);
    1.59 +	    _rmatching.set(_graph->bNode(jt), jt);
    1.60 +	    _reached.set(it, -1);
    1.61 +            ++_size;
    1.62              break;
    1.63            }
    1.64          }
    1.65        }
    1.66 +      _phase = 0;
    1.67      }
    1.68  
    1.69      /// \brief Initalize the data structures with an initial matching.
    1.70 @@ -124,20 +129,23 @@
    1.71      /// It initalizes the data structures with an initial matching.
    1.72      template <typename MatchingMap>
    1.73      void matchingInit(const MatchingMap& mm) {
    1.74 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
    1.75 -        anode_matching[it] = INVALID;
    1.76 +      for (ANodeIt it(*_graph); it != INVALID; ++it) {
    1.77 +        _matching.set(it, INVALID);
    1.78        }
    1.79 -      for (BNodeIt it(*graph); it != INVALID; ++it) {
    1.80 -        bnode_matching[it] = INVALID;
    1.81 +      for (BNodeIt it(*_graph); it != INVALID; ++it) {
    1.82 +        _rmatching.set(it, INVALID);
    1.83 +	_reached.set(it, 0);
    1.84        }
    1.85 -      matching_size = 0;
    1.86 -      for (UEdgeIt it(*graph); it != INVALID; ++it) {
    1.87 +      _size = 0;
    1.88 +      for (UEdgeIt it(*_graph); it != INVALID; ++it) {
    1.89          if (mm[it]) {
    1.90 -          ++matching_size;
    1.91 -          anode_matching[graph->aNode(it)] = it;
    1.92 -          bnode_matching[graph->bNode(it)] = it;
    1.93 +          ++_size;
    1.94 +          _matching.set(_graph->aNode(it), it);
    1.95 +          _rmatching.set(_graph->bNode(it), it);
    1.96 +	  _reached.set(it, 0);
    1.97          }
    1.98        }
    1.99 +      _phase = 0;
   1.100      }
   1.101  
   1.102      /// \brief Initalize the data structures with an initial matching.
   1.103 @@ -145,185 +153,188 @@
   1.104      /// It initalizes the data structures with an initial matching.
   1.105      /// \return %True when the given map contains really a matching.
   1.106      template <typename MatchingMap>
   1.107 -    void checkedMatchingInit(const MatchingMap& mm) {
   1.108 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
   1.109 -        anode_matching[it] = INVALID;
   1.110 +    bool checkedMatchingInit(const MatchingMap& mm) {
   1.111 +      for (ANodeIt it(*_graph); it != INVALID; ++it) {
   1.112 +        _matching.set(it, INVALID);
   1.113        }
   1.114 -      for (BNodeIt it(*graph); it != INVALID; ++it) {
   1.115 -        bnode_matching[it] = INVALID;
   1.116 +      for (BNodeIt it(*_graph); it != INVALID; ++it) {
   1.117 +        _rmatching.set(it, INVALID);
   1.118 +	_reached.set(it, 0);
   1.119        }
   1.120 -      matching_size = 0;
   1.121 -      for (UEdgeIt it(*graph); it != INVALID; ++it) {
   1.122 +      _size = 0;
   1.123 +      for (UEdgeIt it(*_graph); it != INVALID; ++it) {
   1.124          if (mm[it]) {
   1.125 -          ++matching_size;
   1.126 -          if (anode_matching[graph->aNode(it)] != INVALID) {
   1.127 +          ++_size;
   1.128 +          if (_matching[_graph->aNode(it)] != INVALID) {
   1.129              return false;
   1.130            }
   1.131 -          anode_matching[graph->aNode(it)] = it;
   1.132 -          if (bnode_matching[graph->aNode(it)] != INVALID) {
   1.133 +          _matching.set(_graph->aNode(it), it);
   1.134 +          if (_matching[_graph->bNode(it)] != INVALID) {
   1.135              return false;
   1.136            }
   1.137 -          bnode_matching[graph->bNode(it)] = it;
   1.138 +          _matching.set(_graph->bNode(it), it);
   1.139 +	  _reached.set(_graph->bNode(it), -1);
   1.140          }
   1.141        }
   1.142 +      _phase = 0;
   1.143 +      return true;
   1.144 +    }
   1.145 +
   1.146 +  private:
   1.147 +    
   1.148 +    bool _find_path(Node anode, int maxlevel,
   1.149 +		    typename Graph::template BNodeMap<int>& level) {
   1.150 +      for (IncEdgeIt it(*_graph, anode); it != INVALID; ++it) {
   1.151 +	Node bnode = _graph->bNode(it); 
   1.152 +	if (level[bnode] == maxlevel) {
   1.153 +	  level.set(bnode, -1);
   1.154 +	  if (maxlevel == 0) {
   1.155 +	    _matching.set(anode, it);
   1.156 +	    _rmatching.set(bnode, it);
   1.157 +	    return true;
   1.158 +	  } else {
   1.159 +	    Node nnode = _graph->aNode(_rmatching[bnode]);
   1.160 +	    if (_find_path(nnode, maxlevel - 1, level)) {
   1.161 +	      _matching.set(anode, it);
   1.162 +	      _rmatching.set(bnode, it);
   1.163 +	      return true;
   1.164 +	    }
   1.165 +	  }
   1.166 +	}
   1.167 +      }
   1.168        return false;
   1.169      }
   1.170  
   1.171 +  public:
   1.172 +
   1.173      /// \brief An augmenting phase of the Hopcroft-Karp algorithm
   1.174      ///
   1.175      /// It runs an augmenting phase of the Hopcroft-Karp
   1.176 -    /// algorithm. This phase finds maximum count of edge disjoint
   1.177 -    /// augmenting paths and augments on these paths. The algorithm
   1.178 -    /// consists at most of \f$ O(\sqrt{n}) \f$ phase and one phase is
   1.179 -    /// \f$ O(e) \f$ long.
   1.180 +    /// algorithm. This phase finds maximal edge disjoint augmenting
   1.181 +    /// paths and augments on these paths. The algorithm consists at
   1.182 +    /// most of \f$ O(\sqrt{n}) \f$ phase and one phase is \f$ O(e)
   1.183 +    /// \f$ long.
   1.184      bool augment() {
   1.185  
   1.186 -      typename Graph::template ANodeMap<bool> areached(*graph, false);
   1.187 -      typename Graph::template BNodeMap<bool> breached(*graph, false);
   1.188 +      ++_phase;
   1.189 +      
   1.190 +      typename Graph::template BNodeMap<int> _level(*_graph, -1);
   1.191 +      typename Graph::template ANodeMap<bool> _found(*_graph, false);
   1.192  
   1.193 -      typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
   1.194 -
   1.195 -      std::vector<Node> queue, bqueue;
   1.196 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
   1.197 -        if (anode_matching[it] == INVALID) {
   1.198 +      std::vector<Node> queue, aqueue;
   1.199 +      for (BNodeIt it(*_graph); it != INVALID; ++it) {
   1.200 +        if (_rmatching[it] == INVALID) {
   1.201            queue.push_back(it);
   1.202 -          areached[it] = true;
   1.203 +          _reached.set(it, _phase);
   1.204 +	  _level.set(it, 0);
   1.205          }
   1.206        }
   1.207  
   1.208        bool success = false;
   1.209  
   1.210 +      int level = 0;
   1.211        while (!success && !queue.empty()) {
   1.212 -        std::vector<Node> newqueue;
   1.213 +        std::vector<Node> nqueue;
   1.214          for (int i = 0; i < int(queue.size()); ++i) {
   1.215 -          Node anode = queue[i];
   1.216 -          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
   1.217 -            Node bnode = graph->bNode(jt);
   1.218 -            if (breached[bnode]) continue;
   1.219 -            breached[bnode] = true;
   1.220 -            bpred[bnode] = jt;
   1.221 -            if (bnode_matching[bnode] == INVALID) {
   1.222 -              bqueue.push_back(bnode);
   1.223 +          Node bnode = queue[i];
   1.224 +          for (IncEdgeIt jt(*_graph, bnode); jt != INVALID; ++jt) {
   1.225 +            Node anode = _graph->aNode(jt);
   1.226 +            if (_matching[anode] == INVALID) {
   1.227 +
   1.228 +	      if (!_found[anode]) {
   1.229 +		if (_find_path(anode, level, _level)) {
   1.230 +		  ++_size;
   1.231 +		}
   1.232 +		_found.set(anode, true);
   1.233 +	      }
   1.234                success = true;
   1.235              } else {           
   1.236 -              Node newanode = graph->aNode(bnode_matching[bnode]);
   1.237 -              if (!areached[newanode]) {
   1.238 -                areached[newanode] = true;
   1.239 -                newqueue.push_back(newanode);
   1.240 +              Node nnode = _graph->bNode(_matching[anode]);
   1.241 +              if (_reached[nnode] != _phase) {
   1.242 +                _reached.set(nnode, _phase);
   1.243 +                nqueue.push_back(nnode);
   1.244 +		_level.set(nnode, level + 1);
   1.245                }
   1.246              }
   1.247            }
   1.248          }
   1.249 -        queue.swap(newqueue);
   1.250 +	++level;
   1.251 +        queue.swap(nqueue);
   1.252        }
   1.253 -
   1.254 -      if (success) {
   1.255 -
   1.256 -        typename Graph::template ANodeMap<bool> aused(*graph, false);
   1.257 -        
   1.258 -        for (int i = 0; i < int(bqueue.size()); ++i) {
   1.259 -          Node bnode = bqueue[i];
   1.260 -
   1.261 -          bool used = false;
   1.262 -
   1.263 -          while (bnode != INVALID) {
   1.264 -            UEdge uedge = bpred[bnode];
   1.265 -            Node anode = graph->aNode(uedge);
   1.266 -            
   1.267 -            if (aused[anode]) {
   1.268 -              used = true;
   1.269 -              break;
   1.270 -            }
   1.271 -            
   1.272 -            bnode = anode_matching[anode] != INVALID ? 
   1.273 -              graph->bNode(anode_matching[anode]) : INVALID;
   1.274 -            
   1.275 -          }
   1.276 -          
   1.277 -          if (used) continue;
   1.278 -
   1.279 -          bnode = bqueue[i];
   1.280 -          while (bnode != INVALID) {
   1.281 -            UEdge uedge = bpred[bnode];
   1.282 -            Node anode = graph->aNode(uedge);
   1.283 -            
   1.284 -            bnode_matching[bnode] = uedge;
   1.285 -
   1.286 -            bnode = anode_matching[anode] != INVALID ? 
   1.287 -              graph->bNode(anode_matching[anode]) : INVALID;
   1.288 -            
   1.289 -            anode_matching[anode] = uedge;
   1.290 -
   1.291 -            aused[anode] = true;
   1.292 -          }
   1.293 -          ++matching_size;
   1.294 -
   1.295 -        }
   1.296 -      }
   1.297 +      
   1.298        return success;
   1.299      }
   1.300 +  private:
   1.301 +    
   1.302 +    void _find_path_bfs(Node anode,
   1.303 +			typename Graph::template ANodeMap<UEdge>& pred) {
   1.304 +      while (true) {
   1.305 +	UEdge uedge = pred[anode];
   1.306 +	Node bnode = _graph->bNode(uedge);
   1.307  
   1.308 -    /// \brief An augmenting phase of the Ford-Fulkerson algorithm
   1.309 +	UEdge nedge = _rmatching[bnode];
   1.310 +
   1.311 +	_matching.set(anode, uedge);
   1.312 +	_rmatching.set(bnode, uedge);
   1.313 +
   1.314 +	if (nedge == INVALID) break;
   1.315 +	anode = _graph->aNode(nedge);
   1.316 +      }
   1.317 +    }
   1.318 +
   1.319 +  public:
   1.320 +
   1.321 +    /// \brief An augmenting phase with single path augementing
   1.322      ///
   1.323 -    /// It runs an augmenting phase of the Ford-Fulkerson
   1.324 -    /// algorithm. This phase finds only one augmenting path and 
   1.325 -    /// augments only on this paths. The algorithm consists at most 
   1.326 -    /// of \f$ O(n) \f$ simple phase and one phase is at most 
   1.327 -    /// \f$ O(e) \f$ long.
   1.328 -    bool simpleAugment() {
   1.329 +    /// This phase finds only one augmenting paths and augments on
   1.330 +    /// these paths. The algorithm consists at most of \f$ O(n) \f$
   1.331 +    /// phase and one phase is \f$ O(e) \f$ long.
   1.332 +    bool simpleAugment() { 
   1.333 +      ++_phase;
   1.334 +      
   1.335 +      typename Graph::template ANodeMap<UEdge> _pred(*_graph);
   1.336  
   1.337 -      typename Graph::template ANodeMap<bool> areached(*graph, false);
   1.338 -      typename Graph::template BNodeMap<bool> breached(*graph, false);
   1.339 -
   1.340 -      typename Graph::template BNodeMap<UEdge> bpred(*graph, INVALID);
   1.341 -
   1.342 -      std::vector<Node> queue;
   1.343 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
   1.344 -        if (anode_matching[it] == INVALID) {
   1.345 +      std::vector<Node> queue, aqueue;
   1.346 +      for (BNodeIt it(*_graph); it != INVALID; ++it) {
   1.347 +        if (_rmatching[it] == INVALID) {
   1.348            queue.push_back(it);
   1.349 -          areached[it] = true;
   1.350 +          _reached.set(it, _phase);
   1.351          }
   1.352        }
   1.353  
   1.354 -      while (!queue.empty()) {
   1.355 -        std::vector<Node> newqueue;
   1.356 +      bool success = false;
   1.357 +
   1.358 +      int level = 0;
   1.359 +      while (!success && !queue.empty()) {
   1.360 +        std::vector<Node> nqueue;
   1.361          for (int i = 0; i < int(queue.size()); ++i) {
   1.362 -          Node anode = queue[i];
   1.363 -          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
   1.364 -            Node bnode = graph->bNode(jt);
   1.365 -            if (breached[bnode]) continue;
   1.366 -            breached[bnode] = true;
   1.367 -            bpred[bnode] = jt;
   1.368 -            if (bnode_matching[bnode] == INVALID) {
   1.369 -              while (bnode != INVALID) {
   1.370 -                UEdge uedge = bpred[bnode];
   1.371 -                anode = graph->aNode(uedge);
   1.372 -                
   1.373 -                bnode_matching[bnode] = uedge;
   1.374 -                
   1.375 -                bnode = anode_matching[anode] != INVALID ? 
   1.376 -                  graph->bNode(anode_matching[anode]) : INVALID;
   1.377 -                
   1.378 -                anode_matching[anode] = uedge;
   1.379 -                
   1.380 -              }
   1.381 -              ++matching_size;
   1.382 -              return true;
   1.383 +          Node bnode = queue[i];
   1.384 +          for (IncEdgeIt jt(*_graph, bnode); jt != INVALID; ++jt) {
   1.385 +            Node anode = _graph->aNode(jt);
   1.386 +            if (_matching[anode] == INVALID) {
   1.387 +	      _pred.set(anode, jt);
   1.388 +	      _find_path_bfs(anode, _pred);
   1.389 +	      ++_size;
   1.390 +	      return true;
   1.391              } else {           
   1.392 -              Node newanode = graph->aNode(bnode_matching[bnode]);
   1.393 -              if (!areached[newanode]) {
   1.394 -                areached[newanode] = true;
   1.395 -                newqueue.push_back(newanode);
   1.396 +              Node nnode = _graph->bNode(_matching[anode]);
   1.397 +              if (_reached[nnode] != _phase) {
   1.398 +		_pred.set(anode, jt);
   1.399 +		_reached.set(nnode, _phase);
   1.400 +                nqueue.push_back(nnode);
   1.401                }
   1.402              }
   1.403            }
   1.404          }
   1.405 -        queue.swap(newqueue);
   1.406 +	++level;
   1.407 +        queue.swap(nqueue);
   1.408        }
   1.409        
   1.410 -      return false;
   1.411 +      return success;
   1.412      }
   1.413  
   1.414 +
   1.415 +
   1.416      /// \brief Starts the algorithm.
   1.417      ///
   1.418      /// Starts the algorithm. It runs augmenting phases until the optimal
   1.419 @@ -350,6 +361,27 @@
   1.420      
   1.421      ///@{
   1.422  
   1.423 +    /// \brief Return true if the given uedge is in the matching.
   1.424 +    /// 
   1.425 +    /// It returns true if the given uedge is in the matching.
   1.426 +    bool matchingEdge(const UEdge& edge) const {
   1.427 +      return _matching[_graph->aNode(edge)] == edge;
   1.428 +    }
   1.429 +
   1.430 +    /// \brief Returns the matching edge from the node.
   1.431 +    /// 
   1.432 +    /// Returns the matching edge from the node. If there is not such
   1.433 +    /// edge it gives back \c INVALID.
   1.434 +    /// \note If the parameter node is a B-node then the running time is
   1.435 +    /// propotional to the degree of the node.
   1.436 +    UEdge matchingEdge(const Node& node) const {
   1.437 +      if (_graph->aNode(node)) {
   1.438 +        return _matching[node];
   1.439 +      } else {
   1.440 +        return _rmatching[node];
   1.441 +      }
   1.442 +    }
   1.443 +
   1.444      /// \brief Set true all matching uedge in the map.
   1.445      /// 
   1.446      /// Set true all matching uedge in the map. It does not change the
   1.447 @@ -357,12 +389,12 @@
   1.448      /// \return The number of the matching edges.
   1.449      template <typename MatchingMap>
   1.450      int quickMatching(MatchingMap& mm) const {
   1.451 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
   1.452 -        if (anode_matching[it] != INVALID) {
   1.453 -          mm.set(anode_matching[it], true);
   1.454 +      for (ANodeIt it(*_graph); it != INVALID; ++it) {
   1.455 +        if (_matching[it] != INVALID) {
   1.456 +          mm.set(_matching[it], true);
   1.457          }
   1.458        }
   1.459 -      return matching_size;
   1.460 +      return _size;
   1.461      }
   1.462  
   1.463      /// \brief Set true all matching uedge in the map and the others to false.
   1.464 @@ -371,10 +403,10 @@
   1.465      /// \return The number of the matching edges.
   1.466      template <typename MatchingMap>
   1.467      int matching(MatchingMap& mm) const {
   1.468 -      for (UEdgeIt it(*graph); it != INVALID; ++it) {
   1.469 -        mm.set(it, it == anode_matching[graph->aNode(it)]);
   1.470 +      for (UEdgeIt it(*_graph); it != INVALID; ++it) {
   1.471 +        mm.set(it, it == _matching[_graph->aNode(it)]);
   1.472        }
   1.473 -      return matching_size;
   1.474 +      return _size;
   1.475      }
   1.476  
   1.477      ///Gives back the matching in an ANodeMap.
   1.478 @@ -384,10 +416,10 @@
   1.479      ///\return The number of the matching edges.
   1.480      template<class MatchingMap>
   1.481      int aMatching(MatchingMap& mm) const {
   1.482 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
   1.483 -        mm.set(it, anode_matching[it]);
   1.484 +      for (ANodeIt it(*_graph); it != INVALID; ++it) {
   1.485 +        mm.set(it, _matching[it]);
   1.486        }
   1.487 -      return matching_size;
   1.488 +      return _size;
   1.489      }
   1.490  
   1.491      ///Gives back the matching in a BNodeMap.
   1.492 @@ -397,36 +429,10 @@
   1.493      ///\return The number of the matching edges.
   1.494      template<class MatchingMap>
   1.495      int bMatching(MatchingMap& mm) const {
   1.496 -      for (BNodeIt it(*graph); it != INVALID; ++it) {
   1.497 -        mm.set(it, bnode_matching[it]);
   1.498 +      for (BNodeIt it(*_graph); it != INVALID; ++it) {
   1.499 +        mm.set(it, _rmatching[it]);
   1.500        }
   1.501 -      return matching_size;
   1.502 -    }
   1.503 -
   1.504 -    /// \brief Return true if the given uedge is in the matching.
   1.505 -    /// 
   1.506 -    /// It returns true if the given uedge is in the matching.
   1.507 -    bool matchingEdge(const UEdge& edge) const {
   1.508 -      return anode_matching[graph->aNode(edge)] == edge;
   1.509 -    }
   1.510 -
   1.511 -    /// \brief Returns the matching edge from the node.
   1.512 -    /// 
   1.513 -    /// Returns the matching edge from the node. If there is not such
   1.514 -    /// edge it gives back \c INVALID.
   1.515 -    UEdge matchingEdge(const Node& node) const {
   1.516 -      if (graph->aNode(node)) {
   1.517 -        return anode_matching[node];
   1.518 -      } else {
   1.519 -        return bnode_matching[node];
   1.520 -      }
   1.521 -    }
   1.522 -
   1.523 -    /// \brief Gives back the number of the matching edges.
   1.524 -    ///
   1.525 -    /// Gives back the number of the matching edges.
   1.526 -    int matchingSize() const {
   1.527 -      return matching_size;
   1.528 +      return _size;
   1.529      }
   1.530  
   1.531      /// \brief Returns a minimum covering of the nodes.
   1.532 @@ -438,54 +444,23 @@
   1.533      template <typename CoverMap>
   1.534      int coverSet(CoverMap& covering) const {
   1.535  
   1.536 -      typename Graph::template ANodeMap<bool> areached(*graph, false);
   1.537 -      typename Graph::template BNodeMap<bool> breached(*graph, false);
   1.538 -      
   1.539 -      std::vector<Node> queue;
   1.540 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
   1.541 -        if (anode_matching[it] == INVALID) {
   1.542 -          queue.push_back(it);
   1.543 -        }
   1.544 +      int size = 0;
   1.545 +      for (ANodeIt it(*_graph); it != INVALID; ++it) {
   1.546 +	bool cn = _matching[it] != INVALID && 
   1.547 +	  _reached[_graph->bNode(_matching[it])] == _phase;
   1.548 +        covering.set(it, cn);
   1.549 +        if (cn) ++size;
   1.550        }
   1.551 -
   1.552 -      while (!queue.empty()) {
   1.553 -        std::vector<Node> newqueue;
   1.554 -        for (int i = 0; i < int(queue.size()); ++i) {
   1.555 -          Node anode = queue[i];
   1.556 -          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
   1.557 -            Node bnode = graph->bNode(jt);
   1.558 -            if (breached[bnode]) continue;
   1.559 -            breached[bnode] = true;
   1.560 -            if (bnode_matching[bnode] != INVALID) {
   1.561 -              Node newanode = graph->aNode(bnode_matching[bnode]);
   1.562 -              if (!areached[newanode]) {
   1.563 -                areached[newanode] = true;
   1.564 -                newqueue.push_back(newanode);
   1.565 -              }
   1.566 -            }
   1.567 -          }
   1.568 -        }
   1.569 -        queue.swap(newqueue);
   1.570 -      }
   1.571 -
   1.572 -      int size = 0;
   1.573 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
   1.574 -        covering.set(it, !areached[it] && anode_matching[it] != INVALID);
   1.575 -        if (!areached[it] && anode_matching[it] != INVALID) {
   1.576 -          ++size;
   1.577 -        }
   1.578 -      }
   1.579 -      for (BNodeIt it(*graph); it != INVALID; ++it) {
   1.580 -        covering.set(it, breached[it]);
   1.581 -        if (breached[it]) {
   1.582 -          ++size;
   1.583 -        }
   1.584 +      for (BNodeIt it(*_graph); it != INVALID; ++it) {
   1.585 +	bool cn = _reached[it] != _phase;
   1.586 +        covering.set(it, cn);
   1.587 +        if (cn) ++size;
   1.588        }
   1.589        return size;
   1.590      }
   1.591  
   1.592      /// \brief Gives back a barrier on the A-nodes
   1.593 -    
   1.594 +    ///    
   1.595      /// The barrier is s subset of the nodes on the same side of the
   1.596      /// graph, which size minus its neighbours is exactly the
   1.597      /// unmatched nodes on the A-side.  
   1.598 @@ -493,43 +468,14 @@
   1.599      template <typename BarrierMap>
   1.600      void aBarrier(BarrierMap& barrier) const {
   1.601  
   1.602 -      typename Graph::template ANodeMap<bool> areached(*graph, false);
   1.603 -      typename Graph::template BNodeMap<bool> breached(*graph, false);
   1.604 -      
   1.605 -      std::vector<Node> queue;
   1.606 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
   1.607 -        if (anode_matching[it] == INVALID) {
   1.608 -          queue.push_back(it);
   1.609 -        }
   1.610 -      }
   1.611 -
   1.612 -      while (!queue.empty()) {
   1.613 -        std::vector<Node> newqueue;
   1.614 -        for (int i = 0; i < int(queue.size()); ++i) {
   1.615 -          Node anode = queue[i];
   1.616 -          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
   1.617 -            Node bnode = graph->bNode(jt);
   1.618 -            if (breached[bnode]) continue;
   1.619 -            breached[bnode] = true;
   1.620 -            if (bnode_matching[bnode] != INVALID) {
   1.621 -              Node newanode = graph->aNode(bnode_matching[bnode]);
   1.622 -              if (!areached[newanode]) {
   1.623 -                areached[newanode] = true;
   1.624 -                newqueue.push_back(newanode);
   1.625 -              }
   1.626 -            }
   1.627 -          }
   1.628 -        }
   1.629 -        queue.swap(newqueue);
   1.630 -      }
   1.631 -
   1.632 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
   1.633 -        barrier.set(it, areached[it] || anode_matching[it] == INVALID);
   1.634 +      for (ANodeIt it(*_graph); it != INVALID; ++it) {
   1.635 +        barrier.set(it, _matching[it] == INVALID || 
   1.636 +		    _reached[_graph->bNode(_matching[it])] != _phase);
   1.637        }
   1.638      }
   1.639  
   1.640      /// \brief Gives back a barrier on the B-nodes
   1.641 -    
   1.642 +    ///    
   1.643      /// The barrier is s subset of the nodes on the same side of the
   1.644      /// graph, which size minus its neighbours is exactly the
   1.645      /// unmatched nodes on the B-side.  
   1.646 @@ -537,50 +483,31 @@
   1.647      template <typename BarrierMap>
   1.648      void bBarrier(BarrierMap& barrier) const {
   1.649  
   1.650 -      typename Graph::template ANodeMap<bool> areached(*graph, false);
   1.651 -      typename Graph::template BNodeMap<bool> breached(*graph, false);
   1.652 -      
   1.653 -      std::vector<Node> queue;
   1.654 -      for (ANodeIt it(*graph); it != INVALID; ++it) {
   1.655 -        if (anode_matching[it] == INVALID) {
   1.656 -          queue.push_back(it);
   1.657 -        }
   1.658 +      for (BNodeIt it(*_graph); it != INVALID; ++it) {
   1.659 +        barrier.set(it, _reached[it] == _phase);
   1.660        }
   1.661 +    }
   1.662  
   1.663 -      while (!queue.empty()) {
   1.664 -        std::vector<Node> newqueue;
   1.665 -        for (int i = 0; i < int(queue.size()); ++i) {
   1.666 -          Node anode = queue[i];
   1.667 -          for (IncEdgeIt jt(*graph, anode); jt != INVALID; ++jt) {
   1.668 -            Node bnode = graph->bNode(jt);
   1.669 -            if (breached[bnode]) continue;
   1.670 -            breached[bnode] = true;
   1.671 -            if (bnode_matching[bnode] != INVALID) {
   1.672 -              Node newanode = graph->aNode(bnode_matching[bnode]);
   1.673 -              if (!areached[newanode]) {
   1.674 -                areached[newanode] = true;
   1.675 -                newqueue.push_back(newanode);
   1.676 -              }
   1.677 -            }
   1.678 -          }
   1.679 -        }
   1.680 -        queue.swap(newqueue);
   1.681 -      }
   1.682 -
   1.683 -      for (BNodeIt it(*graph); it != INVALID; ++it) {
   1.684 -        barrier.set(it, !breached[it]);
   1.685 -      }
   1.686 +    /// \brief Gives back the number of the matching edges.
   1.687 +    ///
   1.688 +    /// Gives back the number of the matching edges.
   1.689 +    int matchingSize() const {
   1.690 +      return _size;
   1.691      }
   1.692  
   1.693      /// @}
   1.694  
   1.695    private:
   1.696  
   1.697 -    ANodeMatchingMap anode_matching;
   1.698 -    BNodeMatchingMap bnode_matching;
   1.699 -    const Graph *graph;
   1.700 +    typename BpUGraph::template ANodeMap<UEdge> _matching;
   1.701 +    typename BpUGraph::template BNodeMap<UEdge> _rmatching;
   1.702  
   1.703 -    int matching_size;
   1.704 +    typename BpUGraph::template BNodeMap<int> _reached;
   1.705 +
   1.706 +    int _phase;
   1.707 +    const Graph *_graph;
   1.708 +
   1.709 +    int _size;
   1.710    
   1.711    };
   1.712