parallel implementation of edmonds' algorithm

edmonds.c

#define _GNU_SOURCE
#include <stdint.h>
#include <stdatomic.h>
#include <stdbool.h>
#include <threads.h>
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <time.h>
#include <assert.h>

#ifdef TRACY
  #include "tracy/public/tracy/TracyC.h"
#else
  #define TracyCZone(...)
  #define TracyCZoneS(...)
  #define TracyCZoneEnd(...)
#endif

typedef uint32_t vertex_t;
typedef int32_t cost_t;

#define NIL_VERTEX UINT32_MAX
#define INF_WEIGHT INT32_MAX

struct inbound_edge {
	vertex_t source;
	cost_t cost;
};

struct outbound_edge {
	vertex_t dest;
	cost_t cost;
};

struct full_edge {
	vertex_t source;
	vertex_t dest;
	cost_t cost;
};

struct inbound_star {
	uint32_t num_edges;
	struct inbound_edge *edges;
};

struct outbound_star {
	uint32_t num_edges;
	struct outbound_edge *edges;
};

struct edge_list {
	_Atomic uint32_t num_edges;
	struct full_edge *edges;
};

struct inbound_graph {
	uint32_t num_nodes;
	struct inbound_star *star;
};

struct outbound_graph {
	uint32_t num_nodes;
	struct outbound_star *star;
};

struct thread_context {
	vertex_t root;
	// used as a work queue, increases monotonically
	_Atomic vertex_t next_job;
	struct inbound_graph input_graph;

	struct leftist_heap_node *heap;

	// scc partition
	vertex_t *path;
	vertex_t *heap_root;
	_Atomic vertex_t *parent0;
	_Atomic uint32_t *count0;

	// tree partition
	_Atomic vertex_t *parent1;
	_Atomic uint32_t *count1;
	atomic_flag *lock;

	struct edge_list solution;
};

struct leftist_heap_node {
	cost_t bias;
	vertex_t left;
	vertex_t right;
	uint32_t count;
};

static inline void swap_costs(cost_t *u, cost_t *v) { 
	cost_t t = *u;
	*u = *v;
	*v = t;
}

static inline void swap_vertices(vertex_t *u, vertex_t *v) { 
	vertex_t t = *u;
	*u = *v;
	*v = t;
}

static inline void swap_inbound_edges(struct inbound_edge *u, struct inbound_edge *v) { 
	struct inbound_edge t = *u;
	*u = *v;
	*v = t;
}

static inline void swap_full_edges(struct full_edge *u, struct full_edge *v) { 
	struct full_edge t = *u;
	*u = *v;
	*v = t;
}

vertex_t find_scc(struct thread_context *context, vertex_t u);
vertex_t tighten(struct thread_context *context, vertex_t u);
bool join_tree(struct thread_context *context, vertex_t u, vertex_t v);

void inbound_graph_delete(struct inbound_graph *g);
struct inbound_graph inbound_graph_clone(struct inbound_graph *g);
struct inbound_graph random_graph(uint32_t num_nodes, uint32_t num_edges, cost_t max_cost);
struct inbound_graph adversarial_graph(uint32_t num_nodes);
struct edge_list edmonds_spanning_arborescence(struct inbound_graph in, vertex_t root, size_t num_threads);
int edmonds_worker(void* thread_context);
struct thread_context init_context(struct inbound_graph in, vertex_t root);
struct edge_list final_dfs(struct thread_context *context);


void make_binary_heap(struct inbound_edge*edges, uint32_t size);
struct full_edge leftist_heap_extract_min(struct thread_context *context, vertex_t u, cost_t *bias);
vertex_t leftist_heap_join(struct thread_context *context, vertex_t u, vertex_t v);
void leftist_heap_reconstruct(struct thread_context *context, vertex_t u);
void binary_heap_reconstruct(struct thread_context *context, vertex_t u, vertex_t v);

