diff options
Diffstat (limited to 'glpk-5.0/examples/tsp/mincut.c')
| -rw-r--r-- | glpk-5.0/examples/tsp/mincut.c | 392 |
1 files changed, 392 insertions, 0 deletions
diff --git a/glpk-5.0/examples/tsp/mincut.c b/glpk-5.0/examples/tsp/mincut.c new file mode 100644 index 0000000..225fb7f --- /dev/null +++ b/glpk-5.0/examples/tsp/mincut.c @@ -0,0 +1,392 @@ +/* mincut.c */ + +/* Written by Andrew Makhorin <mao@gnu.org>, October 2015. */ + +#include <limits.h> + +#include "maxflow.h" +#include "mincut.h" +#include "misc.h" + +/*********************************************************************** +* NAME +* +* min_cut - find min cut in undirected capacitated network +* +* SYNOPSIS +* +* #include "mincut.h" +* int min_cut(int nn, int ne, const int beg[], const int end[], +* const cap[], int cut[]); +* +* DESCRIPTION +* +* This routine finds min cut in a given undirected network. +* +* The undirected capacitated network is specified by the parameters +* nn, ne, beg, end, and cap. The parameter nn specifies the number of +* vertices (nodes), nn >= 2, and the parameter ne specifies the number +* of edges, ne >= 0. The network edges are specified by triplets +* (beg[k], end[k], cap[k]) for k = 1, ..., ne, where beg[k] < end[k] +* are numbers of the first and second nodes of k-th edge, and +* cap[k] > 0 is a capacity of k-th edge. Loops and multiple edges are +* not allowed. +* +* Let V be the set of nodes of the network and let W be an arbitrary +* non-empty proper subset of V. A cut associated with the subset W is +* a subset of all the edges, one node of which belongs to W and other +* node belongs to V \ W. The capacity of a cut (W, V \ W) is the sum +* of the capacities of all the edges, which belong to the cut. Minimal +* cut is a cut, whose capacity is minimal. +* +* On exit the routine stores flags of nodes v[i], i = 1, ..., nn, to +* locations cut[i], where cut[i] = 1 means that v[i] belongs to W and +* cut[i] = 0 means that v[i] belongs to V \ W, where W corresponds to +* the minimal cut found. +* +* RETURNS +* +* The routine returns the capacity of the min cut found. */ + +int min_cut(int nn, int ne, const int beg[/*1+ne*/], + const int end[/*1+ne*/], const cap[/*1+ne*/], int cut[/*1+nn*/]) +{ int k; + /* sanity checks */ + xassert(nn >= 2); + xassert(ne >= 0); + for (k = 1; k <= ne; k++) + { xassert(1 <= beg[k] && beg[k] < end[k] && end[k] <= nn); + xassert(cap[k] > 0); + } + /* find min cut */ + return min_cut_sw(nn, ne, beg, end, cap, cut); +} + +/*********************************************************************** +* NAME +* +* min_st_cut - find min (s,t)-cut for known max flow +* +* SYNOPSIS +* +* #include "mincut.h" +* +* DESCRIPTION +* +* This routine finds min (s,t)-cut in a given undirected network that +* corresponds to a known max flow from s to t in the network. +* +* Parameters nn, ne, beg, end, and cap specify the network in the same +* way as for the routine min_cut (see above). +* +* Parameters s and t specify, resp., the number of the source node +* and the number of the sink node, s != t, for which the min (s,t)-cut +* has to be found. +* +* Parameter x specifies the known max flow from s to t in the network, +* where locations x[1], ..., x[ne] contains elementary flow thru edges +* of the network (positive value of x[k] means that the elementary +* flow goes from node beg[k] to node end[k], and negative value means +* that the flow goes in opposite direction). +* +* This routine splits the set of nodes V of the network into two +* non-empty subsets V(s) and V(t) = V \ V(s), where the source node s +* belongs to V(s), the sink node t belongs to V(t), and all edges, one +* node of which belongs to V(s) and other one belongs to V(t), are +* saturated (that is, x[k] = +cap[k] or x[k] = -cap[k]). +* +* On exit the routine stores flags of the nodes v[i], i = 1, ..., nn, +* to locations cut[i], where cut[i] = 1 means that v[i] belongs to V(s) +* and cut[i] = 0 means that v[i] belongs to V(t) = V \ V(s). +* +* RETURNS +* +* The routine returns the capacity of min (s,t)-cut, which is the sum +* of the capacities of all the edges, which belong to the cut. (Note +* that due to theorem by Ford and Fulkerson this value is always equal +* to corresponding max flow.) +* +* ALGORITHM +* +* To determine the set V(s) the routine simply finds all nodes, which +* can be reached from the source node s via non-saturated edges. The +* set V(t) is determined as the complement V \ V(s). */ + +int min_st_cut(int nn, int ne, const int beg[/*1+ne*/], + const int end[/*1+ne*/], const int cap[/*1+ne*/], int s, int t, + const int x[/*1+ne*/], int cut[/*1+nn*/]) +{ int i, j, k, p, q, temp, *head1, *next1, *head2, *next2, *list; + /* head1[i] points to the first edge with beg[k] = i + * next1[k] points to the next edge with the same beg[k] + * head2[i] points to the first edge with end[k] = i + * next2[k] points to the next edge with the same end[k] */ + head1 = xalloc(1+nn, sizeof(int)); + head2 = xalloc(1+nn, sizeof(int)); + next1 = xalloc(1+ne, sizeof(int)); + next2 = xalloc(1+ne, sizeof(int)); + for (i = 1; i <= nn; i++) + head1[i] = head2[i] = 0; + for (k = 1; k <= ne; k++) + { i = beg[k], next1[k] = head1[i], head1[i] = k; + j = end[k], next2[k] = head2[j], head2[j] = k; + } + /* on constructing the set V(s) list[1], ..., list[p-1] contain + * nodes, which can be reached from source node and have been + * visited, and list[p], ..., list[q] contain nodes, which can be + * reached from source node but havn't been visited yet */ + list = xalloc(1+nn, sizeof(int)); + for (i = 1; i <= nn; i++) + cut[i] = 0; + p = q = 1, list[1] = s, cut[s] = 1; + while (p <= q) + { /* pick next node, which is reachable from the source node and + * has not visited yet, and visit it */ + i = list[p++]; + /* walk through edges with beg[k] = i */ + for (k = head1[i]; k != 0; k = next1[k]) + { j = end[k]; + xassert(beg[k] == i); + /* from v[i] we can reach v[j], if the elementary flow from + * v[i] to v[j] is non-saturated */ + if (cut[j] == 0 && x[k] < +cap[k]) + list[++q] = j, cut[j] = 1; + } + /* walk through edges with end[k] = i */ + for (k = head2[i]; k != 0; k = next2[k]) + { j = beg[k]; + xassert(end[k] == i); + /* from v[i] we can reach v[j], if the elementary flow from + * v[i] to v[j] is non-saturated */ + if (cut[j] == 0 && x[k] > -cap[k]) + list[++q] = j, cut[j] = 1; + } + } + /* sink cannot belong to V(s) */ + xassert(!cut[t]); + /* free working arrays */ + xfree(head1); + xfree(head2); + xfree(next1); + xfree(next2); + xfree(list); + /* compute capacity of the minimal (s,t)-cut found */ + temp = 0; + for (k = 1; k <= ne; k++) + { i = beg[k], j = end[k]; + if (cut[i] && !cut[j] || !cut[i] && cut[j]) + temp += cap[k]; + } + /* return to the calling program */ + return temp; +} + +/*********************************************************************** +* NAME +* +* min_cut_sw - find min cut with Stoer and Wagner algorithm +* +* SYNOPSIS +* +* #include "mincut.h" +* int min_cut_sw(int nn, int ne, const int beg[], const int end[], +* const cap[], int cut[]); +* +* DESCRIPTION +* +* This routine find min cut in a given undirected network with the +* algorithm proposed by Stoer and Wagner (see references below). +* +* Parameters of this routine have the same meaning as for the routine +* min_cut (see above). +* +* RETURNS +* +* The routine returns the capacity of the min cut found. +* +* ALGORITHM +* +* The basic idea of Stoer&Wagner algorithm is the following. Let G be +* a capacitated network, and G(s,t) be a network, in which the nodes s +* and t are merged into one new node, loops are deleted, but multiple +* edges are retained. It is obvious that a minimum cut in G is the +* minimum of two quantities: the minimum cut in G(s,t) and a minimum +* cut that separates s and t. This allows to find a minimum cut in the +* original network solving at most nn max flow problems. +* +* REFERENCES +* +* M. Stoer, F. Wagner. A Simple Min Cut Algorithm. Algorithms, ESA'94 +* LNCS 855 (1994), pp. 141-47. +* +* J. Cheriyan, R. Ravi. Approximation Algorithms for Network Problems. +* Univ. of Waterloo (1998), p. 147. */ + +int min_cut_sw(int nn, int ne, const int beg[/*1+ne*/], + const int end[/*1+ne*/], const cap[/*1+ne*/], int cut[/*1+nn*/]) +{ int i, j, k, min_cut, flow, temp, *head1, *next1, *head2, *next2; + int I, J, K, S, T, DEG, NV, NE, *HEAD, *NEXT, *NUMB, *BEG, *END, + *CAP, *X, *ADJ, *SUM, *CUT; + /* head1[i] points to the first edge with beg[k] = i + * next1[k] points to the next edge with the same beg[k] + * head2[i] points to the first edge with end[k] = i + * next2[k] points to the next edge with the same end[k] */ + head1 = xalloc(1+nn, sizeof(int)); + head2 = xalloc(1+nn, sizeof(int)); + next1 = xalloc(1+ne, sizeof(int)); + next2 = xalloc(1+ne, sizeof(int)); + for (i = 1; i <= nn; i++) + head1[i] = head2[i] = 0; + for (k = 1; k <= ne; k++) + { i = beg[k], next1[k] = head1[i], head1[i] = k; + j = end[k], next2[k] = head2[j], head2[j] = k; + } + /* an auxiliary network used in the algorithm is resulted from + * the original network by merging some nodes into one supernode; + * all variables and arrays related to this auxiliary network are + * denoted in CAPS */ + /* HEAD[I] points to the first node of the original network that + * belongs to the I-th supernode + * NEXT[i] points to the next node of the original network that + * belongs to the same supernode as the i-th node + * NUMB[i] is a supernode, which the i-th node belongs to */ + /* initially the auxiliary network is equivalent to the original + * network, i.e. each supernode consists of one node */ + NV = nn; + HEAD = xalloc(1+nn, sizeof(int)); + NEXT = xalloc(1+nn, sizeof(int)); + NUMB = xalloc(1+nn, sizeof(int)); + for (i = 1; i <= nn; i++) + HEAD[i] = i, NEXT[i] = 0, NUMB[i] = i; + /* number of edges in the auxiliary network is never greater than + * in the original one */ + BEG = xalloc(1+ne, sizeof(int)); + END = xalloc(1+ne, sizeof(int)); + CAP = xalloc(1+ne, sizeof(int)); + X = xalloc(1+ne, sizeof(int)); + /* allocate some auxiliary arrays */ + ADJ = xalloc(1+nn, sizeof(int)); + SUM = xalloc(1+nn, sizeof(int)); + CUT = xalloc(1+nn, sizeof(int)); + /* currently no min cut is found so far */ + min_cut = INT_MAX; + /* main loop starts here */ + while (NV > 1) + { /* build the set of edges of the auxiliary network */ + NE = 0; + /* multiple edges are not allowed in the max flow algorithm, + * so we can replace each multiple edge, which is the result + * of merging nodes into supernodes, by a single edge, whose + * capacity is the sum of capacities of particular edges; + * these summary capacities will be stored in the array SUM */ + for (I = 1; I <= NV; I++) + SUM[I] = 0.0; + for (I = 1; I <= NV; I++) + { /* DEG is number of single edges, which connects I-th + * supernode and some J-th supernode, where I < J */ + DEG = 0; + /* walk thru nodes that belong to I-th supernode */ + for (i = HEAD[I]; i != 0; i = NEXT[i]) + { /* i-th node belongs to I-th supernode */ + /* walk through edges with beg[k] = i */ + for (k = head1[i]; k != 0; k = next1[k]) + { j = end[k]; + /* j-th node belongs to J-th supernode */ + J = NUMB[j]; + /* ignore loops and edges with I > J */ + if (I >= J) + continue; + /* add an edge that connects I-th and J-th supernodes + * (if not added yet) */ + if (SUM[J] == 0.0) + ADJ[++DEG] = J; + /* sum up the capacity of the original edge */ + xassert(cap[k] > 0.0); + SUM[J] += cap[k]; + } + /* walk through edges with end[k] = i */ + for (k = head2[i]; k != 0; k = next2[k]) + { j = beg[k]; + /* j-th node belongs to J-th supernode */ + J = NUMB[j]; + /* ignore loops and edges with I > J */ + if (I >= J) + continue; + /* add an edge that connects I-th and J-th supernodes + * (if not added yet) */ + if (SUM[J] == 0.0) + ADJ[++DEG] = J; + /* sum up the capacity of the original edge */ + xassert(cap[k] > 0.0); + SUM[J] += cap[k]; + } + } + /* add single edges connecting I-th supernode with other + * supernodes to the auxiliary network; restore the array + * SUM for subsequent use */ + for (K = 1; K <= DEG; K++) + { NE++; + xassert(NE <= ne); + J = ADJ[K]; + BEG[NE] = I, END[NE] = J, CAP[NE] = SUM[J]; + SUM[J] = 0.0; + } + } + /* choose two arbitrary supernodes of the auxiliary network, + * one of which is the source and other is the sink */ + S = 1, T = NV; + /* determine max flow from S to T */ + flow = max_flow(NV, NE, BEG, END, CAP, S, T, X); + /* if the min cut that separates supernodes S and T is less + * than the currently known, remember it */ + if (min_cut > flow) + { min_cut = flow; + /* find min (s,t)-cut in the auxiliary network */ + temp = min_st_cut(NV, NE, BEG, END, CAP, S, T, X, CUT); + /* (Ford and Fulkerson insist on this) */ + xassert(flow == temp); + /* build corresponding min cut in the original network */ + for (i = 1; i <= nn; i++) cut[i] = CUT[NUMB[i]]; + /* if the min cut capacity is zero (i.e. the network has + * unconnected components), the search can be prematurely + * terminated */ + if (min_cut == 0) + break; + } + /* now merge all nodes of the original network, which belong + * to the supernodes S and T, into one new supernode; this is + * attained by carrying all nodes from T to S (for the sake of + * convenience T should be the last supernode) */ + xassert(T == NV); + /* assign new references to nodes from T */ + for (i = HEAD[T]; i != 0; i = NEXT[i]) + NUMB[i] = S; + /* find last entry in the node list of S */ + i = HEAD[S]; + xassert(i != 0); + while (NEXT[i] != 0) + i = NEXT[i]; + /* and attach to it the node list of T */ + NEXT[i] = HEAD[T]; + /* decrease number of nodes in the auxiliary network */ + NV--; + } + /* free working arrays */ + xfree(HEAD); + xfree(NEXT); + xfree(NUMB); + xfree(BEG); + xfree(END); + xfree(CAP); + xfree(X); + xfree(ADJ); + xfree(SUM); + xfree(CUT); + xfree(head1); + xfree(head2); + xfree(next1); + xfree(next2); + /* return to the calling program */ + return min_cut; +} + +/* eof */ |