Index: lemon/matching.h
===================================================================
--- lemon/matching.h	(revision 956)
+++ lemon/matching.h	(revision 1081)
@@ -3,5 +3,5 @@
  * This file is a part of LEMON, a generic C++ optimization library.
  *
- * Copyright (C) 2003-2010
+ * Copyright (C) 2003-2011
  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
  * (Egervary Research Group on Combinatorial Optimization, EGRES).
@@ -17,6 +17,6 @@
  */
 
-#ifndef LEMON_MATCHING_H
-#define LEMON_MATCHING_H
+#ifndef LEMON_MAX_MATCHING_H
+#define LEMON_MAX_MATCHING_H
 
 #include <vector>
@@ -29,5 +29,4 @@
 #include <lemon/bin_heap.h>
 #include <lemon/maps.h>
-#include <lemon/fractional_matching.h>
 
 ///\ingroup matching
@@ -747,5 +746,5 @@
 
     enum Status {
-      EVEN = -1, MATCHED = 0, ODD = 1
+      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
     };
 
@@ -799,8 +798,4 @@
 
     Value _delta_sum;
-    int _unmatched;
-
-    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
-    FractionalMatching *_fractional;
 
     void createStructures() {
@@ -869,4 +864,7 @@
 
     void destroyStructures() {
+      _node_num = countNodes(_graph);
+      _blossom_num = _node_num * 3 / 2;
+
       if (_matching) {
         delete _matching;
@@ -944,5 +942,205 @@
               _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, Arc>::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 (InArcIt e(_graph, n); e != INVALID; ++e) {
+          Node v = _graph.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) {
+
+              Arc r = _graph.oppositeArc(e);
+
+              typename std::map<int, Arc>::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, Arc>::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 (InArcIt e(_graph, n); e != INVALID; ++e) {
+          Node v = _graph.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, Arc>::iterator it =
               (*_node_data)[vi].heap_index.find(tree);
@@ -979,18 +1177,8 @@
     }
 
-    void matchedToOdd(int blossom) {
+
+    void matchedToUnmatched(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;
       }
 
@@ -998,7 +1186,31 @@
            n != INVALID; ++n) {
         int ni = (*_node_index)[n];
-        (*_node_data)[ni].pot -= _delta_sum;
-
-        _delta1->erase(n);
+
+        _blossom_set->increase(n, std::numeric_limits<Value>::max());
+
+        (*_node_data)[ni].heap.clear();
+        (*_node_data)[ni].heap_index.clear();
+
+        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+          Node v = _graph.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 (InArcIt e(_graph, n); e != INVALID; ++e) {
@@ -1022,146 +1234,40 @@
             int vt = _tree_set->find(vb);
 
-            if (vt != tree) {
-
-              Arc r = _graph.oppositeArc(e);
-
-              typename std::map<int, Arc>::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 {
+            Arc r = _graph.oppositeArc(e);
 
             typename std::map<int, Arc>::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 (InArcIt e(_graph, n); e != INVALID; ++e) {
-          Node v = _graph.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, Arc>::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;
+              (*_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)[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);
-                }
+              (*_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);
+            }
           }
         }
       }
-      (*_blossom_data)[blossom].offset = 0;
     }
 
@@ -1208,7 +1314,9 @@
       destroyTree(tree);
 
+      (*_blossom_data)[blossom].status = UNMATCHED;
       (*_blossom_data)[blossom].base = node;
-      (*_blossom_data)[blossom].next = INVALID;
-    }
+      matchedToUnmatched(blossom);
+    }
+
 
     void augmentOnEdge(const Edge& edge) {
@@ -1217,29 +1325,24 @@
       int right = _blossom_set->find(_graph.v(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);
+      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 = _graph.direct(edge, true);
       (*_blossom_data)[right].next = _graph.direct(edge, false);
-    }
-
-    void augmentOnArc(const Arc& arc) {
-
-      int left = _blossom_set->find(_graph.source(arc));
-      int right = _blossom_set->find(_graph.target(arc));
-
-      (*_blossom_data)[left].status = MATCHED;
-
-      int right_tree = _tree_set->find(right);
-      alternatePath(right, right_tree);
-      destroyTree(right_tree);
-
-      (*_blossom_data)[left].next = arc;
-      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
     }
 
@@ -1446,5 +1549,5 @@
           (*_blossom_data)[sb].pred = pred;
           (*_blossom_data)[sb].next =
-            _graph.oppositeArc((*_blossom_data)[tb].next);
+                           _graph.oppositeArc((*_blossom_data)[tb].next);
 
           pred = (*_blossom_data)[ub].next;
@@ -1546,5 +1649,5 @@
 
       for (int i = 0; i < int(blossoms.size()); ++i) {
-        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
+        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
 
           Value offset = (*_blossom_data)[blossoms[i]].offset;
@@ -1584,14 +1687,8 @@
         _delta4_index(0), _delta4(0),
 
-        _delta_sum(), _unmatched(0),
-
-        _fractional(0)
-    {}
+        _delta_sum() {}
 
     ~MaxWeightedMatching() {
       destroyStructures();
-      if (_fractional) {
-        delete _fractional;
-      }
     }
 
@@ -1624,6 +1721,4 @@
         (*_delta4_index)[i] = _delta4->PRE_HEAP;
       }
-
-      _unmatched = _node_num;
 
       _delta1->clear();
@@ -1670,159 +1765,9 @@
     }
 
-    /// \brief Initialize the algorithm with fractional matching
-    ///
-    /// This function initializes the algorithm with a fractional
-    /// matching. This initialization is also called jumpstart heuristic.
-    void fractionalInit() {
-      createStructures();
-
-      _blossom_node_list.clear();
-      _blossom_potential.clear();
-
-      if (_fractional == 0) {
-        _fractional = new FractionalMatching(_graph, _weight, false);
-      }
-      _fractional->run();
-
-      for (ArcIt e(_graph); e != INVALID; ++e) {
-        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
-      }
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        (*_delta1_index)[n] = _delta1->PRE_HEAP;
-      }
-      for (EdgeIt e(_graph); e != INVALID; ++e) {
-        (*_delta3_index)[e] = _delta3->PRE_HEAP;
-      }
-      for (int i = 0; i < _blossom_num; ++i) {
-        (*_delta2_index)[i] = _delta2->PRE_HEAP;
-        (*_delta4_index)[i] = _delta4->PRE_HEAP;
-      }
-
-      _unmatched = 0;
-
-      _delta1->clear();
-      _delta2->clear();
-      _delta3->clear();
-      _delta4->clear();
-      _blossom_set->clear();
-      _tree_set->clear();
-
-      int index = 0;
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        Value pot = _fractional->nodeValue(n);
-        (*_node_index)[n] = index;
-        (*_node_data)[index].pot = pot;
-        (*_node_data)[index].heap_index.clear();
-        (*_node_data)[index].heap.clear();
-        int blossom =
-          _blossom_set->insert(n, std::numeric_limits<Value>::max());
-
-        (*_blossom_data)[blossom].status = MATCHED;
-        (*_blossom_data)[blossom].pred = INVALID;
-        (*_blossom_data)[blossom].next = _fractional->matching(n);
-        if (_fractional->matching(n) == INVALID) {
-          (*_blossom_data)[blossom].base = n;
-        }
-        (*_blossom_data)[blossom].pot = 0;
-        (*_blossom_data)[blossom].offset = 0;
-        ++index;
-      }
-
-      typename Graph::template NodeMap<bool> processed(_graph, false);
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        if (processed[n]) continue;
-        processed[n] = true;
-        if (_fractional->matching(n) == INVALID) continue;
-        int num = 1;
-        Node v = _graph.target(_fractional->matching(n));
-        while (n != v) {
-          processed[v] = true;
-          v = _graph.target(_fractional->matching(v));
-          ++num;
-        }
-
-        if (num % 2 == 1) {
-          std::vector<int> subblossoms(num);
-
-          subblossoms[--num] = _blossom_set->find(n);
-          _delta1->push(n, _fractional->nodeValue(n));
-          v = _graph.target(_fractional->matching(n));
-          while (n != v) {
-            subblossoms[--num] = _blossom_set->find(v);
-            _delta1->push(v, _fractional->nodeValue(v));
-            v = _graph.target(_fractional->matching(v));
-          }
-
-          int surface =
-            _blossom_set->join(subblossoms.begin(), subblossoms.end());
-          (*_blossom_data)[surface].status = EVEN;
-          (*_blossom_data)[surface].pred = INVALID;
-          (*_blossom_data)[surface].next = INVALID;
-          (*_blossom_data)[surface].pot = 0;
-          (*_blossom_data)[surface].offset = 0;
-
-          _tree_set->insert(surface);
-          ++_unmatched;
-        }
-      }
-
-      for (EdgeIt e(_graph); e != INVALID; ++e) {
-        int si = (*_node_index)[_graph.u(e)];
-        int sb = _blossom_set->find(_graph.u(e));
-        int ti = (*_node_index)[_graph.v(e)];
-        int tb = _blossom_set->find(_graph.v(e));
-        if ((*_blossom_data)[sb].status == EVEN &&
-            (*_blossom_data)[tb].status == EVEN && sb != tb) {
-          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
-                            dualScale * _weight[e]) / 2);
-        }
-      }
-
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        int nb = _blossom_set->find(n);
-        if ((*_blossom_data)[nb].status != MATCHED) continue;
-        int ni = (*_node_index)[n];
-
-        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
-          Node v = _graph.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) {
-
-            int vt = _tree_set->find(vb);
-
-            typename std::map<int, Arc>::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, e);
-                (*_node_data)[ni].heap.decrease(e, rw);
-                it->second = e;
-              }
-            } else {
-              (*_node_data)[ni].heap.push(e, rw);
-              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
-            }
-          }
-        }
-
-        if (!(*_node_data)[ni].heap.empty()) {
-          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
-          _delta2->push(nb, _blossom_set->classPrio(nb));
-        }
-      }
-    }
-
     /// \brief Start the algorithm
     ///
     /// This function starts the algorithm.
     ///