int main (int argc, char **argv) {
	TracyCZone(ctx, 1);
	assert(argc == 4);
	uint32_t n = atoi(argv[2]);
	uint32_t num_threads = atoi(argv[3]);

	struct timespec begin, end;
	clock_gettime(CLOCK_MONOTONIC, &begin);
	struct inbound_graph g;
	if (!strcmp(argv[1], "random")) {
		g = random_graph(n, n*(n-1), INF_WEIGHT / n);
	} else if (!strcmp(argv[1], "adversarial")) {
		g = adversarial_graph(n);
	} else {
		assert(false && "unknown CLI option");
	}
	clock_gettime(CLOCK_MONOTONIC, &end);
	printf("gen input: %gμs\n", (double)(end.tv_sec - begin.tv_sec) * 1e6 + (double)(end.tv_nsec - begin.tv_nsec) / 1e3);

	for (size_t runs = 0; runs < 10; ++runs) {
		for (size_t i = 0; i < num_threads; ++i) {
			struct inbound_graph gg = inbound_graph_clone(&g);
			clock_gettime(CLOCK_MONOTONIC, &begin);
			struct edge_list solution = edmonds_spanning_arborescence(gg, 0, i+1);
			clock_gettime(CLOCK_MONOTONIC, &end);
			printf("num_threads: %zu, time: %gus\n", i+1, (double)(end.tv_sec - begin.tv_sec) * 1e6 + (double)(end.tv_nsec - begin.tv_nsec) / 1e3);
			inbound_graph_delete(&gg);
			free(solution.edges);
		}
	}

	inbound_graph_delete(&g);
	TracyCZoneEnd(ctx);
}


struct edge_list edmonds_spanning_arborescence(struct inbound_graph in, vertex_t root, size_t num_threads) {
	TracyCZoneS(ctx, 10, 1);
	thrd_t *pool = calloc(num_threads, sizeof *pool);

	// the result graph starts out empty; while the input graph
	// uses inbound stars, the output graph has outbound stars, which are
	// more useful for the final traversal step.
	struct thread_context context = init_context(in, root);

	for (size_t i = 0; i < num_threads; ++i)
		thrd_create(&pool[i], &edmonds_worker, &context);

	for (size_t i = 0; i < num_threads; ++i)
		thrd_join(pool[i], NULL);

	free(pool);

	struct edge_list result = final_dfs(&context);
	TracyCZoneEnd(ctx);
	return result;
}

vertex_t next_job(struct thread_context *context) {
	vertex_t job = atomic_fetch_add(&context->next_job, 1);
	if (job >= context->input_graph.num_nodes) return NIL_VERTEX;
	return job;
}

struct thread_context init_context(struct inbound_graph in, vertex_t root) {
	uint32_t size = in.num_nodes;
	struct thread_context context = {
		.root = root,
		.next_job = 0,

		// scc partition
		.path = malloc(size * sizeof *context.path),
		.heap_root = malloc(size * sizeof *context.heap_root),
		.heap = malloc(size * sizeof *context.heap),
		.parent0 = malloc(size * sizeof *context.parent0),
		.count0 = malloc(size * sizeof *context.count0),

		// tree partition
		.parent1 = malloc(size * sizeof *context.parent1),
		.count1 = malloc(size * sizeof *context.count1),

		.lock = malloc(size * sizeof *context.lock),

		.input_graph = in,
		.solution = {
			.num_edges = 0,
			.edges = malloc(2UL * (size - 1) * sizeof *context.solution.edges),
		},
	};

	for (vertex_t i = 0; i < size; ++i) {
		context.path[i] = i;
		context.heap_root[i] = i;
		context.heap[i] = (struct leftist_heap_node) {
			.bias = 0,
			.left = NIL_VERTEX,
			.right = NIL_VERTEX, 
			.count = 1,
		};
		context.parent0[i] = i;
		context.count0[i] = 1;

		context.parent1[i] = i;
		context.count1[i] = 1;

		context.lock[i] = (atomic_flag)ATOMIC_FLAG_INIT;
	}

