diff -r f393a8162688 -r a3ba22ebccc6 lemon/max_matching.h
--- a/lemon/max_matching.h	Thu Dec 27 13:40:16 2007 +0000
+++ b/lemon/max_matching.h	Fri Dec 28 11:00:51 2007 +0000
@@ -19,10 +19,14 @@
 #ifndef LEMON_MAX_MATCHING_H
 #define LEMON_MAX_MATCHING_H
 
+#include <vector>
 #include <queue>
+#include <set>
+
 #include <lemon/bits/invalid.h>
 #include <lemon/unionfind.h>
 #include <lemon/graph_utils.h>
+#include <lemon/bin_heap.h>
 
 ///\ingroup matching
 ///\file
@@ -31,7 +35,7 @@
 namespace lemon {
 
   ///\ingroup matching
-
+  ///
   ///\brief Edmonds' alternating forest maximum matching algorithm.
   ///
   ///This class provides Edmonds' alternating forest matching
@@ -618,6 +622,2500 @@
     }
 
   };
+
+  /// \ingroup matching
+  ///
+  /// \brief Weighted matching in general undirected graphs
+  ///
+  /// This class provides an efficient implementation of Edmond's
+  /// maximum weighted matching algorithm. The implementation is based
+  /// on extensive use of priority queues and provides
+  /// \f$O(nm\log(n))\f$ time complexity.
+  ///
+  /// The maximum weighted matching problem is to find undirected
+  /// edges in the graph with maximum overall weight and no two of
+  /// them shares their endpoints. The problem can be formulated with
+  /// the next linear program: 
+  /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
+  ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f] 
+  /// \f[x_e \ge 0\quad \forall e\in E\f]
+  /// \f[\max \sum_{e\in E}x_ew_e\f]
+  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
+  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in
+  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
+  /// the nodes.
+  ///
+  /// The algorithm calculates an optimal matching and a proof of the
+  /// optimality. The solution of the dual problem can be used to check
+  /// the result of the algorithm. The dual linear problem is the next:
+  /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
+  /// \f[y_u \ge 0 \quad \forall u \in V\f]
+  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
+  /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
+  ///
+  /// The algorithm can be executed with \c run() or the \c init() and
+  /// then the \c start() member functions. After it the matching can
+  /// be asked with \c matching() or mate() functions. The dual
+  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
+  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
+  /// "BlossomIt" nested class which is able to iterate on the nodes
+  /// of a blossom. If the value type is integral then the dual
+  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
+  ///
+  /// \author Balazs Dezso
+  template <typename _UGraph, 
+	    typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
+  class MaxWeightedMatching {
+  public:
+
+    typedef _UGraph UGraph;
+    typedef _WeightMap WeightMap;
+    typedef typename WeightMap::Value Value;
+
+    /// \brief Scaling factor for dual solution
+    ///
+    /// Scaling factor for dual solution, it is equal to 4 or 1
+    /// according to the value type.
+    static const int dualScale = 
+      std::numeric_limits<Value>::is_integer ? 4 : 1;
+
+    typedef typename UGraph::template NodeMap<typename UGraph::Edge> 
+    MatchingMap;
+    
+  private:
+
+    UGRAPH_TYPEDEFS(typename UGraph);
+
+    typedef typename UGraph::template NodeMap<Value> NodePotential;
+    typedef std::vector<Node> BlossomNodeList;
+
+    struct BlossomVariable {
+      int begin, end;
+      Value value;
+      
+      BlossomVariable(int _begin, int _end, Value _value)
+        : begin(_begin), end(_end), value(_value) {}
+
+    };
+
+    typedef std::vector<BlossomVariable> BlossomPotential;
+
+    const UGraph& _ugraph;
+    const WeightMap& _weight;
+
+    MatchingMap* _matching;
+
+    NodePotential* _node_potential;
+
+    BlossomPotential _blossom_potential;
+    BlossomNodeList _blossom_node_list;
+
+    int _node_num;
+    int _blossom_num;
+
+    typedef typename UGraph::template NodeMap<int> NodeIntMap;
+    typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
+    typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
+    typedef IntegerMap<int> IntIntMap;
+
+    enum Status {
+      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
+    };
+
+    typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
+    struct BlossomData {
+      int tree;
+      Status status;
+      Edge pred, next;
+      Value pot, offset;
+      Node base;
+    };
+
+    NodeIntMap *_blossom_index;
+    BlossomSet *_blossom_set;
+    IntegerMap<BlossomData>* _blossom_data;
+
+    NodeIntMap *_node_index;
+    EdgeIntMap *_node_heap_index;
+
+    struct NodeData {
+      
+      NodeData(EdgeIntMap& node_heap_index) 
+	: heap(node_heap_index) {}
+      
+      int blossom;
+      Value pot;
+      BinHeap<Value, EdgeIntMap> heap;
+      std::map<int, Edge> heap_index;
+      
+      int tree;
+    };
+
+    IntegerMap<NodeData>* _node_data;
+
+    typedef ExtendFindEnum<IntIntMap> TreeSet;
+
+    IntIntMap *_tree_set_index;
+    TreeSet *_tree_set;
+
+    NodeIntMap *_delta1_index;
+    BinHeap<Value, NodeIntMap> *_delta1;
+
+    IntIntMap *_delta2_index;
+    BinHeap<Value, IntIntMap> *_delta2;
+    
+    UEdgeIntMap *_delta3_index;
+    BinHeap<Value, UEdgeIntMap> *_delta3;
+
+    IntIntMap *_delta4_index;
+    BinHeap<Value, IntIntMap> *_delta4;
+
+    Value _delta_sum;
+
+    void createStructures() {
+      _node_num = countNodes(_ugraph); 
+      _blossom_num = _node_num * 3 / 2;
+
+      if (!_matching) {
+	_matching = new MatchingMap(_ugraph);
+      }
+      if (!_node_potential) {
+	_node_potential = new NodePotential(_ugraph);
+      }
+      if (!_blossom_set) {
+	_blossom_index = new NodeIntMap(_ugraph);
+	_blossom_set = new BlossomSet(*_blossom_index);
+	_blossom_data = new IntegerMap<BlossomData>(_blossom_num);
+      }
+
+      if (!_node_index) {
+	_node_index = new NodeIntMap(_ugraph);
+	_node_heap_index = new EdgeIntMap(_ugraph);
+	_node_data = new IntegerMap<NodeData>(_node_num, 
+					      NodeData(*_node_heap_index));
+      }
+
+      if (!_tree_set) {
+	_tree_set_index = new IntIntMap(_blossom_num);
+	_tree_set = new TreeSet(*_tree_set_index);
+      }
+      if (!_delta1) {
+	_delta1_index = new NodeIntMap(_ugraph);
+	_delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
+      }
+      if (!_delta2) {
+	_delta2_index = new IntIntMap(_blossom_num);
+	_delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
+      }
+      if (!_delta3) {
+	_delta3_index = new UEdgeIntMap(_ugraph);
+	_delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
+      }
+      if (!_delta4) {
+	_delta4_index = new IntIntMap(_blossom_num);
+	_delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
+      }
+    }
+
+    void destroyStructures() {
+      _node_num = countNodes(_ugraph); 
+      _blossom_num = _node_num * 3 / 2;
+
+      if (_matching) {
+	delete _matching;
+      }
+      if (_node_potential) {
+	delete _node_potential;
+      }
+      if (_blossom_set) {
+	delete _blossom_index;
+	delete _blossom_set;
+	delete _blossom_data;
+      }
+
+      if (_node_index) {
+	delete _node_index;
+	delete _node_heap_index;
+	delete _node_data;			      
+      }
+
+      if (_tree_set) {
+	delete _tree_set_index;
+	delete _tree_set;
+      }
+      if (_delta1) {
+	delete _delta1_index;
+	delete _delta1;
+      }
+      if (_delta2) {
+	delete _delta2_index;
+	delete _delta2;
+      }
+      if (_delta3) {
+	delete _delta3_index;
+	delete _delta3;
+      }
+      if (_delta4) {
+	delete _delta4_index;
+	delete _delta4;
+      }
+    }
+
+    void matchedToEven(int blossom, int tree) {
+      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+	_delta2->erase(blossom);
+      }
+
+      if (!_blossom_set->trivial(blossom)) {
+	(*_blossom_data)[blossom].pot -= 
+	  2 * (_delta_sum - (*_blossom_data)[blossom].offset);
+      }
+
+      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 
+	   n != INVALID; ++n) {
+
+	_blossom_set->increase(n, std::numeric_limits<Value>::max());
+	int ni = (*_node_index)[n];
+
+	(*_node_data)[ni].heap.clear();
+	(*_node_data)[ni].heap_index.clear();
+
+	(*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
+
+	_delta1->push(n, (*_node_data)[ni].pot);
+
+	for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  Node v = _ugraph.source(e);
+	  int vb = _blossom_set->find(v);
+	  int vi = (*_node_index)[v];
+
+	  Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 
+	    dualScale * _weight[e];
+
+	  if ((*_blossom_data)[vb].status == EVEN) {
+	    if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
+	      _delta3->push(e, rw / 2);
+	    }
+	  } else if ((*_blossom_data)[vb].status == UNMATCHED) {
+	    if (_delta3->state(e) != _delta3->IN_HEAP) {
+	      _delta3->push(e, rw);
+	    }	    
+	  } else {
+	    typename std::map<int, Edge>::iterator it = 
+	      (*_node_data)[vi].heap_index.find(tree);	  
+
+	    if (it != (*_node_data)[vi].heap_index.end()) {
+	      if ((*_node_data)[vi].heap[it->second] > rw) {
+		(*_node_data)[vi].heap.replace(it->second, e);
+		(*_node_data)[vi].heap.decrease(e, rw);
+		it->second = e;
+	      }
+	    } else {
+	      (*_node_data)[vi].heap.push(e, rw);
+	      (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+	    }
+
+	    if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
+	      _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
+
+	      if ((*_blossom_data)[vb].status == MATCHED) {
+		if (_delta2->state(vb) != _delta2->IN_HEAP) {
+		  _delta2->push(vb, _blossom_set->classPrio(vb) -
+			       (*_blossom_data)[vb].offset);
+		} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
+			   (*_blossom_data)[vb].offset){
+		  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
+				   (*_blossom_data)[vb].offset);
+		}
+	      }
+	    }
+	  }
+	}
+      }
+      (*_blossom_data)[blossom].offset = 0;
+    }
+
+    void matchedToOdd(int blossom) {
+      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+	_delta2->erase(blossom);
+      }
+      (*_blossom_data)[blossom].offset += _delta_sum;
+      if (!_blossom_set->trivial(blossom)) {
+	_delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + 
+		     (*_blossom_data)[blossom].offset);
+      }
+    }
+
+    void evenToMatched(int blossom, int tree) {
+      if (!_blossom_set->trivial(blossom)) {
+	(*_blossom_data)[blossom].pot += 2 * _delta_sum;
+      }
+
+      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+	   n != INVALID; ++n) {
+	int ni = (*_node_index)[n];
+	(*_node_data)[ni].pot -= _delta_sum;
+
+	_delta1->erase(n);
+
+	for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  Node v = _ugraph.source(e);
+	  int vb = _blossom_set->find(v);
+	  int vi = (*_node_index)[v];
+
+	  Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 
+	    dualScale * _weight[e];
+
+	  if (vb == blossom) {
+	    if (_delta3->state(e) == _delta3->IN_HEAP) {
+	      _delta3->erase(e);
+	    }
+	  } else if ((*_blossom_data)[vb].status == EVEN) {
+
+	    if (_delta3->state(e) == _delta3->IN_HEAP) {
+	      _delta3->erase(e);
+	    }
+
+	    int vt = _tree_set->find(vb);
+	    
+	    if (vt != tree) {
+
+	      Edge r = _ugraph.oppositeEdge(e);
+
+	      typename std::map<int, Edge>::iterator it = 
+		(*_node_data)[ni].heap_index.find(vt);	  
+
+	      if (it != (*_node_data)[ni].heap_index.end()) {
+		if ((*_node_data)[ni].heap[it->second] > rw) {
+		  (*_node_data)[ni].heap.replace(it->second, r);
+		  (*_node_data)[ni].heap.decrease(r, rw);
+		  it->second = r;
+		}
+	      } else {
+		(*_node_data)[ni].heap.push(r, rw);
+		(*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
+	      }
+	      
+	      if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
+		_blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
+		
+		if (_delta2->state(blossom) != _delta2->IN_HEAP) {
+		  _delta2->push(blossom, _blossom_set->classPrio(blossom) - 
+			       (*_blossom_data)[blossom].offset);
+		} else if ((*_delta2)[blossom] > 
+			   _blossom_set->classPrio(blossom) - 
+			   (*_blossom_data)[blossom].offset){
+		  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
+				   (*_blossom_data)[blossom].offset);
+		}
+	      }
+	    }
+	    
+	  } else if ((*_blossom_data)[vb].status == UNMATCHED) {
+	    if (_delta3->state(e) == _delta3->IN_HEAP) {
+	      _delta3->erase(e);
+	    }
+	  } else {
+	  
+	    typename std::map<int, Edge>::iterator it = 
+	      (*_node_data)[vi].heap_index.find(tree);
+
+	    if (it != (*_node_data)[vi].heap_index.end()) {
+	      (*_node_data)[vi].heap.erase(it->second);
+	      (*_node_data)[vi].heap_index.erase(it);
+	      if ((*_node_data)[vi].heap.empty()) {
+		_blossom_set->increase(v, std::numeric_limits<Value>::max());
+	      } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
+		_blossom_set->increase(v, (*_node_data)[vi].heap.prio());
+	      }
+	      
+	      if ((*_blossom_data)[vb].status == MATCHED) {
+		if (_blossom_set->classPrio(vb) == 
+		    std::numeric_limits<Value>::max()) {
+		  _delta2->erase(vb);
+		} else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
+			   (*_blossom_data)[vb].offset) {
+		  _delta2->increase(vb, _blossom_set->classPrio(vb) -
+				   (*_blossom_data)[vb].offset);
+		}
+	      }	
+	    }	
+	  }
+	}
+      }
+    }
+
+    void oddToMatched(int blossom) {
+      (*_blossom_data)[blossom].offset -= _delta_sum;
+
+      if (_blossom_set->classPrio(blossom) != 
+	  std::numeric_limits<Value>::max()) {
+	_delta2->push(blossom, _blossom_set->classPrio(blossom) - 
+		       (*_blossom_data)[blossom].offset);
+      }
+
+      if (!_blossom_set->trivial(blossom)) {
+	_delta4->erase(blossom);
+      }
+    }
+
+    void oddToEven(int blossom, int tree) {
+      if (!_blossom_set->trivial(blossom)) {
+	_delta4->erase(blossom);
+	(*_blossom_data)[blossom].pot -= 
+	  2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
+      }
+
+      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 
+	   n != INVALID; ++n) {
+	int ni = (*_node_index)[n];
+
+	_blossom_set->increase(n, std::numeric_limits<Value>::max());
+
+	(*_node_data)[ni].heap.clear();
+	(*_node_data)[ni].heap_index.clear();
+	(*_node_data)[ni].pot += 
+	  2 * _delta_sum - (*_blossom_data)[blossom].offset;
+
+	_delta1->push(n, (*_node_data)[ni].pot);
+
+	for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  Node v = _ugraph.source(e);
+	  int vb = _blossom_set->find(v);
+	  int vi = (*_node_index)[v];
+
+	  Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 
+	    dualScale * _weight[e];
+
+	  if ((*_blossom_data)[vb].status == EVEN) {
+	    if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
+	      _delta3->push(e, rw / 2);
+	    }
+	  } else if ((*_blossom_data)[vb].status == UNMATCHED) {
+	    if (_delta3->state(e) != _delta3->IN_HEAP) {
+	      _delta3->push(e, rw);
+	    }
+	  } else {
+	  
+	    typename std::map<int, Edge>::iterator it = 
+	      (*_node_data)[vi].heap_index.find(tree);	  
+
+	    if (it != (*_node_data)[vi].heap_index.end()) {
+	      if ((*_node_data)[vi].heap[it->second] > rw) {
+		(*_node_data)[vi].heap.replace(it->second, e);
+		(*_node_data)[vi].heap.decrease(e, rw);
+		it->second = e;
+	      }
+	    } else {
+	      (*_node_data)[vi].heap.push(e, rw);
+	      (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+	    }
+
+	    if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
+	      _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
+
+	      if ((*_blossom_data)[vb].status == MATCHED) {
+		if (_delta2->state(vb) != _delta2->IN_HEAP) {
+		  _delta2->push(vb, _blossom_set->classPrio(vb) - 
+			       (*_blossom_data)[vb].offset);
+		} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - 
+			   (*_blossom_data)[vb].offset) {
+		  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
+				   (*_blossom_data)[vb].offset);
+		}
+	      }
+	    }
+	  }
+	}
+      }
+      (*_blossom_data)[blossom].offset = 0;
+    }
+
+
+    void matchedToUnmatched(int blossom) {
+      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+	_delta2->erase(blossom);
+      }
+
+      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 
+	   n != INVALID; ++n) {
+	int ni = (*_node_index)[n];
+
+	_blossom_set->increase(n, std::numeric_limits<Value>::max());
+
+	(*_node_data)[ni].heap.clear();
+	(*_node_data)[ni].heap_index.clear();
+
+	for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  Node v = _ugraph.target(e);
+	  int vb = _blossom_set->find(v);
+	  int vi = (*_node_index)[v];
+
+	  Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 
+	    dualScale * _weight[e];
+
+	  if ((*_blossom_data)[vb].status == EVEN) {
+	    if (_delta3->state(e) != _delta3->IN_HEAP) {
+	      _delta3->push(e, rw);
+	    }
+	  }
+	}
+      }
+    }
+
+    void unmatchedToMatched(int blossom) {
+      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+	   n != INVALID; ++n) {
+	int ni = (*_node_index)[n];
+
+	for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  Node v = _ugraph.source(e);
+	  int vb = _blossom_set->find(v);
+	  int vi = (*_node_index)[v];
+
+	  Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 
+	    dualScale * _weight[e];
+
+	  if (vb == blossom) {
+	    if (_delta3->state(e) == _delta3->IN_HEAP) {
+	      _delta3->erase(e);
+	    }
+	  } else if ((*_blossom_data)[vb].status == EVEN) {
+
+	    if (_delta3->state(e) == _delta3->IN_HEAP) {
+	      _delta3->erase(e);
+	    }
+
+	    int vt = _tree_set->find(vb);
+	    
+	    Edge r = _ugraph.oppositeEdge(e);
+	    
+	    typename std::map<int, Edge>::iterator it = 
+	      (*_node_data)[ni].heap_index.find(vt);	  
+	    
+	    if (it != (*_node_data)[ni].heap_index.end()) {
+	      if ((*_node_data)[ni].heap[it->second] > rw) {
+		(*_node_data)[ni].heap.replace(it->second, r);
+		(*_node_data)[ni].heap.decrease(r, rw);
+		it->second = r;
+	      }
+	    } else {
+	      (*_node_data)[ni].heap.push(r, rw);
+	      (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
+	    }
+	      
+	    if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
+	      _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
+	      
+	      if (_delta2->state(blossom) != _delta2->IN_HEAP) {
+		_delta2->push(blossom, _blossom_set->classPrio(blossom) - 
+			     (*_blossom_data)[blossom].offset);
+	      } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
+			 (*_blossom_data)[blossom].offset){
+		_delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
+				 (*_blossom_data)[blossom].offset);
+	      }
+	    }
+	    
+	  } else if ((*_blossom_data)[vb].status == UNMATCHED) {
+	    if (_delta3->state(e) == _delta3->IN_HEAP) {
+	      _delta3->erase(e);
+	    }
+	  }
+	}
+      }
+    }
+
+    void alternatePath(int even, int tree) {
+      int odd;
+
+      evenToMatched(even, tree);
+      (*_blossom_data)[even].status = MATCHED;
+
+      while ((*_blossom_data)[even].pred != INVALID) {
+	odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred));
+	(*_blossom_data)[odd].status = MATCHED;
+	oddToMatched(odd);
+	(*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
+      
+	even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
+	(*_blossom_data)[even].status = MATCHED;
+	evenToMatched(even, tree);
+	(*_blossom_data)[even].next = 
+	  _ugraph.oppositeEdge((*_blossom_data)[odd].pred);
+      }
+      
+    }
+
+    void destroyTree(int tree) {
+      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
+	if ((*_blossom_data)[b].status == EVEN) {
+	  (*_blossom_data)[b].status = MATCHED;
+	  evenToMatched(b, tree);
+	} else if ((*_blossom_data)[b].status == ODD) {
+	  (*_blossom_data)[b].status = MATCHED;
+	  oddToMatched(b);
+	}
+      }
+      _tree_set->eraseClass(tree);
+    }
+    
+
+    void unmatchNode(const Node& node) {
+      int blossom = _blossom_set->find(node);
+      int tree = _tree_set->find(blossom);
+
+      alternatePath(blossom, tree);
+      destroyTree(tree);
+
+      (*_blossom_data)[blossom].status = UNMATCHED;
+      (*_blossom_data)[blossom].base = node;
+      matchedToUnmatched(blossom);
+    }
+
+
+    void augmentOnEdge(const UEdge& edge) {
+      
+      int left = _blossom_set->find(_ugraph.source(edge));
+      int right = _blossom_set->find(_ugraph.target(edge));
+
+      if ((*_blossom_data)[left].status == EVEN) {
+	int left_tree = _tree_set->find(left);
+	alternatePath(left, left_tree);
+	destroyTree(left_tree);
+      } else {
+	(*_blossom_data)[left].status = MATCHED;
+	unmatchedToMatched(left);
+      }
+
+      if ((*_blossom_data)[right].status == EVEN) {
+	int right_tree = _tree_set->find(right);
+	alternatePath(right, right_tree);
+	destroyTree(right_tree);
+      } else {
+	(*_blossom_data)[right].status = MATCHED;
+	unmatchedToMatched(right);
+      }
+
+      (*_blossom_data)[left].next = _ugraph.direct(edge, true);
+      (*_blossom_data)[right].next = _ugraph.direct(edge, false);
+    }
+
+    void extendOnEdge(const Edge& edge) {
+      int base = _blossom_set->find(_ugraph.target(edge));
+      int tree = _tree_set->find(base);
+      
+      int odd = _blossom_set->find(_ugraph.source(edge));
+      _tree_set->insert(odd, tree);
+      (*_blossom_data)[odd].status = ODD;
+      matchedToOdd(odd);
+      (*_blossom_data)[odd].pred = edge;
+
+      int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next));
+      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
+      _tree_set->insert(even, tree);
+      (*_blossom_data)[even].status = EVEN;
+      matchedToEven(even, tree);
+    }
+    
+    void shrinkOnEdge(const UEdge& uedge, int tree) {
+      int nca = -1;
+      std::vector<int> left_path, right_path;
+
+      {
+	std::set<int> left_set, right_set;
+	int left = _blossom_set->find(_ugraph.source(uedge));
+	left_path.push_back(left);
+	left_set.insert(left);
+
+	int right = _blossom_set->find(_ugraph.target(uedge));
+	right_path.push_back(right);
+	right_set.insert(right);
+
+	while (true) {
+
+	  if ((*_blossom_data)[left].pred == INVALID) break;
+
+	  left = 
+	    _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
+	  left_path.push_back(left);
+	  left = 
+	    _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
+	  left_path.push_back(left);
+
+	  left_set.insert(left);
+
+	  if (right_set.find(left) != right_set.end()) {
+	    nca = left;
+	    break;
+	  }
+
+	  if ((*_blossom_data)[right].pred == INVALID) break;
+
+	  right = 
+	    _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
+	  right_path.push_back(right);
+	  right = 
+	    _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
+	  right_path.push_back(right);
+
+	  right_set.insert(right);
+ 
+	  if (left_set.find(right) != left_set.end()) {
+	    nca = right;
+	    break;
+	  }
+
+	}
+
+	if (nca == -1) {
+	  if ((*_blossom_data)[left].pred == INVALID) {
+	    nca = right;
+	    while (left_set.find(nca) == left_set.end()) {
+	      nca = 
+		_blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
+	      right_path.push_back(nca);
+	      nca = 
+		_blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
+	      right_path.push_back(nca);
+	    }
+	  } else {
+	    nca = left;
+	    while (right_set.find(nca) == right_set.end()) {
+	      nca = 
+		_blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
+	      left_path.push_back(nca);
+	      nca = 
+		_blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
+	      left_path.push_back(nca);
+	    }
+	  }
+	}
+      }
+
+      std::vector<int> subblossoms;
+      Edge prev;
+
+      prev = _ugraph.direct(uedge, true);
+      for (int i = 0; left_path[i] != nca; i += 2) {
+	subblossoms.push_back(left_path[i]);
+	(*_blossom_data)[left_path[i]].next = prev;
+	_tree_set->erase(left_path[i]);
+
+	subblossoms.push_back(left_path[i + 1]);
+	(*_blossom_data)[left_path[i + 1]].status = EVEN;
+	oddToEven(left_path[i + 1], tree);
+	_tree_set->erase(left_path[i + 1]);
+	prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred);
+      }
+
+      int k = 0;
+      while (right_path[k] != nca) ++k;
+
+      subblossoms.push_back(nca);
+      (*_blossom_data)[nca].next = prev;
+      
+      for (int i = k - 2; i >= 0; i -= 2) {
+	subblossoms.push_back(right_path[i + 1]);
+	(*_blossom_data)[right_path[i + 1]].status = EVEN;
+	oddToEven(right_path[i + 1], tree);
+	_tree_set->erase(right_path[i + 1]);
+
+	(*_blossom_data)[right_path[i + 1]].next = 
+	  (*_blossom_data)[right_path[i + 1]].pred;
+
+	subblossoms.push_back(right_path[i]);
+	_tree_set->erase(right_path[i]);
+      }
+
+      int surface = 
+	_blossom_set->join(subblossoms.begin(), subblossoms.end());
+
+      for (int i = 0; i < int(subblossoms.size()); ++i) {
+	if (!_blossom_set->trivial(subblossoms[i])) {
+	  (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
+	}
+	(*_blossom_data)[subblossoms[i]].status = MATCHED;
+      }
+
+      (*_blossom_data)[surface].pot = -2 * _delta_sum;
+      (*_blossom_data)[surface].offset = 0;
+      (*_blossom_data)[surface].status = EVEN;
+      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
+      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
+
+      _tree_set->insert(surface, tree);
+      _tree_set->erase(nca);
+    }
+
+    void splitBlossom(int blossom) {
+      Edge next = (*_blossom_data)[blossom].next; 
+      Edge pred = (*_blossom_data)[blossom].pred;
+
+      int tree = _tree_set->find(blossom);
+
+      (*_blossom_data)[blossom].status = MATCHED;
+      oddToMatched(blossom);
+      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+	_delta2->erase(blossom);
+      }
+
+      std::vector<int> subblossoms;
+      _blossom_set->split(blossom, std::back_inserter(subblossoms));
+
+      Value offset = (*_blossom_data)[blossom].offset;
+      int b = _blossom_set->find(_ugraph.source(pred));
+      int d = _blossom_set->find(_ugraph.source(next));
+      
+      int ib, id;
+      for (int i = 0; i < int(subblossoms.size()); ++i) {
+	if (subblossoms[i] == b) ib = i;
+	if (subblossoms[i] == d) id = i;
+
+	(*_blossom_data)[subblossoms[i]].offset = offset;
+	if (!_blossom_set->trivial(subblossoms[i])) {
+	  (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
+	}
+	if (_blossom_set->classPrio(subblossoms[i]) != 
+	    std::numeric_limits<Value>::max()) {
+	  _delta2->push(subblossoms[i], 
+			_blossom_set->classPrio(subblossoms[i]) - 
+			(*_blossom_data)[subblossoms[i]].offset);
+	}
+      }
+      
+      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
+	for (int i = (id + 1) % subblossoms.size(); 
+	     i != ib; i = (i + 2) % subblossoms.size()) {
+	  int sb = subblossoms[i];
+	  int tb = subblossoms[(i + 1) % subblossoms.size()];
+	  (*_blossom_data)[sb].next = 
+	    _ugraph.oppositeEdge((*_blossom_data)[tb].next);
+	}
+
+	for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
+	  int sb = subblossoms[i];
+	  int tb = subblossoms[(i + 1) % subblossoms.size()];
+	  int ub = subblossoms[(i + 2) % subblossoms.size()];
+
+	  (*_blossom_data)[sb].status = ODD;
+	  matchedToOdd(sb);
+	  _tree_set->insert(sb, tree);
+	  (*_blossom_data)[sb].pred = pred;
+	  (*_blossom_data)[sb].next = 
+			   _ugraph.oppositeEdge((*_blossom_data)[tb].next);
+	  
+	  pred = (*_blossom_data)[ub].next;
+      
+	  (*_blossom_data)[tb].status = EVEN;
+	  matchedToEven(tb, tree);
+	  _tree_set->insert(tb, tree);
+	  (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
+	}
+      
+	(*_blossom_data)[subblossoms[id]].status = ODD;
+	matchedToOdd(subblossoms[id]);
+	_tree_set->insert(subblossoms[id], tree);
+	(*_blossom_data)[subblossoms[id]].next = next;
+	(*_blossom_data)[subblossoms[id]].pred = pred;
+      
+      } else {
+
+	for (int i = (ib + 1) % subblossoms.size(); 
+	     i != id; i = (i + 2) % subblossoms.size()) {
+	  int sb = subblossoms[i];
+	  int tb = subblossoms[(i + 1) % subblossoms.size()];
+	  (*_blossom_data)[sb].next = 
+	    _ugraph.oppositeEdge((*_blossom_data)[tb].next);
+	}	
+
+	for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
+	  int sb = subblossoms[i];
+	  int tb = subblossoms[(i + 1) % subblossoms.size()];
+	  int ub = subblossoms[(i + 2) % subblossoms.size()];
+
+	  (*_blossom_data)[sb].status = ODD;
+	  matchedToOdd(sb);
+	  _tree_set->insert(sb, tree);
+	  (*_blossom_data)[sb].next = next; 
+	  (*_blossom_data)[sb].pred =
+	    _ugraph.oppositeEdge((*_blossom_data)[tb].next);
+
+	  (*_blossom_data)[tb].status = EVEN;
+	  matchedToEven(tb, tree);
+	  _tree_set->insert(tb, tree);
+	  (*_blossom_data)[tb].pred =
+	    (*_blossom_data)[tb].next = 
+	    _ugraph.oppositeEdge((*_blossom_data)[ub].next);
+	  next = (*_blossom_data)[ub].next;
+	}
+
+	(*_blossom_data)[subblossoms[ib]].status = ODD;
+	matchedToOdd(subblossoms[ib]);
+	_tree_set->insert(subblossoms[ib], tree);
+	(*_blossom_data)[subblossoms[ib]].next = next;
+	(*_blossom_data)[subblossoms[ib]].pred = pred;
+      }
+      _tree_set->erase(blossom);
+    }
+
+    void extractBlossom(int blossom, const Node& base, const Edge& matching) {
+      if (_blossom_set->trivial(blossom)) {
+	int bi = (*_node_index)[base];
+	Value pot = (*_node_data)[bi].pot;
+
+	_matching->set(base, matching);
+	_blossom_node_list.push_back(base);
+	_node_potential->set(base, pot);
+      } else {
+
+	Value pot = (*_blossom_data)[blossom].pot;
+	int bn = _blossom_node_list.size();
+	
+	std::vector<int> subblossoms;
+	_blossom_set->split(blossom, std::back_inserter(subblossoms));
+	int b = _blossom_set->find(base);
+	int ib = -1;
+	for (int i = 0; i < int(subblossoms.size()); ++i) {
+	  if (subblossoms[i] == b) { ib = i; break; }
+	}
+	
+	for (int i = 1; i < int(subblossoms.size()); i += 2) {
+	  int sb = subblossoms[(ib + i) % subblossoms.size()];
+	  int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
+	  
+	  Edge m = (*_blossom_data)[tb].next;
+	  extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
+	  extractBlossom(tb, _ugraph.source(m), m);
+	}
+	extractBlossom(subblossoms[ib], base, matching);      
+	
+	int en = _blossom_node_list.size();
+	
+	_blossom_potential.push_back(BlossomVariable(bn, en, pot));
+      }
+    }
+
+    void extractMatching() {
+      std::vector<int> blossoms;
+      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
+	blossoms.push_back(c);
+      }
+
+      for (int i = 0; i < int(blossoms.size()); ++i) {
+	if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
+
+	  Value offset = (*_blossom_data)[blossoms[i]].offset;
+	  (*_blossom_data)[blossoms[i]].pot += 2 * offset;
+	  for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); 
+	       n != INVALID; ++n) {
+	    (*_node_data)[(*_node_index)[n]].pot -= offset;
+	  }
+
+	  Edge matching = (*_blossom_data)[blossoms[i]].next;
+	  Node base = _ugraph.source(matching);
+	  extractBlossom(blossoms[i], base, matching);
+	} else {
+	  Node base = (*_blossom_data)[blossoms[i]].base;	  
+	  extractBlossom(blossoms[i], base, INVALID);
+	}
+      }
+    }
+    
+  public:
+
+    /// \brief Constructor
+    ///
+    /// Constructor.
+    MaxWeightedMatching(const UGraph& ugraph, const WeightMap& weight)
+      : _ugraph(ugraph), _weight(weight), _matching(0),
+	_node_potential(0), _blossom_potential(), _blossom_node_list(),
+	_node_num(0), _blossom_num(0),
+
+	_blossom_index(0), _blossom_set(0), _blossom_data(0),
+	_node_index(0), _node_heap_index(0), _node_data(0),
+	_tree_set_index(0), _tree_set(0),
+
+	_delta1_index(0), _delta1(0),
+	_delta2_index(0), _delta2(0),
+	_delta3_index(0), _delta3(0), 
+	_delta4_index(0), _delta4(0),
+
+	_delta_sum() {}
+
+    ~MaxWeightedMatching() {
+      destroyStructures();
+    }
+
+    /// \name Execution control 
+    /// The simplest way to execute the algorithm is to use the member
+    /// \c run() member function.
+
+    ///@{
+
+    /// \brief Initialize the algorithm
+    ///
+    /// Initialize the algorithm
+    void init() {
+      createStructures();
+
+      for (EdgeIt e(_ugraph); e != INVALID; ++e) {
+	_node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
+      }
+      for (NodeIt n(_ugraph); n != INVALID; ++n) {
+	_delta1_index->set(n, _delta1->PRE_HEAP);
+      }
+      for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
+	_delta3_index->set(e, _delta3->PRE_HEAP);
+      }
+      for (int i = 0; i < _blossom_num; ++i) {
+	_delta2_index->set(i, _delta2->PRE_HEAP);
+	_delta4_index->set(i, _delta4->PRE_HEAP);
+      }
+
+      int index = 0;
+      for (NodeIt n(_ugraph); n != INVALID; ++n) {
+	Value max = 0;
+	for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  if (_ugraph.target(e) == n) continue;
+	  if ((dualScale * _weight[e]) / 2 > max) {
+	    max = (dualScale * _weight[e]) / 2;
+	  }
+	}
+	_node_index->set(n, index);
+	(*_node_data)[index].pot = max;
+	_delta1->push(n, max);
+	int blossom = 
+	  _blossom_set->insert(n, std::numeric_limits<Value>::max());
+
+	_tree_set->insert(blossom);
+
+	(*_blossom_data)[blossom].status = EVEN;
+	(*_blossom_data)[blossom].pred = INVALID;
+	(*_blossom_data)[blossom].next = INVALID;
+	(*_blossom_data)[blossom].pot = 0;
+	(*_blossom_data)[blossom].offset = 0;
+	++index;
+      }
+      for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
+	int si = (*_node_index)[_ugraph.source(e)];
+	int ti = (*_node_index)[_ugraph.target(e)];
+	if (_ugraph.source(e) != _ugraph.target(e)) {
+	  _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - 
+			    dualScale * _weight[e]) / 2);
+	}
+      }
+    }
+
+    /// \brief Starts the algorithm
+    ///
+    /// Starts the algorithm
+    void start() {
+      enum OpType {
+	D1, D2, D3, D4
+      };
+
+      int unmatched = _node_num;
+      while (unmatched > 0) {
+	Value d1 = !_delta1->empty() ? 
+	  _delta1->prio() : std::numeric_limits<Value>::max();
+
+	Value d2 = !_delta2->empty() ? 
+	  _delta2->prio() : std::numeric_limits<Value>::max();
+
+	Value d3 = !_delta3->empty() ? 
+	  _delta3->prio() : std::numeric_limits<Value>::max();
+
+	Value d4 = !_delta4->empty() ? 
+	  _delta4->prio() : std::numeric_limits<Value>::max();
+
+	_delta_sum = d1; OpType ot = D1;
+	if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
+	if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
+	if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
+
+	
+	switch (ot) {
+	case D1:
+	  {
+	    Node n = _delta1->top();
+	    unmatchNode(n);
+	    --unmatched;
+	  }
+	  break;
+	case D2:
+	  {
+	    int blossom = _delta2->top();
+	    Node n = _blossom_set->classTop(blossom);
+	    Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
+	    extendOnEdge(e);
+	  }
+	  break;
+	case D3:
+	  {
+	    UEdge e = _delta3->top();
+	    
+	    int left_blossom = _blossom_set->find(_ugraph.source(e));
+	    int right_blossom = _blossom_set->find(_ugraph.target(e));
+	  
+	    if (left_blossom == right_blossom) {
+	      _delta3->pop();
+	    } else {
+	      int left_tree;
+	      if ((*_blossom_data)[left_blossom].status == EVEN) {
+		left_tree = _tree_set->find(left_blossom);
+	      } else {
+		left_tree = -1;
+		++unmatched;
+	      }
+	      int right_tree;
+	      if ((*_blossom_data)[right_blossom].status == EVEN) {
+		right_tree = _tree_set->find(right_blossom);
+	      } else {
+		right_tree = -1;
+		++unmatched;
+	      }
+	    
+	      if (left_tree == right_tree) {
+		shrinkOnEdge(e, left_tree);
+	      } else {
+		augmentOnEdge(e);
+		unmatched -= 2;
+	      }
+	    }
+	  } break;
+	case D4:
+	  splitBlossom(_delta4->top());
+	  break;
+	}
+      }
+      extractMatching();
+    }
+
+    /// \brief Runs %MaxWeightedMatching algorithm.
+    ///
+    /// This method runs the %MaxWeightedMatching algorithm.
+    ///
+    /// \note mwm.run() is just a shortcut of the following code.
+    /// \code
+    ///   mwm.init();
+    ///   mwm.start();
+    /// \endcode
+    void run() {
+      init();
+      start();
+    }
+
+    /// @}
+
+    /// \name Primal solution
+    /// Functions for get the primal solution, ie. the matching.
+    
+    /// @{
+
+    /// \brief Returns the matching value.
+    ///
+    /// Returns the matching value.
+    Value matchingValue() const {
+      Value sum = 0;
+      for (NodeIt n(_ugraph); n != INVALID; ++n) {
+	if ((*_matching)[n] != INVALID) {
+	  sum += _weight[(*_matching)[n]];
+	}
+      }
+      return sum /= 2;
+    }
+
+    /// \brief Returns true when the edge is in the matching.
+    ///
+    /// Returns true when the edge is in the matching.
+    bool matching(const UEdge& edge) const {
+      return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
+    }
+
+    /// \brief Returns the incident matching edge.
+    ///
+    /// Returns the incident matching edge from given node. If the
+    /// node is not matched then it gives back \c INVALID.
+    Edge matching(const Node& node) const {
+      return (*_matching)[node];
+    }
+
+    /// \brief Returns the mate of the node.
+    ///
+    /// Returns the adjancent node in a mathcing edge. If the node is
+    /// not matched then it gives back \c INVALID.
+    Node mate(const Node& node) const {
+      return (*_matching)[node] != INVALID ?
+	_ugraph.target((*_matching)[node]) : INVALID;
+    }
+
+    /// @}
+
+    /// \name Dual solution
+    /// Functions for get the dual solution.
+    
+    /// @{
+
+    /// \brief Returns the value of the dual solution.
+    ///
+    /// Returns the value of the dual solution. It should be equal to
+    /// the primal value scaled by \ref dualScale "dual scale".
+    Value dualValue() const {
+      Value sum = 0;
+      for (NodeIt n(_ugraph); n != INVALID; ++n) {
+	sum += nodeValue(n);
+      }
+      for (int i = 0; i < blossomNum(); ++i) {
+        sum += blossomValue(i) * (blossomSize(i) / 2);
+      }
+      return sum;
+    }
+
+    /// \brief Returns the value of the node.
+    ///
+    /// Returns the the value of the node.
+    Value nodeValue(const Node& n) const {
+      return (*_node_potential)[n];
+    }
+
+    /// \brief Returns the number of the blossoms in the basis.
+    ///
+    /// Returns the number of the blossoms in the basis.
+    /// \see BlossomIt
+    int blossomNum() const {
+      return _blossom_potential.size();
+    }
+
+
+    /// \brief Returns the number of the nodes in the blossom.
+    ///
+    /// Returns the number of the nodes in the blossom.
+    int blossomSize(int k) const {
+      return _blossom_potential[k].end - _blossom_potential[k].begin;
+    }
+
+    /// \brief Returns the value of the blossom.
+    ///
+    /// Returns the the value of the blossom.
+    /// \see BlossomIt
+    Value blossomValue(int k) const {
+      return _blossom_potential[k].value;
+    }
+
+    /// \brief Lemon iterator for get the items of the blossom.
+    ///
+    /// Lemon iterator for get the nodes of the blossom. This class
+    /// provides a common style lemon iterator which gives back a
+    /// subset of the nodes.
+    class BlossomIt {
+    public:
+
+      /// \brief Constructor.
+      ///
+      /// Constructor for get the nodes of the variable. 
+      BlossomIt(const MaxWeightedMatching& algorithm, int variable) 
+        : _algorithm(&algorithm)
+      {
+        _index = _algorithm->_blossom_potential[variable].begin;
+        _last = _algorithm->_blossom_potential[variable].end;
+      }
+
+      /// \brief Invalid constructor.
+      ///
+      /// Invalid constructor.
+      BlossomIt(Invalid) : _index(-1) {}
+
+      /// \brief Conversion to node.
+      ///
+      /// Conversion to node.
+      operator Node() const { 
+        return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
+      }
+
+      /// \brief Increment operator.
+      ///
+      /// Increment operator.
+      BlossomIt& operator++() {
+        ++_index;
+        if (_index == _last) {
+          _index = -1;
+        }
+        return *this; 
+      }
+
+      bool operator==(const BlossomIt& it) const { 
+        return _index == it._index; 
+      }
+      bool operator!=(const BlossomIt& it) const { 
+        return _index != it._index;
+      }
+
+    private:
+      const MaxWeightedMatching* _algorithm;
+      int _last;
+      int _index;
+    };
+
+    /// @}
+    
+  };
+
+  /// \ingroup matching
+  ///
+  /// \brief Weighted perfect matching in general undirected graphs
+  ///
+  /// This class provides an efficient implementation of Edmond's
+  /// maximum weighted perfecr matching algorithm. The implementation
+  /// is based on extensive use of priority queues and provides
+  /// \f$O(nm\log(n))\f$ time complexity.
+  ///
+  /// The maximum weighted matching problem is to find undirected
+  /// edges in the graph with maximum overall weight and no two of
+  /// them shares their endpoints and covers all nodes. The problem
+  /// can be formulated with the next linear program:
+  /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
+  ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f] 
+  /// \f[x_e \ge 0\quad \forall e\in E\f]
+  /// \f[\max \sum_{e\in E}x_ew_e\f]
+  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
+  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both endpoints in
+  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
+  /// the nodes.
+  ///
+  /// The algorithm calculates an optimal matching and a proof of the
+  /// optimality. The solution of the dual problem can be used to check
+  /// the result of the algorithm. The dual linear problem is the next:
+  /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
+  /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
+  /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
+  ///
+  /// The algorithm can be executed with \c run() or the \c init() and
+  /// then the \c start() member functions. After it the matching can
+  /// be asked with \c matching() or mate() functions. The dual
+  /// solution can be get with \c nodeValue(), \c blossomNum() and \c
+  /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
+  /// "BlossomIt" nested class which is able to iterate on the nodes
+  /// of a blossom. If the value type is integral then the dual
+  /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
+  ///
+  /// \author Balazs Dezso
+  template <typename _UGraph, 
+	    typename _WeightMap = typename _UGraph::template UEdgeMap<int> >
+  class MaxWeightedPerfectMatching {
+  public:
+
+    typedef _UGraph UGraph;
+    typedef _WeightMap WeightMap;
+    typedef typename WeightMap::Value Value;
+
+    /// \brief Scaling factor for dual solution
+    ///
+    /// Scaling factor for dual solution, it is equal to 4 or 1
+    /// according to the value type.
+    static const int dualScale = 
+      std::numeric_limits<Value>::is_integer ? 4 : 1;
+
+    typedef typename UGraph::template NodeMap<typename UGraph::Edge> 
+    MatchingMap;
+    
+  private:
+
+    UGRAPH_TYPEDEFS(typename UGraph);
+
+    typedef typename UGraph::template NodeMap<Value> NodePotential;
+    typedef std::vector<Node> BlossomNodeList;
+
+    struct BlossomVariable {
+      int begin, end;
+      Value value;
+      
+      BlossomVariable(int _begin, int _end, Value _value)
+        : begin(_begin), end(_end), value(_value) {}
+
+    };
+
+    typedef std::vector<BlossomVariable> BlossomPotential;
+
+    const UGraph& _ugraph;
+    const WeightMap& _weight;
+
+    MatchingMap* _matching;
+
+    NodePotential* _node_potential;
+
+    BlossomPotential _blossom_potential;
+    BlossomNodeList _blossom_node_list;
+
+    int _node_num;
+    int _blossom_num;
+
+    typedef typename UGraph::template NodeMap<int> NodeIntMap;
+    typedef typename UGraph::template EdgeMap<int> EdgeIntMap;
+    typedef typename UGraph::template UEdgeMap<int> UEdgeIntMap;
+    typedef IntegerMap<int> IntIntMap;
+
+    enum Status {
+      EVEN = -1, MATCHED = 0, ODD = 1
+    };
+
+    typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
+    struct BlossomData {
+      int tree;
+      Status status;
+      Edge pred, next;
+      Value pot, offset;
+    };
+
+    NodeIntMap *_blossom_index;
+    BlossomSet *_blossom_set;
+    IntegerMap<BlossomData>* _blossom_data;
+
+    NodeIntMap *_node_index;
+    EdgeIntMap *_node_heap_index;
+
+    struct NodeData {
+      
+      NodeData(EdgeIntMap& node_heap_index) 
+	: heap(node_heap_index) {}
+      
+      int blossom;
+      Value pot;
+      BinHeap<Value, EdgeIntMap> heap;
+      std::map<int, Edge> heap_index;
+      
+      int tree;
+    };
+
+    IntegerMap<NodeData>* _node_data;
+
+    typedef ExtendFindEnum<IntIntMap> TreeSet;
+
+    IntIntMap *_tree_set_index;
+    TreeSet *_tree_set;
+
+    IntIntMap *_delta2_index;
+    BinHeap<Value, IntIntMap> *_delta2;
+    
+    UEdgeIntMap *_delta3_index;
+    BinHeap<Value, UEdgeIntMap> *_delta3;
+
+    IntIntMap *_delta4_index;
+    BinHeap<Value, IntIntMap> *_delta4;
+
+    Value _delta_sum;
+
+    void createStructures() {
+      _node_num = countNodes(_ugraph); 
+      _blossom_num = _node_num * 3 / 2;
+
+      if (!_matching) {
+	_matching = new MatchingMap(_ugraph);
+      }
+      if (!_node_potential) {
+	_node_potential = new NodePotential(_ugraph);
+      }
+      if (!_blossom_set) {
+	_blossom_index = new NodeIntMap(_ugraph);
+	_blossom_set = new BlossomSet(*_blossom_index);
+	_blossom_data = new IntegerMap<BlossomData>(_blossom_num);
+      }
+
+      if (!_node_index) {
+	_node_index = new NodeIntMap(_ugraph);
+	_node_heap_index = new EdgeIntMap(_ugraph);
+	_node_data = new IntegerMap<NodeData>(_node_num, 
+					      NodeData(*_node_heap_index));
+      }
+
+      if (!_tree_set) {
+	_tree_set_index = new IntIntMap(_blossom_num);
+	_tree_set = new TreeSet(*_tree_set_index);
+      }
+      if (!_delta2) {
+	_delta2_index = new IntIntMap(_blossom_num);
+	_delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
+      }
+      if (!_delta3) {
+	_delta3_index = new UEdgeIntMap(_ugraph);
+	_delta3 = new BinHeap<Value, UEdgeIntMap>(*_delta3_index);
+      }
+      if (!_delta4) {
+	_delta4_index = new IntIntMap(_blossom_num);
+	_delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
+      }
+    }
+
+    void destroyStructures() {
+      _node_num = countNodes(_ugraph); 
+      _blossom_num = _node_num * 3 / 2;
+
+      if (_matching) {
+	delete _matching;
+      }
+      if (_node_potential) {
+	delete _node_potential;
+      }
+      if (_blossom_set) {
+	delete _blossom_index;
+	delete _blossom_set;
+	delete _blossom_data;
+      }
+
+      if (_node_index) {
+	delete _node_index;
+	delete _node_heap_index;
+	delete _node_data;			      
+      }
+
+      if (_tree_set) {
+	delete _tree_set_index;
+	delete _tree_set;
+      }
+      if (_delta2) {
+	delete _delta2_index;
+	delete _delta2;
+      }
+      if (_delta3) {
+	delete _delta3_index;
+	delete _delta3;
+      }
+      if (_delta4) {
+	delete _delta4_index;
+	delete _delta4;
+      }
+    }
+
+    void matchedToEven(int blossom, int tree) {
+      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+	_delta2->erase(blossom);
+      }
+
+      if (!_blossom_set->trivial(blossom)) {
+	(*_blossom_data)[blossom].pot -= 
+	  2 * (_delta_sum - (*_blossom_data)[blossom].offset);
+      }
+
+      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 
+	   n != INVALID; ++n) {
+
+	_blossom_set->increase(n, std::numeric_limits<Value>::max());
+	int ni = (*_node_index)[n];
+
+	(*_node_data)[ni].heap.clear();
+	(*_node_data)[ni].heap_index.clear();
+
+	(*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
+
+	for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  Node v = _ugraph.source(e);
+	  int vb = _blossom_set->find(v);
+	  int vi = (*_node_index)[v];
+
+	  Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 
+	    dualScale * _weight[e];
+
+	  if ((*_blossom_data)[vb].status == EVEN) {
+	    if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
+	      _delta3->push(e, rw / 2);
+	    }
+	  } else {
+	    typename std::map<int, Edge>::iterator it = 
+	      (*_node_data)[vi].heap_index.find(tree);	  
+
+	    if (it != (*_node_data)[vi].heap_index.end()) {
+	      if ((*_node_data)[vi].heap[it->second] > rw) {
+		(*_node_data)[vi].heap.replace(it->second, e);
+		(*_node_data)[vi].heap.decrease(e, rw);
+		it->second = e;
+	      }
+	    } else {
+	      (*_node_data)[vi].heap.push(e, rw);
+	      (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+	    }
+
+	    if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
+	      _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
+
+	      if ((*_blossom_data)[vb].status == MATCHED) {
+		if (_delta2->state(vb) != _delta2->IN_HEAP) {
+		  _delta2->push(vb, _blossom_set->classPrio(vb) -
+			       (*_blossom_data)[vb].offset);
+		} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
+			   (*_blossom_data)[vb].offset){
+		  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
+				   (*_blossom_data)[vb].offset);
+		}
+	      }
+	    }
+	  }
+	}
+      }
+      (*_blossom_data)[blossom].offset = 0;
+    }
+
+    void matchedToOdd(int blossom) {
+      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+	_delta2->erase(blossom);
+      }
+      (*_blossom_data)[blossom].offset += _delta_sum;
+      if (!_blossom_set->trivial(blossom)) {
+	_delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 + 
+		     (*_blossom_data)[blossom].offset);
+      }
+    }
+
+    void evenToMatched(int blossom, int tree) {
+      if (!_blossom_set->trivial(blossom)) {
+	(*_blossom_data)[blossom].pot += 2 * _delta_sum;
+      }
+
+      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
+	   n != INVALID; ++n) {
+	int ni = (*_node_index)[n];
+	(*_node_data)[ni].pot -= _delta_sum;
+
+	for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  Node v = _ugraph.source(e);
+	  int vb = _blossom_set->find(v);
+	  int vi = (*_node_index)[v];
+
+	  Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 
+	    dualScale * _weight[e];
+
+	  if (vb == blossom) {
+	    if (_delta3->state(e) == _delta3->IN_HEAP) {
+	      _delta3->erase(e);
+	    }
+	  } else if ((*_blossom_data)[vb].status == EVEN) {
+
+	    if (_delta3->state(e) == _delta3->IN_HEAP) {
+	      _delta3->erase(e);
+	    }
+
+	    int vt = _tree_set->find(vb);
+	    
+	    if (vt != tree) {
+
+	      Edge r = _ugraph.oppositeEdge(e);
+
+	      typename std::map<int, Edge>::iterator it = 
+		(*_node_data)[ni].heap_index.find(vt);	  
+
+	      if (it != (*_node_data)[ni].heap_index.end()) {
+		if ((*_node_data)[ni].heap[it->second] > rw) {
+		  (*_node_data)[ni].heap.replace(it->second, r);
+		  (*_node_data)[ni].heap.decrease(r, rw);
+		  it->second = r;
+		}
+	      } else {
+		(*_node_data)[ni].heap.push(r, rw);
+		(*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
+	      }
+	      
+	      if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
+		_blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
+		
+		if (_delta2->state(blossom) != _delta2->IN_HEAP) {
+		  _delta2->push(blossom, _blossom_set->classPrio(blossom) - 
+			       (*_blossom_data)[blossom].offset);
+		} else if ((*_delta2)[blossom] > 
+			   _blossom_set->classPrio(blossom) - 
+			   (*_blossom_data)[blossom].offset){
+		  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
+				   (*_blossom_data)[blossom].offset);
+		}
+	      }
+	    }
+	  } else {
+	  
+	    typename std::map<int, Edge>::iterator it = 
+	      (*_node_data)[vi].heap_index.find(tree);
+
+	    if (it != (*_node_data)[vi].heap_index.end()) {
+	      (*_node_data)[vi].heap.erase(it->second);
+	      (*_node_data)[vi].heap_index.erase(it);
+	      if ((*_node_data)[vi].heap.empty()) {
+		_blossom_set->increase(v, std::numeric_limits<Value>::max());
+	      } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
+		_blossom_set->increase(v, (*_node_data)[vi].heap.prio());
+	      }
+	      
+	      if ((*_blossom_data)[vb].status == MATCHED) {
+		if (_blossom_set->classPrio(vb) == 
+		    std::numeric_limits<Value>::max()) {
+		  _delta2->erase(vb);
+		} else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
+			   (*_blossom_data)[vb].offset) {
+		  _delta2->increase(vb, _blossom_set->classPrio(vb) -
+				   (*_blossom_data)[vb].offset);
+		}
+	      }	
+	    }	
+	  }
+	}
+      }
+    }
+
+    void oddToMatched(int blossom) {
+      (*_blossom_data)[blossom].offset -= _delta_sum;
+
+      if (_blossom_set->classPrio(blossom) != 
+	  std::numeric_limits<Value>::max()) {
+	_delta2->push(blossom, _blossom_set->classPrio(blossom) - 
+		       (*_blossom_data)[blossom].offset);
+      }
+
+      if (!_blossom_set->trivial(blossom)) {
+	_delta4->erase(blossom);
+      }
+    }
+
+    void oddToEven(int blossom, int tree) {
+      if (!_blossom_set->trivial(blossom)) {
+	_delta4->erase(blossom);
+	(*_blossom_data)[blossom].pot -= 
+	  2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
+      }
+
+      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom); 
+	   n != INVALID; ++n) {
+	int ni = (*_node_index)[n];
+
+	_blossom_set->increase(n, std::numeric_limits<Value>::max());
+
+	(*_node_data)[ni].heap.clear();
+	(*_node_data)[ni].heap_index.clear();
+	(*_node_data)[ni].pot += 
+	  2 * _delta_sum - (*_blossom_data)[blossom].offset;
+
+	for (InEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  Node v = _ugraph.source(e);
+	  int vb = _blossom_set->find(v);
+	  int vi = (*_node_index)[v];
+
+	  Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot - 
+	    dualScale * _weight[e];
+
+	  if ((*_blossom_data)[vb].status == EVEN) {
+	    if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
+	      _delta3->push(e, rw / 2);
+	    }
+	  } else {
+	  
+	    typename std::map<int, Edge>::iterator it = 
+	      (*_node_data)[vi].heap_index.find(tree);	  
+
+	    if (it != (*_node_data)[vi].heap_index.end()) {
+	      if ((*_node_data)[vi].heap[it->second] > rw) {
+		(*_node_data)[vi].heap.replace(it->second, e);
+		(*_node_data)[vi].heap.decrease(e, rw);
+		it->second = e;
+	      }
+	    } else {
+	      (*_node_data)[vi].heap.push(e, rw);
+	      (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
+	    }
+
+	    if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
+	      _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
+
+	      if ((*_blossom_data)[vb].status == MATCHED) {
+		if (_delta2->state(vb) != _delta2->IN_HEAP) {
+		  _delta2->push(vb, _blossom_set->classPrio(vb) - 
+			       (*_blossom_data)[vb].offset);
+		} else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) - 
+			   (*_blossom_data)[vb].offset) {
+		  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
+				   (*_blossom_data)[vb].offset);
+		}
+	      }
+	    }
+	  }
+	}
+      }
+      (*_blossom_data)[blossom].offset = 0;
+    }
+    
+    void alternatePath(int even, int tree) {
+      int odd;
+
+      evenToMatched(even, tree);
+      (*_blossom_data)[even].status = MATCHED;
+
+      while ((*_blossom_data)[even].pred != INVALID) {
+	odd = _blossom_set->find(_ugraph.target((*_blossom_data)[even].pred));
+	(*_blossom_data)[odd].status = MATCHED;
+	oddToMatched(odd);
+	(*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
+      
+	even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].pred));
+	(*_blossom_data)[even].status = MATCHED;
+	evenToMatched(even, tree);
+	(*_blossom_data)[even].next = 
+	  _ugraph.oppositeEdge((*_blossom_data)[odd].pred);
+      }
+      
+    }
+
+    void destroyTree(int tree) {
+      for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
+	if ((*_blossom_data)[b].status == EVEN) {
+	  (*_blossom_data)[b].status = MATCHED;
+	  evenToMatched(b, tree);
+	} else if ((*_blossom_data)[b].status == ODD) {
+	  (*_blossom_data)[b].status = MATCHED;
+	  oddToMatched(b);
+	}
+      }
+      _tree_set->eraseClass(tree);
+    }
+
+    void augmentOnEdge(const UEdge& edge) {
+      
+      int left = _blossom_set->find(_ugraph.source(edge));
+      int right = _blossom_set->find(_ugraph.target(edge));
+
+      int left_tree = _tree_set->find(left);
+      alternatePath(left, left_tree);
+      destroyTree(left_tree);
+
+      int right_tree = _tree_set->find(right);
+      alternatePath(right, right_tree);
+      destroyTree(right_tree);
+
+      (*_blossom_data)[left].next = _ugraph.direct(edge, true);
+      (*_blossom_data)[right].next = _ugraph.direct(edge, false);
+    }
+
+    void extendOnEdge(const Edge& edge) {
+      int base = _blossom_set->find(_ugraph.target(edge));
+      int tree = _tree_set->find(base);
+      
+      int odd = _blossom_set->find(_ugraph.source(edge));
+      _tree_set->insert(odd, tree);
+      (*_blossom_data)[odd].status = ODD;
+      matchedToOdd(odd);
+      (*_blossom_data)[odd].pred = edge;
+
+      int even = _blossom_set->find(_ugraph.target((*_blossom_data)[odd].next));
+      (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
+      _tree_set->insert(even, tree);
+      (*_blossom_data)[even].status = EVEN;
+      matchedToEven(even, tree);
+    }
+    
+    void shrinkOnEdge(const UEdge& uedge, int tree) {
+      int nca = -1;
+      std::vector<int> left_path, right_path;
+
+      {
+	std::set<int> left_set, right_set;
+	int left = _blossom_set->find(_ugraph.source(uedge));
+	left_path.push_back(left);
+	left_set.insert(left);
+
+	int right = _blossom_set->find(_ugraph.target(uedge));
+	right_path.push_back(right);
+	right_set.insert(right);
+
+	while (true) {
+
+	  if ((*_blossom_data)[left].pred == INVALID) break;
+
+	  left = 
+	    _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
+	  left_path.push_back(left);
+	  left = 
+	    _blossom_set->find(_ugraph.target((*_blossom_data)[left].pred));
+	  left_path.push_back(left);
+
+	  left_set.insert(left);
+
+	  if (right_set.find(left) != right_set.end()) {
+	    nca = left;
+	    break;
+	  }
+
+	  if ((*_blossom_data)[right].pred == INVALID) break;
+
+	  right = 
+	    _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
+	  right_path.push_back(right);
+	  right = 
+	    _blossom_set->find(_ugraph.target((*_blossom_data)[right].pred));
+	  right_path.push_back(right);
+
+	  right_set.insert(right);
+ 
+	  if (left_set.find(right) != left_set.end()) {
+	    nca = right;
+	    break;
+	  }
+
+	}
+
+	if (nca == -1) {
+	  if ((*_blossom_data)[left].pred == INVALID) {
+	    nca = right;
+	    while (left_set.find(nca) == left_set.end()) {
+	      nca = 
+		_blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
+	      right_path.push_back(nca);
+	      nca = 
+		_blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
+	      right_path.push_back(nca);
+	    }
+	  } else {
+	    nca = left;
+	    while (right_set.find(nca) == right_set.end()) {
+	      nca = 
+		_blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
+	      left_path.push_back(nca);
+	      nca = 
+		_blossom_set->find(_ugraph.target((*_blossom_data)[nca].pred));
+	      left_path.push_back(nca);
+	    }
+	  }
+	}
+      }
+
+      std::vector<int> subblossoms;
+      Edge prev;
+
+      prev = _ugraph.direct(uedge, true);
+      for (int i = 0; left_path[i] != nca; i += 2) {
+	subblossoms.push_back(left_path[i]);
+	(*_blossom_data)[left_path[i]].next = prev;
+	_tree_set->erase(left_path[i]);
+
+	subblossoms.push_back(left_path[i + 1]);
+	(*_blossom_data)[left_path[i + 1]].status = EVEN;
+	oddToEven(left_path[i + 1], tree);
+	_tree_set->erase(left_path[i + 1]);
+	prev = _ugraph.oppositeEdge((*_blossom_data)[left_path[i + 1]].pred);
+      }
+
+      int k = 0;
+      while (right_path[k] != nca) ++k;
+
+      subblossoms.push_back(nca);
+      (*_blossom_data)[nca].next = prev;
+      
+      for (int i = k - 2; i >= 0; i -= 2) {
+	subblossoms.push_back(right_path[i + 1]);
+	(*_blossom_data)[right_path[i + 1]].status = EVEN;
+	oddToEven(right_path[i + 1], tree);
+	_tree_set->erase(right_path[i + 1]);
+
+	(*_blossom_data)[right_path[i + 1]].next = 
+	  (*_blossom_data)[right_path[i + 1]].pred;
+
+	subblossoms.push_back(right_path[i]);
+	_tree_set->erase(right_path[i]);
+      }
+
+      int surface = 
+	_blossom_set->join(subblossoms.begin(), subblossoms.end());
+
+      for (int i = 0; i < int(subblossoms.size()); ++i) {
+	if (!_blossom_set->trivial(subblossoms[i])) {
+	  (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
+	}
+	(*_blossom_data)[subblossoms[i]].status = MATCHED;
+      }
+
+      (*_blossom_data)[surface].pot = -2 * _delta_sum;
+      (*_blossom_data)[surface].offset = 0;
+      (*_blossom_data)[surface].status = EVEN;
+      (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
+      (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
+
+      _tree_set->insert(surface, tree);
+      _tree_set->erase(nca);
+    }
+
+    void splitBlossom(int blossom) {
+      Edge next = (*_blossom_data)[blossom].next; 
+      Edge pred = (*_blossom_data)[blossom].pred;
+
+      int tree = _tree_set->find(blossom);
+
+      (*_blossom_data)[blossom].status = MATCHED;
+      oddToMatched(blossom);
+      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
+	_delta2->erase(blossom);
+      }
+
+      std::vector<int> subblossoms;
+      _blossom_set->split(blossom, std::back_inserter(subblossoms));
+
+      Value offset = (*_blossom_data)[blossom].offset;
+      int b = _blossom_set->find(_ugraph.source(pred));
+      int d = _blossom_set->find(_ugraph.source(next));
+      
+      int ib, id;
+      for (int i = 0; i < int(subblossoms.size()); ++i) {
+	if (subblossoms[i] == b) ib = i;
+	if (subblossoms[i] == d) id = i;
+
+	(*_blossom_data)[subblossoms[i]].offset = offset;
+	if (!_blossom_set->trivial(subblossoms[i])) {
+	  (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
+	}
+	if (_blossom_set->classPrio(subblossoms[i]) != 
+	    std::numeric_limits<Value>::max()) {
+	  _delta2->push(subblossoms[i], 
+			_blossom_set->classPrio(subblossoms[i]) - 
+			(*_blossom_data)[subblossoms[i]].offset);
+	}
+      }
+      
+      if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
+	for (int i = (id + 1) % subblossoms.size(); 
+	     i != ib; i = (i + 2) % subblossoms.size()) {
+	  int sb = subblossoms[i];
+	  int tb = subblossoms[(i + 1) % subblossoms.size()];
+	  (*_blossom_data)[sb].next = 
+	    _ugraph.oppositeEdge((*_blossom_data)[tb].next);
+	}
+
+	for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
+	  int sb = subblossoms[i];
+	  int tb = subblossoms[(i + 1) % subblossoms.size()];
+	  int ub = subblossoms[(i + 2) % subblossoms.size()];
+
+	  (*_blossom_data)[sb].status = ODD;
+	  matchedToOdd(sb);
+	  _tree_set->insert(sb, tree);
+	  (*_blossom_data)[sb].pred = pred;
+	  (*_blossom_data)[sb].next = 
+			   _ugraph.oppositeEdge((*_blossom_data)[tb].next);
+	  
+	  pred = (*_blossom_data)[ub].next;
+      
+	  (*_blossom_data)[tb].status = EVEN;
+	  matchedToEven(tb, tree);
+	  _tree_set->insert(tb, tree);
+	  (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
+	}
+      
+	(*_blossom_data)[subblossoms[id]].status = ODD;
+	matchedToOdd(subblossoms[id]);
+	_tree_set->insert(subblossoms[id], tree);
+	(*_blossom_data)[subblossoms[id]].next = next;
+	(*_blossom_data)[subblossoms[id]].pred = pred;
+      
+      } else {
+
+	for (int i = (ib + 1) % subblossoms.size(); 
+	     i != id; i = (i + 2) % subblossoms.size()) {
+	  int sb = subblossoms[i];
+	  int tb = subblossoms[(i + 1) % subblossoms.size()];
+	  (*_blossom_data)[sb].next = 
+	    _ugraph.oppositeEdge((*_blossom_data)[tb].next);
+	}	
+
+	for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
+	  int sb = subblossoms[i];
+	  int tb = subblossoms[(i + 1) % subblossoms.size()];
+	  int ub = subblossoms[(i + 2) % subblossoms.size()];
+
+	  (*_blossom_data)[sb].status = ODD;
+	  matchedToOdd(sb);
+	  _tree_set->insert(sb, tree);
+	  (*_blossom_data)[sb].next = next; 
+	  (*_blossom_data)[sb].pred =
+	    _ugraph.oppositeEdge((*_blossom_data)[tb].next);
+
+	  (*_blossom_data)[tb].status = EVEN;
+	  matchedToEven(tb, tree);
+	  _tree_set->insert(tb, tree);
+	  (*_blossom_data)[tb].pred =
+	    (*_blossom_data)[tb].next = 
+	    _ugraph.oppositeEdge((*_blossom_data)[ub].next);
+	  next = (*_blossom_data)[ub].next;
+	}
+
+	(*_blossom_data)[subblossoms[ib]].status = ODD;
+	matchedToOdd(subblossoms[ib]);
+	_tree_set->insert(subblossoms[ib], tree);
+	(*_blossom_data)[subblossoms[ib]].next = next;
+	(*_blossom_data)[subblossoms[ib]].pred = pred;
+      }
+      _tree_set->erase(blossom);
+    }
+
+    void extractBlossom(int blossom, const Node& base, const Edge& matching) {
+      if (_blossom_set->trivial(blossom)) {
+	int bi = (*_node_index)[base];
+	Value pot = (*_node_data)[bi].pot;
+
+	_matching->set(base, matching);
+	_blossom_node_list.push_back(base);
+	_node_potential->set(base, pot);
+      } else {
+
+	Value pot = (*_blossom_data)[blossom].pot;
+	int bn = _blossom_node_list.size();
+	
+	std::vector<int> subblossoms;
+	_blossom_set->split(blossom, std::back_inserter(subblossoms));
+	int b = _blossom_set->find(base);
+	int ib = -1;
+	for (int i = 0; i < int(subblossoms.size()); ++i) {
+	  if (subblossoms[i] == b) { ib = i; break; }
+	}
+	
+	for (int i = 1; i < int(subblossoms.size()); i += 2) {
+	  int sb = subblossoms[(ib + i) % subblossoms.size()];
+	  int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
+	  
+	  Edge m = (*_blossom_data)[tb].next;
+	  extractBlossom(sb, _ugraph.target(m), _ugraph.oppositeEdge(m));
+	  extractBlossom(tb, _ugraph.source(m), m);
+	}
+	extractBlossom(subblossoms[ib], base, matching);      
+	
+	int en = _blossom_node_list.size();
+	
+	_blossom_potential.push_back(BlossomVariable(bn, en, pot));
+      }
+    }
+
+    void extractMatching() {
+      std::vector<int> blossoms;
+      for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
+	blossoms.push_back(c);
+      }
+
+      for (int i = 0; i < int(blossoms.size()); ++i) {
+
+	Value offset = (*_blossom_data)[blossoms[i]].offset;
+	(*_blossom_data)[blossoms[i]].pot += 2 * offset;
+	for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]); 
+	     n != INVALID; ++n) {
+	  (*_node_data)[(*_node_index)[n]].pot -= offset;
+	}
+	
+	Edge matching = (*_blossom_data)[blossoms[i]].next;
+	Node base = _ugraph.source(matching);
+	extractBlossom(blossoms[i], base, matching);
+      }
+    }
+    
+  public:
+
+    /// \brief Constructor
+    ///
+    /// Constructor.
+    MaxWeightedPerfectMatching(const UGraph& ugraph, const WeightMap& weight)
+      : _ugraph(ugraph), _weight(weight), _matching(0),
+	_node_potential(0), _blossom_potential(), _blossom_node_list(),
+	_node_num(0), _blossom_num(0),
+
+	_blossom_index(0), _blossom_set(0), _blossom_data(0),
+	_node_index(0), _node_heap_index(0), _node_data(0),
+	_tree_set_index(0), _tree_set(0),
+
+	_delta2_index(0), _delta2(0),
+	_delta3_index(0), _delta3(0), 
+	_delta4_index(0), _delta4(0),
+
+	_delta_sum() {}
+
+    ~MaxWeightedPerfectMatching() {
+      destroyStructures();
+    }
+
+    /// \name Execution control 
+    /// The simplest way to execute the algorithm is to use the member
+    /// \c run() member function.
+
+    ///@{
+
+    /// \brief Initialize the algorithm
+    ///
+    /// Initialize the algorithm
+    void init() {
+      createStructures();
+
+      for (EdgeIt e(_ugraph); e != INVALID; ++e) {
+	_node_heap_index->set(e, BinHeap<Value, EdgeIntMap>::PRE_HEAP);
+      }
+      for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
+	_delta3_index->set(e, _delta3->PRE_HEAP);
+      }
+      for (int i = 0; i < _blossom_num; ++i) {
+	_delta2_index->set(i, _delta2->PRE_HEAP);
+	_delta4_index->set(i, _delta4->PRE_HEAP);
+      }
+
+      int index = 0;
+      for (NodeIt n(_ugraph); n != INVALID; ++n) {
+	Value max = std::numeric_limits<Value>::min();
+	for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  if (_ugraph.target(e) == n) continue;
+	  if ((dualScale * _weight[e]) / 2 > max) {
+	    max = (dualScale * _weight[e]) / 2;
+	  }
+	}
+	_node_index->set(n, index);
+	(*_node_data)[index].pot = max;
+	int blossom = 
+	  _blossom_set->insert(n, std::numeric_limits<Value>::max());
+
+	_tree_set->insert(blossom);
+
+	(*_blossom_data)[blossom].status = EVEN;
+	(*_blossom_data)[blossom].pred = INVALID;
+	(*_blossom_data)[blossom].next = INVALID;
+	(*_blossom_data)[blossom].pot = 0;
+	(*_blossom_data)[blossom].offset = 0;
+	++index;
+      }
+      for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
+	int si = (*_node_index)[_ugraph.source(e)];
+	int ti = (*_node_index)[_ugraph.target(e)];
+	if (_ugraph.source(e) != _ugraph.target(e)) {
+	  _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot - 
+			    dualScale * _weight[e]) / 2);
+	}
+      }
+    }
+
+    /// \brief Starts the algorithm
+    ///
+    /// Starts the algorithm
+    bool start() {
+      enum OpType {
+	D2, D3, D4
+      };
+
+      int unmatched = _node_num;
+      while (unmatched > 0) {
+	Value d2 = !_delta2->empty() ? 
+	  _delta2->prio() : std::numeric_limits<Value>::max();
+
+	Value d3 = !_delta3->empty() ? 
+	  _delta3->prio() : std::numeric_limits<Value>::max();
+
+	Value d4 = !_delta4->empty() ? 
+	  _delta4->prio() : std::numeric_limits<Value>::max();
+
+	_delta_sum = d2; OpType ot = D2;
+	if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
+	if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
+
+	if (_delta_sum == std::numeric_limits<Value>::max()) {
+	  return false;
+	}
+	
+	switch (ot) {
+	case D2:
+	  {
+	    int blossom = _delta2->top();
+	    Node n = _blossom_set->classTop(blossom);
+	    Edge e = (*_node_data)[(*_node_index)[n]].heap.top();
+	    extendOnEdge(e);
+	  }
+	  break;
+	case D3:
+	  {
+	    UEdge e = _delta3->top();
+	    
+	    int left_blossom = _blossom_set->find(_ugraph.source(e));
+	    int right_blossom = _blossom_set->find(_ugraph.target(e));
+	  
+	    if (left_blossom == right_blossom) {
+	      _delta3->pop();
+	    } else {
+	      int left_tree = _tree_set->find(left_blossom);
+	      int right_tree = _tree_set->find(right_blossom);
+	    
+	      if (left_tree == right_tree) {
+		shrinkOnEdge(e, left_tree);
+	      } else {
+		augmentOnEdge(e);
+		unmatched -= 2;
+	      }
+	    }
+	  } break;
+	case D4:
+	  splitBlossom(_delta4->top());
+	  break;
+	}
+      }
+      extractMatching();
+      return true;
+    }
+
+    /// \brief Runs %MaxWeightedPerfectMatching algorithm.
+    ///
+    /// This method runs the %MaxWeightedPerfectMatching algorithm.
+    ///
+    /// \note mwm.run() is just a shortcut of the following code.
+    /// \code
+    ///   mwm.init();
+    ///   mwm.start();
+    /// \endcode
+    bool run() {
+      init();
+      return start();
+    }
+
+    /// @}
+
+    /// \name Primal solution
+    /// Functions for get the primal solution, ie. the matching.
+    
+    /// @{
+
+    /// \brief Returns the matching value.
+    ///
+    /// Returns the matching value.
+    Value matchingValue() const {
+      Value sum = 0;
+      for (NodeIt n(_ugraph); n != INVALID; ++n) {
+	if ((*_matching)[n] != INVALID) {
+	  sum += _weight[(*_matching)[n]];
+	}
+      }
+      return sum /= 2;
+    }
+
+    /// \brief Returns true when the edge is in the matching.
+    ///
+    /// Returns true when the edge is in the matching.
+    bool matching(const UEdge& edge) const {
+      return (*_matching)[_ugraph.source(edge)] == _ugraph.direct(edge, true);
+    }
+
+    /// \brief Returns the incident matching edge.
+    ///
+    /// Returns the incident matching edge from given node.
+    Edge matching(const Node& node) const {
+      return (*_matching)[node];
+    }
+
+    /// \brief Returns the mate of the node.
+    ///
+    /// Returns the adjancent node in a mathcing edge.
+    Node mate(const Node& node) const {
+      return _ugraph.target((*_matching)[node]);
+    }
+
+    /// @}
+
+    /// \name Dual solution
+    /// Functions for get the dual solution.
+    
+    /// @{
+
+    /// \brief Returns the value of the dual solution.
+    ///
+    /// Returns the value of the dual solution. It should be equal to
+    /// the primal value scaled by \ref dualScale "dual scale".
+    Value dualValue() const {
+      Value sum = 0;
+      for (NodeIt n(_ugraph); n != INVALID; ++n) {
+	sum += nodeValue(n);
+      }
+      for (int i = 0; i < blossomNum(); ++i) {
+        sum += blossomValue(i) * (blossomSize(i) / 2);
+      }
+      return sum;
+    }
+
+    /// \brief Returns the value of the node.
+    ///
+    /// Returns the the value of the node.
+    Value nodeValue(const Node& n) const {
+      return (*_node_potential)[n];
+    }
+
+    /// \brief Returns the number of the blossoms in the basis.
+    ///
+    /// Returns the number of the blossoms in the basis.
+    /// \see BlossomIt
+    int blossomNum() const {
+      return _blossom_potential.size();
+    }
+
+
+    /// \brief Returns the number of the nodes in the blossom.
+    ///
+    /// Returns the number of the nodes in the blossom.
+    int blossomSize(int k) const {
+      return _blossom_potential[k].end - _blossom_potential[k].begin;
+    }
+
+    /// \brief Returns the value of the blossom.
+    ///
+    /// Returns the the value of the blossom.
+    /// \see BlossomIt
+    Value blossomValue(int k) const {
+      return _blossom_potential[k].value;
+    }
+
+    /// \brief Lemon iterator for get the items of the blossom.
+    ///
+    /// Lemon iterator for get the nodes of the blossom. This class
+    /// provides a common style lemon iterator which gives back a
+    /// subset of the nodes.
+    class BlossomIt {
+    public:
+
+      /// \brief Constructor.
+      ///
+      /// Constructor for get the nodes of the variable. 
+      BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable) 
+        : _algorithm(&algorithm)
+      {
+        _index = _algorithm->_blossom_potential[variable].begin;
+        _last = _algorithm->_blossom_potential[variable].end;
+      }
+
+      /// \brief Invalid constructor.
+      ///
+      /// Invalid constructor.
+      BlossomIt(Invalid) : _index(-1) {}
+
+      /// \brief Conversion to node.
+      ///
+      /// Conversion to node.
+      operator Node() const { 
+        return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
+      }
+
+      /// \brief Increment operator.
+      ///
+      /// Increment operator.
+      BlossomIt& operator++() {
+        ++_index;
+        if (_index == _last) {
+          _index = -1;
+        }
+        return *this; 
+      }
+
+      bool operator==(const BlossomIt& it) const { 
+        return _index == it._index; 
+      }
+      bool operator!=(const BlossomIt& it) const { 
+        return _index != it._index;
+      }
+
+    private:
+      const MaxWeightedPerfectMatching* _algorithm;
+      int _last;
+      int _index;
+    };
+
+    /// @}
+    
+  };
+
   
 } //END OF NAMESPACE LEMON