2 * Edmonds-Karp algorithm for finding a maximum flow and minimum
3 * cut in a network. Almost identical to the Ford-Fulkerson
4 * algorithm, but apparently using breadth-first search to find the
5 * _shortest_ augmenting path is a good way to guarantee
6 * termination and ensure the time complexity is not dependent on
7 * the actual value of the maximum flow. I don't understand why
8 * that should be, but it's claimed on the Internet that it's been
9 * proved, and that's good enough for me. I prefer BFS to DFS
19 #include "puzzles.h" /* for snewn/sfree */
21 int maxflow_with_scratch(void *scratch, int nv, int source, int sink,
22 int ne, const int *edges, const int *backedges,
23 const int *capacity, int *flow, int *cut)
25 int *todo = (int *)scratch;
26 int *prev = todo + nv;
27 int *firstedge = todo + 2*nv;
28 int *firstbackedge = todo + 3*nv;
29 int i, j, head, tail, from, to;
33 * Scan the edges array to find the index of the first edge
37 for (i = 0; i < ne; i++)
38 while (j <= edges[2*i])
45 * Scan the backedges array to find the index of the first edge
49 for (i = 0; i < ne; i++)
50 while (j <= edges[2*backedges[i]+1])
51 firstbackedge[j++] = i;
53 firstbackedge[j++] = ne;
57 * Start the flow off at zero on every edge.
59 for (i = 0; i < ne; i++)
64 * Repeatedly look for an augmenting path, and follow it.
69 * Set up the prev array.
71 for (i = 0; i < nv; i++)
75 * Initialise the to-do list for BFS.
78 todo[tail++] = source;
81 * Now do the BFS loop.
83 while (head < tail && prev[sink] <= 0) {
87 * Try all the forward edges out of node `from'. For a
88 * forward edge to be valid, it must have flow
89 * currently less than its capacity.
91 for (i = firstedge[from]; i < ne && edges[2*i] == from; i++) {
93 if (to == source || prev[to] >= 0)
95 if (capacity[i] >= 0 && flow[i] >= capacity[i])
98 * This is a valid augmenting edge. Visit node `to'.
105 * Try all the backward edges into node `from'. For a
106 * backward edge to be valid, it must have flow
107 * currently greater than zero.
109 for (i = firstbackedge[from];
110 j = backedges[i], i < ne && edges[2*j+1]==from; i++) {
112 if (to == source || prev[to] >= 0)
117 * This is a valid augmenting edge. Visit node `to'.
125 * If prev[sink] is non-null, we have found an augmenting
128 if (prev[sink] >= 0) {
132 * Work backwards along the path figuring out the
133 * maximum flow we can add.
137 while (to != source) {
141 * Find the edge we're currently moving along.
148 * Determine the spare capacity of this edge.
151 spare = flow[i / 2]; /* backward edge */
152 else if (capacity[i / 2] >= 0)
153 spare = capacity[i / 2] - flow[i / 2]; /* forward edge */
155 spare = -1; /* unlimited forward edge */
159 if (max < 0 || (spare >= 0 && spare < max))
165 * Fail an assertion if max is still < 0, i.e. there is
166 * an entirely unlimited path from source to sink. Also
167 * max should not _be_ zero, because by construction
168 * this _should_ be an augmenting path.
173 * Now work backwards along the path again, this time
174 * actually adjusting the flow.
177 while (to != source) {
179 * Find the edge we're currently moving along.
189 flow[i / 2] -= max; /* backward edge */
191 flow[i / 2] += max; /* forward edge */
197 * And adjust the overall flow counter.
205 * If we reach here, we have failed to find an augmenting
206 * path, which means we're done. Output the `cut' array if
207 * required, and leave.
210 for (i = 0; i < nv; i++) {
211 if (i == source || prev[i] >= 0)
221 int maxflow_scratch_size(int nv)
223 return (nv * 4) * sizeof(int);
226 void maxflow_setup_backedges(int ne, const int *edges, int *backedges)
230 for (i = 0; i < ne; i++)
234 * We actually can't use the C qsort() function, because we'd
235 * need to pass `edges' as a context parameter to its
236 * comparator function. So instead I'm forced to implement my
237 * own sorting algorithm internally, which is a pest. I'll use
238 * heapsort, because I like it.
241 #define LESS(i,j) ( (edges[2*(i)+1] < edges[2*(j)+1]) || \
242 (edges[2*(i)+1] == edges[2*(j)+1] && \
243 edges[2*(i)] < edges[2*(j)]) )
244 #define PARENT(n) ( ((n)-1)/2 )
245 #define LCHILD(n) ( 2*(n)+1 )
246 #define RCHILD(n) ( 2*(n)+2 )
247 #define SWAP(i,j) do { int swaptmp = (i); (i) = (j); (j) = swaptmp; } while (0)
250 * Phase 1: build the heap. We want the _largest_ element at
258 * Swap element n with its parent repeatedly to preserve
266 if (LESS(backedges[p], backedges[i])) {
267 SWAP(backedges[p], backedges[i]);
275 * Phase 2: repeatedly remove the largest element and stick it
276 * at the top of the array.
280 * The largest element is at position 0. Put it at the top,
281 * and swap the arbitrary element from that position into
285 SWAP(backedges[0], backedges[n]);
288 * Now repeatedly move that arbitrary element down the heap
289 * by swapping it with the more suitable of its children.
299 break; /* we've hit bottom */
303 * Special case: there is only one child to check.
305 if (LESS(backedges[i], backedges[lc]))
306 SWAP(backedges[i], backedges[lc]);
308 /* _Now_ we've hit bottom. */
312 * The common case: there are two children and we
313 * must check them both.
315 if (LESS(backedges[i], backedges[lc]) ||
316 LESS(backedges[i], backedges[rc])) {
318 * Pick the more appropriate child to swap with
319 * (i.e. the one which would want to be the
320 * parent if one were above the other - as one
323 if (LESS(backedges[lc], backedges[rc])) {
324 SWAP(backedges[i], backedges[rc]);
327 SWAP(backedges[i], backedges[lc]);
331 /* This element is in the right place; we're done. */
346 int maxflow(int nv, int source, int sink,
347 int ne, const int *edges, const int *capacity,
356 * Allocate the space.
358 size = ne * sizeof(int) + maxflow_scratch_size(nv);
359 backedges = smalloc(size);
362 scratch = backedges + ne;
365 * Set up the backedges array.
367 maxflow_setup_backedges(ne, edges, backedges);
370 * Call the main function.
372 ret = maxflow_with_scratch(scratch, nv, source, sink, ne, edges,
373 backedges, capacity, flow, cut);
376 * Free the scratch space.
389 #define MAXVERTICES 128
390 #define ADDEDGE(i,j) do{edges[ne*2] = (i); edges[ne*2+1] = (j); ne++;}while(0)
392 int compare_edge(const void *av, const void *bv)
394 const int *a = (const int *)av;
395 const int *b = (const int *)bv;
399 else if (a[0] > b[0])
401 else if (a[1] < b[1])
403 else if (a[1] > b[1])
411 int edges[MAXEDGES*2], ne, nv;
412 int capacity[MAXEDGES], flow[MAXEDGES], cut[MAXVERTICES];
413 int source, sink, p, q, i, j, ret;
416 * Use this algorithm to find a maximal complete matching in a
427 for (i = 0; i < 5; i++) {
429 ADDEDGE(source, p+i);
431 for (i = 0; i < 5; i++) {
436 capacity[ne] = 1; ADDEDGE(p+0,q+0);
437 capacity[ne] = 1; ADDEDGE(p+1,q+0);
438 capacity[ne] = 1; ADDEDGE(p+1,q+1);
439 capacity[ne] = 1; ADDEDGE(p+2,q+1);
440 capacity[ne] = 1; ADDEDGE(p+2,q+2);
441 capacity[ne] = 1; ADDEDGE(p+3,q+2);
442 capacity[ne] = 1; ADDEDGE(p+3,q+3);
443 capacity[ne] = 1; ADDEDGE(p+4,q+3);
444 /* capacity[ne] = 1; ADDEDGE(p+2,q+4); */
445 qsort(edges, ne, 2*sizeof(int), compare_edge);
447 ret = maxflow(nv, source, sink, ne, edges, capacity, flow, cut);
449 printf("ret = %d\n", ret);
451 for (i = 0; i < ne; i++)
452 printf("flow %d: %d -> %d\n", flow[i], edges[2*i], edges[2*i+1]);
454 for (i = 0; i < nv; i++)
456 printf("difficult set includes %d\n", i);