	return context;
}

static void add_solution_edge(struct thread_context *context, struct full_edge edge) {
	uint32_t index = atomic_fetch_add(&context->solution.num_edges, 1);
	context->solution.edges[index] = edge;
}

int edmonds_worker(void* thrd_context) {
	struct thread_context *context = thrd_context;
	vertex_t u, v;
	while ((u = next_job(context)) != NIL_VERTEX) {
		TracyCZone(ctx, 1);
		if (u == context->root) {
			TracyCZoneEnd(ctx);
			continue;
		}
		struct inbound_star *star = &context->input_graph.star[u];
		make_binary_heap(star->edges, star->num_edges);
		do {
			u = tighten(context, u);

			struct full_edge edge;
			cost_t bias;
			size_t retries = 0;
			do {
				// once we weed out the edges, it should not be possible to fail again.
				assert(retries < 100);
				++retries;
				if (retries == 100) {
					// after 100 consecutive edges extracted without leaving the scc,
					// try weeding out the internal edges all at once.
					// it should be cheaper to do this than to extract them one at a time.
					leftist_heap_reconstruct(context, u);
				}
				edge = leftist_heap_extract_min(context, u, &bias);
			} while (u == find_scc(context, edge.source));
			v = edge.source;

			add_solution_edge(context, edge);
			context->heap[context->heap_root[u]].bias -= edge.cost + bias;
			context->path[u] = v;
		} while (!join_tree(context, u, v));
		TracyCZoneEnd(ctx);
	}
	return 0;
}

void
__attribute__((noinline))
binary_heap_heapify(struct inbound_edge *edges, uint32_t size, uint32_t i) {
	TracyCZone(ctx, 1);
	while (true) {
		uint32_t left = i * 2 + 1, right = i * 2 + 2;
		if (left >= size) 
			break;
		uint32_t min = left;
		if (right < size && edges[right].cost < edges[left].cost)
			min = right;
		if (edges[i].cost <= edges[min].cost)
			break;
		swap_inbound_edges(&edges[i], &edges[min]);
		i = min;
	}
	TracyCZoneEnd(ctx);
}

void make_binary_heap(struct inbound_edge *edges, uint32_t size) {
	TracyCZone(ctx, 1);
	for (uint32_t i = size/2; i--; /**/)
		binary_heap_heapify(edges, size, i);
	TracyCZoneEnd(ctx);
}

struct full_edge binary_heap_extract_min(struct thread_context *context, vertex_t v) {
	TracyCZone(ctx, 1);
	struct inbound_star *star = &context->input_graph.star[v];
	if (star->num_edges) {
		struct full_edge result = {
			.source = star->edges[0].source,
			.cost = star->edges[0].cost,
			.dest = v,
		};
		star->num_edges -= 1;
		swap_inbound_edges(&star->edges[0], &star->edges[star->num_edges]);	
		binary_heap_heapify(star->edges, star->num_edges, 0);
		TracyCZoneEnd(ctx);
		return result;
	} else {
		TracyCZoneEnd(ctx);
		return (struct full_edge) {
			.source = context->root,
			.dest = v,
			.cost = INF_WEIGHT,
		};
	}
}

cost_t peek_cost(struct thread_context *context, vertex_t u, cost_t cumulated_bias) {
	TracyCZone(ctx, 1);
	cost_t result;
	if (u == NIL_VERTEX || context->input_graph.star[u].num_edges == 0) {
		result = INF_WEIGHT;
	} else {
		result = context->input_graph.star[u].edges[0].cost + context->heap[u].bias + cumulated_bias;
	}
	TracyCZoneEnd(ctx);
	return result;
}