-    /// \pre \ref init() or \ref fractionalInit() must be called
-    /// before using this function.
+    /// \pre \ref init() must be called before using this function.
     void start() {
       enum OpType {
@@ -1830,5 +1775,6 @@
       };
 
-      while (_unmatched > 0) {
+      int unmatched = _node_num;
+      while (unmatched > 0) {
         Value d1 = !_delta1->empty() ?
           _delta1->prio() : std::numeric_limits<Value>::max();
@@ -1843,8 +1789,9 @@
           _delta4->prio() : std::numeric_limits<Value>::max();
 
-        _delta_sum = d3; OpType ot = D3;
-        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
+        _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) {
@@ -1853,5 +1800,5 @@
             Node n = _delta1->top();
             unmatchNode(n);
-            --_unmatched;
+            --unmatched;
           }
           break;
@@ -1860,11 +1807,6 @@
             int blossom = _delta2->top();
             Node n = _blossom_set->classTop(blossom);
-            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
-            if ((*_blossom_data)[blossom].next == INVALID) {
-              augmentOnArc(a);
-              --_unmatched;
-            } else {
-              extendOnArc(a);
-            }
+            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
+            extendOnArc(e);
           }
           break;
@@ -1879,6 +1821,18 @@
               _delta3->pop();
             } else {
-              int left_tree = _tree_set->find(left_blossom);
-              int right_tree = _tree_set->find(right_blossom);
+              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) {
@@ -1886,5 +1840,5 @@
               } else {
                 augmentOnEdge(e);
-                _unmatched -= 2;
+                unmatched -= 2;
               }
             }
