1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
|
/* maxflow.c */
/* Written by Andrew Makhorin <mao@gnu.org>, October 2015. */
#include <math.h>
#include <glpk.h>
#include "maxflow.h"
#include "misc.h"
/***********************************************************************
* NAME
*
* max_flow - find max flow in undirected capacitated network
*
* SYNOPSIS
*
* #include "maxflow.h"
* int max_flow(int nn, int ne, const int beg[], const int end[],
* const int cap[], int s, int t, int x[]);
*
* DESCRIPTION
*
* This routine finds max flow 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.
*
* The parameter s is the number of a source node, and the parameter t
* is the number of a sink node, s != t.
*
* On exit the routine computes elementary flows thru edges and stores
* their values to locations x[1], ..., x[ne]. 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.
*
* RETURNS
*
* The routine returns the total maximum flow through the network. */
int max_flow(int nn, int ne, const int beg[/*1+ne*/],
const int end[/*1+ne*/], const int cap[/*1+ne*/], int s, int t,
int x[/*1+ne*/])
{ int k;
/* sanity checks */
xassert(nn >= 2);
xassert(ne >= 0);
xassert(1 <= s && s <= nn);
xassert(1 <= t && t <= nn);
xassert(s != t);
for (k = 1; k <= ne; k++)
{ xassert(1 <= beg[k] && beg[k] < end[k] && end[k] <= nn);
xassert(cap[k] > 0);
}
/* find max flow */
return max_flow_lp(nn, ne, beg, end, cap, s, t, x);
}
/***********************************************************************
* NAME
*
* max_flow_lp - find max flow with simplex method
*
* SYNOPSIS
*
* #include "maxflow.h"
* int max_flow_lp(int nn, int ne, const int beg[], const int end[],
* const int cap[], int s, int t, int x[]);
*
* DESCRIPTION
*
* This routine finds max flow in a given undirected network with the
* simplex method.
*
* Parameters of this routine have the same meaning as for the routine
* max_flow (see above).
*
* RETURNS
*
* The routine returns the total maximum flow through the network. */
int max_flow_lp(int nn, int ne, const int beg[/*1+ne*/],
const int end[/*1+ne*/], const int cap[/*1+ne*/], int s, int t,
int x[/*1+ne*/])
{ glp_prob *lp;
glp_smcp smcp;
int i, k, nz, flow, *rn, *cn;
double temp, *aa;
/* create LP problem instance */
lp = glp_create_prob();
/* create LP rows; i-th row is the conservation condition of the
* flow at i-th node, i = 1, ..., nn */
glp_add_rows(lp, nn);
for (i = 1; i <= nn; i++)
glp_set_row_bnds(lp, i, GLP_FX, 0.0, 0.0);
/* create LP columns; k-th column is the elementary flow thru
* k-th edge, k = 1, ..., ne; the last column with the number
* ne+1 is the total flow through the network, which goes along
* a dummy feedback edge from the sink to the source */
glp_add_cols(lp, ne+1);
for (k = 1; k <= ne; k++)
{ xassert(cap[k] > 0);
glp_set_col_bnds(lp, k, GLP_DB, -cap[k], +cap[k]);
}
glp_set_col_bnds(lp, ne+1, GLP_FR, 0.0, 0.0);
/* build the constraint matrix; structurally this matrix is the
* incidence matrix of the network, so each its column (including
* the last column for the dummy edge) has exactly two non-zero
* entries */
rn = xalloc(1+2*(ne+1), sizeof(int));
cn = xalloc(1+2*(ne+1), sizeof(int));
aa = xalloc(1+2*(ne+1), sizeof(double));
nz = 0;
for (k = 1; k <= ne; k++)
{ /* x[k] > 0 means the elementary flow thru k-th edge goes from
* node beg[k] to node end[k] */
nz++, rn[nz] = beg[k], cn[nz] = k, aa[nz] = -1.0;
nz++, rn[nz] = end[k], cn[nz] = k, aa[nz] = +1.0;
}
/* total flow thru the network goes from the sink to the source
* along the dummy feedback edge */
nz++, rn[nz] = t, cn[nz] = ne+1, aa[nz] = -1.0;
nz++, rn[nz] = s, cn[nz] = ne+1, aa[nz] = +1.0;
/* check the number of non-zero entries */
xassert(nz == 2*(ne+1));
/* load the constraint matrix into the LP problem object */
glp_load_matrix(lp, nz, rn, cn, aa);
xfree(rn);
xfree(cn);
xfree(aa);
/* objective function is the total flow through the network to
* be maximized */
glp_set_obj_dir(lp, GLP_MAX);
glp_set_obj_coef(lp, ne + 1, 1.0);
/* solve LP instance with the (primal) simplex method */
glp_term_out(0);
glp_adv_basis(lp, 0);
glp_term_out(1);
glp_init_smcp(&smcp);
smcp.msg_lev = GLP_MSG_ON;
smcp.out_dly = 5000;
xassert(glp_simplex(lp, &smcp) == 0);
xassert(glp_get_status(lp) == GLP_OPT);
/* obtain optimal elementary flows thru edges of the network */
/* (note that the constraint matrix is unimodular and the data
* are integral, so all elementary flows in basic solution should
* also be integral) */
for (k = 1; k <= ne; k++)
{ temp = glp_get_col_prim(lp, k);
x[k] = (int)floor(temp + .5);
xassert(fabs(x[k] - temp) <= 1e-6);
}
/* obtain the maximum flow thru the original network which is the
* flow thru the dummy feedback edge */
temp = glp_get_col_prim(lp, ne+1);
flow = (int)floor(temp + .5);
xassert(fabs(flow - temp) <= 1e-6);
/* delete LP problem instance */
glp_delete_prob(lp);
/* return to the calling program */
return flow;
}
/* eof */
|