void binary_heap_reconstruct(struct thread_context *context, vertex_t u, vertex_t v) {
	size_t n = 0;
	struct inbound_star *star = &context->input_graph.star[v];

	// copy only the edges that are not internal.
	for (size_t i = 0; i < star->num_edges; ++i) {
		if (u != find_scc(context, star->edges[i].source)) {
			star->edges[n++] = star->edges[i];
		}
	}
	star->num_edges = n;

	make_binary_heap(star->edges, star->num_edges);
}

void leftist_heap_reconstruct(struct thread_context *context, vertex_t u) {
	TracyCZone(ctx, 1);
	vertex_t root = NIL_VERTEX;
	vertex_t unprocessed[32] = {context->heap_root[u]};
	int num_unprocessed = 1;
	while (num_unprocessed--) {
		// detach the root of the subtree to be processed.
		vertex_t next = unprocessed[num_unprocessed]; 
		vertex_t left = context->heap[next].left;
		vertex_t right = context->heap[next].right;

		context->heap[next].left = NIL_VERTEX;
		context->heap[next].right = NIL_VERTEX;
		context->heap[next].count = 1;

		// leave the heavier subtree for recursion last, this guarantees that the number of elements in the work queue
		// is less than log2(tree size), so 32 is enough.
		if (left != NIL_VERTEX && right != NIL_VERTEX && context->heap[left].count < context->heap[right].count) {
			swap_vertices(&left, &right);
		}

		// append the two subtrees to the work queue
		cost_t bias = context->heap[next].bias;
		if (left != NIL_VERTEX) {
			context->heap[left].bias += bias;
			unprocessed[num_unprocessed++] = left;
			assert(num_unprocessed <= 32);
		}

		if (right != NIL_VERTEX) {
			context->heap[right].bias += bias;
			unprocessed[num_unprocessed++] = right;
			assert(num_unprocessed <= 32);
		}
		
		binary_heap_reconstruct(context, u, next);
		root = leftist_heap_join(context, root, next);
	}
	context->heap_root[u] = root;
	TracyCZoneEnd(ctx);
	return;
}

struct full_edge leftist_heap_extract_min(struct thread_context *context, vertex_t u, cost_t *bias) {
	TracyCZone(ctx, 1);
	vertex_t root = context->heap_root[u];
	struct leftist_heap_node *heap = context->heap;
	struct full_edge edge = binary_heap_extract_min(context, context->heap_root[u]);
	*bias = context->heap[root].bias;

	vertex_t left = heap[root].left;
	vertex_t right = heap[root].right;
	if (left != NIL_VERTEX)
		heap[left].bias += heap[root].bias;
	if (right != NIL_VERTEX)
		heap[right].bias += heap[root].bias;

	// detach root from the tree
	heap[root].left = NIL_VERTEX;
	heap[root].right = NIL_VERTEX;
	heap[root].count = 1;

	context->heap_root[u] = leftist_heap_join(context, left, leftist_heap_join(context, root, right));
	TracyCZoneEnd(ctx);
	return edge;
}