@@ -1904,9 +1858,9 @@
     /// \note mwm.run() is just a shortcut of the following code.
     /// \code
-    ///   mwm.fractionalInit();
+    ///   mwm.init();
     ///   mwm.start();
     /// \endcode
     void run() {
-      fractionalInit();
+      init();
       start();
     }
@@ -1934,5 +1888,5 @@
         }
       }
-      return sum / 2;
+      return sum /= 2;
     }
 
@@ -2268,9 +2222,4 @@
 
     Value _delta_sum;
-    int _unmatched;
-
-    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap>
-    FractionalMatching;
-    FractionalMatching *_fractional;
 
     void createStructures() {
@@ -2334,4 +2283,7 @@
 
     void destroyStructures() {
+      _node_num = countNodes(_graph);
+      _blossom_num = _node_num * 3 / 2;
+
       if (_matching) {
         delete _matching;
@@ -3006,14 +2958,8 @@
         _delta4_index(0), _delta4(0),
 
-        _delta_sum(), _unmatched(0),
-
-        _fractional(0)
-    {}
+        _delta_sum() {}
 
     ~MaxWeightedPerfectMatching() {
       destroyStructures();
-      if (_fractional) {
-        delete _fractional;
-      }
     }
 
@@ -3043,6 +2989,4 @@
         (*_delta4_index)[i] = _delta4->PRE_HEAP;
       }
-
-      _unmatched = _node_num;
 
       _delta2->clear();
@@ -3087,153 +3031,9 @@
     }
 
-    /// \brief Initialize the algorithm with fractional matching
-    ///
-    /// This function initializes the algorithm with a fractional
-    /// matching. This initialization is also called jumpstart heuristic.
-    void fractionalInit() {
-      createStructures();
-
-      _blossom_node_list.clear();
-      _blossom_potential.clear();
-
-      if (_fractional == 0) {
-        _fractional = new FractionalMatching(_graph, _weight, false);
-      }
-      if (!_fractional->run()) {
-        _unmatched = -1;
-        return;
-      }
-
-      for (ArcIt e(_graph); e != INVALID; ++e) {
-        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
-      }
-      for (EdgeIt e(_graph); e != INVALID; ++e) {
-        (*_delta3_index)[e] = _delta3->PRE_HEAP;
-      }
-      for (int i = 0; i < _blossom_num; ++i) {
-        (*_delta2_index)[i] = _delta2->PRE_HEAP;
-        (*_delta4_index)[i] = _delta4->PRE_HEAP;
-      }
-
-      _unmatched = 0;
-
-      _delta2->clear();
-      _delta3->clear();
-      _delta4->clear();
-      _blossom_set->clear();
-      _tree_set->clear();
-
-      int index = 0;
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        Value pot = _fractional->nodeValue(n);
-        (*_node_index)[n] = index;
-        (*_node_data)[index].pot = pot;
-        (*_node_data)[index].heap_index.clear();
-        (*_node_data)[index].heap.clear();
-        int blossom =
-          _blossom_set->insert(n, std::numeric_limits<Value>::max());
-
-        (*_blossom_data)[blossom].status = MATCHED;
-        (*_blossom_data)[blossom].pred = INVALID;
-        (*_blossom_data)[blossom].next = _fractional->matching(n);
-        (*_blossom_data)[blossom].pot = 0;
-        (*_blossom_data)[blossom].offset = 0;
-        ++index;
-      }
-
-      typename Graph::template NodeMap<bool> processed(_graph, false);
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        if (processed[n]) continue;
-        processed[n] = true;
-        if (_fractional->matching(n) == INVALID) continue;
-        int num = 1;
-        Node v = _graph.target(_fractional->matching(n));
-        while (n != v) {
-          processed[v] = true;
-          v = _graph.target(_fractional->matching(v));
-          ++num;
-        }
-
-        if (num % 2 == 1) {
-          std::vector<int> subblossoms(num);
-
-          subblossoms[--num] = _blossom_set->find(n);
-          v = _graph.target(_fractional->matching(n));
-          while (n != v) {
-            subblossoms[--num] = _blossom_set->find(v);
-            v = _graph.target(_fractional->matching(v));
-          }
-
-          int surface =
-            _blossom_set->join(subblossoms.begin(), subblossoms.end());
-          (*_blossom_data)[surface].status = EVEN;
-          (*_blossom_data)[surface].pred = INVALID;
-          (*_blossom_data)[surface].next = INVALID;
-          (*_blossom_data)[surface].pot = 0;
-          (*_blossom_data)[surface].offset = 0;
-
-          _tree_set->insert(surface);
-          ++_unmatched;
-        }
-      }
-
-      for (EdgeIt e(_graph); e != INVALID; ++e) {
-        int si = (*_node_index)[_graph.u(e)];
-        int sb = _blossom_set->find(_graph.u(e));
-        int ti = (*_node_index)[_graph.v(e)];
-        int tb = _blossom_set->find(_graph.v(e));
-        if ((*_blossom_data)[sb].status == EVEN &&
-            (*_blossom_data)[tb].status == EVEN && sb != tb) {
-          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
-                            dualScale * _weight[e]) / 2);
-        }
-      }
-
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        int nb = _blossom_set->find(n);
-        if ((*_blossom_data)[nb].status != MATCHED) continue;
-        int ni = (*_node_index)[n];
-
-        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
-          Node v = _graph.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) {
-
-            int vt = _tree_set->find(vb);
-
-            typename std::map<int, Arc>::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, e);
-                (*_node_data)[ni].heap.decrease(e, rw);
-                it->second = e;
-              }
-            } else {
-              (*_node_data)[ni].heap.push(e, rw);
-              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
-            }
-          }
-        }
-
-        if (!(*_node_data)[ni].heap.empty()) {
-          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
-          _delta2->push(nb, _blossom_set->classPrio(nb));
-        }
-      }
-    }
-
     /// \brief Start the algorithm
     ///
     /// This function starts the algorithm.
     ///
