38 files changed, 7972 insertions, 65 deletions
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..dae9299
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+*.o
+*.dat
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..a64e9f3
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,24 @@
+BLOSSOM_DIR := blossom5-v2.05.src
+BLOSSOM_DIRS := $(BLOSSOM_DIR) $(BLOSSOM_DIR)/MinCost $(BLOSSOM_DIR)/GEOM
+
+BLOSSOM_SOURCES := $(filter-out $(BLOSSOM_DIR)/example.cpp, $(foreach dir, $(BLOSSOM_DIRS), $(wildcard $(dir)/*.cpp)))
+BLOSSOM_OBJS := $(patsubst %.cpp, %.o, $(BLOSSOM_SOURCES))
+
+CFLAGS := -flto -O3 -D_NDEBUG -DPERFECT_MATCHING_DOUBLE
+CXX = clang++
+LD = ld.lld
+LIBS := -lrt
+
+all: rbmp
+
+rbmp: rbmp.cpp $(BLOSSOM_DIR)/blossom5.o
+	$(CXX) $(CFLAGS) $(BLOSSOM_DIR)/blossom5.o rbmp.cpp -o $@
+
+$(BLOSSOM_DIR)/blossom5.o: ${BLOSSOM_OBJS}
+	$(LD) -r -o $@ ${BLOSSOM_OBJS}
+
+.cpp.o:
+	$(CXX) $(CFLAGS) $< -c -o $@
+
+clean:
+	rm -f ${BLOSSOM_OBJS} $(BLOSSOM_DIR)/blossom5.o rbmp
diff --git a/blossom5-v2.05.src/CHANGES.TXT b/blossom5-v2.05.src/CHANGES.TXT
new file mode 100644
index 0000000..61320a6
--- /dev/null
+++ b/blossom5-v2.05.src/CHANGES.TXT
@@ -0,0 +1,45 @@
+Changes from version 1.0 to 2.0:
+
+- Replaced Fibonacci heaps with pairing heaps 
+  (M. Fredman, R. Sedgewick, D. Sleator, R. Tarjan, Algorithmica 1(1):111-129, 1986).
+  Pairing heaps take less memory (namely, 2 pointers per edge less) compared to Fibonacci heaps,
+  and seem to be marginally faster.
+
+  Finonacci heaps are still available - replace PQ.h with PQ-Fibonacci.h .
+
+- Changed data structures so that the time in SHRINK operations
+  is O(m log n) per augmentation (I believe). This was not the case
+  in version 1.0 (see footnote 5 in the MPC paper).
+
+I re-ran experiments corresponding to tables 1,3,4,5,9 in the paper.
+The new version was marginally faster (e.g. up to 10% faster) on all examples
+except for lrb744710, where it was about 3 times faster.
+
+Changes from version 2.0 to 2.01 (thanks to Nic Schraudolph and Dmitry Kamenetsky for useful suggestions):
+
+- Fixed bug in block.h (replaced "new char[] ... delete ..." with "new char[] ... delete[] ...";
+in the first case the behavious is not specified, though most compilers would compile it correctly.)
+
+- Removed PQ-Fibonacci.h
+
+- Added disclaimer about using floating point numbers
+
+Changes from version 2.01 to 2.02:
+
+- Tweaks to stop compiler warnings,
+  changed "delete rev_mapping" to "delete [] rev_mapping" in misc.cpp (thanks to Nic Schraudolph for suggestions)
+
+- Added a statement to license conditions
+
+Changes from version 2.02 to 2.03:
+
+- Fixed a bug: if the number of expands was too large (more than node_num/4) then the previous version could have crashed.
+This was more likely to happen when dynamic updates were used (with multiple calls to PerfectMatching::Solve()).
+
+Changes from version 2.03 to 2.04:
+
+- Changes to make it compile under OS X and with latest gcc compilers, such as gcc 4.7.0 or 4.8.0
+(thanks to Adam Whiteside and Michael Cook for suggesting these changes).
+
+
+
diff --git a/blossom5-v2.05.src/GEOM/GPMinit.cpp b/blossom5-v2.05.src/GEOM/GPMinit.cpp
new file mode 100644
index 0000000..0e31231
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GPMinit.cpp
@@ -0,0 +1,200 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "GeomPerfectMatching.h"
+#include "GPMkdtree.h"
+
+
+
+
+// greedy procedure to make sure that a perfect matching exists:
+// 1. construct a matching among existing edges
+//    (greedy procedure: pick a node, check whether there are edges leading
+//     to unmatched nodes, if there are pick the edge with the smallest length).
+// 2. take remaining unmatched nodes, construct kd-tree for them,
+//    assign an ordering to nodes (last visited time during left-most depth-first search),
+//    add edges between consecutive nodes (2*i,2*i+1)
+void GeomPerfectMatching::CompleteInitialMatching()
+{
+	if (options.verbose) printf("adding edges to make sure that a perfect matching exists...");
+	PointId p, q;
+	Edge* e;
+	double len, len_min;
+	int unmatched_num = 0, edge_num0 = edge_num;
+
+	// construct greedy matching
+	for (p=0; p<node_num; p++)
+	{
+		if (nodes[p].is_marked) continue;
+		q = -1;
+		for (e=nodes[p].first[0]; e; e=e->next[0])
+		{
+			if (nodes[e->head[0]].is_marked) continue;
+			len = Dist2(p, e->head[0]);
+			if (q < 0 || len_min > len)
+			{
+				q = e->head[0];
+				len_min = len;
+			}
+		}
+		if (q >= 0)
+		{
+			nodes[p].is_marked = nodes[q].is_marked = 1;
+		}
+		else unmatched_num ++;
+	}
+
+	if (unmatched_num == 0)
+	{
+		for (p=0; p<node_num; p++) nodes[p].is_marked = 0;
+		return;
+	}
+
+	//printf("%d unmatched\n", unmatched_num);
+
+	REAL* unmatched_coords = new REAL[unmatched_num*DIM];
+	int* rev_mapping = new int[unmatched_num];
+	unmatched_num = 0;
+	for (p=0; p<node_num; p++)
+	{
+		if (nodes[p].is_marked) nodes[p].is_marked = 0;
+		else
+		{
+			memcpy(unmatched_coords+unmatched_num*DIM, coords+p*DIM, DIM*sizeof(REAL));
+			rev_mapping[unmatched_num ++] = p;
+		}
+	}
+
+	GPMKDTree* kd_tree = new GPMKDTree(DIM, unmatched_num, unmatched_coords, this);
+	kd_tree->AddPerfectMatching(rev_mapping);
+
+	delete kd_tree;
+	delete [] unmatched_coords;
+	delete [] rev_mapping;
+
+	if (options.verbose) printf("done (%d edges)\n", edge_num-edge_num0);
+}
+
+void GeomPerfectMatching::InitKNN(int K)
+{
+	if (node_num != node_num_max) { printf("InitKNN() cannot be called before all points have been added!\n"); exit(1); }
+	if (options.verbose) printf("adding K nearest neighbors (K=%d)\n", K);
+
+	int dir, k;
+	PointId p;
+	Edge* e;
+
+	if (K > node_num - 1) K = node_num - 1;
+
+	GPMKDTree* kd_tree = new GPMKDTree(DIM, node_num, coords, this);
+	PointId* neighbors = new PointId[K];
+
+	for (p=0; p<node_num; p++)
+	{
+		for (dir=0; dir<2; dir++)
+		for (e=nodes[p].first[dir]; e; e=e->next[dir])
+		{
+			nodes[e->head[dir]].is_marked = 1;
+		}
+
+		kd_tree->ComputeKNN(p, K, neighbors);
+		for (k=0; k<K; k++)
+		{
+			if (nodes[neighbors[k]].is_marked) continue;
+			AddInitialEdge(p, neighbors[k]);
+			nodes[neighbors[k]].is_marked = 1;
+		}
+
+		for (dir=0; dir<2; dir++)
+		for (e=nodes[p].first[dir]; e; e=e->next[dir])
+		{
+			nodes[e->head[dir]].is_marked = 0;
+		}
+	}
+
+	delete kd_tree;
+	delete [] neighbors;
+}
+
+#ifdef DELAUNAY_TRIANGLE
+
+#ifdef _MSC_VER
+#pragma warning(disable: 4311)
+#pragma warning(disable: 4312)
+#endif
+
+extern "C" {
+#define ANSI_DECLARATORS
+#define TRILIBRARY
+#define NO_TIMER
+#define main NO_MAIN_FUNCTION
+#include "../triangle/triangle.c"
+}
+
+void GeomPerfectMatching::InitDelaunay()
+{
+	if (node_num < 16) return;
+	if (options.verbose) printf("adding edges in Delaunay triangulation\n");
+
+	int k;
+
+	struct triangulateio in, out, vorout;
+	in.numberofpoints = node_num;
+	in.numberofpointattributes = 0;
+	in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL));
+	for (k=0; k<2*node_num; k++) in.pointlist[k] = coords[k];
+	in.pointattributelist = NULL;
+	in.pointmarkerlist = NULL;
+	in.numberofsegments = 0;
+	in.numberofholes = 0;
+	in.numberofregions = 0;
+	in.regionlist = 0;
+
+	out.pointlist = (REAL *) NULL;
+	out.pointattributelist = (REAL *) NULL;
+	out.pointmarkerlist = (int *) NULL;
+	out.trianglelist = (int *) NULL;
+	out.triangleattributelist = (REAL *) NULL;
+	out.neighborlist = (int *) NULL;
+	out.segmentlist = (int *) NULL;
+	out.segmentmarkerlist = (int *) NULL;
+	out.edgelist = (int *) NULL;
+	out.edgemarkerlist = (int *) NULL;
+
+	vorout.pointlist = (REAL *) NULL;
+	vorout.pointattributelist = (REAL *) NULL;
+	vorout.edgelist = (int *) NULL;
+	vorout.normlist = (REAL *) NULL;
+
+	triangulate("pczAevn", &in, &out, &vorout);
+
+	free(in.pointlist);
+	free(out.pointlist);
+	free(out.pointmarkerlist);
+	free(out.trianglelist);
+	free(out.neighborlist);
+	free(out.segmentlist);
+	free(out.segmentmarkerlist);
+	free(out.edgemarkerlist);
+	free(vorout.pointlist);
+	free(vorout.pointattributelist);
+	free(vorout.edgelist);
+	free(vorout.normlist);
+
+	for (k=0; k<out.numberofedges; k++) AddInitialEdge(out.edgelist[2*k], out.edgelist[2*k+1]);
+
+	free(out.edgelist);
+}
+
+
+
+#else
+
+void GeomPerfectMatching::InitDelaunay()
+{
+	printf("You need to download the 'Triangle' software from \n\thttp://www.cs.cmu.edu/~quake/triangle.html ,\nextract it to the directory GeomPerfectMatching and define DELAUNARY_TRIANG in GeomPerfectMatching.h\n");
+	exit(1);
+}
+
+#endif
+
diff --git a/blossom5-v2.05.src/GEOM/GPMinterface.cpp b/blossom5-v2.05.src/GEOM/GPMinterface.cpp
new file mode 100644
index 0000000..97a5834
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GPMinterface.cpp
@@ -0,0 +1,84 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "GeomPerfectMatching.h"
+#include "GPMkdtree.h"
+
+
+GeomPerfectMatching::GeomPerfectMatching(int nodeNum, int _DIM)
+	: DIM(_DIM),
+	  node_num(0),
+	  node_num_max(nodeNum),
+	  edge_num(0)
+{
+	if (node_num_max < 1) { printf("too few nodes\n"); exit(1); }
+	if (node_num_max & 1) { printf("# of points is odd: perfect matching cannot exist\n"); exit(1); }
+	nodes = (Node*) malloc(node_num_max*sizeof(Node));
+	memset(nodes, 0, node_num_max*sizeof(Node));
+	edges = new Block<Edge>(512);
+	coords = (REAL*) malloc((DIM+1)*node_num_max*sizeof(REAL));
+	sums = coords + DIM*node_num_max;
+	matching = (int*) malloc(node_num_max*sizeof(int));
+	int i;
+	for (i=0; i<node_num_max; i++) matching[i] = i;
+}
+
+GeomPerfectMatching::~GeomPerfectMatching()
+{
+	free(nodes);
+	delete edges;
+	free(coords);
+	free(matching);
+}
+
+GeomPerfectMatching::PointId GeomPerfectMatching::AddPoint(REAL* coord)
+{
+	if (node_num >= node_num_max)
+	{
+		printf("Error: you are trying to add too many points!\n");
+		exit(1);
+	}
+	memcpy(coords+DIM*node_num, coord, DIM*sizeof(REAL));
+	return node_num ++;
+}
+
+void GeomPerfectMatching::AddInitialEdge(PointId _i, PointId _j)
+{
+	assert(_i>=0 && _i<node_num_max && _j>=0 && _j<node_num_max && _i!=_j);
+	if (_j < _i) { int _k = _i; _i = _j; _j = _k; }
+	Node* i = nodes + _i;
+	Node* j = nodes + _j;
+	Edge* e = edges->New();
+	edge_num ++;
+
+	e->head[1] = _i;
+	e->head[0] = _j;
+	e->next[0] = i->first[0];
+	e->next[1] = j->first[1];
+	i->first[0] = e;
+	j->first[1] = e;
+}
+
+
+GeomPerfectMatching::REAL GeomPerfectMatching::ComputeCost(PointId* matching)
+{
+	if (node_num != node_num_max) { printf("ComputeCost() cannot be called before all points have been added!\n"); exit(1); }
+
+	REAL cost = 0;
+	int i;
+	for (i=0; i<node_num; i++)
+	{
+		if (matching[i]==i || matching[i]<0 || matching[i]>=node_num || matching[matching[i]]!=i)
+		{
+			printf("ComputeCost(): not a valid matching!\n");
+			exit(1);
+		}
+		if (matching[i] > i)
+		{
+			cost += Dist(i, matching[i]);
+		}
+	}
+	return cost;
+}
+
+
diff --git a/blossom5-v2.05.src/GEOM/GPMkdtree.cpp b/blossom5-v2.05.src/GEOM/GPMkdtree.cpp
new file mode 100644
index 0000000..f8cd17b
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GPMkdtree.cpp
@@ -0,0 +1,688 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "GPMkdtree.h"
+#include "../timer.h"
+
+
+// 'mapping' must be of size 2*N. (array[N], ... array[2*N-1] is used as a temporary buffer).
+// After the call  array[mapping[0]] <= array[mapping[1]] <= ... <= array[mapping[N-1]].
+// array is not modified.
+template <typename Type> inline void sort(Type* array, int array_skip, int* mapping, int N)
+{
+	// mergesort
+	int i;
+	int* mappingSrc = mapping;
+	int* mappingDst = mapping + N;
+	int* pSrc1;
+	int* pSrc2;
+	int* pSrc1End;
+	int* pSrc2End;
+	int* pDst;
+
+	for (i=0; i<(N&(~1)); i+=2)
+	{
+		if (array[array_skip*i] < array[array_skip*(i+1)])
+		{
+			mappingSrc[i]   = i;
+			mappingSrc[i+1] = i+1;
+		}
+		else
+		{
+			mappingSrc[i]   = i+1;
+			mappingSrc[i+1] = i;
+		}
+	}
+	if (i != N) mappingSrc[i] = i;
+
+	int step;
+	for (step=2; step<N; step*=2)
+	{
+		pSrc2End = mappingSrc;
+		pDst = mappingDst;
+		while ( 1 )
+		{
+			pSrc1 = pSrc2End;
+			pSrc1End = pSrc1 + step;
+			if (pSrc1End >= mappingSrc + N)
+			{
+				memcpy(pDst, pSrc1, (int)((char*)(mappingSrc + N) - (char*)pSrc1));
+				break;
+			}
+			pSrc2 = pSrc1End;
+			pSrc2End = pSrc2 + step;
+			if (pSrc2End > mappingSrc + N) pSrc2End = mappingSrc + N;
+			while ( 1 )
+			{
+				if (array[(array_skip)*(*pSrc1)] < array[array_skip*(*pSrc2)])
+				{
+					*pDst ++ = *pSrc1 ++;
+					if (pSrc1 == pSrc1End)
+					{
+						memcpy(pDst, pSrc2, (int)((char*)pSrc2End - (char*)pSrc2));
+						pDst = (int*) ((char*)pDst + (int)((char*)pSrc2End - (char*)pSrc2));
+						break;
+					}
+				}
+				else
+				{
+					*pDst ++ = *pSrc2 ++;
+					if (pSrc2 == pSrc2End)
+					{
+						memcpy(pDst, pSrc1, (int)((char*)pSrc1End - (char*)pSrc1));
+						pDst = (int*) ((char*)pDst + (int)((char*)pSrc1End - (char*)pSrc1));
+						break;
+					}
+				}
+			}
+		}
+		pDst = mappingDst;
+		mappingDst = mappingSrc;
+		mappingSrc = pDst;
+	}
+	if (mappingSrc != mapping) memcpy(mapping, mappingSrc, N*sizeof(int));
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+#define NEIGHBOR_PARENT(k) (((k)-1)>>1)
+#define NEIGHBOR_FIRST_CHILD(k) (((k)<<1)+1)
+
+Neighbors::Neighbors()
+{
+	K_max = 0;
+	dist_array = NULL;
+}
+
+Neighbors::~Neighbors()
+{
+	if (dist_array) delete [] dist_array;
+}
+
+void Neighbors::Init(int _K, PointId* _array)
+{
+	K = _K;
+	array = _array;
+	num = 0;
+	if (K > K_max)
+	{
+		if (dist_array) delete [] dist_array;
+		K_max = K;
+		dist_array = new double[K_max];
+	}
+}
+
+inline void Neighbors::Swap(int k1, int k2)
+{
+	PointId p = array[k1]; array[k1] = array[k2]; array[k2] = p;
+	double d = dist_array[k1]; dist_array[k1] = dist_array[k2]; dist_array[k2] = d;
+}
+
+inline void Neighbors::Add(PointId p, double dist)
+{
+	int k;
+	if (num < K)
+	{
+		k = num ++;
+		array[k] = p;
+		dist_array[k] = dist;
+		while ( k > 0 )
+		{
+			int k_parent = NEIGHBOR_PARENT(k);
+			if (dist_array[k] <= dist_array[k_parent]) break;
+			Swap(k, k_parent);
+			k = k_parent;
+		}
+	}
+	else
+	{
+		if (dist_array[0] <= dist) return;
+		array[0] = p;
+		dist_array[0] = dist;
+		k = 0;
+		while ( 1 )
+		{
+			int k_child = NEIGHBOR_FIRST_CHILD(k);
+			if (k_child >= K) break;
+			if (k_child+1 < K && dist_array[k_child+1] > dist_array[k_child]) k_child ++;
+			if (dist_array[k] >= dist_array[k_child]) break;
+			Swap(k, k_child);
+			k = k_child;
+		}
+	}
+	//for (k=1; k<num; k++) assert(dist_array[k] <= dist_array[NEIGHBOR_PARENT(k)]);
+}
+
+inline double Neighbors::GetMax()
+{
+	assert(num > 0);
+	return dist_array[0];
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+GPMKDTree::GPMKDTree(int _D, int _point_num, REAL* coords, GeomPerfectMatching* _GPM)
+	: D(_D), DIM(_GPM->DIM), point_num(_point_num), GPM(_GPM)
+{
+	Node* i;
+	Node* j;
+	int d, d0, k;
+	int* mapping = new int[(D+2)*point_num];
+	int* buf = mapping + D*point_num;
+	int* marking = buf + point_num;
+	memset(marking, 0, point_num*sizeof(int));
+	int* ptr = mapping;
+
+	int node_num_max = 4*point_num/3+2;
+	nodes = (Node*)malloc(node_num_max*sizeof(Node));
+	rev_mapping = (Node**)malloc(point_num*sizeof(Node*));
+	memset(rev_mapping, 0, point_num*sizeof(Node*));
+
+	REAL** coords_array = new REAL*[D];
+	int* skip_array = new int[D];
+	for (d=0; d<DIM; d++) { coords_array[d] = coords + d; skip_array[d] = DIM; }
+	if (d < D)            { coords_array[d] = GPM->sums; skip_array[d] = 1; }
+
+	for (d=0; d<D; d++)
+	{
+		sort<REAL>(coords_array[d], skip_array[d], ptr, point_num);
+		if (d == DIM) sum_max = GPM->sums[ptr[point_num-1]];
+		ptr += point_num;
+	}
+
+	nodes[0].parent = NULL;
+	nodes[0].order = 0;
+	nodes[0].d = point_num;
+	nodes[0].first_child = NULL;
+	node_num = 1;
+	Node* first_unprocessed = &nodes[0];
+	while ( (i=first_unprocessed) )
+	{
+		first_unprocessed = i->first_child;
+
+		int start = i->order;
+		int num0 = i->d, num;
+		if ((DIM==D && num0<=2) || (DIM<D && num0 <= CHILDREN_MAX))
+		{
+			// leaf
+			i->d = -num0;
+			for (k=0; k<num0; k++)
+			{
+				i->points[k] = mapping[start+k];
+				rev_mapping[mapping[start+k]] = i;
+			}
+			continue;
+		}
+
+		// not a leaf.
+		if (node_num + 2 > node_num_max)
+		{
+			node_num_max = 3*node_num_max/2 + 16;
+			Node* nodes_old = nodes;
+			nodes = (Node*)realloc(nodes, node_num_max*sizeof(Node));
+#define UPDATE_NODE_PTR(ptr) ptr = (Node*)((char*)ptr + ((char*)nodes-(char*)nodes_old))
+			UPDATE_NODE_PTR(i);
+			if (first_unprocessed) UPDATE_NODE_PTR(first_unprocessed);
+			for (k=0; k<node_num; k++)
+			{
+				if (nodes[k].parent)                       UPDATE_NODE_PTR(nodes[k].parent);
+				if (nodes[k].d>=0 && nodes[k].first_child) UPDATE_NODE_PTR(nodes[k].first_child);
+			}
+			for (k=0; k<point_num; k++) if (rev_mapping[k]) UPDATE_NODE_PTR(rev_mapping[k]);
+		}
+
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		// choose split dimension d0 and number of items 'num' in the first group.
+
+		//d0 = (i->parent) ? ((i->parent->d + 1) % D) : 0;
+		//num = num0/2;
+
+		const int FRACTION = 20;
+		int num_min = 1;      if (num_min < num0/FRACTION) num_min = num0/FRACTION;
+		int num_max = num0-1; if (num_max > (FRACTION-1)*num0/FRACTION) num_max = (FRACTION-1)*num0/FRACTION;
+		int num_max_DIM = num0-1;
+
+		if (D>DIM && (!i->parent || !i->parent->parent || !i->parent->parent->parent)) d0 = D-1;
+		else
+		{
+			d0 = -1;
+			REAL diff_max = 0;
+			for (d=0; d<D; d++)
+			{
+				ptr = mapping + d*point_num;
+
+				REAL c_min, c_max;
+				c_min = coords_array[d][ptr[start +          num_min             ]*skip_array[d]];
+				c_max = coords_array[d][ptr[start + ((d<DIM)?num_max:num_max_DIM)]*skip_array[d]];
+				REAL diff = c_max - c_min;
+				//if (d == DIM) diff *= 2;
+				if (d0<0 || diff_max < diff) { d0 = d; diff_max = diff; }
+			}
+		}
+
+		ptr = mapping + d0*point_num;
+		REAL c_min, c_max;
+		if (d < DIM)
+		{
+			c_min = coords_array[d0][ptr[start+num_min]*skip_array[d0]];
+			c_max = coords_array[d0][ptr[start+num_max]*skip_array[d0]];
+		}
+		else
+		{
+			c_min = coords_array[d0][ptr[start+1]*skip_array[d0]];
+			c_max = coords_array[d0][ptr[start+num0-1]*skip_array[d0]]; // for this dimension upper tail is important, so don't discard it!
+		}
+		REAL split = (c_min + c_max) / 2;
+		for (num=num_min; num<num_max; num++)
+		{
+			if (coords_array[d0][ptr[start+num]*skip_array[d0]] > split) break;
+		}
+
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+
+		ptr = mapping + d0*point_num;
+		for (k=start; k<start+num; k++) marking[ptr[k]] = 1;
+		for (d=0; d<D; d++)
+		{
+			if (d == d0) continue;
+			ptr = mapping + d*point_num;
+			int* bufa = buf;
+			int* bufb = buf + num;
+			for (k=start; k<start+num0; k++)
+			{
+				if (marking[ptr[k]]) *bufa ++ = ptr[k];
+				else                 *bufb ++ = ptr[k];
+			}
+			memcpy(ptr+start, buf, num0*sizeof(int));
+		}
+		ptr = mapping + d0*point_num;
+		for (k=start; k<start+num; k++) marking[ptr[k]] = 0;
+
+		i->d = d0;
+		PointId p = ptr[start + num - ((num >= num0-num) ? 1 : 0)];
+		i->coord = coords_array[d0][p*skip_array[d0]];
+		i->first_child = j = &nodes[node_num];
+		node_num += 2;
+		j->parent = (j+1)->parent = i;
+		j->order = start;
+		j->d = num;
+		(j+1)->order = start+num;
+		(j+1)->d = num0-num;
+		(j+1)->first_child = first_unprocessed;
+		j->first_child = j+1;
+		first_unprocessed = j;
+	}
+
+	delete [] coords_array;
+	delete [] skip_array;
+	delete [] mapping;
+
+	////////////////////////////////////////////////////////////////////////////////////
+	// set ordering and depth_max
+	int depth = 0, depth_max = 0;
+
+	// set ordering
+	i = &nodes[0];
+	k = 0;
+	if (D > DIM)
+	{
+		// ordering for AddNegativeEdges() - tree preorder
+		while ( 1 )
+		{
+			if (!IS_LEAF(i))
+			{
+				i = i->first_child;
+				depth ++;
+				if (depth_max < depth) depth_max = depth;
+			}
+			else
+			{
+				while ( 1 )
+				{
+					i->order = k ++;
+
+					if (!i->parent) break;
+					if (i->parent->first_child == i) { i ++; break; }
+					i = i->parent;
+					depth --;
+				}
+				if (!i->parent) break;
+			}
+		}
+	}
+	else
+	{
+		// compute tree inorder - useful for nearest neighbor search so that that branches close to the input node are explored first
+		while ( 1 )
+		{
+			if (!IS_LEAF(i))
+			{
+				i = i->first_child;
+				depth ++;
+				if (depth_max < depth) depth_max = depth;
+			}
+			else
+			{
+				i->order = k ++;
+				while ( 1 )
+				{
+					if (!i->parent) break;
+					if (i->parent->first_child == i) { i->parent->order = k ++; i ++; break; }
+					i = i->parent;
+					depth --;
+				}
+				if (!i->parent) break;
+			}
+		}
+	}
+
+	traversing_buf = (REAL*) malloc((D + depth_max + 2)*sizeof(REAL));
+}
+
+GPMKDTree::~GPMKDTree()
+{
+	free(nodes);
+	free(rev_mapping);
+	free(traversing_buf);
+}
+
+void GPMKDTree::AddPerfectMatching(PointId* rev_mapping)
+{
+	Node* i;
+	int k;
+	PointId p, q = -1;
+	i = &nodes[0];
+	do
+	{
+		if (IS_LEAF(i))
+		{
+			for (k=0; k<-i->d; k++)
+			{
+				p = i->points[k];
+				if (q < 0) q = p;
+				else { GPM->AddInitialEdge(rev_mapping[p], rev_mapping[q]); q = -1; }
+			}
+		}
+		else
+		{
+			i = i->first_child;
+			continue;
+		}
+
+		while ( i->parent )
+		{
+			if (i->parent->first_child == i) { i ++; break; }
+			i = i->parent;
+		}
+	} while (i->parent);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+#define MOVE_DOWN_LEFT(i)\
+	{\
+		*stack ++ = current_diff[i->d];\
+		if (current_diff[i->d] <= 0 && (diff = i->coord - coord0[i->d]) < 0) current_diff[i->d] = diff;\
+		i = i->first_child;\
+	}
+
+#define MOVE_DOWN_RIGHT(i)\
+	{\
+		*stack ++ = current_diff[i->d];\
+		if (current_diff[i->d] >= 0 && (diff = i->coord - coord0[i->d]) > 0) current_diff[i->d] = diff;\
+		i = i->first_child+1;\
+	}
+
+#define MOVE_LEFT(i)\
+	{\
+		int d_prev = i->parent->d;\
+		current_diff[d_prev] = stack[-1];\
+		if (current_diff[d_prev] <= 0 && (diff = i->parent->coord - coord0[d_prev]) < 0) current_diff[d_prev] = diff;\
+		i --;\
+	}
+
+#define MOVE_RIGHT(i)\
+	{\
+		int d_prev = i->parent->d;\
+		current_diff[d_prev] = stack[-1];\
+		if (current_diff[d_prev] >= 0 && (diff = i->parent->coord - coord0[d_prev]) > 0) current_diff[d_prev] = diff;\
+		i ++;\
+	}
+
+#define MOVE_UP(i)\
+	{\
+		i = i->parent;\
+		current_diff[i->d] = *(-- stack);\
+	}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+#define MOVE_DOWN_LEFT_X(i)\
+	{\
+		*stack ++ = current_diff[i->d];\
+		if (i->d < D-1) { if (current_diff[i->d] <= 0 && (diff = i->coord - coord0[i->d]) < 0) current_diff[i->d] = diff; }\
+		else current_diff[i->d] = i->coord;\
+		i = i->first_child;\
+	}
+
+#define MOVE_DOWN_RIGHT_X(i)\
+	{\
+		*stack ++ = current_diff[i->d];\
+		if (i->d < D-1) { if (current_diff[i->d] >= 0 && (diff = i->coord - coord0[i->d]) > 0) current_diff[i->d] = diff; }\
+		i = i->first_child+1;\
+	}
+
+#define MOVE_LEFT_X(i)\
+	{\
+		int d_prev = i->parent->d;\
+		current_diff[d_prev] = stack[-1];\
+		if (d_prev < D-1) { if (current_diff[d_prev] <= 0 && (diff = i->parent->coord - coord0[d_prev]) < 0) current_diff[d_prev] = diff; }\
+		else current_diff[d_prev] = i->parent->coord;\
+		i --;\
+	}
+
+#define MOVE_RIGHT_X(i)\
+	{\
+		int d_prev = i->parent->d;\
+		current_diff[d_prev] = stack[-1];\
+		if (d_prev < D-1) { if (current_diff[d_prev] >= 0 && (diff = i->parent->coord - coord0[d_prev]) > 0) current_diff[d_prev] = diff; }\
+		i ++;\
+	}
+
+#define MOVE_UP_X(i)\
+	{\
+		i = i->parent;\
+		current_diff[i->d] = *(-- stack);\
+	}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+/*
+void GPMKDTree::ComputeKNN(PointId p0, int K, PointId* neighbors_array)
+{
+	int p;
+	REAL* coord0 = GPM->coords + p0*DIM;
+
+	neighbors.Init(K, neighbors_array);
+
+	for (p=0; p<point_num; p++)
+	if (p != p0)
+	{
+		neighbors.Add(p, GPM->Dist(coord0, GPM->coords+p*DIM));
+	}
+}
+*/
+
+void GPMKDTree::ComputeKNN(PointId p0, int K, PointId* neighbors_array)
+{
+	int neighbor_num = 0;
+	Node* i = rev_mapping[p0];
+	int k, order0 = i->order;
+	REAL diff;
+	REAL* coords = GPM->coords;
+	REAL* coord0 = GPM->coords + p0*DIM;
+	REAL* current_diff = traversing_buf;
+	REAL* stack = traversing_buf + D;
+
+	neighbors.Init(K, neighbors_array);
+
+	for (k=0; k<D; k++) current_diff[k] = 0;
+	i = &nodes[0];
+	do
+	{
+		if (IS_LEAF(i))
+		{
+			for (k=0; k<-i->d; k++)
+			{
+				PointId p = i->points[k];
+				if (p == p0) continue;
+				double dist2;
+				REAL* coord2 = coords+p*DIM;
+				GPM_GET_DIST2(dist2, coord0, coord2);
+				neighbors.Add(p, dist2);
+			}
+		}
+		else
+		{
+			if (neighbors.GetNum() < K || GPM->Norm2(current_diff) < neighbors.GetMax())
+			{
+				if (i->order > order0)
+				{
+					MOVE_DOWN_LEFT(i);
+				}
+				else
+				{
+					MOVE_DOWN_RIGHT(i);
+				}
+				continue;
+			}
+		}
+
+		while ( i->parent )
+		{
+			if (i->parent->order > order0)
+			{
+				if (i->parent->first_child == i)
+				{
+					MOVE_RIGHT(i);
+					break;
+				}
+			}
+			else
+			{
+				if (i->parent->first_child != i)
+				{
+					MOVE_LEFT(i);
+					break;
+				}
+			}
+			MOVE_UP(i);
+		}
+	} while ( i->parent );
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+/*
+void GPMKDTree::AddNegativeEdges(PointId p, PerfectMatching* pm)
+{
+	PointId q;
+	for (q=p+1; q<node_num; q++)
+	{
+		if (GPM->nodes[q].is_marked) continue;
+		REAL len = GPM->Dist(p, q);
+		if (2*len - GPM->sums[p] - GPM->sums[q] < 0)
+		{
+			if (pm->AddNewEdge(p, q, len, true)>=0)
+			{
+				GPM->AddInitialEdge(p, q);
+			}
+		}
+	}
+}
+*/
+
+
+void GPMKDTree::AddNegativeEdges(PointId p0, PerfectMatching* pm)
+{
+	Node* i = rev_mapping[p0];
+	int k, order0 = i->order;
+	bool check;
+	REAL diff;
+	REAL* coords = GPM->coords;
+	REAL* coord0 = coords + p0*DIM;
+	REAL* sums = GPM->sums;
+	REAL sum0 = sums[p0];
+	REAL* current_diff = traversing_buf;
+	REAL* stack = traversing_buf + D;
+
+	for (k=0; i->points[k]!=p0; k++) {}
+	for (k++; k<-i->d; k++)
+	{
+		PointId p = i->points[k];
+		
+		REAL len = GPM->Dist(coord0, GPM->coords+p*DIM);
+		if (2*len - GPM->sums[p] < GPM->sums[p0])
+		{
+			double start_time = get_time();
+			if (pm->AddNewEdge(p0, p, len, true)>=0) GPM->AddInitialEdge(p0, p);
+			GPM->graph_update_time += get_time() - start_time;
+		}
+	}
+
+	for (k=0; k<DIM; k++) current_diff[k] = 0;
+	current_diff[k] = sum_max;
+	i = &nodes[0];
+	do
+	{
+		if (i->order > order0)
+		{
+			if (IS_LEAF(i))
+			{
+				for (k=0; k<-i->d; k++)
+				{
+					PointId p = i->points[k];
+					if (!GPM->nodes[p].is_marked)
+					{
+						//REAL* coord2 = coords+p*DIM;
+						//REAL threshold = sums[p0]+sums[p];
+						//GPM_CHECK_DIST(check, coord0, coord2, threshold);
+						//if (check)
+						{
+							REAL len = GPM->Dist(coord0, GPM->coords+p*DIM);
+							if (2*len - GPM->sums[p] < GPM->sums[p0])
+							{
+								double start_time = get_time();
+								if (pm->AddNewEdge(p0, p, len, true)>=0) GPM->AddInitialEdge(p0, p);
+								GPM->graph_update_time += get_time() - start_time;
+							}
+						}
+					}
+				}
+			}
+			else
+			{
+				REAL threshold = current_diff[D-1] + sum0;
+				GPM_CHECK_NORM(check, current_diff, threshold);
+				if (check)
+				{
+					MOVE_DOWN_LEFT_X(i);
+					continue;
+				}
+			}
+		}
+
+		while ( i->parent )
+		{
+			if (i->parent->first_child == i)
+			{
+				MOVE_RIGHT_X(i);
+				break;
+			}
+			MOVE_UP_X(i);
+		}
+	} while (i->parent);
+}
diff --git a/blossom5-v2.05.src/GEOM/GPMkdtree.h b/blossom5-v2.05.src/GEOM/GPMkdtree.h
new file mode 100644
index 0000000..67e3d91
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GPMkdtree.h
@@ -0,0 +1,105 @@
+/*
+    GPMkdtree.h - kd-tree data structure for pricing in complete geometric instances
+
+    Copyright 2008 Vladimir Kolmogorov (vnk@adastral.ucl.ac.uk)
+
+    This software can be used for research purposes only. Commercial use is prohibited.
+    Public redistribution of the code or its derivatives is prohibited.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+
+#ifndef NJKASKJTASJNBAJSNRBAJS
+#define NJKASKJTASJNBAJSNRBAJS
+
+#include "GeomPerfectMatching.h"
+
+struct Neighbors
+{
+	typedef GeomPerfectMatching::REAL REAL;
+	typedef GeomPerfectMatching::PointId PointId;
+
+	Neighbors();
+	~Neighbors();
+
+	int GetNum() { return num; }
+
+	void Init(int K, PointId* array);
+	void Add(PointId p, double dist);
+	double GetMax();
+
+private:
+	void Swap(int k1, int k2);
+	PointId* array;
+	double* dist_array;
+	int num, K, K_max;
+};
+
+
+
+struct GPMKDTree
+{
+	typedef GeomPerfectMatching::REAL REAL;
+	typedef GeomPerfectMatching::PointId PointId;
+
+	GPMKDTree(int D, int point_num, REAL* coords, GeomPerfectMatching* GPM);
+	~GPMKDTree();
+
+	// if D == DIM
+	void AddPerfectMatching(PointId* rev_mapping);
+	void ComputeKNN(PointId p, int K, PointId* neighbors);
+
+	// if D == DIM+1
+	void AddNegativeEdges(PointId p, PerfectMatching* pm);
+
+	//////////////////////////////////////////////////////////////////////////
+
+public:
+#define CHILDREN_MAX 2
+	struct Node
+	{
+		Node* parent;
+		int d; // split dimension. d<0 indicates a leaf.
+#define IS_LEAF(i) ((i)->d < 0)
+		union
+		{
+			struct // for non-leaves
+			{
+				REAL coord;
+				Node* first_child; // the second child is first_child+1
+			};
+			struct // for leaves
+			{
+				PointId points[CHILDREN_MAX]; // the number of points is -d
+			};
+		};
+		int order;
+	}* nodes;
+
+
+	//////////////////////////////////////////////////////////////////////////
+	Node** rev_mapping;
+
+	int D, DIM, point_num, node_num;
+	REAL sum_max;
+	REAL* traversing_buf;
+	GeomPerfectMatching* GPM;
+
+	Neighbors neighbors;
+};
+
+
+
+
+#endif
diff --git a/blossom5-v2.05.src/GEOM/GPMmain.cpp b/blossom5-v2.05.src/GEOM/GPMmain.cpp
new file mode 100644
index 0000000..d2bcb3c
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GPMmain.cpp
@@ -0,0 +1,148 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "GeomPerfectMatching.h"
+#include "GPMkdtree.h"
+#include "../timer.h"
+
+
+
+GeomPerfectMatching::REAL GeomPerfectMatching::SolveComplete()
+{
+	if (node_num != node_num_max) { printf("ComputeCost() cannot be called before all points have been added!\n"); exit(1); }
+
+	PointId p, q;
+	int e = 0, E = node_num*(node_num-1)/2;
+	PerfectMatching* pm = new PerfectMatching(node_num, E);
+	for (p=0; p<node_num; p++)
+	{
+		for (q=p+1; q<node_num; q++)
+		{
+			pm->AddEdge(p, q, Dist(p, q));
+		}
+	}
+	pm->options = options;
+	pm->Solve();
+	for (p=0; p<node_num; p++)
+	{
+		for (q=p+1; q<node_num; q++)
+		{
+			if (pm->GetSolution(e++))
+			{
+				matching[p] = q;
+				matching[q] = p;
+			}
+		}
+	}
+	delete pm;
+	return ComputeCost(matching);
+}
+
+GeomPerfectMatching::REAL GeomPerfectMatching::Solve()
+{
+	double start_time = get_time();
+	double perfect_matching_time = 0;
+	double negative_edges_time = 0;
+	if (options.verbose) { printf("starting geometric matching with %d points\n", node_num); fflush(stdout); }
+	PointId p, q;
+	Edge* e;
+	int _e;
+	int iter;
+	bool success = false;
+	PerfectMatching* pm = NULL;
+	GPMKDTree* kd_tree;
+
+	double init_matching_time = get_time();
+
+	if (gpm_options.init_Delaunay) InitDelaunay();
+	if (gpm_options.init_KNN > 0) InitKNN(gpm_options.init_KNN);
+	if (gpm_options.init_greedy) CompleteInitialMatching();
+
+	init_matching_time = get_time() - init_matching_time;
+
+	graph_update_time = 0;
+
+	int iter_max = gpm_options.iter_max;
+	for (iter=0; iter_max<=0 || iter<iter_max; iter++)
+	{
+		if (pm)
+		{
+			double negative_edges_start_time = get_time();
+			int edge_num0 = edge_num;
+			pm->StartUpdate();
+			for (p=0; p<node_num; p++)
+			{
+				PerfectMatching::REAL s = pm->GetTwiceSum(p);
+				if ( ((REAL)1 / 2) == 0 && ((PerfectMatching::REAL)1 / 2) != 0 ) sums[p] = (REAL)ceil((double)s);
+				else                                                             sums[p] = (REAL)s;
+			}
+			if (options.verbose) { printf("building kd_tree..."); fflush(stdout); }
+			{
+				kd_tree = new GPMKDTree(DIM+1, node_num, coords, this);
+			}
+			if (options.verbose) { printf(" done. Now adding negative edges:\n    "); fflush(stdout); }
+
+			for (p=0; p<node_num; p++)
+			{
+				if (options.verbose && (p%(node_num/72)==0)) { printf("+"); fflush(stdout); }
+				for (e=nodes[p].first[0]; e; e=e->next[0]) nodes[e->head[0]].is_marked = 1;
+				for (e=nodes[p].first[1]; e; e=e->next[1]) nodes[e->head[1]].is_marked = 1;
+				kd_tree->AddNegativeEdges(p, pm);
+				for (e=nodes[p].first[0]; e; e=e->next[0]) nodes[e->head[0]].is_marked = 0;
+				for (e=nodes[p].first[1]; e; e=e->next[1]) nodes[e->head[1]].is_marked = 0;
+			}
+			delete kd_tree;
+			//if (edge_num - edge_num0 > node_num / 32)
+			if ( 0 ) // always reuse previous computation
+			{
+				delete pm;
+				pm = NULL;
+			}
+			else
+			{
+				pm->FinishUpdate();
+				if (edge_num0 == edge_num) success = true;
+			}
+			if (options.verbose) { printf("\ndone (%d edges added)\n", edge_num-edge_num0); fflush(stdout); }
+			negative_edges_time += get_time() - negative_edges_start_time;
+		}
+		if (!pm)
+		{
+			int E = 5*node_num;
+			if (E < 5*edge_num/4) E = 5*edge_num/4;
+			pm = new PerfectMatching(node_num, E);
+			for (e=edges->ScanFirst(); e; e=edges->ScanNext())
+			{
+				p = e->head[1]; q = e->head[0];
+				pm->AddEdge(p, q, Dist(p, q));
+			}
+		}
+		if (options.verbose) printf("iter %d: ", iter+1);
+		pm->options = options;
+		double perfect_matching_start = get_time();
+		pm->Solve();
+		perfect_matching_time += get_time() - perfect_matching_start;
+		if (success) break;
+	}
+
+	for (_e=0, e=edges->ScanFirst(); e; _e++, e=edges->ScanNext())
+	{
+		if (pm->GetSolution(_e))
+		{
+			p = e->head[1]; q = e->head[0];
+			matching[p] = q;
+			matching[q] = p;
+		}
+	}
+	delete pm;
+	REAL cost = ComputeCost(matching);
+	if (options.verbose)
+	{
+		printf("geometric matching finished [%.3f secs]. cost=%.1f \n", get_time()-start_time, (double)cost); 
+		printf("    selecting initial edges: [%.3f secs], perfect matching: [%.3f secs]\n", init_matching_time, perfect_matching_time);
+		printf("    pricing: [%.3f secs] including graph updates: [%.3f secs]\n", negative_edges_time, graph_update_time); 
+		fflush(stdout);
+	}
+	return cost;
+}
+
diff --git a/blossom5-v2.05.src/GEOM/GeomPerfectMatching.h b/blossom5-v2.05.src/GEOM/GeomPerfectMatching.h
new file mode 100644
index 0000000..b54d471
--- /dev/null
+++ b/blossom5-v2.05.src/GEOM/GeomPerfectMatching.h
@@ -0,0 +1,177 @@
+/*
+    GeomPerfectMatching.h - computing min cost perfect matching in complete geometric instances (with edge weights equal to Euclidean distances)
+
+    Copyright 2008 Vladimir Kolmogorov (vnk@adastral.ucl.ac.uk)
+
+    This software can be used for research purposes only. Commercial use is prohibited.
+    Public redistribution of the code or its derivatives is prohibited.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef ASKHALKSJBRASMNABFAJSTAS
+#define ASKHALKSJBRASMNABFAJSTAS
+
+#include <assert.h>
+#include <math.h>
+#include "../PerfectMatching.h"
+
+//#define DELAUNAY_TRIANGLE
+
+struct GPMKDTree;
+
+class GeomPerfectMatching
+{
+public:
+	typedef int REAL; // if you change it to double, you should also change PerfectMatching::REAL to double!
+	typedef int PointId;
+
+	// pointNum must be a positive even number.
+	GeomPerfectMatching(int pointNum, int DIM);
+	~GeomPerfectMatching();
+
+	// Must be called exactly point_num times (before anything else can be called).
+	// First call returns 0, second returns 1, and so on.
+	// coord must be an array of size DIM. (This array is read into internal memory.)
+	PointId AddPoint(REAL* coord);
+
+
+	////////////// SOLVING //////////////
+
+	// solves perfect matching in the complete graph. Inefficient, just for testing.
+	REAL SolveComplete(); 
+
+
+	// options for Solve()
+	struct GPMOptions
+	{
+		GPMOptions() : init_Delaunay(true), init_KNN(0), init_greedy(true), iter_max(0) {}
+
+		// three variables below determine the initial subset of edges. 
+		// Delaunay initialization seems to be more robust than K nearest neighbors
+		// (but you need to download the "Triangle" package of Shewchuk
+		// from http://www.cs.cmu.edu/~quake/triangle.html , extract it to the directory 'triangle'
+		// and define DELAUNAY_TRIANGLE above). 
+		bool	init_Delaunay;  // add Delaunay triangulation edges.
+		int		init_KNN;       // use init_KNN nearest neighbors for each point.
+		bool	init_greedy;    // add edges greedily to make sure that a perfect matching exists (see comments before CompleteInitialMatching() in GPMinit.cpp for details)
+
+		int		iter_max;   // If iter_max <= 0 then adds subsets of edges until an optimal solution is found.
+				            // Otherwise runs at most iter_max iterations, so the solution may be suboptimal. 
+				            // (iter_max=1 runs perfect matching just for the initial subset).
+	};
+	struct PerfectMatching::Options options;
+	struct GPMOptions gpm_options;
+
+	REAL Solve(); 
+
+	// You can also specify the initial subset manually
+	void AddInitialEdge(PointId i, PointId j);
+
+	//////////// READING RESULTS //////////
+	// Can be called after Solve() or SolveComplete(). Returns pointer to an array 'matching' 
+	// of size node_num. (matching[i] is the point corresponding to point i).
+	// User must not modify this array.
+	PointId GetMatch(PointId p) { return matching[p]; }
+
+
+
+	REAL Dist(REAL* coord1, REAL* coord2);
+	REAL Dist(PointId p, PointId q);
+
+
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+private:
+
+friend struct GPMKDTree;
+
+	struct Edge
+	{
+		PointId		head[2];
+		Edge*		next[2];
+	};
+	struct Node
+	{
+		Edge*		first[2];
+		int			is_marked;
+	};
+
+	Node*			nodes;
+	Block<Edge>*	edges;
+	REAL*			coords; // array of size DIM*node_num_max
+	REAL*			sums; // array of size node_num_max
+	PointId*		matching; // array of size node_num_max
+
+	int DIM;
+	int node_num, node_num_max;
+	int edge_num;
+
+	double graph_update_time;
+
+// below x denotes REAL, X denotes double, c,c1,c2 are coords (REAL*)
+#define GPM_ROUND(X) (REAL)( ( ((REAL)1 / 2) == 0 ) ? ((X)+0.5) : (X) )
+#define GPM_GET_NORM2(X, c) { int d; X = 0; for (d=0; d<DIM; d++) X += ((double)(c)[d])*(c)[d]; }
+#define GPM_GET_NORM(x, c) { double X; GPM_GET_NORM2(X, c); X = sqrt(X); x = GPM_ROUND(X); }
+#define GPM_GET_DIST2(X, c1, c2) { int d; X = 0; for (d=0; d<DIM; d++) X += ((double)((c1)[d]-(c2)[d]))*((c1)[d]-(c2)[d]); }
+#define GPM_GET_DIST(x, c1, c2) { double X; GPM_GET_DIST2(X, c1, c2); X = sqrt(X); x = GPM_ROUND(X); }
+// if returns false, then 2*||coord||-treshold >= 0
+#define GPM_CHECK_NORM(result, c, x)\
+	{\
+		if (threshold <= 0) result = false;\
+		else\
+		{\
+			double X;\
+			GPM_GET_NORM2(X, c);\
+			result = (4*X < (double)x*threshold);\
+		}\
+	}
+#define GPM_CHECK_DIST(result, c1, c2, x)\
+	{\
+		if (threshold <= 0) result = false;\
+		else\
+		{\
+			double X;\
+			GPM_GET_DIST2(X, c1, c2);\
+			result = (4*X < (double)x*threshold);\
+		}\
+	}
+
+	double Norm2(REAL* coord)                { double norm2; GPM_GET_NORM2(norm2, coord); return norm2; }
+	REAL Norm(REAL* coord)                   { REAL norm;    GPM_GET_NORM (norm,  coord); return norm;  }
+	double Dist2(REAL* coord1, REAL* coord2) { double dist2; GPM_GET_DIST2(dist2, coord1, coord2); return dist2; }
+	double Dist2(PointId p, PointId q) { return Dist2(coords+DIM*p, coords+DIM*q); }
+
+	void CompleteInitialMatching(); // add edges so that a perfect matching is guaranteed to exist
+	void InitKNN(int K);
+	void InitDelaunay();
+
+	REAL ComputeCost(PointId* matching);
+};
+
+inline GeomPerfectMatching::REAL GeomPerfectMatching::Dist(REAL* coord1, REAL* coord2)
+{
+	REAL dist;
+	GPM_GET_DIST (dist,  coord1, coord2);
+	return dist; 
+}
+
+inline GeomPerfectMatching::REAL GeomPerfectMatching::Dist(PointId p, PointId q)
+{
+	return Dist(coords+DIM*p,  coords+DIM*q);
+}
+
+#endif
diff --git a/blossom5-v2.05.src/GRAPH1.TXT b/blossom5-v2.05.src/GRAPH1.TXT
new file mode 100644
index 0000000..13ec532
--- /dev/null
+++ b/blossom5-v2.05.src/GRAPH1.TXT
@@ -0,0 +1,10 @@
+6 9
+0 1 3
+0 3 10
+0 4 7
+1 2 -1
+1 4 5
+1 5 4
+2 5 -7
+3 4 0
+4 5 4
diff --git a/blossom5-v2.05.src/GRAPH2.TXT b/blossom5-v2.05.src/GRAPH2.TXT
new file mode 100644
index 0000000..b3e487f
--- /dev/null
+++ b/blossom5-v2.05.src/GRAPH2.TXT
@@ -0,0 +1,11 @@
+c example of a graph in DIMACS format
+p edge 6 9
+e 1 2 3
+e 1 4 10
+e 1 5 7
+e 2 3 -1
+e 2 5 5
+e 2 6 4
+e 3 6 -7
+e 4 5 0
+e 5 6 4
diff --git a/blossom5-v2.05.src/LCA.h b/blossom5-v2.05.src/LCA.h
new file mode 100644
index 0000000..3426f7e
--- /dev/null
+++ b/blossom5-v2.05.src/LCA.h
@@ -0,0 +1,297 @@
+/*
+
+LCA.h - Imlementation of LCA data structure described in
+
+O. Berkman, U. Vishkin: Recursive Star-Tree Parallel Data Structure. SIAM J. Comput. 22(2): 221-242 (1993)
+
+It takes O(n \log n) space (and O(n \log n) time for construction). The memory requirement can be
+reduced to O(n) by defining LCA_BLOCKS. The implementation then follows the approach described in
+
+J. Fischer, V. Heun: Theoretical and Practical Improvements on the RMQ-Problem, with Applications to LCA and LCE.
+	Proceedings of the 17th Annual Symposium on Combinatorial Pattern Matching
+	(CPM'06), Lecture Notes in Computer Science 4009, 36-48, Springer-Verlag, 2006.
+
+(although in-block queries are computed naively in O(K) worst-case time, where K is the size of the block,
+rather than as described by Fischer and Heun).
+
+Written by Vladimir Kolmogorov, vnk@adastral.ucl.ac.uk.
+
+//  example for 
+//      4
+//     / \
+//    2   3
+//   / \
+//  0   1
+// (Note: the ordering of nodes is in the tree preorder!!!)
+
+	LCATree* lca = new LCATree(5);
+	char A0, A1, A2, A3, A4;
+
+	lca->Add(&A0, &A2);
+	lca->Add(&A1, &A2);
+	lca->Add(&A2, &A4);
+	lca->Add(&A3, &A4);
+	lca->AddRoot(&A4);
+
+	int result = lca->GetLCA(1, 3); // should be 4
+	delete lca;
+
+
+*/
+
+#ifndef GNAKDLATHJSTHAJSRNAKSJDA
+#define GNAKDLATHJSTHAJSRNAKSJDA
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+//#define LCA_BLOCKS
+
+class LCATree
+{
+public:
+	typedef void* NodeId; // can be any type, e.g. int. (The code checks NodeId's only for equalities.)
+	typedef int PreorderId;
+
+	LCATree(int node_num_max);
+	~LCATree();
+
+	// construct tree. Nodes must be added in the tree preorder!!!
+	// First call returns 0, second returns 1, and so on.
+	PreorderId Add(NodeId i, NodeId i_parent);
+	PreorderId AddRoot(NodeId i); // completes tree construction
+
+	PreorderId GetLCA(PreorderId i, PreorderId j);
+	// Let i0=i, j0=j be the input nodes, and let r = LCA(i0,j0).
+	// This function sets i and j to be the immediate children of r 
+	// such that i is a descendant of i0 and j is a descendant of j0.
+	// There must hold i0!=r, j0!=r.
+	void GetPenultimateNodes(PreorderId& i, PreorderId& j);
+
+
+//////////////////////////////////////////////////////////////////////////
+private:
+	int n, n_max, K, k_max;
+	int** array;
+
+	NodeId* buf0;
+	int* buf1;
+	NodeId* parent_current; 
+	int* child_current;
+
+	int* parents;
+
+	int _GetLCA(int i, int j); // same as GetLCA, but assumes that i<j
+	int GetLCADirect(int i, int j);
+};
+
+inline LCATree::LCATree(int _n_max) : n_max(_n_max), array(NULL)
+{
+#ifdef LCA_BLOCKS
+	K = -2;
+	n = n_max;
+	while (n > 0) { K ++; n /= 2; }
+	if (K < 1) K = 1;
+#else
+	K = 1;
+	n = 0;
+#endif
+	parents = new int[n_max];
+	buf0 = new NodeId[n_max];
+	buf1 = new int[n_max];
+	parent_current = buf0;
+	child_current = buf1;
+}
+
+inline LCATree::~LCATree()
+{
+	int k;
+	delete [] parents;
+	if (buf0) delete [] buf0;
+	if (buf1) delete [] buf1;
+	if (array)
+	{
+		for (k=1; k<=k_max; k++) delete [] array[k];
+		delete [] array;
+	}
+}
+
+inline LCATree::PreorderId LCATree::Add(NodeId i, NodeId i_parent)
+{
+	assert(n < n_max);
+
+	if (n == 0)
+	{
+		*parent_current = i;
+		*(++ parent_current) = i_parent;
+		parents[0] = -1;
+	}
+	else
+	{
+		if (i == *parent_current)
+		{
+			int c = *child_current --;
+			while ( 1 )
+			{
+				int c_next = parents[c];
+				parents[c] = n;
+				if (c_next < 0) break;
+				c = c_next;
+			}
+			parent_current --;
+		}
+		if (i_parent == *parent_current) parents[n] = *child_current;
+		else
+		{
+			*(++ parent_current) = i_parent;
+			parents[n] = -1;
+			child_current ++;
+		}
+	}
+	*child_current = n;
+	return n ++;
+}
+
+
+inline LCATree::PreorderId LCATree::AddRoot(NodeId i)
+{
+	assert(n < n_max);
+
+	if (n > 0)
+	{
+		if (i != *parent_current || parent_current != buf0+1)
+		{
+			printf("Error in LCATree construction: wrong sequence of calls!\n");
+			exit(1);
+		}
+		int c = *child_current --;
+		while ( 1 )
+		{
+			int c_next = parents[c];
+			parents[c] = n;
+			if (c_next < 0) break;
+			c = c_next;
+		}
+		child_current ++;
+	}
+	parents[n++] = -1;
+
+	delete [] buf0;
+	buf0 = NULL;
+	delete [] buf1;
+	buf1 = NULL;
+
+	// initialize array
+	int b, k = 1, block_num = (n-1)/K+1;
+	if (block_num < 3) return n-1;
+	int d = (block_num-1)/4;
+	while (d) { k ++; d >>= 1; }
+	k_max = k;
+
+	array = new int*[k_max+1];
+	array[0] = parents;
+	for (k=1, d=2; k<=k_max; k++, d*=2)
+	{
+		array[k] = new int[block_num-d];
+		if (k == 1)
+		{
+			for (b=0; b<block_num-d; b++) array[1][b] = GetLCADirect((b+1)*K-1, (b+d)*K);
+		}
+		else
+		{
+			for (b=0; b<block_num-d; b++)
+			{
+				int i = array[k-1][b];
+				int j = array[k-1][b+d/2];
+				if (i < j) i = j;
+				j = array[1][b+d/2-1];
+				array[k][b] = (i > j) ? i : j;
+			}
+		}
+	}
+
+	return n-1;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+
+inline int LCATree::GetLCADirect(int i, int j)
+{
+	while (i < j) i = parents[i];
+	return i;
+}
+
+
+inline int LCATree::_GetLCA(int i, int j)
+{
+#ifdef LCA_BLOCKS
+
+	int bi = i/K, bj = j/K;
+	if (bi == bj) return GetLCADirect(i, j);
+	int i_last = (bi+1)*K-1, j_first = bj*K;
+	i = GetLCADirect(i, i_last);
+	j = GetLCADirect(j_first, j);
+	if (i < j) i = j;
+	// set j = LCA(i_last, j_first)
+	if (j_first - i_last == 1) j = parents[i_last];
+	else
+	{
+		int k = 1, d = (bj-bi)/4;
+		while (d) { k ++; d >>= 1; }
+		int diff = 1<<k;
+		//assert(bi+diff <= bj && bi+diff>bj-diff);
+		j = (array[k][bi] > array[k][bj-diff]) ? array[k][bi] : array[k][bj-diff];
+	}
+	return (i > j) ? i : j;
+
+#else
+
+	if (j == i) return i;
+
+	int k = 0, d = (j-i)/2;
+	while (d) { k ++; d >>= 1; }
+	int diff = 1<<k;
+	//assert(i+diff <= j && i+diff>j-diff);
+	return (array[k][i] > array[k][j-diff]) ? array[k][i] : array[k][j-diff];
+
+#endif
+}
+
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+
+inline LCATree::PreorderId LCATree::GetLCA(PreorderId i, PreorderId j)
+{
+	if (i > j) { PreorderId k = i; i = j; j = k; }
+	return _GetLCA(i, j);
+}
+
+inline void LCATree::GetPenultimateNodes(PreorderId& _i, PreorderId& _j)
+{
+	int i, j, d, swap;
+	if (_i < _j) { i = _i; j = _j; swap = 0; }
+	else         { i = _j; j = _i; swap = 1; }
+	int r = _GetLCA(i, j);
+	assert(i!=r && j!=r);
+	while (parents[i] != r)
+	{
+		int i0 = parents[i];
+		d = (j - i0)/2;
+		while ( (i=_GetLCA(i0, i0+d)) == r ) d /= 2;
+	}
+	while (parents[j] != r)
+	{
+		int j0 = parents[j];
+		d = (r - j0)/2;
+		while ( (j=_GetLCA(j0, j0+d)) == r ) d /= 2;
+	}
+	if (swap == 0) { _i = i; _j = j; }
+	else           { _j = i; _i = j; }
+}
+
+#endif
diff --git a/blossom5-v2.05.src/LICENSE.TXT b/blossom5-v2.05.src/LICENSE.TXT
new file mode 100644
index 0000000..346e2db
--- /dev/null
+++ b/blossom5-v2.05.src/LICENSE.TXT
@@ -0,0 +1,25 @@
+    Copyright 2008-2009 UCL Business PLC, Author Vladimir Kolmogorov (vnk@ist.ac.at)
+
+    This software can be used for evaluation and non-commercial research purposes only. Commercial use is prohibited.
+    Public redistribution of the code or its derivatives is prohibited.
+    If you use this software for research purposes, you should cite the following paper in any resulting publication:
+
+        Vladimir Kolmogorov. "Blossom V: A new implementation of a minimum cost perfect matching algorithm."
+        In Mathematical Programming Computation (MPC), July 2009, 1(1):43-67.
+
+    For commercial use of the software not covered by this agreement, you may obtain a licence from
+    the copyright holders UCL Business via their licensing site: www.e-lucid.com/i/software/optimisation_software/BlossomV.html.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
diff --git a/blossom5-v2.05.src/Makefile b/blossom5-v2.05.src/Makefile
new file mode 100644
index 0000000..da30258
--- /dev/null
+++ b/blossom5-v2.05.src/Makefile
@@ -0,0 +1,26 @@
+DIRS := . MinCost GEOM 
+
+SOURCES := $(foreach dir, $(DIRS), $(wildcard $(dir)/*.cpp))
+OBJS := $(patsubst %.cpp, %.o, $(SOURCES))
+
+CFLAGS := -O3 -D_NDEBUG
+CXX ?= c++
+LIBS := 
+INCLUDES := 
+LIBDIR := 
+
+# Add librt if the target platform is not Darwin (OS X)
+ifneq ($(shell uname -s),Darwin)
+    LIBS += -lrt
+endif
+ 
+all: blossom5
+
+blossom5: ${OBJS}
+	$(CXX) $(CFLAGS) ${LIBDIR} -o $@ ${OBJS} ${LIBS}
+
+.cpp.o:
+	$(CXX) $(CFLAGS) ${INCLUDES} $< -c -o $@
+
+clean:
+	rm -f ${OBJS} blossom5
diff --git a/blossom5-v2.05.src/MinCost/MinCost.cpp b/blossom5-v2.05.src/MinCost/MinCost.cpp
new file mode 100644
index 0000000..a0a9759
--- /dev/null
+++ b/blossom5-v2.05.src/MinCost/MinCost.cpp
@@ -0,0 +1,287 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "MinCost.h"
+
+
+template <typename FlowType, typename CostType> 
+	MinCost<FlowType, CostType>::MinCost(int _nodeNum, int _edgeNumMax, void (*err_function)(const char *))
+	: nodeNum(_nodeNum),
+	  edgeNum(0),
+	  edgeNumMax(_edgeNumMax),
+	  counter(0),
+	  cost(0),
+	  error_function(err_function)
+{
+	nodes = (Node*) malloc(nodeNum*sizeof(Node));
+	arcs = (Arc*) malloc(2*edgeNumMax*sizeof(Arc));
+	if (!nodes || !arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
+
+	memset(nodes, 0, nodeNum*sizeof(Node));
+	memset(arcs, 0, 2*edgeNumMax*sizeof(Arc));
+	firstActive = &nodes[nodeNum];
+#ifdef MINCOST_DEBUG
+	for (int i=0; i<nodeNum; i++) nodes[i].id = i;
+#endif
+}
+
+template <typename FlowType, typename CostType> 
+	MinCost<FlowType, CostType>::~MinCost()
+{
+	free(nodes);
+	free(arcs);
+}
+
+template <typename FlowType, typename CostType> 
+	void MinCost<FlowType, CostType>::Init()
+{
+	Node* i;
+	Arc* a;
+
+	for (a=arcs; a<arcs+2*edgeNum; a++)
+	{
+		if (a->r_cap > 0 && a->GetRCost() < 0) PushFlow(a, a->r_cap);
+	}
+
+	Node** lastActivePtr = &firstActive;
+	for (i=nodes; i<nodes+nodeNum; i++)
+	{
+		if (i->excess > 0)
+		{
+			*lastActivePtr = i;
+			lastActivePtr = &i->next;
+		}
+		else i->next = NULL;
+	}
+	*lastActivePtr = &nodes[nodeNum];
+}
+
+
+template <typename FlowType, typename CostType> 
+	FlowType MinCost<FlowType, CostType>::Augment(Node* start, Node* end)
+{
+	FlowType delta = (start->excess < -end->excess) ? start->excess : -end->excess;
+	Arc* a;
+
+	for (a=end->parent; a; a=a->sister->head->parent)
+	{
+		if (delta > a->r_cap) delta = a->r_cap;
+	}
+	assert(delta > 0);
+
+	end->excess += delta;
+	for (a=end->parent; a; a=a->head->parent)
+	{
+		DecreaseRCap(a, delta);
+		a = a->sister;
+		IncreaseRCap(a, delta);
+	}
+	start->excess -= delta;
+
+	return delta;
+}
+
+template <typename FlowType, typename CostType> 
+	void MinCost<FlowType, CostType>::Dijkstra(Node* start)
+{
+	assert(start->excess > 0);
+
+	Node* i;
+	Node* j;
+	Arc* a;
+	CostType d;
+	Node* permanentNodes;
+
+	int FLAG0 = ++ counter; // permanently labeled nodes
+	int FLAG1 = ++ counter; // temporarily labeled nodes
+
+	start->parent = NULL;
+	start->flag = FLAG1;
+	queue.Reset();
+	queue.Add(start, 0);
+
+	permanentNodes = NULL;
+
+	while ( (i=queue.RemoveMin(d)) )
+	{
+		if (i->excess < 0)
+		{
+			FlowType delta = Augment(start, i);
+			cost += delta*(d - i->pi + start->pi);
+			for (i=permanentNodes; i; i=i->next_permanent) i->pi += d;
+			break;
+		}
+
+		i->pi -= d;
+		i->flag = FLAG0;
+		i->next_permanent = permanentNodes;
+		permanentNodes = i;
+
+		for (a=i->firstNonsaturated; a; a=a->next)
+		{
+			j = a->head;
+			if (j->flag == FLAG0) continue;
+			d = a->GetRCost();
+			if (j->flag == FLAG1)
+			{
+				if (d >= queue.GetKey(j)) continue;
+				queue.DecreaseKey(j, d);
+			}
+			else
+			{
+				queue.Add(j, d);
+				j->flag = FLAG1;
+			}
+			j->parent = a;
+		}
+
+	}
+}
+
+
+template <typename FlowType, typename CostType> 
+	CostType MinCost<FlowType, CostType>::Solve()
+{
+	Node* i;
+	//Init();
+	while ( 1 )
+	{
+		i = firstActive;
+		if (i == &nodes[nodeNum]) break;
+		firstActive = i->next;
+		i->next = NULL;
+		if (i->excess > 0)
+		{
+			Dijkstra(i);
+			if (i->excess > 0 && !i->next) 
+			{ 
+				i->next = firstActive; 
+				firstActive = i; 
+			}
+		}
+	}
+#ifdef MINCOST_DEBUG
+	TestOptimality();
+	TestCosts();
+#endif
+
+	return cost;
+}
+
+
+template <typename FlowType, typename CostType> 
+	void MinCost<FlowType, CostType>::TestOptimality()
+{
+	Node* i;
+	Arc* a;
+
+	for (i=nodes; i<nodes+nodeNum; i++)
+	{
+		if (i->excess != 0)
+		{
+			assert(0);
+		}
+		for (a=i->firstSaturated; a; a=a->next)
+		{
+			if (a->r_cap != 0)
+			{
+				assert(0);
+			}
+		}
+		for (a=i->firstNonsaturated; a; a=a->next)
+		{
+			CostType c = a->GetRCost();
+			if (a->r_cap <= 0 || a->GetRCost() < -1e-5)
+			{
+				assert(0);
+			}
+		}
+	}
+}
+
+#ifdef MINCOST_DEBUG
+
+template <typename FlowType, typename CostType> 
+	void MinCost<FlowType, CostType>::TestCosts()
+{
+	Arc* a;
+
+	CostType _cost = 0;
+
+	for (a=arcs; a<arcs+2*edgeNum; a+=2)
+	{
+		assert(a->r_cap + a->sister->r_cap == a->cap_orig + a->sister->cap_orig);
+		_cost += a->cost*(a->cap_orig - a->r_cap);
+	}
+
+	CostType delta = cost - _cost;
+	if (delta < 0) delta = -delta;
+	if (delta >= 1e-5)
+	{
+		assert(0);
+	}
+}
+
+#endif
+
+
+
+
+///////////////////////////////////////////////////////////////////////////////////////
+
+#define FLOW_INFTY ((int)0x00fffffff)
+
+template <typename CostType> 
+	DualMinCost<CostType>::DualMinCost(int _nodeNum, int _edgeNumMax)
+	: MinCost<int,CostType>(_nodeNum+1, _edgeNumMax+2*_nodeNum)
+{
+	source = _nodeNum;
+}
+
+template <typename CostType> 
+	DualMinCost<CostType>::~DualMinCost()
+{
+}
+
+template <typename CostType> 
+	void DualMinCost<CostType>::AddUnaryTerm(NodeId i, int objective_coef)
+{
+	MinCost<int, CostType>::AddNodeExcess(i, objective_coef);
+	MinCost<int, CostType>::AddNodeExcess(source, -objective_coef);
+}
+
+template <typename CostType> 
+	void DualMinCost<CostType>::SetLowerBound(NodeId i, CostType cmin)
+{
+	DualMinCost<CostType>::AddEdge(i, source, FLOW_INFTY, 0, -cmin);
+}
+
+template <typename CostType> 
+	void DualMinCost<CostType>::SetUpperBound(NodeId i, CostType cmax)
+{
+	DualMinCost<CostType>::AddEdge(source, i, FLOW_INFTY, 0, cmax);
+}
+
+template <typename CostType> 
+	void DualMinCost<CostType>::AddConstraint(NodeId i, NodeId j, CostType cmax)
+{
+	DualMinCost<CostType>::AddEdge(i, j, FLOW_INFTY, 0, cmax);
+}
+
+template <typename CostType> 
+	void DualMinCost<CostType>::Solve()
+{
+	MinCost<int,CostType>::Solve();
+}
+
+template <typename CostType> 
+	CostType DualMinCost<CostType>::GetSolution(NodeId i)
+{
+	return MinCost<int, CostType>::nodes[source].pi - MinCost<int, CostType>::nodes[i].pi;
+}
+
+
+
+#include "instances.inc"
+
+
diff --git a/blossom5-v2.05.src/MinCost/MinCost.cpp.orig b/blossom5-v2.05.src/MinCost/MinCost.cpp.orig
new file mode 100644
index 0000000..43d94bd
--- /dev/null
+++ b/blossom5-v2.05.src/MinCost/MinCost.cpp.orig
@@ -0,0 +1,287 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "MinCost.h"
+
+
+template <typename FlowType, typename CostType> 
+	MinCost<FlowType, CostType>::MinCost(int _nodeNum, int _edgeNumMax, void (*err_function)(const char *))
+	: nodeNum(_nodeNum),
+	  edgeNum(0),
+	  edgeNumMax(_edgeNumMax),
+	  counter(0),
+	  cost(0),
+	  error_function(err_function)
+{
+	nodes = (Node*) malloc(nodeNum*sizeof(Node));
+	arcs = (Arc*) malloc(2*edgeNumMax*sizeof(Arc));
+	if (!nodes || !arcs) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
+
+	memset(nodes, 0, nodeNum*sizeof(Node));
+	memset(arcs, 0, 2*edgeNumMax*sizeof(Arc));
+	firstActive = &nodes[nodeNum];
+#ifdef MINCOST_DEBUG
+	for (int i=0; i<nodeNum; i++) nodes[i].id = i;
+#endif
+}
+
+template <typename FlowType, typename CostType> 
+	MinCost<FlowType, CostType>::~MinCost()
+{
+	free(nodes);
+	free(arcs);
+}
+
+template <typename FlowType, typename CostType> 
+	void MinCost<FlowType, CostType>::Init()
+{
+	Node* i;
+	Arc* a;
+
+	for (a=arcs; a<arcs+2*edgeNum; a++)
+	{
+		if (a->r_cap > 0 && a->GetRCost() < 0) PushFlow(a, a->r_cap);
+	}
+
+	Node** lastActivePtr = &firstActive;
+	for (i=nodes; i<nodes+nodeNum; i++)
+	{
+		if (i->excess > 0)
+		{
+			*lastActivePtr = i;
+			lastActivePtr = &i->next;
+		}
+		else i->next = NULL;
+	}
+	*lastActivePtr = &nodes[nodeNum];
+}
+
+
+template <typename FlowType, typename CostType> 
+	FlowType MinCost<FlowType, CostType>::Augment(Node* start, Node* end)
+{
+	FlowType delta = (start->excess < -end->excess) ? start->excess : -end->excess;
+	Arc* a;
+
+	for (a=end->parent; a; a=a->sister->head->parent)
+	{
+		if (delta > a->r_cap) delta = a->r_cap;
+	}
+	assert(delta > 0);
+
+	end->excess += delta;
+	for (a=end->parent; a; a=a->head->parent)
+	{
+		DecreaseRCap(a, delta);
+		a = a->sister;
+		IncreaseRCap(a, delta);
+	}
+	start->excess -= delta;
+
+	return delta;
+}
+
+template <typename FlowType, typename CostType> 
+	void MinCost<FlowType, CostType>::Dijkstra(Node* start)
+{
+	assert(start->excess > 0);
+
+	Node* i;
+	Node* j;
+	Arc* a;
+	CostType d;
+	Node* permanentNodes;
+
+	int FLAG0 = ++ counter; // permanently labeled nodes
+	int FLAG1 = ++ counter; // temporarily labeled nodes
+
+	start->parent = NULL;
+	start->flag = FLAG1;
+	queue.Reset();
+	queue.Add(start, 0);
+
+	permanentNodes = NULL;
+
+	while ( (i=queue.RemoveMin(d)) )
+	{
+		if (i->excess < 0)
+		{
+			FlowType delta = Augment(start, i);
+			cost += delta*(d - i->pi + start->pi);
+			for (i=permanentNodes; i; i=i->next_permanent) i->pi += d;
+			break;
+		}
+
+		i->pi -= d;
+		i->flag = FLAG0;
+		i->next_permanent = permanentNodes;
+		permanentNodes = i;
+
+		for (a=i->firstNonsaturated; a; a=a->next)
+		{
+			j = a->head;
+			if (j->flag == FLAG0) continue;
+			d = a->GetRCost();
+			if (j->flag == FLAG1)
+			{
+				if (d >= queue.GetKey(j)) continue;
+				queue.DecreaseKey(j, d);
+			}
+			else
+			{
+				queue.Add(j, d);
+				j->flag = FLAG1;
+			}
+			j->parent = a;
+		}
+
+	}
+}
+
+
+template <typename FlowType, typename CostType> 
+	CostType MinCost<FlowType, CostType>::Solve()
+{
+	Node* i;
+	//Init();
+	while ( 1 )
+	{
+		i = firstActive;
+		if (i == &nodes[nodeNum]) break;
+		firstActive = i->next;
+		i->next = NULL;
+		if (i->excess > 0)
+		{
+			Dijkstra(i);
+			if (i->excess > 0 && !i->next) 
+			{ 
+				i->next = firstActive; 
+				firstActive = i; 
+			}
+		}
+	}
+#ifdef MINCOST_DEBUG
+	TestOptimality();
+	TestCosts();
+#endif
+
+	return cost;
+}
+
+
+template <typename FlowType, typename CostType> 
+	void MinCost<FlowType, CostType>::TestOptimality()
+{
+	Node* i;
+	Arc* a;
+
+	for (i=nodes; i<nodes+nodeNum; i++)
+	{
+		if (i->excess != 0)
+		{
+			assert(0);
+		}
+		for (a=i->firstSaturated; a; a=a->next)
+		{
+			if (a->r_cap != 0)
+			{
+				assert(0);
+			}
+		}
+		for (a=i->firstNonsaturated; a; a=a->next)
+		{
+			CostType c = a->GetRCost();
+			if (a->r_cap <= 0 || a->GetRCost() < -1e-5)
+			{
+				assert(0);
+			}
+		}
+	}
+}
+
+#ifdef MINCOST_DEBUG
+
+template <typename FlowType, typename CostType> 
+	void MinCost<FlowType, CostType>::TestCosts()
+{
+	Arc* a;
+
+	CostType _cost = 0;
+
+	for (a=arcs; a<arcs+2*edgeNum; a+=2)
+	{
+		assert(a->r_cap + a->sister->r_cap == a->cap_orig + a->sister->cap_orig);
+		_cost += a->cost*(a->cap_orig - a->r_cap);
+	}
+
+	CostType delta = cost - _cost;
+	if (delta < 0) delta = -delta;
+	if (delta >= 1e-5)
+	{
+		assert(0);
+	}
+}
+
+#endif
+
+
+
+
+///////////////////////////////////////////////////////////////////////////////////////
+
+#define FLOW_INFTY ((int)0x00fffffff)
+
+template <typename CostType> 
+	DualMinCost<CostType>::DualMinCost(int _nodeNum, int _edgeNumMax)
+	: MinCost<int,CostType>(_nodeNum+1, _edgeNumMax+2*_nodeNum)
+{
+	source = _nodeNum;
+}
+
+template <typename CostType> 
+	DualMinCost<CostType>::~DualMinCost()
+{
+}
+
+template <typename CostType> 
+	void DualMinCost<CostType>::AddUnaryTerm(NodeId i, int objective_coef)
+{
+	MinCost<int, CostType>::AddNodeExcess(i, objective_coef);
+	MinCost<int, CostType>::AddNodeExcess(source, -objective_coef);
+}
+
+template <typename CostType> 
+	void DualMinCost<CostType>::SetLowerBound(NodeId i, CostType cmin)
+{
+	AddEdge(i, source, FLOW_INFTY, 0, -cmin);
+}
+
+template <typename CostType> 
+	void DualMinCost<CostType>::SetUpperBound(NodeId i, CostType cmax)
+{
+	AddEdge(source, i, FLOW_INFTY, 0, cmax);
+}
+
+template <typename CostType> 
+	void DualMinCost<CostType>::AddConstraint(NodeId i, NodeId j, CostType cmax)
+{
+	AddEdge(i, j, FLOW_INFTY, 0, cmax);
+}
+
+template <typename CostType> 
+	void DualMinCost<CostType>::Solve()
+{
+	MinCost<int,CostType>::Solve();
+}
+
+template <typename CostType> 
+	CostType DualMinCost<CostType>::GetSolution(NodeId i)
+{
+	return MinCost<int, CostType>::nodes[source].pi - MinCost<int, CostType>::nodes[i].pi;
+}
+
+
+
+#include "instances.inc"
+
+
diff --git a/blossom5-v2.05.src/MinCost/MinCost.cpp.rej b/blossom5-v2.05.src/MinCost/MinCost.cpp.rej
new file mode 100644
index 0000000..1597c49
--- /dev/null
+++ b/blossom5-v2.05.src/MinCost/MinCost.cpp.rej
@@ -0,0 +1,25 @@
+--- MinCost.cpp
++++ MinCost.cpp
+@@ -253,19 +253,19 @@
+ template <typename CostType> 
+ 	void DualMinCost<CostType>::SetLowerBound(NodeId i, CostType cmin)
+ {
+-	AddEdge(i, source, FLOW_INFTY, 0, -cmin);
++	DualMinCost<CostType>::AddEdge(i, source, FLOW_INFTY, 0, -cmin);
+ }
+ 
+ template <typename CostType> 
+ 	void DualMinCost<CostType>::SetUpperBound(NodeId i, CostType cmax)
+ {
+-	AddEdge(source, i, FLOW_INFTY, 0, cmax);
++	DualMinCost<CostType>::AddEdge(source, i, FLOW_INFTY, 0, cmax);
+ }
+ 
+ template <typename CostType> 
+ 	void DualMinCost<CostType>::AddConstraint(NodeId i, NodeId j, CostType cmax)
+ {
+-	AddEdge(i, j, FLOW_INFTY, 0, cmax);
++	DualMinCost<CostType>::AddEdge(i, j, FLOW_INFTY, 0, cmax);
+ }
+ 
+ template <typename CostType> 
diff --git a/blossom5-v2.05.src/MinCost/MinCost.h b/blossom5-v2.05.src/MinCost/MinCost.h
new file mode 100644
index 0000000..fa3bc91
--- /dev/null
+++ b/blossom5-v2.05.src/MinCost/MinCost.h
@@ -0,0 +1,504 @@
+/*
+    MinCost.h - successive shortest path algorithm of Ford and Fulkerson for solving a minimum cost flow problem
+
+    Copyright 2008 Vladimir Kolmogorov (vnk@adastral.ucl.ac.uk)
+
+    This software can be used for research purposes only. Commercial use is prohibited.
+    Public redistribution of the code or its derivatives is prohibited.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+
+#ifndef __MINCOST_H__
+#define __MINCOST_H__
+
+#include <string.h>
+#include <assert.h>
+
+// if GRAPH_ASSERT is defined then all calls to graph construction functions are assert'ed for correctness
+// (e.g. that node_id's are valid id's and edge capacities are non-negative).
+//#define GRAPH_ASSERT 
+
+//#define MINCOST_DEBUG
+
+
+
+template <typename FlowType, typename CostType> class MinCost
+{
+public:
+	typedef int NodeId;
+	typedef int EdgeId;
+
+	MinCost(int NodeNum, int edgeNumMax, void (*err_function)(const char *) = NULL);
+
+	// Destructor
+	~MinCost();
+
+	void AddNodeExcess(NodeId i, FlowType excess);
+
+	// first call returns 0, second 1, and so on.
+	// cap, rev_cap must be non-negative. 
+	// cost can be negative.
+	EdgeId AddEdge(NodeId i, NodeId j, FlowType cap, FlowType rev_cap, CostType cost);
+
+	CostType Solve();
+
+	///////////////////////////////////////////////////
+
+	FlowType GetRCap(EdgeId e);
+	void SetRCap(EdgeId e, FlowType new_rcap);
+	FlowType GetReverseRCap(EdgeId e);
+	void SetReverseRCap(EdgeId e, FlowType new_rcap);
+	void PushFlow(EdgeId e, FlowType delta);
+	void UpdateCost(EdgeId e, FlowType cap_orig, CostType delta);
+
+	CostType GetDual(NodeId i) { return nodes[i].pi; }
+
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+	
+protected:
+	// internal variables and functions
+
+	struct Node;
+	struct Arc;
+
+	struct Node
+	{
+		Arc			*firstNonsaturated;
+		Arc			*firstSaturated;
+
+		Arc			*parent;
+		Node		*next; // list of nodes with positive excesses
+
+		FlowType	excess;
+		CostType	pi;
+		int			flag;
+		union
+		{
+			int		heap_ptr;
+			Node*	next_permanent;
+		};
+#ifdef MINCOST_DEBUG
+		int			id;
+#endif
+	};
+
+	struct Arc
+	{
+		Node		*head;
+		Arc			*prev;
+		Arc			*next;
+		Arc			*sister;	// reverse arc
+
+		FlowType	r_cap;		// residual capacity
+#ifdef MINCOST_DEBUG
+		FlowType	cap_orig;
+#endif
+		CostType	cost;
+		CostType GetRCost() { return cost + head->pi - sister->head->pi; }
+	};
+
+	int		nodeNum, edgeNum, edgeNumMax;
+	Node	*nodes;
+	Arc		*arcs;
+	Node*	firstActive;
+	int		counter;
+	CostType cost;
+
+
+	void	(*error_function)(const char *);	// this function is called if a error occurs,
+										// with a corresponding error message
+										// (or exit(1) is called if it's NULL)
+
+	/////////////////////////////////////////////////////////////////////////
+
+	struct PriorityQueue
+	{
+		PriorityQueue();
+		~PriorityQueue();
+		void Reset();
+		CostType GetKey(Node* i);
+		void Add(Node* i, CostType key);
+		void DecreaseKey(Node* i, CostType key);
+		Node* RemoveMin(CostType& key);
+
+	private:
+		struct Item
+		{
+			Node*		i;
+			CostType	key;
+		}* array;
+		int N, arraySize;
+		void Swap(int k1, int k2);
+	};
+
+	PriorityQueue queue;
+
+	/////////////////////////////////////////////////////////////////////////
+
+	void SetRCap(Arc* a, FlowType new_rcap);
+	void PushFlow(Arc* a, FlowType delta);
+
+	void Init();
+	void DecreaseRCap(Arc* a, FlowType delta);
+	void IncreaseRCap(Arc* a, FlowType delta);
+	FlowType Augment(Node* start, Node* end);
+	void Dijkstra(Node* start);
+
+	void TestOptimality();
+#ifdef MINCOST_DEBUG
+	void TestCosts();
+#endif
+};
+
+
+template <typename CostType> class DualMinCost : private MinCost<int, CostType>
+{
+public:
+	typedef int NodeId;
+	DualMinCost(int node_num, int constraint_num_max);
+	~DualMinCost();
+
+	void AddUnaryTerm(NodeId i, int objective_coef);
+	void SetLowerBound(NodeId, CostType cmin);
+	void SetUpperBound(NodeId, CostType cmax);
+	void AddConstraint(NodeId i, NodeId j, CostType cmax); // xj - xi <= cmax
+
+	void Solve();
+	CostType GetSolution(NodeId i);
+
+private:
+	NodeId	source;
+};
+
+
+
+
+
+
+
+
+
+///////////////////////////////////////
+// Implementation - inline functions //
+///////////////////////////////////////
+
+
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::AddNodeExcess(NodeId _i, FlowType excess)
+{
+	assert(_i>=0 && _i<nodeNum);
+	nodes[_i].excess += excess;
+	if (nodes[_i].excess > 0 && !nodes[_i].next)
+	{
+		nodes[_i].next = firstActive;
+		firstActive = &nodes[_i];
+	}
+}
+
+template <typename FlowType, typename CostType> 
+	inline typename MinCost<FlowType, CostType>::EdgeId MinCost<FlowType, CostType>::AddEdge(NodeId _i, NodeId _j, FlowType cap, FlowType rev_cap, CostType cost)
+{
+	assert(_i>=0 && _i<nodeNum);
+	assert(_j>=0 && _j<nodeNum);
+	assert(_i!=_j && edgeNum<edgeNumMax);
+	assert(cap >= 0);
+	assert(rev_cap >= 0);
+
+	Arc *a = &arcs[2*edgeNum];
+	Arc *a_rev = a+1;
+	edgeNum ++;
+
+	Node* i = nodes + _i;
+	Node* j = nodes + _j;
+
+	a -> sister = a_rev;
+	a_rev -> sister = a;
+	if (cap > 0)
+	{
+		if (i->firstNonsaturated) i->firstNonsaturated->prev = a;
+		a -> next = i -> firstNonsaturated;
+		i -> firstNonsaturated = a;
+	}
+	else
+	{
+		if (i->firstSaturated) i->firstSaturated->prev = a;
+		a -> next = i -> firstSaturated;
+		i -> firstSaturated = a;
+	}
+	a->prev = NULL;
+	if (rev_cap > 0)
+	{
+		if (j->firstNonsaturated) j->firstNonsaturated->prev = a_rev;
+		a_rev -> next = j -> firstNonsaturated;
+		j -> firstNonsaturated = a_rev;
+	}
+	else
+	{
+		if (j->firstSaturated) j->firstSaturated->prev = a_rev;
+		a_rev -> next = j -> firstSaturated;
+		j -> firstSaturated = a_rev;
+	}
+	a_rev->prev = NULL;
+
+	a -> head = j;
+	a_rev -> head = i;
+	a -> r_cap = cap;
+	a_rev -> r_cap = rev_cap;
+	a -> cost = cost;
+	a_rev -> cost = -cost;
+#ifdef MINCOST_DEBUG
+	a->cap_orig = cap;
+	a_rev->cap_orig = rev_cap;
+#endif
+
+	if (a->r_cap > 0 && a->GetRCost() < 0) PushFlow(a, a->r_cap);
+	if (a_rev->r_cap > 0 && a_rev->GetRCost() < 0) PushFlow(a_rev, a_rev->r_cap);
+
+	return edgeNum-1;
+}
+
+///////////////////////////////////////
+///////////////////////////////////////
+///////////////////////////////////////
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::DecreaseRCap(Arc* a, FlowType delta)
+{
+	a->r_cap -= delta;
+	if (a->r_cap == 0)
+	{
+		Node* i = a->sister->head;
+		if (a->next) a->next->prev = a->prev;
+		if (a->prev) a->prev->next = a->next;
+		else         i->firstNonsaturated = a->next;
+		a->next = i->firstSaturated;
+		if (a->next) a->next->prev = a;
+		a->prev = NULL;
+		i->firstSaturated = a;
+	}
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::IncreaseRCap(Arc* a, FlowType delta)
+{
+	if (a->r_cap == 0)
+	{
+		Node* i = a->sister->head;
+		if (a->next) a->next->prev = a->prev;
+		if (a->prev) a->prev->next = a->next;
+		else         i->firstSaturated = a->next;
+		a->next = i->firstNonsaturated;
+		if (a->next) a->next->prev = a;
+		a->prev = NULL;
+		i->firstNonsaturated = a;
+	}
+	a->r_cap += delta;
+}
+
+template <typename FlowType, typename CostType> 
+	inline FlowType MinCost<FlowType, CostType>::GetRCap(EdgeId e)
+{
+	Arc* a = &arcs[2*e];
+	return a->r_cap;
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::SetRCap(Arc* a, FlowType new_rcap)
+{
+	assert(new_rcap >= 0);
+#ifdef MINCOST_DEBUG
+	a->cap_orig += new_rcap - a->r_cap;
+#endif
+	if (a->r_cap == 0)
+	{
+		Node* i = a->sister->head;
+		if (a->next) a->next->prev = a->prev;
+		if (a->prev) a->prev->next = a->next;
+		else         i->firstSaturated = a->next;
+		a->next = i->firstNonsaturated;
+		if (a->next) a->next->prev = a;
+		a->prev = NULL;
+		i->firstNonsaturated = a;
+	}
+	a->r_cap = new_rcap;
+	if (a->r_cap == 0)
+	{
+		Node* i = a->sister->head;
+		if (a->next) a->next->prev = a->prev;
+		if (a->prev) a->prev->next = a->next;
+		else         i->firstNonsaturated = a->next;
+		a->next = i->firstSaturated;
+		if (a->next) a->next->prev = a;
+		a->prev = NULL;
+		i->firstSaturated = a;
+	}
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::SetRCap(EdgeId e, FlowType new_rcap)
+{
+	SetRCap(&arcs[2*e], new_rcap);
+}
+
+template <typename FlowType, typename CostType> 
+	inline FlowType MinCost<FlowType, CostType>::GetReverseRCap(EdgeId e)
+{
+	Arc* a = &arcs[2*e+1];
+	return a->r_cap;
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::SetReverseRCap(EdgeId e, FlowType new_rcap)
+{
+	SetRCap(&arcs[2*e+1], new_rcap);
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::PushFlow(Arc* a, FlowType delta)
+{
+	if (delta < 0) { a = a->sister; delta = -delta; }
+	DecreaseRCap(a, delta);
+	IncreaseRCap(a->sister, delta);
+	a->head->excess += delta;
+	a->sister->head->excess -= delta;
+	cost += delta*a->cost;
+	if (a->head->excess > 0 && !a->head->next)
+	{
+		a->head->next = firstActive;
+		firstActive = a->head;
+	}
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::PushFlow(EdgeId e, FlowType delta)
+{
+	PushFlow(&arcs[2*e], delta);
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::UpdateCost(EdgeId e, FlowType cap_orig, CostType delta)
+{
+	Arc* a = &arcs[2*e];
+	cost += delta*(cap_orig-a->r_cap);
+	a->cost += delta;
+	a->sister->cost = -a->cost;
+
+	if (a->GetRCost() > 0) a = a->sister;
+	if (a->r_cap > 0 && a->GetRCost() < 0) PushFlow(a, a->r_cap);
+}
+
+///////////////////////////////////////
+///////////////////////////////////////
+///////////////////////////////////////
+
+template <typename FlowType, typename CostType> 
+	inline MinCost<FlowType, CostType>::PriorityQueue::PriorityQueue()
+{
+	N = 0;
+	arraySize = 16;
+	array = (Item*) malloc(arraySize*sizeof(Item));
+}
+
+template <typename FlowType, typename CostType> 
+	inline MinCost<FlowType, CostType>::PriorityQueue::~PriorityQueue()
+{
+	free(array);
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::PriorityQueue::Reset()
+{
+	N = 0;
+}
+
+template <typename FlowType, typename CostType> 
+	inline CostType MinCost<FlowType, CostType>::PriorityQueue::GetKey(Node* i)
+{
+	return array[i->heap_ptr].key;
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::PriorityQueue::Swap(int k1, int k2)
+{
+	Item* a = array+k1;
+	Item* b = array+k2;
+	a->i->heap_ptr = k2;
+	b->i->heap_ptr = k1;
+	Node* i = a->i;   a->i   = b->i;   b->i   = i;
+	CostType key = a->key; a->key = b->key; b->key = key;
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::PriorityQueue::Add(Node* i, CostType key)
+{
+	if (N == arraySize)
+	{
+		arraySize *= 2;
+		array = (Item*) realloc(array, arraySize*sizeof(Item));
+	}
+	int k = i->heap_ptr = N ++;
+	array[k].i = i;
+	array[k].key = key;
+	while (k > 0)
+	{
+		int k2 = (k-1)/2;
+		if (array[k2].key <= array[k].key) break;
+		Swap(k, k2);
+		k = k2;
+	}
+}
+
+template <typename FlowType, typename CostType> 
+	inline void MinCost<FlowType, CostType>::PriorityQueue::DecreaseKey(Node* i, CostType key)
+{
+	int k = i->heap_ptr;
+	array[k].key = key;
+	while (k > 0)
+	{
+		int k2 = (k-1)/2;
+		if (array[k2].key <= array[k].key) break;
+		Swap(k, k2);
+		k = k2;
+	}
+}
+
+template <typename FlowType, typename CostType> 
+	inline typename MinCost<FlowType, CostType>::Node* MinCost<FlowType, CostType>::PriorityQueue::RemoveMin(CostType& key)
+{
+	if (N == 0) return NULL;
+
+	Swap(0, N-1);
+	N --;
+
+	int k = 0;
+	while ( 1 )
+	{
+		int k1 = 2*k + 1, k2 = k1 + 1;
+		if (k1 >= N) break;
+		int k_min = (k2 >= N || array[k1].key <= array[k2].key) ? k1 : k2;
+		if (array[k].key <= array[k_min].key) break;
+		Swap(k, k_min);
+		k = k_min;
+	}
+
+	key = array[N].key;
+	return array[N].i;
+}
+
+
+#endif
diff --git a/blossom5-v2.05.src/MinCost/instances.inc b/blossom5-v2.05.src/MinCost/instances.inc
new file mode 100644
index 0000000..a6d099e
--- /dev/null
+++ b/blossom5-v2.05.src/MinCost/instances.inc
@@ -0,0 +1,12 @@
+#include "MinCost.h"
+
+#ifdef _MSC_VER
+#pragma warning(disable: 4661)
+#endif
+
+
+template class MinCost<int,int>;
+template class MinCost<int,double>;
+
+template class DualMinCost<int>;
+template class DualMinCost<double>;
diff --git a/blossom5-v2.05.src/MinCost/p b/blossom5-v2.05.src/MinCost/p
new file mode 100644
index 0000000..84855b0
--- /dev/null
+++ b/blossom5-v2.05.src/MinCost/p
@@ -0,0 +1,26 @@
+index 3b6bedd..44338f7 100644
+--- Autotune/blossomv/MinCost/MinCost.cpp
++++ Autotune/blossomv/MinCost/MinCost.cpp
+@@ -253,19 +253,19 @@ template <typename CostType>
+ template <typename CostType> 
+ 	void DualMinCost<CostType>::SetLowerBound(NodeId i, CostType cmin)
+ {
+-	AddEdge(i, source, FLOW_INFTY, 0, -cmin);
++	DualMinCost<CostType>::AddEdge(i, source, FLOW_INFTY, 0, -cmin);
+ }
+ 
+ template <typename CostType> 
+ 	void DualMinCost<CostType>::SetUpperBound(NodeId i, CostType cmax)
+ {
+-	AddEdge(source, i, FLOW_INFTY, 0, cmax);
++	DualMinCost<CostType>::AddEdge(source, i, FLOW_INFTY, 0, cmax);
+ }
+ 
+ template <typename CostType> 
+ 	void DualMinCost<CostType>::AddConstraint(NodeId i, NodeId j, CostType cmax)
+ {
+-	AddEdge(i, j, FLOW_INFTY, 0, cmax);
++	DualMinCost<CostType>::AddEdge(i, j, FLOW_INFTY, 0, cmax);
+ }
+ 
+ template <typename CostType> 
diff --git a/blossom5-v2.05.src/PMduals.cpp b/blossom5-v2.05.src/PMduals.cpp
new file mode 100644
index 0000000..ac3fa66
--- /dev/null
+++ b/blossom5-v2.05.src/PMduals.cpp
@@ -0,0 +1,441 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "PMimplementation.h"
+#include "MinCost/MinCost.h"
+
+
+void PerfectMatching::ComputeEpsGlobal()
+{
+	Node* r;
+	PriorityQueue<REAL>::Item* q;
+	Tree* t;
+	Tree* t2;
+	TreeEdge* e;
+	int i, j, k, N = 0, E = 0;
+
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		t = r->tree;
+		t->id = N;
+		N += 2;
+		for (k=0; k<2; k++)
+		for (e=t->first[k]; e; e=e->next[k]) E += 6;
+	}
+	DualMinCost<REAL>* m = new DualMinCost<REAL>(N, E);
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		t = r->tree;
+		i = t->id;
+		m->AddUnaryTerm(i, -1);
+		m->SetLowerBound(i, 0);
+		m->AddUnaryTerm(i+1, 1);
+		m->SetUpperBound(i+1, 0);
+
+		if (t->eps_delta < PM_INFTY)
+		{
+			m->SetUpperBound(i, t->eps_delta);
+			m->SetLowerBound(i+1, -t->eps_delta);
+		}
+		for (e=t->first[0]; e; e=e->next[0])
+		{
+			t2 = e->head[0];
+			if (t2 == NULL) continue;
+			j = e->head[0]->id;
+			if ((q=e->pq01[0].GetMin()))
+			{
+				m->AddConstraint(j, i, q->slack - t->eps + t2->eps);
+				m->AddConstraint(i+1, j+1, q->slack - t->eps + t2->eps);
+			}
+			if ((q=e->pq01[1].GetMin()))
+			{
+				m->AddConstraint(i, j, q->slack - t2->eps + t->eps);
+				m->AddConstraint(j+1, i+1, q->slack - t2->eps + t->eps);
+			}
+			if ((q=e->pq00.GetMin()))
+			{
+				m->AddConstraint(i+1, j, q->slack - t->eps - t2->eps);
+				m->AddConstraint(j+1, i, q->slack - t->eps - t2->eps);
+			}
+		}
+	}
+	m->Solve();
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		t = r->tree;
+		i = t->id;
+		t->eps_delta = (m->GetSolution(i) - m->GetSolution(i+1))/2;
+	}
+	delete m;
+}
+
+void PerfectMatching::ComputeEpsSingle()
+{
+	Node* r;
+	PriorityQueue<REAL>::Item* q;
+	Tree* t;
+	Tree* t2;
+	TreeEdge* e;
+	REAL eps = PM_INFTY;
+
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		t = r->tree;
+		if (eps > t->eps_delta) eps = t->eps_delta;
+		for (e=t->first[0]; e; e=e->next[0])
+		{
+			t2 = e->head[0];
+			if ((q=e->pq00.GetMin()) && 2*eps > q->slack-t->eps-t2->eps)
+			{
+				eps = (q->slack-t->eps-t2->eps)/2;
+			}
+		}
+	}
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		r->tree->eps_delta = eps;
+	}
+}
+
+
+void PerfectMatching::ComputeEpsCC()
+{
+	Node* r;
+	PriorityQueue<REAL>::Item* q;
+	Tree* t0;
+	Tree* t;
+	Tree* t2;
+	Tree* t_next;
+	TreeEdge* e;
+	REAL eps, eps2;
+	Tree* queue_last;
+	int dir;
+	Tree* FIXED_TREE = trees-1;
+	int component_num = 0;
+	TreeEdge** e_ptr;
+
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		t0 = r->tree;
+		t0->next = NULL;
+	}
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		t0 = r->tree;
+		if (t0->next) continue;
+		eps = t0->eps_delta;
+
+		t0->next = queue_last = t = t0;
+		while ( 1 )
+		{
+			for (dir=0; dir<2; dir++)
+			for (e_ptr=&t->first[dir], e=*e_ptr; e; e=*e_ptr)
+			{
+				t2 = e->head[dir];
+				if (t2 == NULL) { *e_ptr = e->next[dir]; tree_edges->Delete(e); continue; }
+				e_ptr = &e->next[dir];
+
+				REAL eps00 = ((q=e->pq00.GetMin())) ? (q->slack - t->eps - t2->eps) : PM_INFTY;
+				if (t2->next && t2->next != FIXED_TREE)
+				{
+					if (2*eps > eps00) eps = eps00/2;
+					continue;
+				}
+
+				REAL eps01[2];
+				eps01[dir]   = ((q=e->pq01[dir].GetMin()))   ? (q->slack - t->eps + t2->eps) : PM_INFTY;
+				eps01[1-dir] = ((q=e->pq01[1-dir].GetMin())) ? (q->slack - t2->eps + t->eps) : PM_INFTY;
+
+				if (t2->next == FIXED_TREE) eps2 = t2->eps_delta;
+				else if (eps01[0] > 0 && eps01[1] > 0) eps2 = 0;
+				else
+				{
+					queue_last->next = t2;
+					queue_last = t2;
+					t2->next = t2;
+					if (eps > eps00) eps = eps00;
+					if (eps > t2->eps_delta) eps = t2->eps_delta;
+					continue;
+				}
+				if (eps > eps00 - eps2)      eps = eps00 - eps2;
+				if (eps > eps2 + eps01[dir]) eps = eps2 + eps01[dir];
+			}
+
+			if (t->next == t) break;
+			t = t->next;
+		}
+		for (t=t0; ; t=t_next)
+		{
+			t->eps_delta = eps;
+			t_next = t->next;
+			t->next = FIXED_TREE;
+			if (t_next == t) break;
+		}
+		component_num ++;
+	}
+	//printf("%d CCs ", component_num);
+}
+
+
+void PerfectMatching::ComputeEpsSCC()
+{
+	PriorityQueue<REAL>::Item* q;
+	Node* r;
+	Tree* t0;
+	Tree* t;
+	Tree* t2;
+	TreeEdge* e;
+	TreeEdge** e_ptr;
+	REAL eps;
+	int c, dir;
+
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		t0 = r->tree;
+		t0->dfs_parent = NULL;
+
+		for (dir=0; dir<2; dir++)
+		for (e_ptr=&t0->first[dir], e=*e_ptr; e; e=*e_ptr)
+		{
+			t2 = e->head[dir];
+			if (t2 == NULL) { *e_ptr = e->next[dir]; tree_edges->Delete(e); continue; }
+			e_ptr = &e->next[dir];
+		}
+	}
+	Tree* stack = NULL;
+
+	// first DFS
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		t0 = r->tree;
+		if (t0->dfs_parent) continue;
+		t = t0;
+		e = (t->first[0]) ? t->first[0] : t->first[1];
+		t->dfs_parent = (TreeEdge*)trees;
+		while ( 1 )
+		{
+			if (e == NULL)
+			{
+				t->next = stack;
+				stack = t;
+
+				if (t == t0) break;
+
+				e = t->dfs_parent;
+				if (t == e->head[0]) { t = e->head[1]; e = (e->next[0]) ? e->next[0] : t->first[1]; }
+				else                 { t = e->head[0]; e = e->next[1]; }
+				continue;
+			}
+
+			if (e->head[1] == t)
+			{
+				if (e->head[0]->dfs_parent || !(q=e->pq01[0].GetMin()) || q->slack - t->eps + e->head[0]->eps > 0) { e = (e->next[0]) ? e->next[0] : t->first[1]; continue; }
+				t = e->head[0];
+			}
+			else
+			{
+				if (e->head[1]->dfs_parent || !(q=e->pq01[1].GetMin()) || q->slack - t->eps + e->head[1]->eps > 0) { e = e->next[1]; continue; }
+				t = e->head[1];
+			}
+			t->dfs_parent = e;
+			e = (t->first[0]) ? t->first[0] : t->first[1];
+		}
+	}
+
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) r->tree->dfs_parent = NULL;
+
+	int component_num = 0;
+	while (stack)
+	{
+		t0 = stack;
+		stack = t0->next;
+		if (t0->dfs_parent) continue;
+		t = t0;
+		e = (t->first[0]) ? t->first[0] : t->first[1];
+		t->dfs_parent = (TreeEdge*)trees;
+		while ( 1 )
+		{
+			if (e == NULL)
+			{
+				e = t->dfs_parent;
+				t->dfs_parent = (TreeEdge*)((char*)trees + component_num);
+				if (t == t0) break;
+				if (t == e->head[0]) { t = e->head[1]; e = (e->next[0]) ? e->next[0] : t->first[1]; }
+				else                 { t = e->head[0]; e = e->next[1]; }
+				continue;
+			}
+
+			if (e->head[1] == t)
+			{
+				if (e->head[0]->dfs_parent || !(q=e->pq01[1].GetMin()) || q->slack - e->head[0]->eps + t->eps > 0) { e = (e->next[0]) ? e->next[0] : t->first[1]; continue; }
+				t = e->head[0];
+			}
+			else
+			{
+				if (e->head[1]->dfs_parent || !(q=e->pq01[0].GetMin()) || q->slack - e->head[1]->eps + t->eps > 0) { e = e->next[1]; continue; }
+				t = e->head[1];
+			}
+			t->dfs_parent = e;
+			e = (t->first[0]) ? t->first[0] : t->first[1];
+		}
+		component_num ++;
+	}
+
+	Tree** array = new Tree*[component_num];
+	memset(array, 0, component_num*sizeof(Tree*));
+
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		t = r->tree;
+		t->id = (int)((char*)t->dfs_parent - (char*)trees);
+		t->next = array[t->id];
+		array[t->id] = t;
+	}
+
+	for (c=component_num-1; c>=0; c--)
+	{
+		eps = PM_INFTY;
+		for (t=array[c]; t; t=t->next)
+		{
+			if (eps > t->eps_delta) eps = t->eps_delta;
+			FOR_ALL_TREE_EDGES(t, e, dir)
+			{
+				t2 = e->head[dir];
+				REAL eps00 = (q=e->pq00.GetMin()) ? (q->slack-t->eps-t2->eps) : PM_INFTY;
+				REAL eps01[2];
+				eps01[dir]   = ((q=e->pq01[dir].GetMin()))   ? (q->slack - t->eps + t2->eps) : PM_INFTY;
+				eps01[1-dir] = ((q=e->pq01[1-dir].GetMin())) ? (q->slack - t2->eps + t->eps) : PM_INFTY;
+				if (t2->id < c)
+				{
+					if (eps > eps01[dir]) eps = eps01[dir];
+					if (eps > eps00) eps = eps00;
+				}
+				else if (t2->id == c)
+				{
+					if (2*eps > eps00) eps = eps00 / 2;
+				}
+				else
+				{
+					if (eps > eps01[dir] + t2->eps_delta) eps = eps01[dir] + t2->eps_delta;
+					if (eps > eps00 - t2->eps_delta)  eps = eps00 - t2->eps_delta;
+				}
+			}
+		}
+		for (t=array[c]; t; t=t->next) t->eps_delta = eps;
+	}
+
+	delete [] array;
+	//printf("%d SCCs ", component_num);
+}
+
+void PerfectMatching::CommitEps()
+{
+	printf("CommitEps()\n");
+	Node* i;
+	Node* j;
+	Node* r;
+	int dir;
+	Edge* a;
+	EdgeIterator I;
+	Tree* t;
+	TreeEdge* e;
+	TreeEdge** e_ptr;
+	REAL eps, eps2;
+	PriorityQueue<REAL>::Item* q;
+
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		t = r->tree;
+		eps = t->eps;
+
+		i = r;
+		while ( 1 )
+		{
+			i->y += eps;
+			if (!i->is_tree_root)
+			{
+				Node* i0 = i;
+				i = ARC_HEAD(i0->match);
+				if (i->is_blossom) ARC_TO_EDGE_PTR(i0->match)->slack -= eps;
+				else               i->y                              -= eps;
+				FOR_ALL_EDGES(i, a, dir, I)
+				{
+					GET_OUTER_HEAD(a, dir, j);
+
+					a->slack += eps;
+					if (j->flag == 0) a->slack -= j->tree->eps;
+				}
+				i = i0;
+			}
+
+			MOVE_NODE_IN_TREE(i);
+		}
+
+		t->pq0.Update(-eps);
+
+		PriorityQueue<REAL> pq00 = t->pq00;
+		t->pq00.Reset();
+		for (q=pq00.GetAndResetFirst(); q; q=pq00.GetAndResetNext())
+		{
+			a = (Edge*)q;
+			if (ProcessEdge00(a)) t->pq00.Add(a);
+		}
+
+		for (e_ptr=&t->first[0], e=*e_ptr; e; e=*e_ptr)
+		{
+			if (e->head[0] == NULL) { *e_ptr = e->next[0]; tree_edges->Delete(e); continue; }
+			e_ptr = &e->next[0];
+
+			eps2 = e->head[0]->eps;
+			e->pq00.Update( - eps - eps2 );
+		}
+	}
+
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next) r->tree->eps = 0;
+}
+
+bool PerfectMatching::UpdateDuals()
+{
+	Node* r;
+
+	double start_time = get_time();
+
+	////////////////////////////////////////////////////////////////////////////////////
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		Tree* t = r->tree;
+		PriorityQueue<REAL>::Item* q;
+		REAL eps = PM_INFTY;
+		if ((q=t->pq0.GetMin())) eps = q->slack;
+		if ((q=t->pq_blossoms.GetMin()) && eps > q->slack) eps = q->slack;
+		while ((q=t->pq00.GetMin()))
+		{
+			if (ProcessEdge00((Edge*)q, false)) break;
+			t->pq00.Remove(q, pq_buf);
+		}
+		if (q && 2*eps > q->slack) eps = q->slack/2;
+		t->eps_delta = eps - t->eps;
+	}
+
+	if (tree_num >= options.dual_LP_threshold*node_num)
+	{
+		if      (options.dual_greedy_update_option == 0) ComputeEpsCC();
+		else if (options.dual_greedy_update_option == 1) ComputeEpsSCC();
+		else                                             ComputeEpsSingle();
+	}
+	else ComputeEpsGlobal();
+
+	REAL delta = 0;
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		if (r->tree->eps_delta > 0)
+		{
+			delta += r->tree->eps_delta;
+			r->tree->eps += r->tree->eps_delta;
+		}
+	}
+
+	stat.dual_time += get_time() - start_time;
+
+	return (delta > PM_THRESHOLD);
+}
diff --git a/blossom5-v2.05.src/PMexpand.cpp b/blossom5-v2.05.src/PMexpand.cpp
new file mode 100644
index 0000000..79ba684
--- /dev/null
+++ b/blossom5-v2.05.src/PMexpand.cpp
@@ -0,0 +1,304 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "PMimplementation.h"
+
+
+
+
+inline void PerfectMatching::ProcessSelfloop(Node* b, Edge* a)
+{
+	int dir;
+	Node* j;
+	Node* prev[2];
+	for (dir=0; dir<2; dir++)
+	{
+		j = a->head[dir];
+		GET_PENULTIMATE_BLOSSOM(j);
+		prev[dir] = j;
+	}
+	if (prev[0] != prev[1])
+	{
+		ADD_EDGE(prev[0], a, 1);
+		ADD_EDGE(prev[1], a, 0);
+		a->slack -= 2*prev[0]->blossom_eps;
+	}
+	else
+	{
+		a->next[0] = prev[0]->blossom_selfloops;
+		prev[0]->blossom_selfloops = a;
+	}
+
+}
+
+
+
+void PerfectMatching::Expand(Node* b)
+{
+	assert(b->is_blossom);
+	assert(b->is_outer);
+	assert(b->flag == 1);
+
+	double start_time = get_time();
+
+	Node* i;
+	Node* j;
+	Node* k;
+	Edge* a;
+	EdgeIterator I;
+	int dir;
+	ExpandTmpItem* tmp_item;
+	Tree* t = b->tree;
+	REAL eps = t->eps;
+	Edge* a_augment = NULL;
+
+	GET_TREE_PARENT(b, i);
+	a = ARC_TO_EDGE_PTR(b->tree_parent);
+	dir = ARC_TO_EDGE_DIR(b->tree_parent);
+
+	j = a->head0[1-dir];
+	GET_PENULTIMATE_BLOSSOM(j);
+	MOVE_EDGE(b, j, a, dir);
+
+	a = ARC_TO_EDGE_PTR(b->match);
+	dir = ARC_TO_EDGE_DIR(b->match);
+	k = a->head0[1-dir];
+	GET_PENULTIMATE_BLOSSOM(k);
+	MOVE_EDGE(b, k, a, dir);
+
+	i = ARC_HEAD(k->blossom_sibling);
+	while ( 1 )
+	{
+		tmp_item = expand_tmp_list->New();
+		tmp_item->i = i; tmp_item->blossom_parent = i->blossom_parent; tmp_item->blossom_grandparent = i->blossom_grandparent;
+		i->flag = 2;
+
+		// blossom_selfloops
+		i->is_outer = 1;
+		while ((a=i->blossom_selfloops))
+		{
+			i->blossom_selfloops = a->next[0];
+			ProcessSelfloop(i, a);
+		}
+		i->is_outer = 0;
+
+		if (i == k) break;
+		i->match = i->blossom_sibling;
+		j = ARC_HEAD(i->match);
+		tmp_item = expand_tmp_list->New();
+		tmp_item->i = j; tmp_item->blossom_parent = j->blossom_parent; tmp_item->blossom_grandparent = j->blossom_grandparent;
+		j->flag = 2;
+
+		// blossom_selfloops
+		j->is_outer = 1;
+		while ((a=j->blossom_selfloops))
+		{
+			j->blossom_selfloops = a->next[0];
+			ProcessSelfloop(j, a);
+		}
+		j->is_outer = 0;
+
+		j->match = ARC_REV(i->match);
+		i = ARC_HEAD(j->blossom_sibling);
+	}
+	k->match = b->match;
+	i = ARC_TAIL(b->tree_parent);
+	Arc* aa = i->blossom_sibling;
+	i->flag = 1; i->tree = b->tree; i->y += b->tree->eps;
+	i->tree_parent = b->tree_parent;
+	if (i != k)
+	{
+		Node** i_ptr;
+		if (i->match == aa)
+		{
+			i = ARC_HEAD(i->match);
+			i_ptr = &j;
+			while ( 1 )
+			{
+				aa = i->blossom_sibling;
+				i->flag = 0; i->tree = b->tree; i->y -= t->eps;
+				*i_ptr = i;
+				i_ptr = &i->first_tree_child;
+				i->tree_sibling_prev = i;
+				i->tree_sibling_next = NULL;
+				i = ARC_HEAD(aa);
+				i->flag = 1; i->tree = b->tree; i->y += t->eps;
+				i->tree_parent = ARC_REV(aa);
+				if (i == k) break;
+				i = ARC_HEAD(i->match);
+			}
+			*i_ptr = ARC_HEAD(k->match);
+		}
+		else
+		{
+			i = k;
+			j = ARC_HEAD(k->match);
+			do
+			{
+				i->tree_parent = i->blossom_sibling;
+				i->flag = 1; i->tree = b->tree; i->y += b->tree->eps;
+				i = ARC_HEAD(i->tree_parent);
+				i->flag = 0; i->tree = b->tree; i->y -= b->tree->eps;
+				i->first_tree_child = j;
+				j = i;
+				i->tree_sibling_prev = i;
+				i->tree_sibling_next = NULL;
+				i = ARC_HEAD(i->match);
+			} while ( i->flag != 1 );
+		}
+		i = ARC_HEAD(k->match);
+
+		j->tree_sibling_prev = i->tree_sibling_prev;
+		j->tree_sibling_next = i->tree_sibling_next;
+		if (i->tree_sibling_prev->tree_sibling_next) i->tree_sibling_prev->tree_sibling_next = j;
+		else ARC_HEAD(b->tree_parent)->first_tree_child = j;
+		if (i->tree_sibling_next) i->tree_sibling_next->tree_sibling_prev = j;
+		else ARC_HEAD(b->tree_parent)->first_tree_child->tree_sibling_prev = j;
+
+		i->tree_sibling_prev = i;
+		i->tree_sibling_next = NULL;
+	}
+
+	// go through inner arcs
+	i = k;
+	while ( 1 )
+	{
+		// "-" node
+		if (i->is_blossom)
+		{
+			a = ARC_TO_EDGE_PTR(i->match);
+			REAL tmp = a->slack; a->slack = i->y; i->y = tmp;
+			t->pq_blossoms.Add(a);
+		}
+		FOR_ALL_EDGES(i, a, dir, I)
+		{
+			j = a->head[dir];
+			if (j->flag != 0) a->slack -= eps;
+		}
+		i->is_processed = 1;
+		if (i->tree_parent == b->tree_parent) break;
+		i = ARC_HEAD(i->tree_parent);
+		// "+" node
+		FOR_ALL_EDGES(i, a, dir, I)
+		{
+			j = a->head[dir];
+			if (j->flag == 2)
+			{
+				a->slack += eps;
+				t->pq0.Add(a);
+			}
+			else if (j->flag == 0 && i < j)
+			{
+				a->slack += 2*eps;
+				t->pq00.Add(a);
+			}
+		}
+		i->is_processed = 1;
+		i = ARC_HEAD(i->match);
+	}
+
+	// go through boundary arcs
+	for (tmp_item=expand_tmp_list->ScanFirst(); tmp_item; tmp_item=expand_tmp_list->ScanNext())
+	{
+		i = tmp_item->i;
+		j = tmp_item->blossom_parent; tmp_item->blossom_parent = i->blossom_parent; i->blossom_parent = j;
+		j = tmp_item->blossom_grandparent; tmp_item->blossom_grandparent = i->blossom_grandparent; i->blossom_grandparent = j;
+	}
+	for (dir=0; dir<2; dir++)
+	{
+		if (!b->first[dir]) continue;
+		b->first[dir]->prev[dir]->next[dir] = NULL;
+
+		Edge* a_next;
+		for (a=b->first[dir]; a; a=a_next)
+		{
+			a_next = a->next[dir];
+			i = a->head0[1-dir];
+			GET_PENULTIMATE_BLOSSOM2(i);
+			ADD_EDGE(i, a, dir);
+			GET_OUTER_HEAD(a, dir, j);
+
+			if (i->flag == 1) continue;
+
+			if (j->flag == 0 && j->tree != t) j->tree->pq_current->pq01[1-j->tree->dir_current].Remove(a, pq_buf);
+
+			if (i->flag == 2)
+			{
+				a->slack += eps;
+				if (j->flag == 0) j->tree->pq0.Add(a);
+			}
+			else
+			{
+				a->slack += 2*eps;
+				if      (j->flag == 2) t->pq0.Add(a);
+				else if (j->flag == 0)
+				{
+					if (j->tree != t)
+					{
+						if (!j->tree->pq_current) AddTreeEdge(t, j->tree);
+						if (a->slack <= j->tree->eps + eps) a_augment = a;
+					}
+					j->tree->pq_current->pq00.Add(a);
+				}
+				else if (j->tree != t)
+				{
+					if (!j->tree->pq_current) AddTreeEdge(t, j->tree);
+					j->tree->pq_current->pq01[j->tree->dir_current].Add(a); 
+				}
+					
+			}
+		}
+	}
+	for (tmp_item=expand_tmp_list->ScanFirst(); tmp_item; tmp_item=expand_tmp_list->ScanNext())
+	{
+		i = tmp_item->i;
+		i->blossom_parent = tmp_item->blossom_parent;
+		i->blossom_grandparent = tmp_item->blossom_grandparent;
+		i->is_outer = 1;
+	}
+	expand_tmp_list->Reset();
+
+	b->is_removed = 1;
+	b->tree_sibling_next = removed_first;
+	removed_first = b;
+	removed_num ++;
+	if (4*removed_num > node_num) FreeRemoved();
+
+	blossom_num --;
+	stat.expand_count ++;
+
+	stat.expand_time += get_time() - start_time;
+
+	if (a_augment) Augment(a_augment);
+}
+
+void PerfectMatching::FreeRemoved()
+{
+	Node* i0;
+	Node* i;
+	for (i0=nodes; i0<nodes+node_num; i0++)
+	{
+		for (i=i0; !i->is_outer && !i->is_marked; i=i->blossom_parent)
+		{
+			i->is_marked = 1;
+			if (i->blossom_grandparent->is_removed) i->blossom_grandparent = i->blossom_parent;
+		}
+	}
+	for (i0=nodes; i0<nodes+node_num; i0++)
+	{
+		for (i=i0; !i->is_outer && i->is_marked; i=i->blossom_parent)
+		{
+			i->is_marked = 0;
+		}
+	}
+
+	while ((i=removed_first))
+	{
+		removed_first = i->tree_sibling_next;
+		blossoms->Delete(i);
+		removed_num --;
+	}
+
+	assert(removed_num == 0);
+}
+
diff --git a/blossom5-v2.05.src/PMimplementation.h b/blossom5-v2.05.src/PMimplementation.h
new file mode 100644
index 0000000..b949c14
--- /dev/null
+++ b/blossom5-v2.05.src/PMimplementation.h
@@ -0,0 +1,356 @@
+/*
+    PMimplementation.h
+
+    Copyright 2008 Vladimir Kolmogorov (vnk@ist.ac.at)
+
+    This software can be used for research purposes only. Commercial use is prohibited.
+    Public redistribution of the code or its derivatives is prohibited.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef ASKHAKJSNTJAKSNBBAVASRA
+#define ASKHAKJSNTJAKSNBBAVASRA
+
+#ifdef _MSC_VER
+#pragma warning(disable: 4311)
+#pragma warning(disable: 4312)
+#endif
+
+#include "PerfectMatching.h"
+#include "timer.h"
+#include "PQ.h"
+#include "LCA.h"
+
+
+#define LCA_REPAIRS
+
+
+
+
+#define IS_INT ( ((REAL)1 / 2) == 0 )
+#define COST_FACTOR 2
+
+#define PM_THRESHOLD ((REAL)1e-12)
+
+struct PerfectMatching::Node
+{
+	unsigned int	is_outer : 1; // 0 - the node is contained in another blossom, 1 - it's the outermost node
+	unsigned int	flag : 2; // 0 corresponds to +, 1 corresponds to -, 2 corresponds to a free node
+	unsigned int	is_tree_root : 1;
+	unsigned int	is_processed : 1;
+	unsigned int	is_blossom : 1;
+	unsigned int	is_marked : 1;
+	unsigned int	is_removed : 1;
+
+	Edge*		first[2];
+	union
+	{
+		Arc*	match; // used if not a tree root (is_tree_root = 0) or it's an inner node (is_outer == 0)
+		Node*	blossom_grandparent;
+	};
+	REAL		y;
+
+	union
+	{
+		struct // used when is_outer = 0
+		{
+			Arc*	blossom_sibling;
+			Node*	blossom_parent;
+			union
+			{
+				Edge*	blossom_selfloops;
+				Node*	blossom_ptr; // used in repairs
+#ifdef LCA_REPAIRS
+				int		lca_preorder; // used in repairs
+#endif
+			};
+			REAL blossom_eps; // stores 'eps' of the tree at the moment when the node was shrunk into a blossom.
+			                  // (it is used for determining slacks of self-loops)
+		};
+		struct // used when is_outer = 1
+		{
+			union
+			{
+				struct // used for "+" nodes (flag = 0)
+				{
+					Node*	first_tree_child;
+					Node*	tree_sibling_prev; // circular list (with one exception: parent->first_tree_child->tree_sibling_prev->tree_sibling_next is NULL)
+					Node*	tree_sibling_next;
+				};
+				Arc*	tree_parent; // used for "-" nodes (flag = 1)
+			};
+			union
+			{
+				Tree*	tree;
+				Edge*	best_edge;  // used during InitGlobal() for non-tree nodes
+#ifdef LCA_REPAIRS
+				int		lca_size; // used in repairs
+				LCATreeX*	lca;  // used in repairs
+#endif
+			};
+		};
+	};
+};
+
+struct PerfectMatching::Edge : PriorityQueue<REAL>::Item
+{
+	Node*	head[2];
+	Node*	head0[2];
+	Edge*	next[2];
+	Edge*	prev[2];
+};
+
+typedef unsigned long POINTER_TYPE;
+// if the declaration below fails, set POINTER_TYPE to be the appropriate integer type of the same length as (void*)
+extern char dummy_array[2*(sizeof(void*)==sizeof(POINTER_TYPE))-1];
+
+#define ARC_TO_EDGE_PTR(a)       ( (Edge*) ( ((POINTER_TYPE)(a)) & (~1)      ) )
+#define ARC_TO_EDGE_DIR(a)       ( (int)   ( ((POINTER_TYPE)(a)) & 1         ) )
+#define EDGE_DIR_TO_ARC(a, dir)  ( (Arc*)  ( (char*)(a) + (dir)) )
+
+#define ARC_REV(a) ( (Arc*) ( ((POINTER_TYPE)(a)) ^ 1 ) )
+
+#define ARC_TAIL(a)  (ARC_TO_EDGE_PTR(a)->head [1-ARC_TO_EDGE_DIR(a)])
+#define ARC_TAIL0(a) (ARC_TO_EDGE_PTR(a)->head0[1-ARC_TO_EDGE_DIR(a)])
+#define ARC_HEAD(a)  (ARC_TO_EDGE_PTR(a)->head [ARC_TO_EDGE_DIR(a)])
+#define ARC_HEAD0(a) (ARC_TO_EDGE_PTR(a)->head0[ARC_TO_EDGE_DIR(a)])
+
+struct PerfectMatching::PQPointers
+{
+	PriorityQueue<REAL> pq00; // plus-plus edges
+	union
+	{
+		PriorityQueue<REAL> pq01[2]; // plus-minus, minus-plus edges. Used for tree edges.
+		struct // used for trees.
+		{
+			PriorityQueue<REAL> pq0; // plus-free edges
+			PriorityQueue<REAL> pq_blossoms;
+		};
+	};
+};
+
+struct PerfectMatching::Tree : PQPointers
+{
+	REAL		eps;
+	TreeEdge*	first[2];
+	Node*		root;
+
+	PQPointers*	pq_current;
+	int			dir_current;
+
+	/////////////////////////////////////////
+	// used while computing dual updates
+	REAL		eps_delta;
+	Tree*		next;
+	union
+	{
+		int			id;
+		TreeEdge*	dfs_parent;
+	};
+};
+
+struct PerfectMatching::TreeEdge : PQPointers
+{
+	Tree*		head[2];
+	TreeEdge*	next[2];
+};
+
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////////////////
+
+
+#define GET_PENULTIMATE_BLOSSOM(j)\
+	{\
+		Node* jtmp1 = j;\
+		while ( 1 )\
+		{\
+			if (!j->blossom_grandparent->is_outer) j = j->blossom_grandparent;\
+			else if (j->blossom_grandparent != j->blossom_parent) j->blossom_grandparent = j->blossom_parent;\
+			else break;\
+		}\
+		Node* jtmp2;\
+		for ( ; jtmp1!=j; jtmp1=jtmp2)\
+		{\
+			jtmp2 = jtmp1->blossom_grandparent;\
+			jtmp1->blossom_grandparent = j;\
+		}\
+	}
+#define GET_PENULTIMATE_BLOSSOM2(j)\
+	{\
+		Node* jtmp1 = j;\
+		Node* jtmp_prev = NULL;\
+		while ( 1 )\
+		{\
+			if (!j->blossom_grandparent->is_outer) { jtmp_prev = j; j = j->blossom_grandparent; }\
+			else if (j->blossom_grandparent != j->blossom_parent) j->blossom_grandparent = j->blossom_parent;\
+			else break;\
+		}\
+		if (jtmp_prev)\
+		{\
+			Node* jtmp2;\
+			for ( ; jtmp1!=jtmp_prev; jtmp1=jtmp2)\
+			{\
+				jtmp2 = jtmp1->blossom_grandparent;\
+				jtmp1->blossom_grandparent = jtmp_prev;\
+			}\
+		}\
+	}
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////////////////
+
+
+struct PerfectMatching::EdgeIterator
+{
+	Edge*	a_last;
+	int		start_flag;
+};
+
+#define FOR_ALL_EDGES(i, a, dir, I)\
+	for ( dir = (i->first[0]) ? 0 : 1, I.a_last = a = i->first[dir], I.start_flag = (a) ? 0 : 1;\
+		a != I.a_last || (I.start_flag ++ == 0) || (dir ++ == 0 && (I.a_last = a = i->first[1]));\
+		a = a->next[dir] )
+
+#define CONTINUE_FOR_ALL_EDGES(i, a, dir, I)\
+	for ( a = a->next[dir];\
+		a != I.a_last || (I.start_flag ++ == 0) || (dir ++ == 0 && (I.a_last = a = i->first[1]));\
+		a = a->next[dir] )
+
+#define REMOVE_EDGE(i, a, dir)\
+	{\
+		if ((a)->prev[dir]==(a)) (i)->first[dir] = NULL;\
+		else\
+		{\
+			(a)->prev[dir]->next[dir] = (a)->next[dir];\
+			(a)->next[dir]->prev[dir] = (a)->prev[dir];\
+			(i)->first[dir] = (a)->next[dir];\
+		}\
+	}
+
+#define ADD_EDGE(i, a, dir)\
+	{\
+		if ((i)->first[dir])\
+		{\
+			(a)->prev[dir] = (i)->first[dir]->prev[dir];\
+			(a)->next[dir] = (i)->first[dir];\
+			(i)->first[dir]->prev[dir]->next[dir] = (a);\
+			(i)->first[dir]->prev[dir] = (a);\
+		}\
+		else (i)->first[dir] = (a)->prev[dir] = (a)->next[dir] = (a);\
+		(a)->head[1-(dir)] = (i);\
+	}
+
+#define MOVE_EDGE(i_old, i_new, a, dir)\
+	{\
+		REMOVE_EDGE(i_old, a, dir);\
+		ADD_EDGE(i_new, a, dir);\
+	}
+
+#define GET_OUTER_HEAD(a, dir, j)\
+	{\
+		j = (a)->head[dir];\
+		if (!j->is_outer)\
+		{\
+			Node* j_orig = j;\
+			GET_PENULTIMATE_BLOSSOM(j);\
+			j = j->blossom_parent;\
+			int dir_rev = 1 - (dir);\
+			MOVE_EDGE(j_orig, j, a, dir_rev);\
+		}\
+	}
+
+#define GET_TREE_PARENT(child, parent)\
+	{\
+		Arc* a = (child)->tree_parent;\
+		Edge* e = ARC_TO_EDGE_PTR(a);\
+		int dir = ARC_TO_EDGE_DIR(a);\
+		GET_OUTER_HEAD(e, dir, parent);\
+	}
+
+///////////////////////////////////////////////////////////////////////////////////////////////
+
+struct PerfectMatching::TreeEdgeIterator
+{
+	TreeEdge**	e_ptr;
+};
+
+#define FOR_ALL_TREE_EDGES(t, e, dir)\
+	for ( dir = (t->first[0]) ? 0 : 1, e = t->first[dir];\
+	      e || (dir ++ == 0 && (e = t->first[1]));\
+	      e = e->next[dir] )
+
+#define FOR_ALL_TREE_EDGES_X(t, e, dir, T)\
+	for ( dir = (t->first[0]) ? 0 : 1, T.e_ptr = &t->first[dir], e = *T.e_ptr;\
+	      e || (dir ++ == 0 && (e = *(T.e_ptr = &t->first[1])));\
+	      e = *T.e_ptr )\
+	if (e->head[dir] == NULL) { *T.e_ptr = e->next[dir]; tree_edges->Delete(e); }\
+	else if ((T.e_ptr = &e->next[dir]))
+
+///////////////////////////////////////////////////////////////////////////////////////////////
+
+#define MOVE_NODE_IN_TREE(i)\
+	{\
+		if ((i)->first_tree_child) (i) = (i)->first_tree_child;\
+		else\
+		{\
+			while (!(i)->is_tree_root && !(i)->tree_sibling_next) { (i) = ARC_HEAD((i)->match); GET_TREE_PARENT(i, i); }\
+			if ((i)->is_tree_root) break;\
+			(i) = (i)->tree_sibling_next;\
+		}\
+	}
+
+// i=parent, j=child
+#define ADD_TREE_CHILD(i, j)\
+	{\
+		(j)->flag = 0;\
+		(j)->tree = (i)->tree;\
+		(j)->first_tree_child = NULL;\
+		(j)->tree_sibling_next = (i)->first_tree_child;\
+		if ((i)->first_tree_child)\
+		{\
+			(j)->tree_sibling_prev = (i)->first_tree_child->tree_sibling_prev;\
+			(i)->first_tree_child->tree_sibling_prev = j;\
+		}\
+		else\
+		{\
+			(j)->tree_sibling_prev = j;\
+		}\
+		(i)->first_tree_child = j;\
+	}
+
+
+#define REMOVE_FROM_TREE(i)\
+	{\
+		if ((i)->tree_sibling_next) (i)->tree_sibling_next->tree_sibling_prev = (i)->tree_sibling_prev;\
+		else\
+		{\
+			Node* i_NEXT = ARC_HEAD((i)->match); i_NEXT = ARC_HEAD(i_NEXT->tree_parent); i_NEXT = i_NEXT->first_tree_child;\
+			i_NEXT->tree_sibling_prev = (i)->tree_sibling_prev;\
+		}\
+		if ((i)->tree_sibling_prev->tree_sibling_next) (i)->tree_sibling_prev->tree_sibling_next = (i)->tree_sibling_next;\
+		else\
+		{\
+			Node* i_PARENT = ARC_HEAD((i)->match); i_PARENT = ARC_HEAD(i_PARENT->tree_parent);\
+			i_PARENT->first_tree_child = (i)->tree_sibling_next;\
+		}\
+	}
+
+
+#endif
+
diff --git a/blossom5-v2.05.src/PMinit.cpp b/blossom5-v2.05.src/PMinit.cpp
new file mode 100644
index 0000000..cbcbb99
--- /dev/null
+++ b/blossom5-v2.05.src/PMinit.cpp
@@ -0,0 +1,544 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "PMimplementation.h"
+
+
+void PerfectMatching::InitGreedy(bool allocate_trees)
+{
+	Node* i;
+	int dir;
+	Edge* a;
+	EdgeIterator I;
+	Tree* t = NULL;
+	Node* last_root = &nodes[node_num];
+	REAL slack_min;
+
+	for (i=nodes; i<nodes+node_num; i++) i->y = PM_INFTY;
+	for (a=edges; a<edges+edge_num; a++)
+	{
+		if (a->head[0]->y > a->slack) a->head[0]->y = a->slack;
+		if (a->head[1]->y > a->slack) a->head[1]->y = a->slack;
+	}
+	for (a=edges; a<edges+edge_num; a++)
+	{
+		i = a->head[0];
+		if (!i->is_outer)
+		{
+			i->is_outer = 1;
+			i->y /= 2;
+		}
+		a->slack -= i->y;
+		i = a->head[1];
+		if (!i->is_outer)
+		{
+			i->is_outer = 1;
+			i->y /= 2;
+		}
+		a->slack -= i->y;
+	}
+
+	tree_num = node_num;
+	for (i=nodes; i<nodes+node_num; i++)
+	{
+		if (i->flag == 2) continue;
+		slack_min = PM_INFTY;
+		FOR_ALL_EDGES(i, a, dir, I) if (slack_min > a->slack) slack_min = a->slack;
+		i->y += slack_min;
+		FOR_ALL_EDGES(i, a, dir, I)
+		{
+			if (a->slack <= slack_min && i->flag == 0 && a->head[dir]->flag == 0)
+			{
+				i->flag = 2;
+				a->head[dir]->flag = 2;
+				i->match = EDGE_DIR_TO_ARC(a, dir);
+				a->head[dir]->match = EDGE_DIR_TO_ARC(a, 1-dir);
+				tree_num -= 2;
+			}
+			a->slack -= slack_min;
+		}
+	}
+	if (allocate_trees)
+	{
+		if (tree_num > tree_num_max)
+		{
+			if (trees) free(trees);
+			tree_num_max = tree_num;
+			trees = (Tree*) malloc(tree_num_max*sizeof(Tree));
+		}
+		t = trees;
+	}
+	for (i=nodes; i<nodes+node_num; i++)
+	{
+		if (i->flag != 0) continue;
+		i->is_tree_root = 1;
+		i->first_tree_child = NULL;
+		i->tree_sibling_prev = last_root;
+		last_root->tree_sibling_next = i;
+		last_root = i;
+		if (allocate_trees)
+		{
+			i->tree = t;
+			t->root = i;
+			t->eps = 0;
+			t->first[0] = t->first[1] = NULL;
+			t->pq_current = NULL;
+			t->pq00.Reset();
+			t->pq0.Reset();
+			t->pq_blossoms.Reset();
+			t ++;
+		}
+	}
+	last_root->tree_sibling_next = NULL;
+}
+
+////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////
+
+PerfectMatching::Node* PerfectMatching::FindBlossomRootInit(Edge* a0)
+{
+	Node* i;
+	Node* j;
+	Node* _i[2];
+	Node* r;
+	int branch;
+
+	_i[0] = ARC_HEAD(a0);
+	_i[1] = ARC_TAIL(a0);
+	branch = 0;
+	while ( 1 )
+	{
+		if (!_i[branch]->is_outer) 
+		{
+			r = _i[branch]; 
+			j = _i[1-branch];
+			break; 
+		}
+		_i[branch]->is_outer = 0;
+		if (_i[branch]->is_tree_root)
+		{
+			j = _i[branch];
+			i = _i[1-branch];
+			while (i->is_outer)
+			{
+				i->is_outer = 0;
+				i = ARC_HEAD(i->match);
+				i->is_outer = 0;
+				i = ARC_HEAD(i->tree_parent);
+			}
+			r = i;
+			break;
+		}
+		i = ARC_HEAD(_i[branch]->match);
+		i->is_outer = 0;
+		_i[branch] = ARC_HEAD(i->tree_parent);
+		branch = 1 - branch;
+	}
+	i = r;
+	while ( i != j )
+	{
+		i = ARC_HEAD(i->match);
+		i->is_outer = 1;
+		i = ARC_HEAD(i->tree_parent);
+		i->is_outer = 1;
+	}
+	return r;
+}
+
+void PerfectMatching::ShrinkInit(Edge* a0, Node* tree_root)
+{
+	int branch, flag;
+	Node* i;
+	Node* j;
+	Node* r;
+	Arc* a_prev;
+	Arc* aa;
+
+	tree_root->flag = 2;
+	i = tree_root->first_tree_child;
+	if ( i )
+	while ( 1 )
+	{
+		ARC_HEAD(i->match)->flag = 2;
+		i->flag = 2;
+
+		MOVE_NODE_IN_TREE(i);
+	}
+
+	r = FindBlossomRootInit(a0);
+
+	if ( !r->is_tree_root )
+	{
+		j = ARC_HEAD(r->match);
+		j->match = aa = j->tree_parent;
+		i = ARC_HEAD(aa);
+		while ( !i->is_tree_root )
+		{
+			j = ARC_HEAD(i->match);
+			i->match = ARC_REV(aa);
+			j->match = aa = j->tree_parent;
+			i = ARC_HEAD(aa);
+		}
+		i->match = ARC_REV(aa);
+	}
+
+	tree_root->is_tree_root = 0;
+
+	branch = 0;
+	flag = 0;
+	a_prev = EDGE_DIR_TO_ARC(a0, 0);
+	i = ARC_HEAD(a_prev);
+	while ( 1 )
+	{
+		Arc* a_next = (flag == 0) ? i->match : i->tree_parent;
+		flag = 1 - flag;
+		i->flag = 0;
+		i->match = NULL;
+		if (branch == 0)
+		{
+			i->blossom_sibling = a_next;
+			if (i == r)
+			{
+				branch = 1;
+				flag = 0;
+				a_prev = ARC_REV(a0);
+				i = ARC_HEAD(a_prev);
+				if (i == r) break;
+			}
+			else
+			{
+				a_prev = i->blossom_sibling;
+				i = ARC_HEAD(a_prev);
+			}
+		}
+		else
+		{
+			i->blossom_sibling = ARC_REV(a_prev);
+			a_prev = a_next;
+			i = ARC_HEAD(a_prev);
+			if (i == r) break;
+		}
+	}
+	i->blossom_sibling = ARC_REV(a_prev);
+}
+
+void PerfectMatching::ExpandInit(Node* k)
+{
+	Node* i = ARC_HEAD(k->blossom_sibling);
+	Node* j;
+
+	while ( 1 )
+	{
+		i->flag = 2; i->is_outer = 1;
+		if (i == k) break;
+		i->match = i->blossom_sibling;
+		j = ARC_HEAD(i->match);
+		j->flag = 2; j->is_outer = 1;
+		j->match = ARC_REV(i->match);
+		i = ARC_HEAD(j->blossom_sibling);
+	}
+}
+
+void PerfectMatching::AugmentBranchInit(Node* i0, Node* r)
+{
+	Node* tree_root_prev = r->tree_sibling_prev;
+	Node* i;
+	Node* j;
+	Arc* aa;
+
+	r->flag = 2;
+	i = r->first_tree_child;
+	if ( i )
+	while ( 1 )
+	{
+		ARC_HEAD(i->match)->flag = 2;
+		i->flag = 2;
+
+		MOVE_NODE_IN_TREE(i);
+	}
+	i = i0;
+	if ( !i0->is_tree_root )
+	{
+		j = ARC_HEAD(i0->match);
+		j->match = aa = j->tree_parent;
+		i = ARC_HEAD(aa);
+		while ( !i->is_tree_root )
+		{
+			j = ARC_HEAD(i->match);
+			i->match = ARC_REV(aa);
+			j->match = aa = j->tree_parent;
+			i = ARC_HEAD(aa);
+		}
+		i->match = ARC_REV(aa);
+	}
+	r->is_tree_root = 0;
+	tree_root_prev->tree_sibling_next = r->tree_sibling_next;
+	if (r->tree_sibling_next) r->tree_sibling_next->tree_sibling_prev = tree_root_prev;
+	tree_num --;
+}
+
+////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////
+
+// true_slack(a) = slack(a) + ...
+
+// i->flag=0, i->is_processed=1:              true_slack -= eps
+// i->flag=1, i->match->head->is_processed=1: true_slack += eps - slack(i->match)
+
+void PerfectMatching::InitGlobal()
+{
+	Node* i;
+	Node* j;
+	Node* r;
+	Node* r2;
+	Node* r3 = NULL; // initialize to prevent compiler warning
+	Edge* a;
+	EdgeIterator I;
+	int dir;
+	Tree TREE;
+	enum { NONE, AUGMENT, SHRINK } flag;
+
+	InitGreedy();
+
+	for (i=nodes; i<nodes+node_num; i++) i->best_edge = NULL;
+
+	PriorityQueue<REAL> pq;
+
+	for (r=nodes[node_num].tree_sibling_next; r; )
+	{
+		r2 = r->tree_sibling_next;
+		if (r2) r3 = r2->tree_sibling_next;
+		i = r;
+
+		pq.Reset();
+
+		r->tree = &TREE;
+
+		REAL eps = 0;
+		Arc* critical_arc = NULL;
+		REAL critical_eps = PM_INFTY; 
+		flag = NONE;
+		Node* branch_root = i;
+	
+		while ( 1 )
+		{
+			i->is_processed = 1;
+			i->y -= eps;
+			if (!i->is_tree_root) ARC_HEAD(i->match)->y += eps;
+
+			FOR_ALL_EDGES(i, a, dir, I)
+			{
+				a->slack += eps;
+				j = a->head[dir];
+
+				if (j->tree == &TREE)
+				{
+					// same tree
+					if (j->flag == 0)
+					{
+						REAL slack = a->slack;
+						if (!j->is_processed) slack += eps;
+						if (2*critical_eps > slack || critical_arc == NULL)
+						{
+							flag = SHRINK;
+							critical_eps = slack/2;
+							critical_arc = EDGE_DIR_TO_ARC(a, dir);
+							if (critical_eps <= eps) break;
+							//pq.DecreaseUpperBound(critical_eps);
+						}
+					}
+				}
+				else if (j->flag == 0)
+				{
+					// different tree
+					if (critical_eps >= a->slack || critical_arc == NULL)
+					{
+						flag = AUGMENT;
+						critical_eps = a->slack;
+						critical_arc = EDGE_DIR_TO_ARC(a, dir);
+						if (critical_eps <= eps) break;
+						//pq.DecreaseUpperBound(critical_eps);
+					}
+				}
+				else
+				{
+					// free node
+					if (a->slack > eps)
+					{
+						if (a->slack < critical_eps)
+						{
+							if (j->best_edge == NULL)
+							{
+								j->best_edge = a;
+								pq.Add(a);
+							}
+							else
+							{
+								if (a->slack < j->best_edge->slack)
+								{
+									pq.Decrease(j->best_edge, a, pq_buf);
+									j->best_edge = a;
+								}
+							}
+						}
+					}
+					else
+					{
+						assert(j->flag == 2 && !j->is_blossom && !ARC_HEAD(j->match)->is_blossom);
+						if (j->best_edge) pq.Remove(j->best_edge, pq_buf);
+						j->flag = 1;
+						j->tree = i->tree;
+						j->tree_parent = EDGE_DIR_TO_ARC(a, 1-dir);
+						j = ARC_HEAD(j->match);
+						if (j->best_edge) pq.Remove(j->best_edge, pq_buf);
+						ADD_TREE_CHILD(i, j);
+					}
+				}
+			}
+
+			if (dir < 2 && a)
+			{
+				Edge* atmp = a;
+				int dirtmp = dir;
+				CONTINUE_FOR_ALL_EDGES(i, atmp, dirtmp, I) atmp->slack += eps;
+				break;
+			}
+
+			// move i
+			if (i->first_tree_child) i = i->first_tree_child;
+			else
+			{
+				while (i != branch_root && !i->tree_sibling_next) { i = ARC_HEAD(i->match); i = ARC_HEAD(i->tree_parent); }
+				if (i == branch_root)
+				{
+					PriorityQueue<REAL>::Item* q = pq.GetMin();
+					if (q == NULL || q->slack >= critical_eps)
+					{
+						eps = critical_eps;
+						break;
+					}
+					pq.Remove(q, pq_buf);
+					a = (Edge*)q;
+					dir = (a->head[0]->flag == 2) ? 0 : 1;
+					j = a->head[0]; 
+					Arc* aa = EDGE_DIR_TO_ARC(a, dir);
+					eps = a->slack;
+					assert(eps < critical_eps);
+
+					// continue growth
+					i = ARC_TAIL(aa);
+					j = ARC_HEAD(aa);
+
+					assert(j->flag == 2 && !j->is_blossom && !ARC_HEAD(j->match)->is_blossom);
+					j->flag = 1;
+					j->tree = i->tree;
+					j->tree_parent = ARC_REV(aa);
+					j = ARC_HEAD(j->match);
+					if (j->best_edge) pq.Remove(j->best_edge, pq_buf);
+					ADD_TREE_CHILD(i, j);
+					i = branch_root = j;
+					continue;
+				}
+				i = i->tree_sibling_next;
+			}
+		}
+
+		// update slacks
+		i = r;
+		while ( 1 )
+		{
+			if (i->is_processed)
+			{
+				i->y += eps;
+				if (!i->is_tree_root)
+				{
+					j = ARC_HEAD(i->match);
+					j->y -= eps;
+					REAL delta = eps - ARC_TO_EDGE_PTR(i->match)->slack;
+					FOR_ALL_EDGES(j, a, dir, I) a->slack += delta;
+					j->best_edge = NULL;
+				}
+				FOR_ALL_EDGES(i, a, dir, I)
+				{
+					if (!PriorityQueue<REAL>::isReset(a))
+					{
+						assert(a->head[dir]->flag == 2 && a->head[dir]->best_edge == a);
+						a->head[dir]->best_edge = NULL;
+						PriorityQueue<REAL>::ResetItem(a);
+					}
+					a->slack -= eps;
+				}
+
+				i->is_processed = 0;
+			}
+			else
+			{
+				if (!i->is_tree_root) ARC_HEAD(i->match)->best_edge = NULL;
+			}
+			i->best_edge = NULL;
+
+			MOVE_NODE_IN_TREE(i);
+		}
+
+		i = ARC_TAIL(critical_arc);
+		j = ARC_HEAD(critical_arc);
+		if (flag == SHRINK)
+		{
+			// shrink
+			ShrinkInit(ARC_TO_EDGE_PTR(critical_arc), r);
+		}
+		else
+		{
+			// augment
+			AugmentBranchInit(i, r);
+			if (j->is_outer)
+			{
+				AugmentBranchInit(j, j);
+			}
+			else
+			{
+				ExpandInit(j);
+				tree_num --;
+			}
+			i->match = critical_arc;
+			j->match = ARC_REV(critical_arc);
+		}
+
+		r = r2;
+		if (r && !r->is_tree_root) r = r3;
+	}
+
+	if (tree_num > tree_num_max)
+	{
+		if (trees) free(trees);
+		tree_num_max = tree_num;
+		trees = (Tree*) malloc(tree_num_max*sizeof(Tree));
+	}
+	Tree* t = trees;
+	for (r=nodes; r<nodes+node_num; r++)
+	{
+		if (!r->is_outer)
+		{
+			ExpandInit(r);
+			r->is_tree_root = 1;
+			r->flag = 0;
+			r->first_tree_child = NULL;
+			if (t == trees) { nodes[node_num].tree_sibling_next = r; r->tree_sibling_prev = &nodes[node_num]; }
+			else            { (t-1)->root->tree_sibling_next = r;    r->tree_sibling_prev = (t-1)->root; }
+			r->tree = t;
+			t->root = r;
+			t->eps = 0;
+			t->first[0] = t->first[1] = NULL;
+			t->pq_current = NULL;
+			t->pq00.Reset();
+			t->pq0.Reset();
+			t->pq_blossoms.Reset();
+			t ++;
+		}
+	}
+	assert(t == trees+tree_num);
+	if (t == trees) nodes[node_num].tree_sibling_next = NULL;
+	else            (t-1)->root->tree_sibling_next = NULL;
+}
diff --git a/blossom5-v2.05.src/PMinterface.cpp b/blossom5-v2.05.src/PMinterface.cpp
new file mode 100644
index 0000000..baab2d6
--- /dev/null
+++ b/blossom5-v2.05.src/PMinterface.cpp
@@ -0,0 +1,270 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "PMimplementation.h"
+
+
+PerfectMatching::PerfectMatching(int nodeNum, int edgeNumMax)
+	: node_num(nodeNum),
+	  edge_num(0),
+	  edge_num_max(edgeNumMax),
+	  trees(NULL),
+	  tree_num_max(0),
+	  removed_first(NULL),
+	  blossom_num(0),
+	  removed_num(0),
+	  first_solve(true)
+{
+	if (node_num & 1) { printf("# of nodes is odd: perfect matching cannot exist\n"); exit(1); }
+	nodes = (Node*) malloc((node_num+1)*sizeof(Node));
+	edges_orig = (char*) malloc(edge_num_max*sizeof(Edge)+1);
+	edges = (Edge*) ( ( ((POINTER_TYPE)edges_orig) & 1 ) ? (edges_orig + 1) : edges_orig );
+	memset(nodes, 0, (node_num+1)*sizeof(Node));
+
+	blossoms   = new DBlock<Node>(256);
+	tree_edges = new DBlock<TreeEdge>(256);
+	expand_tmp_list = new Block<ExpandTmpItem>(256);
+	pq_buf = PriorityQueue<REAL>::AllocateBuf();
+}
+
+
+void PerfectMatching::Save(char* filename, int format)
+{
+	if (!first_solve) { printf("Save() cannot be called after Solve()!\n"); exit(1); }
+	int e;
+	FILE* fp = fopen(filename, "w");
+	if (!fp) { printf("Can't open %s\n", filename); exit(1); }
+	if (format == 0)
+	{
+		fprintf(fp, "p edge %d %d\n", node_num, edge_num);
+		for (e=0; e<edge_num; e++)
+		{
+			fprintf(fp, "e %d %d %d\n", 1+(int)(edges[e].head0[1]-nodes), 1+(int)(edges[e].head0[0]-nodes), (int)edges[e].slack/COST_FACTOR);
+		}
+	}
+	else
+	{
+		fprintf(fp, "%d %d\n", node_num, edge_num);
+		for (e=0; e<edge_num; e++)
+		{
+			fprintf(fp, "%d %d %d\n", (int)(edges[e].head0[1]-nodes), (int)(edges[e].head0[0]-nodes), (int)edges[e].slack/COST_FACTOR);
+		}
+	}
+	fclose(fp);
+}
+
+
+PerfectMatching::~PerfectMatching()
+{
+	free(nodes);
+	free(edges_orig);
+	delete blossoms;
+	delete tree_edges;
+	delete expand_tmp_list;
+	if (trees) free(trees);
+	PriorityQueue<REAL>::DeallocateBuf(pq_buf);
+}
+
+
+PerfectMatching::EdgeId PerfectMatching::AddEdge(NodeId _i, NodeId _j, REAL cost)
+{
+	if (_i<0 || _i>=node_num || _j<0 || _j>node_num || _i==_j)
+	{
+		printf("wrong node id's! (%d,%d)\n", _i, _j); exit(1);
+	}
+	if (edge_num >= edge_num_max) ReallocateEdges();
+	Node* i = nodes + _i;
+	Node* j = nodes + _j;
+	Edge* a = edges + edge_num;
+
+	ADD_EDGE(i, a, 0);
+	ADD_EDGE(j, a, 1);
+	a->head0[0] = j;
+	a->head0[1] = i;
+
+	a->slack = cost*COST_FACTOR;
+	PriorityQueue<REAL>::ResetItem(a);
+
+	return edge_num ++;
+}
+
+int PerfectMatching::GetSolution(EdgeId e)
+{
+	assert(e>=0 && e<edge_num);
+	Edge* a = edges + e;
+	return (a->head0[1]->match == EDGE_DIR_TO_ARC(a, 0)) ? 1 : 0;
+}
+
+PerfectMatching::NodeId PerfectMatching::GetMatch(NodeId i)
+{
+	assert(i>=0 && i<node_num);
+	return (int)(ARC_HEAD0(nodes[i].match)-nodes);
+}
+
+void PerfectMatching::GetRealEndpoints(Edge* a, Node*& tail, Node*& head)
+{
+	Node* i;
+	Node* j;
+	int delta = 0;
+
+	for (i=a->head0[1]; !i->is_outer; i=i->blossom_parent, delta--) {}
+	for (j=a->head0[0]; !j->is_outer; j=j->blossom_parent, delta++) {}
+	if ( i == j )
+	{
+		i = a->head0[1];
+		j = a->head0[0];
+		while ( delta < 0 ) { i = i->blossom_parent; delta ++; }
+		while ( delta > 0 ) { j = j->blossom_parent; delta --; }
+		while ( i->blossom_parent != j->blossom_parent )
+		{
+			i = i->blossom_parent;
+			j = j->blossom_parent; 
+		}
+	}
+	tail = i;
+	head = j;
+	assert((i->is_outer && j->is_outer) || (i->blossom_parent==j->blossom_parent && !i->is_outer && !j->is_outer));
+}
+
+void PerfectMatching::ReallocateEdges()
+{
+	printf("Warning: reallocating edges. Increasing edge_num_max in the constructor may improve memory efficiency!\n");
+	edge_num_max = edge_num_max*3/2 + 16;
+	char* edges_orig_old = edges_orig;
+	Edge* edges_old = edges;
+	edges_orig = (char*) realloc(edges_orig_old, edge_num_max*sizeof(Edge)+1);
+	edges = (Edge*) ( ( ((POINTER_TYPE)edges_orig_old) & 1 ) ? (edges_orig + 1) : edges_orig );
+	if ( ((POINTER_TYPE)edges) & 1 )
+	{
+		char* edges_orig_old2 = edges_orig;
+		Edge* edges_old2 = edges;
+
+		edges_orig = (char*) malloc(edge_num_max*sizeof(Edge)+1);
+		edges = (Edge*) ( ( ((POINTER_TYPE)edges_orig_old) & 1 ) ? (edges_orig + 1) : edges_orig );
+		memcpy(edges, edges_old2, edge_num*sizeof(Edge));
+		free(edges_orig_old2);
+	}
+
+#define UPDATE_EDGE_PTR(ptr) ptr = (Edge*)((char*)(ptr) + ((char*)edges - (char*)edges_old))
+#define UPDATE_ARC_PTR(ptr) ptr = (Arc*)((char*)(ptr) + ((char*)edges - (char*)edges_old))
+
+	Node* i;
+	Edge* a;
+	for (a=edges; a<edges+edge_num; a++)
+	{
+		if (a->next[0]) UPDATE_EDGE_PTR(a->next[0]);
+		if (a->next[1]) UPDATE_EDGE_PTR(a->next[1]);
+		if (a->prev[0]) UPDATE_EDGE_PTR(a->prev[0]);
+		if (a->prev[1]) UPDATE_EDGE_PTR(a->prev[1]);
+	}
+	if (first_solve)
+	{
+		for (i=nodes; i<nodes+node_num; i++)
+		{
+			if (i->first[0]) UPDATE_EDGE_PTR(i->first[0]);
+			if (i->first[1]) UPDATE_EDGE_PTR(i->first[1]);
+		}
+	}
+	else
+	{
+		Node* i0;
+		for (i0=nodes; i0<nodes+node_num; i0++)
+		{
+			i = i0;
+			while ( 1 )
+			{
+				if (i->is_outer)
+				{
+					UPDATE_ARC_PTR(i->match);
+					if (i->first[0]) UPDATE_EDGE_PTR(i->first[0]);
+					if (i->first[1]) UPDATE_EDGE_PTR(i->first[1]);
+					break;
+				}
+				UPDATE_ARC_PTR(i->blossom_sibling);
+				if (i->first[0]) UPDATE_EDGE_PTR(i->first[0]);
+				if (i->first[1]) UPDATE_EDGE_PTR(i->first[1]);
+
+				i = i->blossom_parent;
+				if (i->is_outer) { if ( i->is_marked) break; i->is_marked = 1; }
+				else             { if (!i->is_marked) break; i->is_marked = 0; }
+			}
+		}
+		for (i0=nodes; i0<nodes+node_num; i0++)
+		{
+			i = i0;
+			while ( 1 )
+			{
+				if (i->is_outer) break;
+
+				i = i->blossom_parent;
+				if (i->is_outer) { if (!i->is_marked) break; i->is_marked = 0; }
+				else             { if ( i->is_marked) break; i->is_marked = 1; }
+			}
+		}
+	}
+}
+
+int PerfectMatching::GetBlossomNum()
+{
+	return blossom_num;
+}
+
+void PerfectMatching::GetDualSolution(int* blossom_parents, REAL* twice_y)
+{
+	int _i0, id = node_num;
+	int* child_ptr;
+	Node* i0;
+	Node* i;
+	int* tmp_array = new int[blossom_num];
+
+	int* tmp_array_ptr = tmp_array;
+	for (_i0=0, i0=nodes; _i0<node_num; _i0++, i0++)
+	{
+		twice_y[_i0] = i0->y;
+		if (i0->is_outer)
+		{
+			blossom_parents[_i0] = -1;
+			continue;
+		}
+		child_ptr = &blossom_parents[_i0];
+		i = i0->blossom_parent;
+		while ( 1 )
+		{
+			if (i->is_marked)
+			{
+				*child_ptr = i->lca_preorder;
+				break;
+			}
+			i->is_marked = 1;
+			*tmp_array_ptr ++ = i->lca_preorder;
+			*child_ptr = i->lca_preorder = id ++;
+			child_ptr = &blossom_parents[i->lca_preorder];
+			twice_y[i->lca_preorder] = i->y;
+			if (i->is_outer)
+			{
+				*child_ptr = -1;
+				break;
+			}
+			i = i->blossom_parent;
+		}
+	}
+
+	assert(id == node_num+blossom_num && tmp_array_ptr == tmp_array + blossom_num);
+
+	tmp_array_ptr = tmp_array;
+	for (_i0=0, i0=nodes; _i0<node_num; _i0++, i0++)
+	{
+		if (i0->is_outer) continue;
+		i = i0->blossom_parent;
+		while ( 1 )
+		{
+			if (!i->is_marked) break;
+			i->is_marked = 0;
+			i->lca_preorder = *tmp_array_ptr ++;
+			if (i->is_outer) break;
+			i = i->blossom_parent;
+		}
+	}
+
+	delete [] tmp_array;
+}
diff --git a/blossom5-v2.05.src/PMmain.cpp b/blossom5-v2.05.src/PMmain.cpp
new file mode 100644
index 0000000..0914fc6
--- /dev/null
+++ b/blossom5-v2.05.src/PMmain.cpp
@@ -0,0 +1,694 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "PMimplementation.h"
+
+
+void PerfectMatching::Finish()
+{
+
+#define IS_VALID_MATCH(i) ((Edge*)(i->match) >= edges && (Edge*)(i->match) < edges + edge_num)
+
+	Node* i0;
+	Node* i;
+	Node* j;
+	Node* k;
+	Node* b;
+	Node* b_prev;
+	Node* b_prev_prev;
+
+	for (i0=nodes; i0<nodes+node_num; i0++)
+	{
+		if (IS_VALID_MATCH(i0)) continue;
+		b_prev = NULL;
+		b = i0;
+		do
+		{
+			b->blossom_grandparent = b_prev;
+			b_prev = b;
+			b = b->blossom_parent;
+		} while (!IS_VALID_MATCH(b));
+
+		b_prev_prev = b_prev->blossom_grandparent;
+		while ( 1 )
+		{
+			for (k=ARC_TAIL0(b->match); k->blossom_parent!=b; k=k->blossom_parent) {}
+			k->match = b->match;
+			i = ARC_HEAD(k->blossom_sibling);
+			while ( i != k )
+			{
+				i->match = i->blossom_sibling;
+				j = ARC_HEAD(i->match);
+				j->match = ARC_REV(i->match);
+				i = ARC_HEAD(j->blossom_sibling);
+			}
+
+			b = b_prev;
+			if (!b->is_blossom) break;
+			b_prev = b_prev_prev;
+			b_prev_prev = b_prev->blossom_grandparent;
+		}
+	}
+}
+
+
+void PerfectMatching::AddTreeEdge(Tree* t0, Tree* t1)
+{
+	TreeEdge* e = tree_edges->New();
+	e->head[0] = t1;
+	e->head[1] = t0;
+	e->next[0] = t0->first[0];
+	t0->first[0] = e;
+	e->next[1] = t1->first[1];
+	t1->first[1] = e;
+
+	e->pq00.Reset();
+	e->pq01[0].Reset();
+	e->pq01[1].Reset();
+
+	t1->pq_current = e;
+	t1->dir_current = 0;
+}
+
+bool PerfectMatching::ProcessEdge00(Edge* a, bool update_boundary_edge)
+{
+	int dir;
+	Node* j;
+	Node* prev[2];
+	Node* last[2];
+	for (dir=0; dir<2; dir++)
+	{
+		if (a->head[dir]->is_outer) { prev[dir] = NULL; last[dir] = a->head[dir]; }
+		else
+		{
+			j = a->head[dir];
+			GET_PENULTIMATE_BLOSSOM(j);
+			prev[dir] = j;
+			last[dir] = prev[dir]->blossom_parent;
+			//assert(last[dir]->is_outer);
+		}
+	}
+
+	if (last[0] != last[1])
+	{
+		for (dir=0; dir<2; dir++)
+		{
+			j = a->head[dir];
+			if (j != last[dir]) { int dir_rev = 1 - dir; MOVE_EDGE(j, last[dir], a, dir_rev); }
+		}
+		if (update_boundary_edge) a->slack -= 2*a->head[0]->tree->eps;
+		return true;
+	}
+
+	if (prev[0] != prev[1])
+	{
+		for (dir=0; dir<2; dir++)
+		{
+			j = a->head[dir];
+			if (j != prev[dir]) { int dir_rev = 1 - dir; MOVE_EDGE(j, prev[dir], a, dir_rev); }
+		}
+		a->slack -= 2*prev[0]->blossom_eps;
+		return false;
+	}
+
+	for (dir=0; dir<2; dir++)
+	{
+		j = a->head[1-dir];
+		REMOVE_EDGE(j, a, dir);
+	}
+	a->next[0] = prev[0]->blossom_selfloops;
+	prev[0]->blossom_selfloops = a;
+	return false;
+}
+
+
+inline void PerfectMatching::AugmentBranch(Node* i0)
+{
+	int dir;
+	Tree* t = i0->tree;
+	Node* r = t->root;
+	Node* tree_root_prev = r->tree_sibling_prev;
+	Node* i;
+	Node* j;
+	Edge* a;
+	EdgeIterator I;
+	Arc* aa;
+	REAL eps = t->eps;
+	PriorityQueue<REAL>::Item* q;
+	TreeEdge* e;
+	TreeEdgeIterator T;
+	Tree* t2;
+
+	t = r->tree;
+	t->pq_current = t;
+
+	FOR_ALL_TREE_EDGES_X(t, e, dir, T)
+	{
+		t2 = e->head[dir];
+		e->head[1-dir] = NULL; // mark it for deletion
+
+		t2->pq_current = e;
+		t2->dir_current = dir;
+	}
+
+	i = r->first_tree_child;
+	if ( i )
+	while ( 1 )
+	{
+		Node* i0 = i;
+		i = ARC_HEAD(i->match);
+		if (i->is_processed)
+		{
+			if (i->is_blossom)
+			{
+				a = ARC_TO_EDGE_PTR(i->match);
+				REAL tmp = a->slack; a->slack = i->y; i->y = tmp;
+				PriorityQueue<REAL>::ResetItem(a);
+			}
+			FOR_ALL_EDGES(i, a, dir, I)
+			{
+				GET_OUTER_HEAD(a, dir, j);
+
+				if (j->flag == 0 && j->is_processed)
+				{
+					if (j->tree != t)
+					{
+						a->slack += eps;
+						if (PriorityQueue<REAL>::isReset(a)) j->tree->pq0.Add(a);
+					}
+				}
+				else a->slack += eps;
+			}
+		}
+
+		i = i0;
+		MOVE_NODE_IN_TREE(i);
+	}
+
+	///////////////////////////////////////////////////////////////////
+
+	FOR_ALL_TREE_EDGES(t, e, dir)
+	{
+		t2 = e->head[dir];
+		t2->pq_current = NULL;
+
+		e->pq01[1-dir].Merge(t2->pq0);
+		for (q=e->pq00.GetFirst(); q; q=e->pq00.GetNext(q))
+		{
+			q->slack -= eps;
+			int dir2;
+			for (dir2=0; dir2<2; dir2++) GET_OUTER_HEAD((Edge*)q, dir2, j);
+		}
+		e->pq00.Merge(t2->pq0);
+		for (q=e->pq01[dir].GetAndResetFirst(); q; q=e->pq01[dir].GetAndResetNext())
+		{
+			q->slack -= eps;
+			int dir2;
+			for (dir2=0; dir2<2; dir2++) GET_OUTER_HEAD((Edge*)q, dir2, j);
+		}
+	}
+	for (q=t->pq0.GetAndResetFirst(); q; q=t->pq0.GetAndResetNext())
+	{
+		q->slack -= eps;
+		int dir2;
+		for (dir2=0; dir2<2; dir2++) GET_OUTER_HEAD((Edge*)q, dir2, j);
+	}
+	for (q=t->pq00.GetAndResetFirst(); q; q=t->pq00.GetAndResetNext())
+	{
+		ProcessEdge00((Edge*)q);
+	}
+
+	///////////////////////////////////////////////////////////////////
+
+	r->flag = 2;
+	r->is_processed = 0;
+	i = r->first_tree_child;
+	r->y += eps;
+	if ( i )
+	while ( 1 )
+	{
+		j = ARC_HEAD(i->match);
+		j->flag = 2;
+		i->flag = 2;
+		j->is_processed = 0;
+		i->is_processed = 0;
+		j->y -= eps;
+		i->y += eps;
+
+		MOVE_NODE_IN_TREE(i);
+	}
+
+	///////////////////////////////////////////////////////////////////
+
+	i = i0;
+	if ( !i0->is_tree_root )
+	{
+		j = ARC_HEAD(i0->match);
+		GET_TREE_PARENT(j, i);
+		j->match = aa = j->tree_parent;
+		while ( !i->is_tree_root )
+		{
+			j = ARC_HEAD(i->match);
+			i->match = ARC_REV(aa);
+			GET_TREE_PARENT(j, i);
+			j->match = aa = j->tree_parent;
+		}
+		i->match = ARC_REV(aa);
+	}
+	r->is_tree_root = 0;
+	tree_root_prev->tree_sibling_next = r->tree_sibling_next;
+	if (r->tree_sibling_next) r->tree_sibling_next->tree_sibling_prev = tree_root_prev;
+	tree_num --;
+}
+
+
+void PerfectMatching::Augment(Edge* a)
+{
+	Node* j;
+	int dir;
+
+	for (dir=0; dir<2; dir++)
+	{
+		GET_OUTER_HEAD(a, dir, j);
+		AugmentBranch(j);
+		j->match = EDGE_DIR_TO_ARC(a, 1-dir);
+	}
+	if (options.verbose)
+	{
+		int k = 1;
+		while (k < tree_num) k *= 2;
+		if (k == tree_num || tree_num<=8 || (tree_num<=64 && (tree_num%8)==0)) { printf("%d.", tree_num); fflush(stdout); }
+	}
+}
+
+inline void PerfectMatching::GrowNode(Node* i)
+{
+	//assert(i->is_outer);
+	//assert(i->flag == 0);
+
+	Edge* a;
+	EdgeIterator I;
+	int dir;
+	Node* j;
+	Tree* t = i->tree;
+	REAL eps = t->eps;
+	Edge* a_augment = NULL;
+
+	FOR_ALL_EDGES(i, a, dir, I)
+	{
+		GET_OUTER_HEAD(a, dir, j);
+
+		if (j->flag == 2)
+		{
+			a->slack += eps; 
+			if (a->slack > 0) 
+			{ 
+				t->pq0.Add(a); 
+			}
+			else
+			{
+				j->flag = 1;
+				j->tree = i->tree;
+				j->tree_parent = EDGE_DIR_TO_ARC(a, 1-dir);
+				j->y += eps;
+				j = ARC_HEAD(j->match);
+				j->y -= eps;
+				ADD_TREE_CHILD(i, j);
+			}
+		}
+		else
+		{
+			if (j->flag == 0 && j->is_processed)
+			{
+				if (!PriorityQueue<REAL>::isReset(a)) j->tree->pq0.Remove(a, pq_buf);
+				if (a->slack <= j->tree->eps && j->tree != t) a_augment = a;
+				a->slack += eps;
+				if (!j->tree->pq_current) AddTreeEdge(t, j->tree);
+				j->tree->pq_current->pq00.Add(a);
+			}
+			else
+			{
+				a->slack += eps;
+				if (j->flag == 1 && j->tree != t)
+				{
+					if (!j->tree->pq_current) AddTreeEdge(t, j->tree);
+					j->tree->pq_current->pq01[j->tree->dir_current].Add(a);
+				}
+			}
+		}
+	}
+
+	//assert(!i->is_processed);
+	i->is_processed = 1;
+
+	if (!i->is_tree_root)
+	{
+		j = ARC_HEAD(i->match);
+		//assert(!j->is_processed);
+		j->is_processed = 1;
+		if (j->is_blossom)
+		{
+			a = ARC_TO_EDGE_PTR(i->match);
+			REAL tmp = a->slack; a->slack = j->y; j->y = tmp;
+			t->pq_blossoms.Add(a);
+		}
+	}
+
+	if (a_augment) Augment(a_augment);
+
+	stat.grow_count ++;
+}
+
+
+
+void PerfectMatching::GrowTree(Node* r, bool new_subtree)
+{
+	//assert(r->flag == 0);
+
+	Node* i = r;
+	Node* j;
+	Node* stop = r->tree_sibling_next;
+	if (new_subtree && r->first_tree_child) stop = r->first_tree_child;
+	Edge* a;
+	EdgeIterator I;
+	int dir;
+	Tree* t = r->tree;
+	REAL eps = t->eps;
+	int tree_num0 = tree_num;
+
+	while ( 1 )
+	{
+		if (!i->is_tree_root)
+		{
+			// process "-" node
+			i = ARC_HEAD(i->match);
+			FOR_ALL_EDGES(i, a, dir, I)
+			{
+				GET_OUTER_HEAD(a, dir, j);
+
+				if (j->flag == 2) a->slack -= eps;
+				else
+				{
+					if (j->flag == 0 && j->is_processed)
+					{
+						if (!PriorityQueue<REAL>::isReset(a)) j->tree->pq0.Remove(a, pq_buf);
+						a->slack -= eps;
+						if (j->tree != t)
+						{
+							if (!j->tree->pq_current) AddTreeEdge(t, j->tree);
+							j->tree->pq_current->pq01[1-j->tree->dir_current].Add(a);
+						}
+					}
+					else a->slack -= eps;
+				}
+			}
+			i = ARC_HEAD(i->match);
+		}
+		// process "+" node
+		GrowNode(i);
+		if (tree_num != tree_num0) break;
+
+		if (i->first_tree_child) i = i->first_tree_child;
+		else
+		{
+			while (i != r && !i->tree_sibling_next) { i = ARC_HEAD(i->match); GET_TREE_PARENT(i, i); }
+			i = i->tree_sibling_next;
+		}
+		if (i == stop) break;
+	}
+}
+
+void PerfectMatching::Solve(bool finish)
+{
+	Node* i;
+	Node* j;
+	Node* r;
+	Node* r2;
+	Node* r3 = NULL; // initialize to prevent compiler warning
+	PriorityQueue<REAL>::Item* q;
+	Edge* a;
+	Tree* t;
+	Tree* t2;
+	TreeEdge* e;
+	TreeEdgeIterator T;
+	int dir;
+	REAL eps;
+
+	double start_time = get_time();
+
+	if (IS_INT)
+	{
+		if (options.dual_greedy_update_option == 2)
+		{
+			printf("Fixed eps approach can only be used with floating point REAL!\n");
+			printf("Change REAL to double in PerfectMatching.h and recompile\n");
+			exit(1);
+		}
+		if (options.dual_LP_threshold > 0)
+		{
+			printf("LP approach can only be used with floating point REAL!\n");
+			printf("Change REAL to double in PerfectMatching.h and recompile\n");
+			exit(1);
+		}
+	}
+	if (options.verbose) { printf("perfect matching with %d nodes and %d edges\n", node_num, edge_num); fflush(stdout); }
+
+	if (first_solve)
+	{
+		if (options.verbose) { printf("    starting init..."); fflush(stdout); }
+		if (options.fractional_jumpstart) InitGlobal();
+		else                              InitGreedy();
+		if (options.verbose) printf("done [%.3f secs]. ", get_time() - start_time);
+		first_solve = false;
+	}
+	else if (options.verbose) printf("    solving updated problem. ");
+
+	if (options.verbose) { printf("%d trees\n    .", tree_num); fflush(stdout); }
+
+	memset(&stat, 0, sizeof(Stat));
+
+	///////////////////////////////////////////////////////
+	//       first pass - initialize auxiliary graph     //
+	///////////////////////////////////////////////////////
+
+	for (r=nodes[node_num].tree_sibling_next; r; r=r->tree_sibling_next)
+	{
+		//assert(!r->is_processed);
+		t = r->tree;
+		//assert(!t->first[0] && !t->first[1]);
+
+		EdgeIterator I;
+		FOR_ALL_EDGES(r, a, dir, I)
+		{
+			j = a->head[dir];
+			if (j->flag == 2) t->pq0.Add(a);
+			else if (j->is_processed)
+			{
+				//assert(j->flag == 0);
+				if (!j->tree->pq_current) AddTreeEdge(t, j->tree);
+				j->tree->pq_current->pq00.Add(a);
+			}
+		}
+		r->is_processed = 1;
+		FOR_ALL_TREE_EDGES(t, e, dir) e->head[dir]->pq_current = NULL;
+	}
+
+	///////////////////////////////////////////////////////
+	//                  main loop                        //
+	///////////////////////////////////////////////////////
+
+	while ( 1 )
+	{
+		int tree_num0 = tree_num;
+		Stat stat0 = stat;
+		REAL delta = 0;
+
+		for (r=nodes[node_num].tree_sibling_next; r; )
+		{
+			r2 = r->tree_sibling_next;
+			if (r2) r3 = r2->tree_sibling_next;
+			t = r->tree;
+
+			int tree_num1 = tree_num;
+
+			//////////////////////////////////////////////////////////////////////
+			// step 1 - traversing auxiliary graph, setting pq_current pointers //
+			//////////////////////////////////////////////////////////////////////
+			t->pq_current = t;
+			if (options.update_duals_before)
+			{
+				eps = PM_INFTY;
+				Edge* a_augment = NULL;
+				REAL eps_augment = PM_INFTY;
+				if ((q=t->pq0.GetMin())) eps = q->slack;
+				if ((q=t->pq_blossoms.GetMin()) && eps > q->slack) eps = q->slack;
+				while ((q=t->pq00.GetMin()))
+				{
+					if (ProcessEdge00((Edge*)q, false)) break;
+					t->pq00.Remove(q, pq_buf);
+				}
+				if (q && 2*eps > q->slack) eps = q->slack/2;
+				FOR_ALL_TREE_EDGES_X(t, e, dir, T)
+				{
+					t2 = e->head[dir];
+					t2->pq_current = e;
+					t2->dir_current = dir;
+					if ((q=e->pq00.GetMin()) && (!a_augment || eps_augment > q->slack-t2->eps)) { a_augment = (Edge*)q; eps_augment = q->slack-t2->eps; }
+					if ((q=e->pq01[dir].GetMin()) && eps > q->slack+t2->eps) eps = q->slack+t2->eps;
+				}
+				if (eps > eps_augment) eps = eps_augment;
+				if (eps > t->eps)
+				{
+					delta += eps - t->eps;
+					t->eps = eps;
+				}
+				if (a_augment && eps_augment <= t->eps) Augment(a_augment);
+			}
+			else
+			{
+				FOR_ALL_TREE_EDGES_X(t, e, dir, T)
+				{
+					t2 = e->head[dir];
+					t2->pq_current = e;
+					t2->dir_current = dir;
+
+					if ((q=e->pq00.GetMin()) && (q->slack - t->eps <= t2->eps))
+					{
+						Augment((Edge*)q);
+						break;
+					}
+				}
+			}
+
+			/////////////////////////////////
+			//   step 2 - growing tree     //
+			/////////////////////////////////
+			eps = t->eps;
+			REAL twice_eps = 2*eps;
+
+			while ( tree_num1 == tree_num )
+			{
+				if ((q=t->pq0.GetMin()) && q->slack <= t->eps)
+				{
+					a = (Edge*)q;
+					dir = (a->head[1]->flag == 2 && a->head[1]->is_outer) ? 1 : 0;
+					GET_OUTER_HEAD(a, 1-dir, i);
+					j = a->head[dir];
+					//assert(i->flag==0 && j->flag==2 && i->is_outer && j->is_outer && i->tree==t);
+
+					j->flag = 1;
+					j->tree = i->tree;
+					j->tree_parent = EDGE_DIR_TO_ARC(a, 1-dir);
+					j->y += eps;
+					j = ARC_HEAD(j->match);
+					j->y -= eps;
+					ADD_TREE_CHILD(i, j);
+
+					GrowTree(j, true);
+				}
+				else if ((q=t->pq00.GetMin()) && q->slack <= twice_eps)
+				{
+					t->pq00.Remove(q, pq_buf);
+					a = (Edge*)q;
+					if (ProcessEdge00(a)) Shrink(a);
+				}
+				else if ((q=t->pq_blossoms.GetMin()) && q->slack <= eps)
+				{
+					t->pq_blossoms.Remove(q, pq_buf);
+					a = (Edge*)q;
+					j = (a->head[0]->flag == 1) ? a->head[0] : a->head[1];
+					REAL tmp = a->slack; a->slack = j->y; j->y = tmp;
+					Expand(j);
+				}
+				else break;
+			}
+
+			///////////////////////////////////////////////////////////////////////
+			// step 3 - traversing auxiliary graph, clearing pq_current pointers //
+			///////////////////////////////////////////////////////////////////////
+			if ( tree_num1 == tree_num )
+			{
+				t->pq_current = NULL;
+				if (options.update_duals_after)
+				{
+					eps = PM_INFTY;
+					Edge* a_augment = NULL;
+					REAL eps_augment = PM_INFTY;
+					if ((q=t->pq0.GetMin())) eps = q->slack;
+					if ((q=t->pq_blossoms.GetMin()) && eps > q->slack) eps = q->slack;
+					while ((q=t->pq00.GetMin()))
+					{
+						if (ProcessEdge00((Edge*)q, false)) break;
+						t->pq00.Remove(q, pq_buf);
+					}
+					if (q && 2*eps > q->slack) eps = q->slack/2;
+					FOR_ALL_TREE_EDGES(t, e, dir)
+					{
+						t2 = e->head[dir];
+						e->head[dir]->pq_current = NULL;
+						if ((q=e->pq00.GetMin()) && (!a_augment || eps_augment > q->slack-t2->eps)) { a_augment = (Edge*)q; eps_augment = q->slack-t2->eps; }
+						if ((q=e->pq01[dir].GetMin()) && eps > q->slack+t2->eps) eps = q->slack+t2->eps;
+					}
+					if (eps > eps_augment) eps = eps_augment;
+					bool progress = false;
+					if (eps > t->eps)
+					{
+						delta += eps - t->eps;
+						t->eps = eps;
+						progress = true;
+					}
+					if (a_augment && eps_augment <= t->eps) Augment(a_augment);
+					else if (progress && tree_num >= options.single_tree_threshold*node_num)
+					{
+						// continue with the same tree
+						r = t->root;
+						continue;
+					}
+				}
+				else
+				{
+					FOR_ALL_TREE_EDGES(t, e, dir) e->head[dir]->pq_current = NULL;
+				}
+			}
+
+			///////////////////////////////////////////////////////////////////////
+			///////////////////////////////////////////////////////////////////////
+			///////////////////////////////////////////////////////////////////////
+
+			r = r2;
+			if (r && !r->is_tree_root) r = r3;
+		}
+
+		if (tree_num == 0) break;
+
+		if ( tree_num == tree_num0 )
+		  //&& stat.grow_count == stat0.grow_count 
+		  //&& stat.shrink_count == stat0.shrink_count 
+		  //&& stat.expand_count == stat0.expand_count )
+		{
+			if (!UpdateDuals())
+			{
+				if (!IS_INT && delta <= PM_THRESHOLD) // for numerical stability
+				{
+					//CommitEps();
+					int dual_greedy_update_option = options.dual_greedy_update_option;
+					options.dual_greedy_update_option = 2;
+					UpdateDuals();
+					options.dual_greedy_update_option = dual_greedy_update_option;
+				}
+			}
+		}
+	}
+
+	if (finish) Finish();
+
+	if (options.verbose)
+	{
+		printf("\ndone [%.3f secs]. %d grows, %d expands, %d shrinks\n", get_time()-start_time, stat.grow_count, stat.expand_count, stat.shrink_count); 
+		printf("    expands: [%.3f secs], shrinks: [%.3f secs], dual updates: [%.3f secs]\n", stat.expand_time, stat.shrink_time, stat.dual_time); 
+		fflush(stdout); 
+	}
+}
+
+
+
diff --git a/blossom5-v2.05.src/PMrepair.cpp b/blossom5-v2.05.src/PMrepair.cpp
new file mode 100644
index 0000000..978b38f
--- /dev/null
+++ b/blossom5-v2.05.src/PMrepair.cpp
@@ -0,0 +1,431 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "PMimplementation.h"
+
+
+struct PerfectMatching::LCATreeX : LCATree
+{
+	LCATreeX(int size) : LCATree(size) { rev_mapping = new Node*[size]; }
+	~LCATreeX() { delete [] rev_mapping; }
+	Node**	rev_mapping;
+};
+
+void PerfectMatching::StartUpdate()
+{
+	Node* i0;
+	Node* i;
+	Node* j;
+	Node* b;
+
+	while ((i=removed_first))
+	{
+		removed_first = i->tree_sibling_next;
+		blossoms->Delete(i);
+		removed_num --;
+	}
+
+	Edge* a;
+	Edge* selfloop_first = NULL;
+	Edge* selfloop_last = NULL;
+
+	for (i0=nodes; i0<nodes+node_num; i0++)
+	{
+		i0->is_processed = 0;
+		if (i0->is_outer) continue;
+
+		i0->is_tree_root = 0;
+		i0->blossom_ptr = NULL;
+		i = i0;
+		while ( 1 )
+		{
+			j = i->blossom_parent;
+			j->is_processed = 0;
+			if (j->is_outer) { j->first_tree_child = i; break; }
+			if (j->is_marked) break;
+			if ((a=j->blossom_selfloops))
+			{
+				if (selfloop_last) selfloop_last->next[1] = a;
+				else               selfloop_first         = a;
+				selfloop_last = a;
+				a->next[1] = NULL;
+			}
+			j->blossom_ptr = i;
+			i = j;
+		}
+		b = (i->blossom_parent->is_outer) ? i->blossom_parent : i->blossom_parent->blossom_grandparent;
+#ifdef LCA_REPAIRS
+		if (!b->is_marked)
+		{
+			b->lca_size = 1;
+			b->is_marked = 1;
+		}
+#endif
+		while ( 1 )
+		{
+#ifdef LCA_REPAIRS
+			b->lca_size ++;
+#endif
+			ARC_TO_EDGE_PTR(i->blossom_sibling)->y_saved = i->y;
+			i->y += i->blossom_parent->y;
+			if (!i->is_blossom) break;
+			i->is_marked = 1;
+			j = i;
+			i = i->blossom_ptr;
+			j->blossom_grandparent = b;
+		}
+		i->blossom_grandparent = b;
+	}
+
+#ifdef LCA_REPAIRS
+	for (i0=nodes; i0<nodes+node_num; i0++)
+	{
+		if (i0->is_outer) continue;
+		b = i0->blossom_grandparent;
+		if (!b->is_marked) continue;
+		b->is_marked = 0;
+		LCATreeX* lca = new LCATreeX(b->lca_size);
+		b->blossom_ptr = b->first_tree_child;
+		i = b;
+		while ( 1 )
+		{
+			if (i->blossom_ptr) i = i->blossom_ptr;
+			else
+			{
+				while ( 1 )
+				{
+					if (i->is_outer) break;
+					i->lca_preorder = lca->Add(i, i->blossom_parent);
+					lca->rev_mapping[i->lca_preorder] = i;
+					i = ARC_HEAD(i->blossom_sibling);
+					if (i != i->blossom_parent->blossom_ptr) break;
+					i = i->blossom_parent;
+				}
+				if (i->is_outer)
+				{
+					lca->AddRoot(i);
+					break;
+				}
+			}
+		}
+		b->lca = lca;
+	}
+#endif
+
+	while ((a=selfloop_first))
+	{
+		selfloop_first = a->next[1];
+		do
+		{
+			Edge* a_next = a->next[0];
+
+#ifdef LCA_REPAIRS
+			int _i = a->head0[1]->lca_preorder;
+			int _j = a->head0[0]->lca_preorder;
+			Node* b = a->head0[1]->blossom_grandparent;
+			b->lca->GetPenultimateNodes(_i, _j);
+			i = b->lca->rev_mapping[_i];
+			j = b->lca->rev_mapping[_j];
+#else
+			GetRealEndpoints(a, i, j);
+#endif
+			ADD_EDGE(i, a, 0);
+			ADD_EDGE(j, a, 1);
+			a->slack -= 2*i->blossom_eps;
+			a = a_next;
+		} while ( a );
+	}
+
+	/*
+	for (i0=nodes; i0<nodes+node_num; i0++)
+	{
+		if (i0->is_outer) continue;
+		b = i0->blossom_grandparent;
+		if (b->lca)
+		{
+			delete b->lca;
+			b->lca = NULL;
+		}
+	}
+	*/
+
+	nodes[node_num].first_tree_child = NULL;
+}
+
+void PerfectMatching::FinishUpdate()
+{
+	Node* i0;
+	Node* i;
+	Node* j;
+	Edge* a;
+	EdgeIterator I;
+	int dir;
+	Tree* t;
+
+	for (i0=nodes; i0<nodes+node_num; i0++)
+	{
+		if (i0->is_outer) continue;
+
+#ifdef LCA_REPAIRS
+		if (i0->blossom_grandparent->lca)
+		{
+			delete i0->blossom_grandparent->lca;
+			i0->blossom_grandparent->lca = NULL;
+		}
+#endif
+
+		//////////////////////////////////////////////////////////////
+		if (!i0->blossom_grandparent->is_removed)
+		{
+			i = i0;
+			do
+			{
+				i->y = ARC_TO_EDGE_PTR(i->blossom_sibling)->y_saved;
+				i->is_marked = 0;
+				i->blossom_selfloops = NULL;
+				i = i->blossom_parent;
+			} while (i->is_marked);
+			continue;
+		}
+		//////////////////////////////////////////////////////////////
+
+		i = i0->blossom_parent;
+		while ( 1 )
+		{
+			if (i->is_removed && !i->is_outer) break;
+			REAL y_parent = (i->is_outer) ? 0 : i->blossom_parent->y;
+			for (dir=0; dir<2; dir++)
+			{
+				if (!i->first[dir]) continue;
+				i->first[dir]->prev[dir]->next[dir] = NULL;
+				Edge* a_next;
+				for (a=i->first[dir]; a; a=a_next)
+				{
+					a_next = a->next[dir];
+					j = a->head0[1-dir];
+					ADD_EDGE(j, a, dir);
+					a->slack += j->blossom_parent->y - y_parent;
+				}
+				i->first[dir] = NULL;
+			}
+			if (i->is_removed) break;
+
+			j = i->blossom_parent;
+			i->is_removed = 1;
+			i->tree_sibling_next = removed_first;
+			removed_first = i;
+			i = j;
+		}
+		i0->y = ARC_TO_EDGE_PTR(i0->blossom_sibling)->y_saved;
+		i0->is_outer = 1;
+		i0->flag = 2;
+		i0->is_tree_root = 1;
+	}
+
+	Node* blossom_list = nodes[node_num].first_tree_child;
+
+
+
+	for (i=nodes; i<nodes+node_num; i++)
+	{
+		if (!i->is_tree_root) continue;
+		i->first_tree_child = nodes[node_num].first_tree_child;
+		nodes[node_num].first_tree_child = i;
+		REAL slack_min = PM_INFTY;
+		FOR_ALL_EDGES(i, a, dir, I)
+		{
+			if (slack_min > a->slack) slack_min = a->slack;
+		}
+		i->y += slack_min;
+		FOR_ALL_EDGES(i, a, dir, I) a->slack -= slack_min;
+	}
+
+	tree_num = 0;
+	for (i=nodes[node_num].first_tree_child; i!=blossom_list; i=i->first_tree_child)
+	{
+		tree_num ++;
+		if (!i->is_tree_root) continue;
+		FOR_ALL_EDGES(i, a, dir, I)
+		{
+			j = a->head[dir];
+			if (a->slack <= 0 && j->is_tree_root)
+			{
+				i->is_tree_root = j->is_tree_root = 0;
+				i->match = EDGE_DIR_TO_ARC(a, dir);
+				j->match = EDGE_DIR_TO_ARC(a, 1-dir);
+				tree_num -= 2;
+				break;
+			}
+		}
+	}
+	for ( ; i; i=i->first_tree_child)
+	{
+		if (i->is_removed) { i->is_tree_root = 0; continue; }
+		tree_num ++;
+	}
+
+	if (tree_num > tree_num_max)
+	{
+		if (trees) free(trees);
+		tree_num_max = tree_num;
+		trees = (Tree*) malloc(tree_num_max*sizeof(Tree));
+	}
+	t = trees;
+
+	Node* last_root = &nodes[node_num];
+	Node* i_next;
+	for (i=nodes; i; i=i_next)
+	{
+		if (!i->is_blossom) i_next = (i<nodes+node_num) ? (i + 1) : blossom_list;
+		else                i_next = i->first_tree_child;
+		if (!i->is_tree_root) continue;
+
+		i->flag = 0;
+		i->first_tree_child = NULL;
+		i->tree_sibling_prev = last_root;
+		last_root->tree_sibling_next = i;
+		last_root = i;
+		i->tree = t;
+		t->root = i;
+		t->eps = 0;
+		t->first[0] = t->first[1] = NULL;
+		t->pq_current = NULL;
+		t->pq00.Reset();
+		t->pq0.Reset();
+		t->pq_blossoms.Reset();
+		t ++;
+	}
+
+	assert(t == trees + tree_num);
+	last_root->tree_sibling_next = NULL;
+
+	while ((i=removed_first))
+	{
+		removed_first = i->tree_sibling_next;
+		blossoms->Delete(i);
+		blossom_num --;
+	}
+}
+
+PerfectMatching::REAL PerfectMatching::GetTwiceSum(NodeId i)
+{
+	assert(i>=0 && i<node_num);
+	return nodes[i].y;
+}
+
+inline void PerfectMatching::ProcessNegativeEdge(Edge* a)
+{
+	int dir;
+	Node* i;
+	for (dir=0; dir<2; dir++)
+	{
+		i = a->head0[dir];
+		if (i->is_outer)
+		{
+			if (!i->is_tree_root)
+			{
+				i->is_tree_root = 1;
+				i = ARC_HEAD(i->match);
+				assert(!i->is_tree_root && i->is_outer);
+				i->is_tree_root = 1;
+				if (i->is_blossom)
+				{
+					i->first_tree_child = nodes[node_num].first_tree_child;
+					nodes[node_num].first_tree_child = i;
+				}
+			}
+            return; 
+		}
+		if (i->blossom_grandparent->is_removed) return;
+	}
+
+	Node* b = i->blossom_grandparent;
+	assert(b->is_outer);
+
+	if (!b->is_tree_root)
+	{
+		b->is_tree_root = 1;
+		i = ARC_HEAD(b->match);
+		assert(!i->is_tree_root && i->is_outer);
+		i->is_tree_root = 1;
+		if (i->is_blossom)
+		{
+			i->first_tree_child = nodes[node_num].first_tree_child;
+			nodes[node_num].first_tree_child = i;
+		}
+	}
+
+	b->is_removed = 1;
+	b->tree_sibling_next = removed_first;
+	removed_first = b;
+}
+
+PerfectMatching::EdgeId PerfectMatching::AddNewEdge(NodeId _i, NodeId _j, REAL cost, bool do_not_add_if_positive_slack)
+{
+	assert(_i>=0 && _i<node_num && _j>=0 && _j<node_num && _i!=_j);
+	if (edge_num >= edge_num_max) ReallocateEdges();
+	Node* i = nodes + _i;
+	Node* j = nodes + _j;
+	Edge* a = edges + edge_num;
+
+	a->slack = cost*COST_FACTOR;
+	a->head0[0] = j;
+	a->head0[1] = i;
+	Node* bi = (i->is_outer) ? i : i->blossom_grandparent;
+	Node* bj = (j->is_outer) ? j : j->blossom_grandparent;
+	if (bi == bj)
+	{
+#ifdef LCA_REPAIRS
+		int _i = i->lca_preorder;
+		int _j = j->lca_preorder;
+		bi->lca->GetPenultimateNodes(_i, _j);
+		i = bi->lca->rev_mapping[_i];
+		j = bi->lca->rev_mapping[_j];
+#else
+		GetRealEndpoints(a, i, j);
+#endif
+		a->slack += i->blossom_parent->y + j->blossom_parent->y;
+	}
+	else
+	{
+		i = bi;
+		j = bj;
+	}
+	a->slack -= a->head0[0]->y + a->head0[1]->y;
+
+	if (a->slack >= 0 && do_not_add_if_positive_slack) return -1;
+
+	ADD_EDGE(i, a, 0);
+	ADD_EDGE(j, a, 1);
+	PriorityQueue<REAL>::ResetItem(a);
+
+	if (a->slack < 0)
+	{
+		ProcessNegativeEdge(a);
+	}
+
+	return edge_num ++;
+}
+
+void PerfectMatching::UpdateCost(EdgeId e, REAL delta_cost)
+{
+	assert(e>=0 && e<edge_num);
+	Edge* a = edges + e;
+	a->slack += delta_cost*COST_FACTOR;
+	if (a->slack == 0) return;
+	if (a->slack > 0)
+	{
+		Node* i = a->head[1];
+		Node* j = a->head[0];
+		if (i->is_outer)
+		{
+			if (ARC_TO_EDGE_PTR(i->match) != a && ARC_TO_EDGE_PTR(j->match) != a) return;
+		}
+		else
+		{
+			if (ARC_TO_EDGE_PTR(i->blossom_sibling) != a && ARC_TO_EDGE_PTR(j->blossom_sibling) != a) return;
+		}
+	}
+	ProcessNegativeEdge(a);
+}
+
diff --git a/blossom5-v2.05.src/PMshrink.cpp b/blossom5-v2.05.src/PMshrink.cpp
new file mode 100644
index 0000000..eabad1a
--- /dev/null
+++ b/blossom5-v2.05.src/PMshrink.cpp
@@ -0,0 +1,318 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "PMimplementation.h"
+
+
+PerfectMatching::Node* PerfectMatching::FindBlossomRoot(Edge* a0)
+{
+	Node* i;
+	Node* j;
+	Node* _i[2];
+	Node* r;
+	int branch;
+
+	_i[0] = ARC_HEAD(a0);
+	_i[1] = ARC_TAIL(a0);
+	branch = 0;
+	while ( 1 )
+	{
+		if (_i[branch]->is_marked)
+		{
+			r = _i[branch]; 
+			j = _i[1-branch];
+			break; 
+		}
+		_i[branch]->is_marked = 1;
+		if (_i[branch]->is_tree_root)
+		{
+			j = _i[branch];
+			i = _i[1-branch];
+			while (!i->is_marked)
+			{
+				i->is_marked = 1;
+				i = ARC_HEAD(i->match);
+				GET_TREE_PARENT(i, i);
+			}
+			r = i;
+			break;
+		}
+		i = ARC_HEAD(_i[branch]->match);
+		GET_TREE_PARENT(i, _i[branch]);
+		branch = 1 - branch;
+	}
+	i = r;
+	while ( i != j )
+	{
+		i = ARC_HEAD(i->match);
+		i = ARC_HEAD(i->tree_parent);
+		i->is_marked = 0;
+	}
+	// clear is_marked and is_outer
+	i = ARC_HEAD(a0);
+	while (i != r)
+	{
+		i->is_marked = 0;
+		i->is_outer = 0;
+		i = ARC_HEAD(i->match);
+		i->is_outer = 0;
+		i = ARC_HEAD(i->tree_parent);
+	}
+	i = ARC_TAIL(a0);
+	while (i != r)
+	{
+		i->is_marked = 0;
+		i->is_outer = 0;
+		i = ARC_HEAD(i->match);
+		i->is_outer = 0;
+		i = ARC_HEAD(i->tree_parent);
+	}
+	r->is_marked = 0;
+	r->is_outer = 0;
+
+	return r;
+}
+
+
+void PerfectMatching::Shrink(Edge* a0)
+{
+	//assert(a0->head[0]->is_outer && a0->head[1]->is_outer);
+	//assert(a0->head[0]->flag == 0 && a0->head[1]->flag == 0);
+
+	double start_time = get_time();
+
+	int branch, dir;
+	Node* r;
+	Node* i;
+	Node* j;
+	Edge* a;
+	Edge** a_inner_ptr;
+	Arc* a_prev;
+	Node* b = blossoms->New();
+	Edge* a_augment = NULL;
+	Edge* b_match;
+
+	b->first[0] = b->first[1] = NULL;
+
+	// set is_outer=0 for all nodes in the blossom
+	r = FindBlossomRoot(a0);
+	Tree* t = r->tree;
+	REAL eps = t->eps;
+
+	b->first_tree_child = NULL;
+	i = ARC_HEAD(a0);
+	branch = 0;
+	while ( 1 )
+	{
+		if (i == r && branch) break;
+		i->is_marked = 1;
+		if (i == r)
+		{
+			branch = 1;
+			i = ARC_TAIL(a0);
+			continue;
+		}
+
+		// remove i from the list of children
+		REMOVE_FROM_TREE(i);
+
+		// move children of i to the list of children of b
+		if (i->first_tree_child)
+		{
+			j = i->first_tree_child;
+			if (!b->first_tree_child) b->first_tree_child = j;
+			else
+			{
+				Node* j_last = j->tree_sibling_prev;
+				j->tree_sibling_prev = b->first_tree_child->tree_sibling_prev;
+				b->first_tree_child->tree_sibling_prev->tree_sibling_next = j;
+				b->first_tree_child->tree_sibling_prev = j_last;
+			}
+		}
+
+		// go to parent
+		i = ARC_HEAD(i->match);
+		i->is_marked = 1;
+		if (i->is_blossom)
+		{
+			a = ARC_TO_EDGE_PTR(i->match);
+			t->pq_blossoms.Remove(a, pq_buf);
+			REAL tmp = a->slack; a->slack = i->y; i->y = tmp;
+		}
+		i = ARC_HEAD(i->tree_parent); 
+	}
+
+	// move children of r to the list of children of b
+	if (i->first_tree_child)
+	{
+		j = i->first_tree_child;
+		if (!b->first_tree_child) b->first_tree_child = j;
+		else
+		{
+			Node* j_last = j->tree_sibling_prev;
+			j->tree_sibling_prev = b->first_tree_child->tree_sibling_prev;
+			b->first_tree_child->tree_sibling_prev->tree_sibling_next = j;
+			b->first_tree_child->tree_sibling_prev = j_last;
+		}
+	}
+
+	// init b
+	b->is_removed = 0;
+	b->is_outer = 1;
+	b->flag = 0;
+	b->is_blossom = 1;
+	b->is_tree_root = r->is_tree_root;
+	b->is_processed = 1;
+	b->tree = t;
+	b->y = -eps;
+	b->is_marked = 0;
+
+	// replace r with b in the tree
+	b->tree_sibling_prev = r->tree_sibling_prev;
+	b->tree_sibling_next = r->tree_sibling_next;
+	Node* b_parent = NULL;
+	if (!b->is_tree_root)
+	{
+		b_parent = ARC_HEAD(r->match); GET_TREE_PARENT(b_parent, b_parent); 
+	}
+	if (b->tree_sibling_prev->tree_sibling_next) b->tree_sibling_prev->tree_sibling_next = b;
+	else b_parent->first_tree_child = b;
+	if (b->tree_sibling_next) b->tree_sibling_next->tree_sibling_prev = b;
+	else if (b_parent) b_parent->first_tree_child->tree_sibling_prev = b;
+
+	if (b->is_tree_root)
+	{
+		b->tree->root = b;
+		b_match = NULL;
+	}
+	else
+	{
+		b->match = r->match;
+		b_match = ARC_TO_EDGE_PTR(b->match);
+	}
+	REAL b_match_slack = 0; // initialize to prevent compiler warning
+	if (b_match && ARC_HEAD(b->match)->is_blossom)
+	{
+		b_match_slack = b_match->slack;
+		b_match->slack = ARC_HEAD(b->match)->y;
+	}
+
+	// second pass over nodes in the blossom
+	branch = 0;
+	a_prev = EDGE_DIR_TO_ARC(a0, 0);
+	i = ARC_HEAD(a_prev);
+	while ( 1 )
+	{
+		// update Arc::next and Arc::head pointers
+		if (i->flag == 0) i->y += eps;
+		else              i->y -= eps;
+		i->is_processed = 0;
+
+		if (i->flag == 1)
+		{
+			Edge* a_prev;
+			for (dir=0; dir<2; dir++)
+			if (i->first[dir])
+			{
+				for (a_inner_ptr=&i->first[dir], a=*a_inner_ptr, a_prev=a->prev[dir], a_prev->next[dir]=NULL; a; a=*a_inner_ptr)
+				{
+					Node* j0 = a->head[dir];
+					for (j=j0; !j->is_outer && !j->is_marked; j = j->blossom_parent) {}
+					if (j != j0) { /*assert(j->flag == 0);*/ int dir_rev = 1 - dir; MOVE_EDGE(j0, j, a, dir_rev); }
+					if (j->is_marked) // "inner" arc
+					{
+						a_inner_ptr = &a->next[dir];
+						a->prev[dir] = a_prev;
+						a_prev = a;
+
+						if (j->flag == 1) a->slack += eps;
+					}
+					else // "boundary" arc
+					{
+						*a_inner_ptr = a->next[dir];
+						ADD_EDGE(b, a, dir);
+
+						if (j->flag == 0 && j->tree != t) 
+						{
+							j->tree->pq_current->pq01[1-j->tree->dir_current].Remove(a, pq_buf);
+							if (a->slack + eps <= j->tree->eps) a_augment = a;
+						}
+						a->slack += 2*eps;
+						if (j->flag == 2) t->pq0.Add(a);
+						else if (j->flag == 0)
+						{
+							if (!j->tree->pq_current) AddTreeEdge(t, j->tree);
+							j->tree->pq_current->pq00.Add(a);
+						}
+						else if (j->tree != t)
+						{
+							if (!j->tree->pq_current) AddTreeEdge(t, j->tree);
+							j->tree->pq_current->pq01[j->tree->dir_current].Add(a);
+						}
+					}
+				}
+				if (i->first[dir])
+				{
+					a_prev->next[dir] = i->first[dir];
+					i->first[dir]->prev[dir] = a_prev;
+				}
+			}
+		}
+
+		Arc* a_next = (i->flag == 0) ? i->match : i->tree_parent;
+		i->blossom_parent = b;
+		i->match = NULL;
+		i->blossom_grandparent = b;
+		i->blossom_selfloops = NULL;
+		if (branch == 0)
+		{
+			i->blossom_sibling = a_next;
+			if (i == r)
+			{
+				branch = 1;
+				a_prev = ARC_REV(a0);
+				i = ARC_HEAD(a_prev);
+				if (i == r) break;
+			}
+			else
+			{
+				a_prev = i->blossom_sibling;
+				i = ARC_HEAD(a_prev);
+			}
+		}
+		else
+		{
+			i->blossom_sibling = ARC_REV(a_prev);
+			a_prev = a_next;
+			i = ARC_HEAD(a_prev);
+			if (i == r) break;
+		}
+	}
+	i->blossom_sibling = ARC_REV(a_prev);
+	r->is_tree_root = 0;
+
+	for (i=ARC_HEAD(r->blossom_sibling); ; i = ARC_HEAD(i->blossom_sibling))
+	{
+		i->is_marked = 0;
+		i->blossom_eps = eps;
+		if (i == r) break;
+	}
+
+	if (b_match)
+	{
+		if (ARC_HEAD(b->match)->is_blossom)
+		{
+			b_match->slack = b_match_slack;
+		}
+		dir = ARC_TO_EDGE_DIR(b->match);
+		//assert(b_match->head[1-dir] == r);
+		MOVE_EDGE(r, b, b_match, dir);
+	}
+
+	stat.shrink_count ++;
+	blossom_num ++;
+	stat.shrink_time += get_time() - start_time;
+
+	if (a_augment) Augment(a_augment);
+}
+
diff --git a/blossom5-v2.05.src/PQ.h b/blossom5-v2.05.src/PQ.h
new file mode 100644
index 0000000..8fb37b7
--- /dev/null
+++ b/blossom5-v2.05.src/PQ.h
@@ -0,0 +1,396 @@
+/*
+    PQ.h - implements pairing heaps priority queue (with multipass 'delete min')
+
+    Copyright 2008 Vladimir Kolmogorov (vnk@ist.ac.at)
+
+    This software can be used for research purposes only. Commercial use is prohibited.
+    Public redistribution of the code or its derivatives is prohibited.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef HFKSJHFKJHARBABDAKFAF
+#define HFKSJHFKJHARBABDAKFAF
+
+// exactly one flag must be defined
+//#define PQ_MULTIPASS
+#define PQ_INTERLEAVED_MULTIPASS
+
+#include <string.h>
+
+template <typename REAL> class PriorityQueue
+{
+public:
+	struct Item
+	{
+		REAL	slack;
+
+		Item*	parentPQ;
+		union
+		{
+			struct
+			{
+				Item*	leftPQ;
+				Item*	rightPQ;
+			};
+			REAL	y_saved; // used in repairs
+		};
+	};
+	static void* AllocateBuf();
+	static void DeallocateBuf(void* buf);
+
+	static void ResetItem(Item* i);
+	static bool isReset(Item* i);
+
+	//////////////////////////////////////////////////////////
+
+	void Reset();
+	void Add(Item* i);
+#define Remove(i, buf) _Remove(i)
+	void _Remove(Item* i);
+	void Decrease(Item* i_old, Item* i_new, void* buf);
+	Item* GetMin();
+
+	//////////////////////////////////////////////////////////
+
+	void Update(REAL delta);
+	void Merge(PriorityQueue<REAL>& dest);
+
+	// traversing items in the order they are stored (irrespective of slack).
+	// The caller must go through all items, no other member functions can be called during the scan.
+	Item* GetAndResetFirst();
+	Item* GetAndResetNext();
+
+	Item* GetFirst();
+	Item* GetNext(Item* i);
+
+	//////////////////////////////////////////////////////////
+
+private:
+	struct Buf
+	{
+	};
+	Item*	rootPQ;
+	void RemoveRoot();
+};
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+template <typename REAL> inline void* PriorityQueue<REAL>::AllocateBuf()
+{
+	return NULL;
+}
+
+template <typename REAL> inline void PriorityQueue<REAL>::DeallocateBuf(void* _buf)
+{
+}
+
+template <typename REAL> inline void PriorityQueue<REAL>::ResetItem(Item* i) 
+{ 
+	i->parentPQ = NULL;
+}
+
+template <typename REAL> inline bool PriorityQueue<REAL>::isReset(Item* i) 
+{ 
+	return (i->parentPQ == NULL);
+}
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+template <typename REAL> inline void PriorityQueue<REAL>::Reset() 
+{ 
+	rootPQ = NULL; 
+}
+
+/*
+template <typename REAL> inline void PriorityQueue<REAL>::RemoveRoot()
+{
+	Item* r = rootPQ;
+	PriorityQueue<REAL> pq;
+	pq.rootPQ = rootPQ;
+	rootPQ = NULL;
+	Item* i;
+	for (i=pq.GetAndResetFirst(); i; i=pq.GetAndResetNext())
+	{
+		if (i != r) Add(i);
+	}
+	r->parentPQ = NULL;
+}
+*/
+
+// sets i = merge(i, j). Ignores parentPQ and rightPQ for i and j.
+#define MERGE_PQ(i, j)\
+	{\
+		if (i->slack <= j->slack)\
+		{\
+			j->rightPQ = i->leftPQ;\
+			if (j->rightPQ) j->rightPQ->parentPQ = j;\
+			j->parentPQ = i;\
+			i->leftPQ = j;\
+		}\
+		else\
+		{\
+			i->rightPQ = j->leftPQ;\
+			if (i->rightPQ) i->rightPQ->parentPQ = i;\
+			i->parentPQ = j;\
+			j->leftPQ = i;\
+			i = j;\
+		}\
+	}
+
+template <typename REAL> inline void PriorityQueue<REAL>::RemoveRoot()
+{
+	Item* i = rootPQ->leftPQ;
+	rootPQ->parentPQ = NULL;
+	if (i)
+	{
+#ifdef PQ_MULTIPASS
+		while ( i->rightPQ )
+		{
+			Item** prev_ptr = &rootPQ;
+			while ( 1 )
+			{
+				if (i->rightPQ)
+				{
+					Item* j = i->rightPQ;
+					Item* next = j->rightPQ;
+					MERGE_PQ(i, j);
+					*prev_ptr = i;
+					if (!next) { i->rightPQ = NULL; break; }
+					prev_ptr = &i->rightPQ;
+					i = next;
+				}
+				else
+				{
+					*prev_ptr = i;
+					i->rightPQ = NULL;
+					break;
+				}
+			}
+			i = rootPQ;
+		}
+#endif
+
+#ifdef PQ_INTERLEAVED_MULTIPASS
+		while ( i->rightPQ )
+		{
+			Item* prev = NULL;
+			while ( i )
+			{
+				Item* next;
+				if (i->rightPQ)
+				{
+					Item* j = i->rightPQ;
+					next = j->rightPQ;
+					MERGE_PQ(i, j);
+				}
+				else next = NULL;
+				i->rightPQ = prev;
+				prev = i;
+				i = next;
+			}
+			i = prev;
+		}
+#endif
+		i->parentPQ = i;
+	}
+	rootPQ = i;
+}
+
+template <typename REAL> inline void PriorityQueue<REAL>::Add(Item* i)
+{
+	if (!rootPQ)
+	{
+		rootPQ = i;
+		i->parentPQ = i;
+		i->leftPQ = i->rightPQ = NULL;
+	}
+	else if (i->slack <= rootPQ->slack)
+	{
+		rootPQ->parentPQ = i;
+		i->leftPQ = rootPQ;
+		i->rightPQ = NULL;
+		rootPQ = i;
+		i->parentPQ = i;
+	}
+	else
+	{
+		i->leftPQ = NULL;
+		i->rightPQ = rootPQ->leftPQ;
+		if (i->rightPQ) i->rightPQ->parentPQ = i;
+		rootPQ->leftPQ = i;
+		i->parentPQ = rootPQ;
+	}
+}
+
+
+template <typename REAL> inline void PriorityQueue<REAL>::_Remove(Item* i)
+{
+	Item* p = i->parentPQ;
+	if (p == i) RemoveRoot();
+	else
+	{
+		if (i->rightPQ) i->rightPQ->parentPQ = p;
+		if (p->leftPQ == i) p->leftPQ  = i->rightPQ;
+		else                p->rightPQ = i->rightPQ;
+		if (i->leftPQ)
+		{
+			i->parentPQ = i;
+			i->rightPQ = NULL;
+			PriorityQueue<REAL> pq;
+			pq.rootPQ = i;
+			pq.RemoveRoot();
+			pq.Merge(*this);
+		}
+		else i->parentPQ = NULL;
+	}
+}
+
+template <typename REAL> inline void PriorityQueue<REAL>::Decrease(Item* i_old, Item* i_new, void* _buf)
+{
+	if (i_old->parentPQ == i_old)
+	{
+		if (i_old != i_new)
+		{
+			rootPQ = i_new;
+			i_new->parentPQ = i_new;
+			i_new->leftPQ = i_old->leftPQ;
+			i_new->rightPQ = NULL;
+			if (i_new->leftPQ) i_new->leftPQ->parentPQ = i_new;
+			i_old->parentPQ = NULL;
+		}
+	}
+	else
+	{
+		Remove(i_old, _buf);
+		Add(i_new);
+	}
+}
+
+template <typename REAL> inline typename PriorityQueue<REAL>::Item* PriorityQueue<REAL>::GetMin()
+{
+	return rootPQ;
+}
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+
+
+template <typename REAL> inline void PriorityQueue<REAL>::Merge(PriorityQueue<REAL>& dest)
+{
+	if (!rootPQ) return;
+	if (!dest.rootPQ) dest.rootPQ = rootPQ;
+	else
+	{
+		if (rootPQ->slack < dest.rootPQ->slack)
+		{
+			Item* j = rootPQ; rootPQ = dest.rootPQ; dest.rootPQ = j;
+		}
+		rootPQ->rightPQ = dest.rootPQ->leftPQ;
+		if (rootPQ->rightPQ) rootPQ->rightPQ->parentPQ = rootPQ;
+		rootPQ->parentPQ = dest.rootPQ;
+		dest.rootPQ->leftPQ = rootPQ;
+	}
+	rootPQ = NULL;
+}
+
+
+
+template <typename REAL> inline void PriorityQueue<REAL>::Update(REAL delta)
+{
+	if (!rootPQ) return;
+
+	Item* i = rootPQ;
+	while (i->leftPQ) i = i->leftPQ;
+
+	while ( 1 )
+	{
+		i->slack += delta;
+
+		if (i->rightPQ)
+		{
+			i = i->rightPQ;
+			while (i->leftPQ) i = i->leftPQ;
+		}
+		else
+		{
+			while ( 1 )
+			{
+				Item* j = i;
+				i = i->parentPQ;
+				if (i == j) return;
+				if (i->leftPQ == j) break;
+			}
+		}
+	}
+}
+
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+template <typename REAL> inline typename PriorityQueue<REAL>::Item* PriorityQueue<REAL>::GetAndResetFirst()
+{
+	if (!rootPQ) return NULL;
+	return GetAndResetNext();
+}
+
+template <typename REAL> inline typename PriorityQueue<REAL>::Item* PriorityQueue<REAL>::GetAndResetNext()
+{
+	if (!rootPQ) return NULL;
+	Item* result = rootPQ;
+	result->parentPQ = NULL;
+	Item* i = rootPQ->leftPQ;
+	if (!i) rootPQ = result->rightPQ;
+	else
+	{
+		rootPQ = i;
+		while (i->rightPQ) i = i->rightPQ;
+		i->rightPQ = result->rightPQ;
+	}
+	return result;
+}
+
+template <typename REAL> inline typename PriorityQueue<REAL>::Item* PriorityQueue<REAL>::GetFirst()
+{
+	if (!rootPQ) return NULL;
+	Item* i = rootPQ;
+	while (i->leftPQ) i = i->leftPQ;
+	return i;
+}
+
+template <typename REAL> inline typename PriorityQueue<REAL>::Item* PriorityQueue<REAL>::GetNext(Item* i)
+{
+	if (i->rightPQ)
+	{
+		i = i->rightPQ;
+		while (i->leftPQ) i = i->leftPQ;
+		return i;
+	}
+	while ( 1 )
+	{
+		Item* j = i;
+		i = i->parentPQ;
+		if (i == j) return NULL;
+		if (i->leftPQ == j) return i;
+	}
+}
+
+#endif
diff --git a/blossom5-v2.05.src/PerfectMatching.h b/blossom5-v2.05.src/PerfectMatching.h
new file mode 100644
index 0000000..7b1ae57
--- /dev/null
+++ b/blossom5-v2.05.src/PerfectMatching.h
@@ -0,0 +1,246 @@
+/*
+    PerfectMatching.h - interface to min cost perfect matching code
+
+
+
+    Copyright 2008-2009 UCL Business PLC, Author Vladimir Kolmogorov (vnk@ist.ac.at)
+
+    This software can be used for evaluation and non-commercial research purposes only. Commercial use is prohibited.
+    Public redistribution of the code or its derivatives is prohibited.
+    If you use this software for research purposes, you should cite the following paper in any resulting publication:
+
+        Vladimir Kolmogorov. "Blossom V: A new implementation of a minimum cost perfect matching algorithm."
+        In Mathematical Programming Computation (MPC), July 2009, 1(1):43-67.
+
+    For commercial use of the software not covered by this agreement, you may obtain a licence from
+    the copyright holders UCL Business via their licensing site: www.e-lucid.com/i/software/optimisation_software/BlossomV.html.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef HALSKDJDFHALSJASFDFASJGLA
+#define HALSKDJDFHALSJASFDFASJGLA
+
+#include <assert.h>
+#include <limits.h>
+#include "block.h"
+
+
+// if defined, edge costs are of type 'double', otherwise 'int'
+//#define PERFECT_MATCHING_DOUBLE
+
+// Note: with floating point numbers polynomial complexity is not guaranteed;
+// the code may even get stuck due to rounding errors. If the code terminates,
+// the solution may not be optimal. It may be worth calling CheckPerfectMatchingOptimality()
+// to see whether complementary slackness conditions are satisfied.
+//
+// Using single precision floating point numbers (float) is really not recommended.
+
+
+class PerfectMatching
+{
+public:
+
+#ifdef PERFECT_MATCHING_DOUBLE
+	typedef double REAL; 
+	#define PM_INFTY ((REAL)1e100)
+#else
+	typedef int REAL;
+	#define PM_INFTY (INT_MAX/2)
+#endif
+
+	typedef int NodeId;
+	typedef int EdgeId;
+
+	PerfectMatching(int nodeNum, int edgeNumMax);
+	~PerfectMatching();
+
+	// first call returns 0, second 1, and so on. 
+	EdgeId AddEdge(NodeId i, NodeId j, REAL cost);
+
+	// Computes a perfect matching of minimum cost. 
+	// NOTE: a perfect matching of finite cost must exist (otherwise the behaviour is not specified).
+	// If finish is false, then the final matching is not computed (so GetSolution() cannot be called afterwards).
+	void Solve(bool finish=true);
+
+	///////////////////////////////////////////////////////////////
+	// Read primal solution (can be called after Solve()).
+	int GetSolution(EdgeId e); // returns 1 if e is in the matching, 0 otherwise
+	NodeId GetMatch(NodeId i); // alternative way to get the result
+
+	///////////////////////////////////////////////////////////////
+	// Read dual solution (can be called after Solve()).
+	// 'blossom_parents' and 'twice_y' must be arrays of size node_num+GetBlossomNum().
+	// The function sets blossom_parent[i] to the parent of i (or to -1 for exterior nodes).
+	void GetDualSolution(int* blossom_parents, REAL* twice_y);
+	int GetBlossomNum();
+
+	///////////////////////////////////////////////////////////////
+	// Dynamic graph updates. After calling Solve() you may call //
+	// StartUpdate(), ..., FinishUpdate() and then Solve() again //
+	///////////////////////////////////////////////////////////////
+	void StartUpdate();
+	void FinishUpdate();
+
+	// 3 functions below can be called only between StartUpdate() and FinishUpdate().
+	REAL GetTwiceSum(NodeId i); // if 2*cost(i,j)>=GetTwiceSum(i)+GetTwiceSum(j) then adding new edge (i,j) is not necessary - optimal solution will not change
+	EdgeId AddNewEdge(NodeId i, NodeId j, REAL cost, bool do_not_add_if_positive_slack=true); // if do_not_add_if_positive_slack is true and the slack of the edge turns out to be non-negative, then the edge will not be added and -1 will be returned
+	void UpdateCost(EdgeId e, REAL delta_cost);
+
+
+	// NOTE: with default options the dual vector is guaranteed to be half-integral
+	// (if all input weights are integral). However, with some options there is
+	// no such guarantee, in particular, if dual_greedy_update_option=2 or dual_LP_threshold>0.
+	// These options can be used only if the code is compiled with REAL=double.
+	struct Options
+	{
+		Options() : fractional_jumpstart(true),
+		            dual_greedy_update_option(0),
+		            dual_LP_threshold(0.00),
+		            update_duals_before(false),
+		            update_duals_after(false),
+		            single_tree_threshold(1.00),
+		            verbose(true)
+		{}
+
+		bool	fractional_jumpstart; // false: greedy, true: compute fractional matching
+
+		int 	dual_greedy_update_option; // 0: compute connected components (as in Blossom IV)
+				                           // 1: compute strongly connected components (discussed by Cook-Rohe, but not implemented)
+				                           // 2: single eps for all trees (fixed eps approach)
+
+		double	dual_LP_threshold; // if tree_num => dual_updates_threshold*node_num: greedy updates
+		                           // if tree_num <  dual_updates_threshold*node_num: global updates (solve LP)
+
+		bool	update_duals_before; // before tree growth
+		bool	update_duals_after;  // after tree growth
+
+		double	single_tree_threshold; // if (tree_num => single_tree_threshold*node_num && update_duals_after): try to grow a single tree as long as possible
+
+		bool	verbose;
+	} options;
+
+
+	// save problem to a file. format=0 corresponds to DIMACS format,
+	// format=1 corresponds to the format used by blossom4.
+	// CANNOT BE CALLED AFTER Solve()!!!
+	void Save(char* filename, int format=0); 
+
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+private:
+	struct Node;
+	struct Arc; // no such struct, only Arc* is used (the pointer can be odd or even)
+	struct Edge; // pointer Edge* is always even
+	struct Tree;
+	struct TreeEdge;
+	struct PQPointers;
+	struct EdgeIterator;
+	struct TreeEdgeIterator;
+	struct LCATreeX;
+
+	Node*	nodes;
+	Edge*	edges;
+	char*	edges_orig;
+	DBlock<Node>* blossoms;
+	Tree*	trees;
+	DBlock<TreeEdge>* tree_edges;
+	struct ExpandTmpItem
+	{
+		Node*	i;
+		Node*	blossom_parent;
+		Node*	blossom_grandparent;
+	};
+	Block<ExpandTmpItem>* expand_tmp_list; // used inside Expand()
+
+	int		node_num;
+	int		edge_num, edge_num_max;
+	int		tree_num, tree_num_max;
+
+	Node*	removed_first;
+	int		blossom_num;
+	int		removed_num;
+
+	void*	pq_buf;
+
+	bool	first_solve;
+
+	// stat
+	struct Stat
+	{
+		int		shrink_count;
+		int		expand_count;
+		int		grow_count;
+		double	shrink_time;
+		double	expand_time;
+		double	dual_time;
+	} stat;
+
+	////////////////////////////////////////////////////////////////////
+
+	void InitGreedy(bool allocate_trees=true);
+
+	void InitGlobal(); // compute fractional matching
+	Node* FindBlossomRootInit(Edge* a0);
+	void ShrinkInit(Edge* a0, Node* tree_root);
+	void ExpandInit(Node* b);
+	void AugmentBranchInit(Node* i0, Node* tree_root);
+
+	void Finish(); // sets matching for inner nodes
+
+	void ProcessNegativeEdge(Edge* a);
+
+	void GetRealEndpoints(Edge* a, Node*& tail, Node*& head);
+	Node* FindBlossomRoot(Edge* a0);
+	void Shrink(Edge* a0);
+	void Expand(Node* b);
+	void Augment(Edge* a0);
+	void AugmentBranch(Node* i0);
+	void GrowNode(Node* i);
+	void GrowTree(Node* root, bool new_subtree);
+	bool ProcessEdge00(Edge* a, bool update_boundary_edge=true); // returns true if boundary edge, false otherwise
+	void ProcessSelfloop(Node* b, Edge* a);
+
+	void AddTreeEdge(Tree* t0, Tree* t1);
+
+	void ComputeEpsSingle(); // called from UpdateDuals()
+	void ComputeEpsCC(); // called from UpdateDuals()
+	void ComputeEpsSCC(); // called from UpdateDuals()
+	void ComputeEpsGlobal(); // called from UpdateDuals()
+	bool UpdateDuals();
+
+	void FreeRemoved();
+	void CommitEps();
+
+	void ReallocateEdges();
+
+	void PrintAll();
+};
+
+
+// in the functions below, 'edges' is an array of size 2*edge_num (edge e = (edges[2*e],edges[2*e+1])), 
+//                         'weights' is an array of size edge_num
+
+// checks complementary slackness conditions. 
+// returns 0 if success.
+// returns 1 if complementary slackness conditions are violated (then the amount of violation is printed - could potentially happen for double's)
+// returns 2 if the blossom tree structure is incorrect (or inconsistent with primal solution)
+int CheckPerfectMatchingOptimality(int node_num, int edge_num, int* edges, int* weights, PerfectMatching* pm, PerfectMatching::REAL threshold=(PerfectMatching::REAL)(1e-10));
+
+
+double ComputePerfectMatchingCost(int node_num, int edge_num, int* edges, int* weights, PerfectMatching* pm);
+
+#endif
diff --git a/blossom5-v2.05.src/README.TXT b/blossom5-v2.05.src/README.TXT
new file mode 100644
index 0000000..2e65742
--- /dev/null
+++ b/blossom5-v2.05.src/README.TXT
@@ -0,0 +1,80 @@
+BLOSSOM V - implementation of Edmonds' algorithm for computing a minimum cost perfect matching in a graph
+Version 2.05
+http://pub.ist.ac.at/~vnk/software.html
+
+Details of the implementation are described in
+
+	Vladimir Kolmogorov. "Blossom V: A new implementation of a minimum cost perfect matching algorithm."
+	In Mathematical Programming Computation (MPC), July 2009, 1(1):43-67.
+
+Please send comments to vnk@ist.ac.at.
+If you use this software for research purposes, you should cite the aforementioned paper in any resulting publication.
+
+##################################################################
+
+License & disclaimer:
+
+    Copyright 2008-2009 UCL Business PLC, Author Vladimir Kolmogorov (vnk@ist.ac.at)
+
+    This software can be used for evaluation and non-commercial research purposes only. Commercial use is prohibited.
+    Public redistribution of the code or its derivatives is prohibited.
+    If you use this software for research purposes, you should cite the following paper in any resulting publication:
+
+        Vladimir Kolmogorov. "Blossom V: A new implementation of a minimum cost perfect matching algorithm."
+        In Mathematical Programming Computation (MPC), July 2009, 1(1):43-67.
+
+    For commercial use of the software not covered by this agreement, you may obtain a licence from
+    the copyright holders UCL Business via their licensing site: www.e-lucid.com/i/software/optimisation_software/BlossomV.html.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+##################################################################
+
+COMPILATION:
+
+In unix, type "make". The code should also compile in Windows with Microsoft Visual C++ compiler.
+Tested on 32-bit machines.
+
+##################################################################
+
+USAGE:
+
+See PerfectMatching.h for interface functions. Alternatively,
+compile the code and run ./blossom5 [options] (see USAGE.TXT for details).
+
+The code also allows solving complete geometric instances; see GEOM/GeomPerfectMatching.h
+for interface functions.
+
+##################################################################
+
+PARAMETERS:
+
+For many types of problems, default parameters should be ok.
+But for certain problems (such as structured geometric instances)
+you may consider setting 
+dual_LP_threshold=0.005 (for example). This corresponds to calling
+./blossom5 -m0.005
+Type PerfectMatching::REAL in PerfectMatcing.h should then be set to double.
+
+##################################################################
+
+EXTERNAL PACKAGE:
+
+When solving complete geometric instances you need to provide the initial subset of edges.
+It may be desirable to use Delaunay triangulation for this purpose.
+Then you need to download the "Triangle" package of J. R. Shewchuk from 
+	http://www.cs.cmu.edu/~quake/triangle.html
+and extract it to the directory "triange".
+Alternatively, you can use nearest neighbours initialization (this does not require external packages).
+
diff --git a/blossom5-v2.05.src/USAGE.TXT b/blossom5-v2.05.src/USAGE.TXT
new file mode 100644
index 0000000..afe93cd
--- /dev/null
+++ b/blossom5-v2.05.src/USAGE.TXT
@@ -0,0 +1,31 @@
+USAGE: ./blossom5 [options]
+
+At least one of the two flags below must be provided. 
+If both are provided, then edges in <f1> are added to the initial subset of edges for the geometric problem.
+
+  -e <f1> read problem (list of edges) from file <f1> in DIMACS or blossom4 format
+          (see GRAPH1.TXT and GRAPH2.TXT for examples)
+  -g <f2> read geometric problem (list of 2D points) from file <f2> in TSPLIB format 
+
+perfect matching parameters:
+
+  -j     do not use fractional jumpstart
+  -m#    update duals by solving LP if the number of trees is smaller than # times node_num. 
+         (# should belong to [0,1], default is 0.)
+  -d1    use SCC dual updates instead of CC
+  -d2    use dual updates with fixed delta instead of CC
+  -b     update duals before processing a tree
+  -a     update duals after processing a tree
+
+geometric matching parameters:
+
+  -D     do NOT add Delaunay triangulation to the initial subset of edges
+  -K#    for each point, add # nearest neighbors to the initial subset of edges
+  -I     do NOT add edges greedily to the initial subset of edges to make sure that a perfect matching exists
+  -T#    use at most # iterations (or run until convergence, if #=0 - default)
+
+other options:
+
+  -w <f> save result to file <f>
+  -c     check complementary slackness conditions for non-geometric problem
+  -V     no verbose
diff --git a/blossom5-v2.05.src/block.h b/blossom5-v2.05.src/block.h
new file mode 100644
index 0000000..3bed9d4
--- /dev/null
+++ b/blossom5-v2.05.src/block.h
@@ -0,0 +1,268 @@
+/* block.h */
+/*
+	Template classes Block and DBlock
+	Implement adding and deleting items of the same type in blocks.
+
+	If there there are many items then using Block or DBlock
+	is more efficient than using 'new' and 'delete' both in terms
+	of memory and time since
+	(1) On some systems there is some minimum amount of memory
+	    that 'new' can allocate (e.g., 64), so if items are
+	    small that a lot of memory is wasted.
+	(2) 'new' and 'delete' are designed for items of varying size.
+	    If all items has the same size, then an algorithm for
+	    adding and deleting can be made more efficient.
+	(3) All Block and DBlock functions are inline, so there are
+	    no extra function calls.
+
+	Differences between Block and DBlock:
+	(1) DBlock allows both adding and deleting items,
+	    whereas Block allows only adding items.
+	(2) Block has an additional operation of scanning
+	    items added so far (in the order in which they were added).
+	(3) Block allows to allocate several consecutive
+	    items at a time, whereas DBlock can add only a single item.
+
+	Note that no constructors or destructors are called for items.
+
+	Example usage for items of type 'MyType':
+
+	///////////////////////////////////////////////////
+	#include "block.h"
+	#define BLOCK_SIZE 1024
+	typedef struct { int a, b; } MyType;
+	MyType *ptr, *array[10000];
+
+	...
+
+	Block<MyType> *block = new Block<MyType>(BLOCK_SIZE);
+
+	// adding items
+	for (int i=0; i<sizeof(array); i++)
+	{
+		ptr = block -> New();
+		ptr -> a = ptr -> b = rand();
+	}
+
+	// reading items
+	for (ptr=block->ScanFirst(); ptr; ptr=block->ScanNext())
+	{
+		printf("%d %d\n", ptr->a, ptr->b);
+	}
+
+	delete block;
+
+	...
+
+	DBlock<MyType> *dblock = new DBlock<MyType>(BLOCK_SIZE);
+	
+	// adding items
+	for (int i=0; i<sizeof(array); i++)
+	{
+		array[i] = dblock -> New();
+	}
+
+	// deleting items
+	for (int i=0; i<sizeof(array); i+=2)
+	{
+		dblock -> Delete(array[i]);
+	}
+
+	// adding items
+	for (int i=0; i<sizeof(array); i++)
+	{
+		array[i] = dblock -> New();
+	}
+
+	delete dblock;
+
+	///////////////////////////////////////////////////
+
+	Note that DBlock deletes items by marking them as
+	empty (i.e., by adding them to the list of free items),
+	so that this memory could be used for subsequently
+	added items. Thus, at each moment the memory allocated
+	is determined by the maximum number of items allocated
+	simultaneously at earlier moments. All memory is
+	deallocated only when the destructor is called.
+*/
+
+#ifndef __BLOCK_H__
+#define __BLOCK_H__
+
+#include <stdlib.h>
+
+/***********************************************************************/
+/***********************************************************************/
+/***********************************************************************/
+
+template <class Type> class Block
+{
+public:
+	/* Constructor. Arguments are the block size and
+	   (optionally) the pointer to the function which
+	   will be called if allocation failed; the message
+	   passed to this function is "Not enough memory!" */
+	Block(int size, void (*err_function)(const char *) = NULL) { first = last = NULL; block_size = size; error_function = err_function; }
+
+	/* Destructor. Deallocates all items added so far */
+	~Block() { while (first) { block *next = first -> next; delete[] ((char*)first); first = next; } }
+
+	/* Allocates 'num' consecutive items; returns pointer
+	   to the first item. 'num' cannot be greater than the
+	   block size since items must fit in one block */
+	Type *New(int num = 1)
+	{
+		Type *t;
+
+		if (!last || last->current + num > last->last)
+		{
+			if (last && last->next) last = last -> next;
+			else
+			{
+				block *next = (block *) new char [sizeof(block) + (block_size-1)*sizeof(Type)];
+				if (!next) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
+				if (last) last -> next = next;
+				else first = next;
+				last = next;
+				last -> current = & ( last -> data[0] );
+				last -> last = last -> current + block_size;
+				last -> next = NULL;
+			}
+		}
+
+		t = last -> current;
+		last -> current += num;
+		return t;
+	}
+
+	/* Returns the first item (or NULL, if no items were added) */
+	Type *ScanFirst()
+	{
+		for (scan_current_block=first; scan_current_block; scan_current_block = scan_current_block->next)
+		{
+			scan_current_data = & ( scan_current_block -> data[0] );
+			if (scan_current_data < scan_current_block -> current) return scan_current_data ++;
+		}
+		return NULL;
+	}
+
+	/* Returns the next item (or NULL, if all items have been read)
+	   Can be called only if previous ScanFirst() or ScanNext()
+	   call returned not NULL. */
+	Type *ScanNext()
+	{
+		while (scan_current_data >= scan_current_block -> current)
+		{
+			scan_current_block = scan_current_block -> next;
+			if (!scan_current_block) return NULL;
+			scan_current_data = & ( scan_current_block -> data[0] );
+		}
+		return scan_current_data ++;
+	}
+
+	/* Marks all elements as empty */
+	void Reset()
+	{
+		block *b;
+		if (!first) return;
+		for (b=first; ; b=b->next)
+		{
+			b -> current = & ( b -> data[0] );
+			if (b == last) break;
+		}
+		last = first;
+	}
+
+/***********************************************************************/
+
+private:
+
+	typedef struct block_st
+	{
+		Type					*current, *last;
+		struct block_st			*next;
+		Type					data[1];
+	} block;
+
+	int		block_size;
+	block	*first;
+	block	*last;
+
+	block	*scan_current_block;
+	Type	*scan_current_data;
+
+	void	(*error_function)(const char *);
+};
+
+/***********************************************************************/
+/***********************************************************************/
+/***********************************************************************/
+
+template <class Type> class DBlock
+{
+public:
+	/* Constructor. Arguments are the block size and
+	   (optionally) the pointer to the function which
+	   will be called if allocation failed; the message
+	   passed to this function is "Not enough memory!" */
+	DBlock(int size, void (*err_function)(const char *) = NULL) { first = NULL; first_free = NULL; block_size = size; error_function = err_function; }
+
+	/* Destructor. Deallocates all items added so far */
+	~DBlock() { while (first) { block *next = first -> next; delete[] ((char*)first); first = next; } }
+
+	/* Allocates one item */
+	Type *New()
+	{
+		block_item *item;
+
+		if (!first_free)
+		{
+			block *next = first;
+			first = (block *) new char [sizeof(block) + (block_size-1)*sizeof(block_item)];
+			if (!first) { if (error_function) (*error_function)("Not enough memory!"); exit(1); }
+			first_free = & (first -> data[0] );
+			for (item=first_free; item<first_free+block_size-1; item++)
+				item -> next_free = item + 1;
+			item -> next_free = NULL;
+			first -> next = next;
+		}
+
+		item = first_free;
+		first_free = item -> next_free;
+		return (Type *) item;
+	}
+
+	/* Deletes an item allocated previously */
+	void Delete(Type *t)
+	{
+		((block_item *) t) -> next_free = first_free;
+		first_free = (block_item *) t;
+	}
+
+/***********************************************************************/
+
+private:
+
+	typedef union block_item_st
+	{
+		Type			t;
+		block_item_st	*next_free;
+	} block_item;
+
+	typedef struct block_st
+	{
+		struct block_st			*next;
+		block_item				data[1];
+	} block;
+
+	int			block_size;
+	block		*first;
+	block_item	*first_free;
+
+	void	(*error_function)(const char *);
+};
+
+
+#endif
+
diff --git a/blossom5-v2.05.src/example.cpp b/blossom5-v2.05.src/example.cpp
new file mode 100644
index 0000000..16c962e
--- /dev/null
+++ b/blossom5-v2.05.src/example.cpp
@@ -0,0 +1,285 @@
+#include <stdio.h>
+#include "PerfectMatching.h"
+#include "GEOM/GeomPerfectMatching.h"
+
+void LoadFile(int& node_num, int& edge_num, int*& edges, int*& weights, char* filename)
+{
+	int e = 0;
+	char LINE[1000];
+	FILE* fp = fopen(filename, "r");
+	if (!fp) { printf("Can't open %s\n", filename); exit(1); }
+
+	int format = -1; // 0: DIMACS format. node id's start from 1
+	                 // 1: simpler format (without "p" and "e"). node id's start from 0
+
+	edge_num = -1;
+	while (fgets(LINE, sizeof(LINE)-1, fp))
+	{
+		if (LINE[0] == 'c') continue;
+		if (format < 0)
+		{
+			if (LINE[0] == 'p')
+			{
+				format = 0;
+				if (sscanf(LINE, "p edge %d %d\n", &node_num, &edge_num) != 2) { printf("%s: wrong format #1\n", filename); exit(1); }
+			}
+			else
+			{
+				format = 1;
+				if (sscanf(LINE, "%d %d\n", &node_num, &edge_num) != 2) { printf("%s: wrong format #1\n", filename); exit(1); }
+			}
+
+			//////////////////////////////////////////////////////////////////////////////////
+			if (node_num <= 0 || edge_num < 0) { printf("# of nodes and edges should be positive\n"); exit(1); }
+			if (node_num & 1) { printf("# of nodes is odd: perfect matching cannot exist\n"); exit(1); }
+			edges = new int[2*edge_num];
+			weights = new int[edge_num];
+			//////////////////////////////////////////////////////////////////////////////////
+		}
+		else
+		{
+			int i, j;
+			char* ptr = LINE;
+			if (format == 0) { if (LINE[0] != 'e') continue; ptr = &LINE[1]; }
+			else             ptr = &LINE[0];
+
+			int len;
+			if (sscanf(ptr, "%d %d %d\n", &i, &j, &len) != 3) continue;
+			if (format == 0) { i --; j --; }
+			edges[2*e] = i;
+			edges[2*e+1] = j;
+			weights[e] = len;
+			e ++;
+		}
+	}
+
+	if (e != edge_num) { printf("%s: wrong format #3\n", filename); exit(1); }
+	fclose(fp);
+}
+
+void SaveMatching(int node_num, PerfectMatching* pm, char* filename)
+{
+	FILE* fp = fopen(filename, "w");
+	if (!fp) { printf("Can't open %s\n", filename); exit(1); }
+	fprintf(fp, "%d %d\n", node_num, node_num/2);
+	int i, j;
+	for (i=0; i<node_num; i++)
+	{
+		j = pm->GetMatch(i);
+		if (i < j) fprintf(fp, "%d %d\n", i, j);
+	}
+	fclose(fp);
+}
+
+void LoadGeomFile(int& node_num, int*& x_array, int*& y_array, char* filename)
+{
+	int i = 0, i_tmp, x, y, DIM = 0;
+	char LINE[1000];
+
+	x_array = y_array = NULL;
+
+	FILE* fp = fopen(filename, "r");
+	if (!fp) { printf("Can't open %s\n", filename); exit(1); }
+
+	while (fgets(LINE, sizeof(LINE)-1, fp))
+	{
+		if (sscanf(LINE, "DIMENSION : %d", &node_num) == 1)
+		{
+			if (node_num < 1) { printf("too few nodes\n"); exit(1); }
+			if (node_num & 1) { printf("# of points is odd: perfect matching cannot exist\n"); exit(1); }
+			if (x_array) { printf("wrong format\n"); exit(1); }
+			x_array = new int[node_num];
+			y_array = new int[node_num];
+			continue;
+		}
+		if (sscanf(LINE, "EDGE_WEIGHT_TYPE : EUC_%dD", &DIM) == 1)
+		{
+			if (DIM != 2) { printf("only EUC_2D is supported"); exit(1); }
+			continue;
+		}
+		if (sscanf(LINE, "%d %d %d", &i_tmp, &x, &y) == 3)
+		{
+			i_tmp --;
+			if (i_tmp != i ++ || i > node_num) { printf("wrong number of points\n"); exit(1); }
+			x_array[i_tmp] = x;
+			y_array[i_tmp] = y;
+			continue;
+		}
+		printf("%s", LINE);
+	}
+	fclose(fp);
+	if (i != node_num || !x_array || DIM != 2) { printf("wrong format\n"); exit(1); }
+}
+
+
+void SaveMatching(int node_num, GeomPerfectMatching* gpm, char* filename)
+{
+	FILE* fp = fopen(filename, "w");
+	if (!fp) { printf("Can't open %s\n", filename); exit(1); }
+	fprintf(fp, "%d %d\n", node_num, node_num/2);
+	int i, j;
+	for (i=0; i<node_num; i++)
+	{
+		j = gpm->GetMatch(i);
+		if (i < j)
+		{
+			GeomPerfectMatching::REAL len = gpm->Dist(i, j);
+			if ( ((GeomPerfectMatching::REAL)1 / 2) == 0 ) fprintf(fp, "%d %d %d\n", i, j, len);
+			else                                           fprintf(fp, "%d %d %f\n", i, j, (double)len);
+		}
+	}
+	fclose(fp);
+}
+
+
+void ShowUsage()
+{
+	printf("Usage: see USAGE.TXT\n");
+	exit(1);
+}
+
+
+int main(int argc, char* argv[])
+{
+	struct PerfectMatching::Options options;
+	struct GeomPerfectMatching::GPMOptions gpm_options;
+	char* filename = NULL;
+	char* geom_filename = NULL;
+	char* save_filename = NULL;
+	bool check_perfect_matching = false;
+	int i, e, node_num, edge_num;
+	int* edges;
+	int* weights;
+
+	for (i=1; i<argc; i++)
+	{
+		if (argv[i][0] != '-') { printf("Unknown option: %s\n", argv[i]); ShowUsage(); }
+		switch (argv[i][1])
+		{
+			case 'e':
+				if (filename || argv[i][2] || ++i == argc) ShowUsage();
+				filename = argv[i];
+				break;
+			case 'g':
+				if (geom_filename || argv[i][2] || ++i == argc) ShowUsage();
+				geom_filename = argv[i];
+				break;
+			case 'j':
+				if (argv[i][2]) ShowUsage();
+				options.fractional_jumpstart = false;
+				break;
+			case 'm':
+				options.dual_LP_threshold = atof(&argv[i][2]);
+				if (options.dual_LP_threshold<0 || options.dual_LP_threshold>1) ShowUsage();
+				break;
+			case 'd':
+				options.dual_greedy_update_option = atoi(&argv[i][2]);
+				if (options.dual_greedy_update_option<1 || options.dual_greedy_update_option>2) ShowUsage();
+				break;
+			case 'b':
+				if (argv[i][2]) ShowUsage();
+				options.update_duals_before = true;
+				break;
+			case 'a':
+				if (argv[i][2]) ShowUsage();
+				options.update_duals_after = true;
+				break;
+			case 'D':
+				if (argv[i][2]) ShowUsage();
+				gpm_options.init_Delaunay = false;
+				break;
+			case 'K':
+				gpm_options.init_KNN = atoi(&argv[i][2]);
+				if (gpm_options.init_KNN<0) ShowUsage();
+				break;
+			case 'I':
+				if (argv[i][2]) ShowUsage();
+				gpm_options.init_greedy = false;
+				break;
+			case 'T':
+				gpm_options.iter_max = atoi(&argv[i][2]);
+				if (gpm_options.iter_max<0) ShowUsage();
+				break;
+			case 'w':
+				if (save_filename || argv[i][2] || ++i == argc) ShowUsage();
+				save_filename = argv[i];
+				break;
+			case 'c':
+				if (argv[i][2]) ShowUsage();
+				check_perfect_matching = true;
+				break;
+			case 'V':
+				if (argv[i][2]) ShowUsage();
+				options.verbose = false;
+				break;
+			default:
+				printf("Unknown option: %s\n", argv[i]);
+				ShowUsage();
+				break;
+		}
+		
+	}
+
+	if (!filename && !geom_filename) ShowUsage();
+
+	if (filename) LoadFile(node_num, edge_num, edges, weights, filename);
+
+	if (!geom_filename)
+	{
+		PerfectMatching *pm = new PerfectMatching(node_num, edge_num);
+		for (e=0; e<edge_num; e++) pm->AddEdge(edges[2*e], edges[2*e+1], weights[e]);
+		pm->options = options;
+		pm->Solve();
+		if (check_perfect_matching)
+		{
+			int res = CheckPerfectMatchingOptimality(node_num, edge_num, edges, weights, pm);
+			printf("check optimality: res=%d (%s)\n", res, (res==0) ? "ok" : ((res==1) ? "error" : "fatal error"));
+		}
+		double cost = ComputePerfectMatchingCost(node_num, edge_num, edges, weights, pm);
+		printf("cost = %.1f\n", cost);
+		if (save_filename) SaveMatching(node_num, pm, save_filename);
+		delete pm;
+	}
+	else
+	{
+		int geom_node_num;
+		int* x_array;
+		int* y_array;
+		LoadGeomFile(geom_node_num, x_array, y_array, geom_filename);
+		GeomPerfectMatching *gpm = new GeomPerfectMatching(geom_node_num, 2);
+		for (i=0; i<geom_node_num; i++)
+		{
+			GeomPerfectMatching::REAL coord[2];
+			coord[0] = x_array[i];
+			coord[1] = y_array[i];
+			gpm->AddPoint(coord);
+		}
+		delete [] x_array;
+		delete [] y_array;
+		gpm->options = options;
+		gpm->gpm_options = gpm_options;
+		if (filename)
+		{
+			if (node_num != geom_node_num) { printf("%s and %s don't match!\n", geom_filename, filename); exit(1); }
+			for (e=0; e<edge_num; e++)
+			{
+				if (weights[e] != gpm->Dist(edges[2*e], edges[2*e+1]))
+				{ printf("edge lengths in %s and %s don't match!\n", geom_filename, filename); exit(1); }
+				gpm->AddInitialEdge(edges[2*e], edges[2*e+1]);
+			}
+		}
+		gpm->Solve();
+		if (save_filename) SaveMatching(geom_node_num, gpm, save_filename);
+		delete gpm;
+	}
+
+	if (filename)
+	{
+		delete [] edges;
+		delete [] weights;
+	}
+
+	return 0;
+}
+
+
diff --git a/blossom5-v2.05.src/misc.cpp b/blossom5-v2.05.src/misc.cpp
new file mode 100644
index 0000000..5bdfd3e
--- /dev/null
+++ b/blossom5-v2.05.src/misc.cpp
@@ -0,0 +1,172 @@
+#include <stdio.h>
+#include "PerfectMatching.h"
+#include "LCA.h"
+
+struct Node
+{
+	PerfectMatching::REAL sum; // = twice_y[i] + twice_y[i->parent] + twice_y[i->parent->parent] + ...
+	Node* match;
+	Node* parent;
+	Node* child;
+	Node* sibling;
+	int lca_preorder;
+};
+
+
+int CheckPerfectMatchingOptimality(int node_num, int edge_num, int* edges, int* weights, PerfectMatching* pm, PerfectMatching::REAL threshold)
+{
+	int _i, _j, _e;
+	Node* i;
+	Node* j;
+	int blossom_num = pm->GetBlossomNum();
+	int* blossom_parents = new int[node_num+blossom_num];
+	PerfectMatching::REAL* twice_y = new PerfectMatching::REAL[node_num+blossom_num];
+
+	PerfectMatching::REAL y_blossom_min = 0;
+	PerfectMatching::REAL slack_min = 0;
+	PerfectMatching::REAL active_slack_max = 0;
+
+	// step 1 - read dual solution and construct tree
+	pm->GetDualSolution(blossom_parents, twice_y);
+	Node* nodes = new Node[node_num+blossom_num+1];
+	memset(nodes, 0, (node_num+blossom_num+1)*sizeof(Node));
+	Node* ROOT = nodes+node_num+blossom_num;
+	for (_i=0, i=nodes; _i<node_num+blossom_num; _i++, i++)
+	{
+		i->sum = twice_y[_i];
+		if (_i >= node_num && y_blossom_min > i->sum) y_blossom_min = i->sum;
+		if (blossom_parents[_i] >= 0)
+		{
+			if (blossom_parents[_i]<node_num || blossom_parents[_i]>=node_num+blossom_num)
+			{
+				delete [] nodes;
+				delete [] blossom_parents;
+				delete [] twice_y;
+				return 2;
+			}
+			i->parent = nodes + blossom_parents[_i];
+			i->sibling = i->parent->child;
+			i->parent->child = i;
+		}
+	}
+	delete [] blossom_parents;
+	delete [] twice_y;
+
+	for (i=nodes; i<nodes+node_num+blossom_num; i++)
+	{
+		if (!i->parent)
+		{
+			i->parent = ROOT;
+			i->sibling = ROOT->child;
+			ROOT->child = i;
+		}
+	}
+
+	LCATree* lca_tree = new LCATree(node_num+blossom_num+1);
+	Node** rev_mapping = new Node*[node_num+blossom_num];
+
+	i = ROOT;
+	while ( 1 )
+	{
+		if (i->child)
+		{
+			if (i < nodes+node_num) { delete [] nodes; delete lca_tree; delete [] rev_mapping; return 2; }
+			i->child->sum += i->sum;
+			i = i->child;
+		}
+		else
+		{
+			if (i >= nodes+node_num) { delete [] nodes; delete lca_tree; delete [] rev_mapping; return 2; }
+			while ( 1 )
+			{
+				i->lca_preorder = lca_tree->Add(i, i->parent);
+				rev_mapping[i->lca_preorder] = i;
+				if (i->sibling) break;
+				i = i->parent;
+				if (i == ROOT)
+				{
+					i->lca_preorder = lca_tree->AddRoot(i);
+					break;
+				}
+			}
+			if (i == ROOT) break;
+			i = i->sibling;
+			i->sum += i->parent->sum;
+		}
+	}
+
+	int matched_num = 0;
+	for (_e=0; _e<edge_num; _e++)
+	{
+		_i = edges[2*_e];
+		_j = edges[2*_e+1];
+		if (_i<0 || _j<0 || _i>=node_num || _j>=node_num || _i==_j) { delete [] nodes; delete lca_tree; delete [] rev_mapping; return 2; }
+
+		int lca_i = nodes[_i].lca_preorder;
+		int lca_j = nodes[_j].lca_preorder;
+		lca_tree->GetPenultimateNodes(lca_i, lca_j);
+		i = rev_mapping[lca_i];
+		j = rev_mapping[lca_j];
+		PerfectMatching::REAL twice_slack = 2*weights[_e] - (nodes[_i].sum - i->parent->sum) - (nodes[_j].sum - j->parent->sum);
+		if (slack_min > twice_slack) slack_min = twice_slack;
+		if (pm->GetSolution(_e))
+		{
+			if (pm->GetMatch(_i)!=_j || pm->GetMatch(_j)!=_i || i->match || j->match) { delete [] nodes; delete lca_tree; delete [] rev_mapping; return 2; }
+			i->match = j;
+			j->match = i;
+			if (active_slack_max < twice_slack) active_slack_max = twice_slack;
+			matched_num += 2;
+		}
+	}
+
+	delete [] nodes;
+	delete lca_tree;
+	delete [] rev_mapping;
+
+	if (matched_num != node_num) return 2;
+
+	if (y_blossom_min < -threshold || slack_min < -threshold || active_slack_max > threshold)
+	{
+		printf("ERROR in CheckPerfectMatchingOptimality():\n");
+		if ( ((PerfectMatching::REAL)1 / 2) == 0 )
+			printf("\ty_blossom_min=%d\n\tslack_min=%d\n\tactive_slack_max=%d\n", (int)y_blossom_min, (int)slack_min, (int)active_slack_max);
+		else
+			printf("\ty_blossom_min=%.15f\n\tslack_min=%.15f\n\tactive_slack_max=%.15f\n", (double)y_blossom_min, (double)slack_min, (double)active_slack_max);
+		return 1;
+	}
+
+	return 0;
+}
+
+double ComputePerfectMatchingCost(int node_num, int edge_num, int* edges, int* weights, PerfectMatching* pm)
+{
+	int i;
+	int j;
+	int e;
+	double cost = 0;
+
+	int* nodes = new int[node_num];
+	memset(nodes, 0, node_num*sizeof(int));
+	for (e=0; e<edge_num; e++)
+	{
+		if (pm->GetSolution(e))
+		{
+			i = edges[2*e];
+			j = edges[2*e+1];
+			nodes[i] ++;
+			nodes[j] ++;
+			cost += weights[e];
+		}
+	}
+	for (i=0; i<node_num; i++)
+	{
+		if (nodes[i] != 1)
+		{
+			printf("ComputeCost(): degree = %d!\n", nodes[i]);
+			exit(1);
+		}
+	}
+	delete [] nodes;
+	return cost;
+}
+
diff --git a/blossom5-v2.05.src/timer.h b/blossom5-v2.05.src/timer.h
new file mode 100644
index 0000000..a0b9373
--- /dev/null
+++ b/blossom5-v2.05.src/timer.h
@@ -0,0 +1,98 @@
+#ifndef NJAKSDTHASKJERAXJGFBZJDLAGZ
+#define NJAKSDTHASKJERAXJGFBZJDLAGZ
+
+// At most one of the flags
+//    PM_TIMER_MSVC
+//    PM_TIMER_CLOCK_GETTIME
+//    PM_TIMER_GETRUSAGE
+//    PM_TIMER_EXTERNAL
+//    PM_TIMER_NONE
+// can be defined externally. If PM_TIMER_EXTERNAL is defined,
+// then there must exist a definition of function "double get_time()" elsewhere.
+
+#if defined (PM_TIMER_MSVC) || defined (PM_TIMER_CLOCK_GETTIME) || defined (PM_TIMER_GETRUSAGE) || defined (PM_TIMER_EXTERNAL) || defined (PM_TIMER_NONE)
+#else
+	// default option
+	#ifdef _MSC_VER
+		#define PM_TIMER_MSVC
+	#elif defined(__APPLE_CC__)
+		#define PM_TIMER_GETRUSAGE
+	#else
+		#define PM_TIMER_CLOCK_GETTIME
+	#endif
+#endif
+
+//////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+
+#ifdef PM_TIMER_MSVC
+
+	#include <windows.h>
+
+	inline double get_time()
+	{
+		LARGE_INTEGER t, frequency;
+		QueryPerformanceCounter(&t);
+		QueryPerformanceFrequency(&frequency);
+		return (double)t.QuadPart/(double)frequency.QuadPart;
+	}
+
+#endif
+
+//////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+
+#ifdef PM_TIMER_CLOCK_GETTIME
+
+	#include <time.h>
+
+	inline double get_time()
+	{
+		struct timespec t;
+		clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &t);
+		return (double)t.tv_nsec*1.00E-9 + (double)t.tv_sec;
+	}
+
+#endif
+
+//////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+
+#ifdef PM_TIMER_GETRUSAGE
+
+	#include <sys/resource.h>
+
+	inline double get_time()
+	{
+		struct rusage t;
+		getrusage (RUSAGE_SELF, &t);
+		return (double)t.ru_utime.tv_usec*1.00E-6 + (double)t.ru_utime.tv_sec;
+	}
+
+#endif
+
+//////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+
+#ifdef PM_TIMER_EXTERNAL
+
+	extern double get_time();
+
+#endif
+
+//////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
+
+#ifdef PM_TIMER_NONE
+
+	inline double get_time() { return 0; }
+
+#endif
+
+#endif
+
diff --git a/hungarian.cpp b/hungarian.cpp
index 00cc91d..b8cb7ab 100644
--- a/hungarian.cpp
+++ b/hungarian.cpp
@@ -75,7 +75,7 @@ public:
     }
   }
   bool isTight() const {
-    return std::abs(s->y + t->y - weight) < 1e-10;
+    return std::abs(s->y + t->y - weight) < 1e-13;
   }
 };
 
@@ -166,10 +166,10 @@ public:
             if (e.isTight() && (e.head().index == i && !e.inPath)) {
               if (!e.tail().inMatching()) {
                 break;
+              } else {
+                e.inPath= true;
+                path.push({e.tail().neighbors.begin(), e.tail().neighbors.end(), e.tail().index});
               }
-
-              e.inPath= true;
-              path.push({e.tail().neighbors.begin(), e.tail().neighbors.end(), e.tail().index});
             } else {
               it++;
               path.top() = {it, itEnd, i};
@@ -177,12 +177,6 @@ public:
           }
         }
 
-
-        if (path.empty()) {
-          exit(1);
-        }
-
-
         while (!path.empty()) {
           auto [it, itEnd, i] = path.top();
           Edge& e = it->get();
diff --git a/rbmp.cpp b/rbmp.cpp
index ffc7c44..ddd76cb 100644
--- a/rbmp.cpp
+++ b/rbmp.cpp
@@ -8,6 +8,8 @@
 #include "randutils/randutils.hpp"
 #include "pcg-cpp/include/pcg_random.hpp"
 
+#include "blossom5-v2.05.src/PerfectMatching.h"
+
 using Rng = randutils::random_generator<pcg32>;
 
 class Edge;
@@ -31,9 +33,8 @@ private:
 
 public:
   const Edge& edge;
-  double X;
 
-  HalfEdge(const Edge& e) : edge(e), X(0) {}
+  HalfEdge(const Edge& e) : edge(e) {}
   void setVertices(Vertex& vt, Vertex& vh) {
     tail = &vt;
     head = &vh;
@@ -58,12 +59,6 @@ public:
     red.addEdge(halfedges[0]);
     blue.addEdge(halfedges[1]);
   }
-  double belief() const {
-    return halfedges[0].X + halfedges[1].X;
-  }
-  bool active() const {
-    return belief() >= 0;
-  }
 };
 
 class Graph {
@@ -86,51 +81,18 @@ public:
       edges[i].weight = r.variate<double, std::exponential_distribution>(1);
     }
   }