vertex_t leftist_heap_join(struct thread_context *context, vertex_t u, vertex_t v) {
	TracyCZone(ctx, 1);
	if (u == NIL_VERTEX) {
		TracyCZoneEnd(ctx);
		return v;
	}
	if (v == NIL_VERTEX) {
		TracyCZoneEnd(ctx);
		return u;
	}

	struct leftist_heap_node *heap = context->heap;
	vertex_t result = NIL_VERTEX;
	vertex_t *hole = &result; // helps make this procedure iterative, like tail recursion modulo CONS.
	// NOTE: I think this could be done in two variables without the reference variable bias_h,
	// but this way the variables have a clear interpretation.
	cost_t bias_u = 0; // the cumulated bias along the path from u's original root to u.
	cost_t bias_v = 0; // the cumulated bias along the path from v's original root to v.
	cost_t bias_h = 0; // the cumulated bias along the path to the hole in the newly constructed tree.
	cost_t cost_u = peek_cost(context, u, bias_u);
	cost_t cost_v = peek_cost(context, v, bias_v);

	while (true) {
		assert(u != NIL_VERTEX);
		assert(v != NIL_VERTEX);

		if (cost_u > cost_v) {
			swap_vertices(&u, &v);
			swap_costs(&bias_u, &bias_v);
			swap_costs(&cost_u, &cost_v);
		}

		*hole = u;

		uint32_t count_left = heap[u].left  == NIL_VERTEX ? 0 : heap[heap[u].left ].count;
		uint32_t count_right = heap[u].right == NIL_VERTEX ? 0 : heap[heap[u].right].count;

		vertex_t next_u = heap[u].right;
		bool swap = count_left < count_right + heap[v].count;
		if (swap) swap_vertices(&heap[u].right, &heap[u].left);
		vertex_t *next_hole = swap ? &heap[u].left : &heap[u].right;

		// adjust biases
		bias_u += heap[u].bias;
		bias_h += heap[u].bias;
		heap[u].count += heap[v].count;
		heap[u].bias += bias_u - bias_h; // we actually care about the previous value of this difference,
										 // but it's unchanged here because the two variables are
										 // incremented together.
										 // could be cleaned up with an intermediate value?
		if (next_u == NIL_VERTEX) {
			assert(v != NIL_VERTEX);
			*next_hole = v;
			heap[v].bias += bias_v - bias_h;
			TracyCZoneEnd(ctx);
			return result;
		}

		// next iteration
		u = next_u;
		hole = next_hole;
		cost_u = peek_cost(context, u, bias_u);
	}
}

static bool try_lock(atomic_flag *lock) {
	return !atomic_flag_test_and_set_explicit(lock, memory_order_acquire);
}

// SAFETY: must be locked by us!
static void unlock(atomic_flag *lock) {
	atomic_flag_clear_explicit(lock, memory_order_release);
}

uint32_t find_scc(struct thread_context *context, vertex_t u) {
	assert(u < context -> input_graph.num_nodes);
	for (vertex_t up; u != (up = context->parent0[u]); u = up);
	return u;
}

vertex_t tighten(struct thread_context *context, vertex_t u) {
	TracyCZone(ctx, 1);
	if (context->path[u] == NIL_VERTEX)  {
		TracyCZoneEnd(ctx);
		return u;
	}
	vertex_t v = find_scc(context, context->path[u]);
	while (u != v) {
		vertex_t w = context->path[v];
		if (context->count0[u] < context->count0[v]) swap_vertices(&u, &v);
		context->parent0[v] = u;
		context->count0[u] += context->count0[v];
		context->heap_root[u] = leftist_heap_join(context, context->heap_root[u], context->heap_root[v]);
		v = find_scc(context, w);
	}
	context->path[u] = u;
	TracyCZoneEnd(ctx);
	return u;
}