-    /// \pre \ref init() or \ref fractionalInit() must be called before
-    /// using this function.
+    /// \pre \ref init() must be called before using this function.
     bool start() {
       enum OpType {
@@ -3241,7 +3041,6 @@
       };
 
-      if (_unmatched == -1) return false;
-
-      while (_unmatched > 0) {
+      int unmatched = _node_num;
+      while (unmatched > 0) {
         Value d2 = !_delta2->empty() ?
           _delta2->prio() : std::numeric_limits<Value>::max();
@@ -3253,6 +3052,6 @@
           _delta4->prio() : std::numeric_limits<Value>::max();
 
-        _delta_sum = d3; OpType ot = D3;
-        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
+        _delta_sum = d2; OpType ot = D2;
+        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
 
@@ -3287,5 +3086,5 @@
               } else {
                 augmentOnEdge(e);
-                _unmatched -= 2;
+                unmatched -= 2;
               }
             }
@@ -3306,9 +3105,9 @@
     /// \note mwpm.run() is just a shortcut of the following code.
     /// \code
-    ///   mwpm.fractionalInit();
+    ///   mwpm.init();
     ///   mwpm.start();
     /// \endcode
     bool run() {
-      fractionalInit();
+      init();
       return start();
     }
@@ -3336,5 +3135,5 @@
         }
       }
-      return sum / 2;
+      return sum /= 2;
     }
 
@@ -3503,3 +3302,3 @@
 } //END OF NAMESPACE LEMON
 
-#endif //LEMON_MATCHING_H
+#endif //LEMON_MAX_MATCHING_H