-
-  void propagateBeliefs() {
-    for (unsigned i : {0, 1}) {
-#pragma omp parallel for
-      for (Edge& e : edges) {
-        HalfEdge& h = e.halfedges[i];
-        h.X = std::numeric_limits<double>::infinity();
-        for (const HalfEdge& hn : h.getHead().neighbors) {
-          if (h.getTail().index != hn.getHead().index) {
-            h.X = std::min(hn.edge.weight - hn.X, h.X);
-          }
-        }
-      }
-    }
-  }
-
-  double minBelief() const {
-    double minBelief = std::numeric_limits<double>::infinity();
-
-    for (const Edge& e : edges) {
-      minBelief = std::min(std::abs(e.belief()), minBelief);
-    }
-
-    return minBelief;
-  }
 };
 
 int main(int argc, char* argv[]) {
   unsigned n = 100;
-  unsigned maxSteps = 1e8;
-  double beliefThreshold = 1;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "n:m:t:")) != -1) {
+  while ((opt = getopt(argc, argv, "n:")) != -1) {
     switch (opt) {
     case 'n':
       n = atoi(optarg);
       break;
-    case 'm':
-      maxSteps = (unsigned)atof(optarg);
-      break;
-    case 't':
-      beliefThreshold = atof(optarg);
-      break;
     default:
       exit(1);
     }
@@ -139,26 +101,60 @@ int main(int argc, char* argv[]) {
   Rng r;
   Graph G(n, r);
 
-  unsigned steps = 0;
-  double minBelief = 0;
+  PerfectMatching pm(G.vertices.size(), G.edges.size());
 
-  while (minBelief < beliefThreshold && steps < maxSteps) {
-    G.propagateBeliefs();
-    if (steps % 100 == 0) {
-      std::cerr << (minBelief = G.minBelief()) << std::endl;
-    }
-    steps++;
+  for (const Edge& e : G.edges) {
+    pm.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, e.weight);
   }
 
+  pm.options.verbose = false;
+  pm.Solve();
+
+  std::reference_wrapper<Edge> eFlip = r.pick(G.edges);
+
+  while (pm.GetMatch(eFlip.get().halfedges[0].getTail().index) != eFlip.get().halfedges[0].getHead().index) {
+    eFlip = r.pick(G.edges);
+  }
+
+  eFlip.get().weight = 1e10;
+
+  PerfectMatching pm2(G.vertices.size(), G.edges.size());
+
   for (const Edge& e : G.edges) {
-    if (e.active()) {
+    pm2.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, e.weight);
+  }
+
+  pm2.options.verbose = false;
+  pm2.Solve();
+
+  std::cout << n << std::endl;
+
+  for (unsigned i = 0; i < G.vertices.size() / 2; i++) {
+    unsigned j1 = pm.GetMatch(i);
+
+    std::cout
+      << G.vertices[i].coordinate[0] << " "
+      << G.vertices[i].coordinate[1] << " "
+      << G.vertices[j1].coordinate[0] << " "
+      << G.vertices[j1].coordinate[1] << std::endl;
+  }
+
+  unsigned m = 0;
+  for (unsigned i = 0; i < G.vertices.size() / 2; i++) {
+    unsigned j1 = pm.GetMatch(i);
+    unsigned j2 = pm2.GetMatch(i);
+
+    if (j1 != j2) {
+      m++;
       std::cout
-        << e.halfedges[0].getTail().coordinate[0] << " "
-        << e.halfedges[0].getTail().coordinate[1] << " "
-        << e.halfedges[0].getHead().coordinate[0] << " "
-        << e.halfedges[0].getHead().coordinate[1] << std::endl;
+        << G.vertices[i].coordinate[0] << " "
+        << G.vertices[i].coordinate[1] << " "
+        << G.vertices[j2].coordinate[0] << " "
+        << G.vertices[j2].coordinate[1] << std::endl;
     }
   }
 
+  std::cerr << "Size of excitation: " << m << std::endl;
+
   return 0;
 }