bool join_tree(struct thread_context *context, uint32_t u, uint32_t v) {
	TracyCZone(ctx, 1);
	while (true) {
		for (vertex_t up; u != (up = context->parent1[u]); u = up);
		for (vertex_t vp; v != (vp = context->parent1[v]); v = vp);

		// join didn't merge two separate partitions
		if (u == v) {
			TracyCZoneEnd(ctx);
			return false; 
		}

		// we always take the smaller lock first; this ensures deadlock freedom
		// because whenever thread t0 can't take one of the locks (l0), that lock
		// is taken by another thread t1, and, if t1 is also waiting on a lock (l1),
		// then l1 > l0 (since t1 acquired l0 before l1).
		// Such an ascending chain of locks can't be infinite, therefore some 
		// thread is able to take both locks and make progress.
		if (u > v) swap_vertices(&u, &v);
		if (try_lock(&context->lock[u])) {
			if (try_lock(&context->lock[v])) {
				if (context->parent1[u] == u && context->parent1[v] == v) {
					// no one made progress on the roots while we were waiting for locks.
					// other threads may start and complete find operations, but they 
					// won't be able to complete joins on these specific roots.
					// if they have any join operations in flight, they will either
					// resume after we release the locks or fail due to us having 
					// joined their partitions before them.

					// attach the lighter part below the heavier one.
					if (context->count1[u] < context->count1[v]) swap_vertices(&u, &v);
					context->parent1[v] = u;
					context->count1[u] += context->count1[v];

					unlock(&context->lock[u]);
					unlock(&context->lock[v]);

					// we joined separate partitions
					TracyCZoneEnd(ctx);
					return true; 
				}
				// Someone has made progress and one of the nodes we hold is not a root anymore;
				// we have to descend and retry.
				//
				// It's not safe to hold on to the other lock, because descending may swap the
				// node order and we would have the greater lock before the smaller one, which
				// violates the invariant we rely on to avoid deadlock.
				unlock(&context->lock[v]);
			}
			unlock(&context->lock[u]);
		}
	}
}

struct edge_list final_dfs(struct thread_context *context) {
	TracyCZone(ctx,1);
	// TODO: could multi-thread this as well, but probably no use.

	free(context->heap);

	// path is used in the DFS to backtrack
	memset(context->path, -1, context->input_graph.num_nodes * sizeof *context->path);
	free(context->heap_root);
	free(context->parent0);
	free(context->count0);

	free(context->parent1);
	free(context->count1);
	free(context->lock);

	// build up a graph from the nodes that we just got.
	// it would have been awkward to do this in the parallel part
	// because of race conditions on growing the vectors,
	// so we just do it once after.
	struct outbound_graph result = {
		.num_nodes = context->input_graph.num_nodes,
		.star = calloc(context->input_graph.num_nodes, sizeof *result.star),
	};

	for (size_t i = 0; i < context->solution.num_edges; ++i)
		result.star[context->solution.edges[i].source].num_edges++;

	for (vertex_t u = 0; u < context->input_graph.num_nodes; ++u) {
		result.star[u].edges = calloc(result.star[u].num_edges, sizeof *result.star[u].edges);
		result.star[u].num_edges = 0;
	}

	for (size_t i = 0; i < context->solution.num_edges; ++i) {
		struct full_edge edge = context->solution.edges[i];
		struct outbound_star *star = &result.star[edge.source];
		star->edges[star->num_edges++] = (struct outbound_edge) {
			.cost = edge.cost,
			.dest = edge.dest,
		};
	}

	context->solution.num_edges = 0;
	vertex_t u = context->root;
	vertex_t *path = context->path;
	struct edge_list solution = context->solution;
	while (u != NIL_VERTEX) {
		if (!result.star[u].num_edges) {
			u = context->path[u];
			continue;
		}
		--result.star[u].num_edges;
		struct outbound_edge edge = result.star[u].edges[result.star[u].num_edges];
		vertex_t v = edge.dest;
		if (path[v] == NIL_VERTEX) {
			path[v] = u;
			solution.edges[solution.num_edges++] = (struct full_edge) {
				.source = u,
				.dest = v,
				.cost = edge.cost,
			};
			u = v;
		}
	}

	for (vertex_t u = 0; u < result.num_nodes; ++u)
		free(result.star[u].edges);
	free(result.star);

	free(context->path);
	TracyCZoneEnd(ctx);
	return solution;
}

struct inbound_graph inbound_graph_clone(struct inbound_graph *g) {
	struct inbound_graph result = {
		.num_nodes = g->num_nodes,
		.star = malloc(g->num_nodes * sizeof *result.star),
	};

	for (uint32_t i = 0; i < result.num_nodes; ++i) {
		struct inbound_star *star = &g->star[i];
		struct inbound_edge *edges = malloc(star->num_edges * sizeof *edges);
		memcpy(edges, star->edges, star->num_edges * sizeof *edges);
		result.star[i].num_edges = star->num_edges;
		result.star[i].edges = edges;
	}

	return result;
}

void inbound_graph_delete(struct inbound_graph *g) {
	for (uint32_t i = 0; i < g->num_nodes; ++i) {
		free(g->star[i].edges);
	}
	free(g->star);
}

struct inbound_graph random_graph(uint32_t num_nodes, uint32_t num_edges, cost_t max_cost) {
	TracyCZoneS(ctx, 10, 1);
	struct inbound_graph result = {
		.num_nodes = num_nodes,
		.star = calloc(num_nodes, sizeof *result.star),
	};

	struct full_edge *edges = malloc((size_t)num_nodes * (num_nodes - 1) * sizeof *edges);

	for (uint32_t i = 0; i < num_nodes * (num_nodes - 1); ++i) {
		edges[i].source = i % num_nodes;
		edges[i].dest = (i / num_nodes) + (i / num_nodes >= i % num_nodes);
		edges[i].cost = rand() % max_cost;
	}

	for (uint32_t i = 0; i < num_edges; ++i) {
		uint32_t j = i + (rand() % (num_nodes * (num_nodes - 1) - i));
		swap_full_edges(&edges[i], &edges[j]);
		vertex_t u = edges[i].dest;

		uint32_t n = result.star[u].num_edges;
		if ((n&(n-1)) == 0) 
			result.star[u].edges = realloc(result.star[u].edges, (n ? 2 * n : 1) * sizeof *result.star[u].edges);

		result.star[u].edges[result.star[u].num_edges++] = (struct inbound_edge) {
			.source = edges[i].source,
			.cost = edges[i].cost,
		};
	}

	free(edges);
	TracyCZoneEnd(ctx);
	return result;
}


struct inbound_graph adversarial_graph(uint32_t num_nodes) {
	TracyCZoneS(ctx, 10, 1);
	struct inbound_graph result = {
		.num_nodes = num_nodes,
		.star = calloc(num_nodes, sizeof *result.star),
	};

	uint32_t loop_size = num_nodes / 2;
	uint32_t rest_size = num_nodes - loop_size - 1;

	for (vertex_t u = 1; u < loop_size + 1; ++u) {
		result.star[u].num_edges = rest_size + 1 + (u == 1);
		result.star[u].edges = malloc((rest_size + 1 + (u == 1)) * sizeof *result.star[u].edges);

		// each element of the loop has inbound edges from each of the rest,
		// with increasing costs. these are there to saturate the priority queues.
		for (uint32_t i = 0; i < rest_size ; ++i) {
			vertex_t v = loop_size + 1 + i;
			result.star[u].edges[i] = (struct inbound_edge) {
				.source = v,
				.cost = (cost_t)(u + ((u + i) % loop_size)),
			};
		}

		// the loop is closed with cheap edges
		result.star[u].edges[rest_size] = (struct inbound_edge) {
			.source = (u == 1) ? loop_size : u - 1,
			.cost = 0,
		};
		if (u == 1) {
			// a very expensive edge connects to the root
			result.star[u].edges[rest_size + 1] = (struct inbound_edge) {
				.source = 0,
				.cost = (cost_t)num_nodes,
			};
		}
	}

	for (vertex_t u = loop_size + 1; u < num_nodes; ++u) {
		result.star[u].num_edges = loop_size;
		result.star[u].edges = malloc(loop_size * sizeof *result.star[u].edges);

		// each element of the rest is connected to each of the loop.
		for (uint32_t i = 0; i < loop_size; ++i) {
			vertex_t v = 1 + i;
			result.star[u].edges[i] = (struct inbound_edge) {
				.source = v,
				.cost = (cost_t)((u + i) % loop_size),
			};
		}
	}

	assert(result.star[0].num_edges == 0);

	TracyCZoneEnd(ctx);
	return result;
}