1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpnet07.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,222 @@
1.4 +/* glpnet07.c (Ford-Fulkerson algorithm) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpenv.h"
1.29 +#include "glpnet.h"
1.30 +
1.31 +/***********************************************************************
1.32 +* NAME
1.33 +*
1.34 +* ffalg - Ford-Fulkerson algorithm
1.35 +*
1.36 +* SYNOPSIS
1.37 +*
1.38 +* #include "glpnet.h"
1.39 +* void ffalg(int nv, int na, const int tail[], const int head[],
1.40 +* int s, int t, const int cap[], int x[], char cut[]);
1.41 +*
1.42 +* DESCRIPTION
1.43 +*
1.44 +* The routine ffalg implements the Ford-Fulkerson algorithm to find a
1.45 +* maximal flow in the specified flow network.
1.46 +*
1.47 +* INPUT PARAMETERS
1.48 +*
1.49 +* nv is the number of nodes, nv >= 2.
1.50 +*
1.51 +* na is the number of arcs, na >= 0.
1.52 +*
1.53 +* tail[a], a = 1,...,na, is the index of tail node of arc a.
1.54 +*
1.55 +* head[a], a = 1,...,na, is the index of head node of arc a.
1.56 +*
1.57 +* s is the source node index, 1 <= s <= nv.
1.58 +*
1.59 +* t is the sink node index, 1 <= t <= nv, t != s.
1.60 +*
1.61 +* cap[a], a = 1,...,na, is the capacity of arc a, cap[a] >= 0.
1.62 +*
1.63 +* NOTE: Multiple arcs are allowed, but self-loops are not allowed.
1.64 +*
1.65 +* OUTPUT PARAMETERS
1.66 +*
1.67 +* x[a], a = 1,...,na, is optimal value of the flow through arc a.
1.68 +*
1.69 +* cut[i], i = 1,...,nv, is 1 if node i is labelled, and 0 otherwise.
1.70 +* The set of arcs, whose one endpoint is labelled and other is not,
1.71 +* defines the minimal cut corresponding to the maximal flow found.
1.72 +* If the parameter cut is NULL, the cut information are not stored.
1.73 +*
1.74 +* REFERENCES
1.75 +*
1.76 +* L.R.Ford, Jr., and D.R.Fulkerson, "Flows in Networks," The RAND
1.77 +* Corp., Report R-375-PR (August 1962), Chap. I "Static Maximal Flow,"
1.78 +* pp.30-33. */
1.79 +
1.80 +void ffalg(int nv, int na, const int tail[], const int head[],
1.81 + int s, int t, const int cap[], int x[], char cut[])
1.82 +{ int a, delta, i, j, k, pos1, pos2, temp,
1.83 + *ptr, *arc, *link, *list;
1.84 + /* sanity checks */
1.85 + xassert(nv >= 2);
1.86 + xassert(na >= 0);
1.87 + xassert(1 <= s && s <= nv);
1.88 + xassert(1 <= t && t <= nv);
1.89 + xassert(s != t);
1.90 + for (a = 1; a <= na; a++)
1.91 + { i = tail[a], j = head[a];
1.92 + xassert(1 <= i && i <= nv);
1.93 + xassert(1 <= j && j <= nv);
1.94 + xassert(i != j);
1.95 + xassert(cap[a] >= 0);
1.96 + }
1.97 + /* allocate working arrays */
1.98 + ptr = xcalloc(1+nv+1, sizeof(int));
1.99 + arc = xcalloc(1+na+na, sizeof(int));
1.100 + link = xcalloc(1+nv, sizeof(int));
1.101 + list = xcalloc(1+nv, sizeof(int));
1.102 + /* ptr[i] := (degree of node i) */
1.103 + for (i = 1; i <= nv; i++)
1.104 + ptr[i] = 0;
1.105 + for (a = 1; a <= na; a++)
1.106 + { ptr[tail[a]]++;
1.107 + ptr[head[a]]++;
1.108 + }
1.109 + /* initialize arc pointers */
1.110 + ptr[1]++;
1.111 + for (i = 1; i < nv; i++)
1.112 + ptr[i+1] += ptr[i];
1.113 + ptr[nv+1] = ptr[nv];
1.114 + /* build arc lists */
1.115 + for (a = 1; a <= na; a++)
1.116 + { arc[--ptr[tail[a]]] = a;
1.117 + arc[--ptr[head[a]]] = a;
1.118 + }
1.119 + xassert(ptr[1] == 1);
1.120 + xassert(ptr[nv+1] == na+na+1);
1.121 + /* now the indices of arcs incident to node i are stored in
1.122 + locations arc[ptr[i]], arc[ptr[i]+1], ..., arc[ptr[i+1]-1] */
1.123 + /* initialize arc flows */
1.124 + for (a = 1; a <= na; a++)
1.125 + x[a] = 0;
1.126 +loop: /* main loop starts here */
1.127 + /* build augmenting tree rooted at s */
1.128 + /* link[i] = 0 means that node i is not labelled yet;
1.129 + link[i] = a means that arc a immediately precedes node i */
1.130 + /* initially node s is labelled as the root */
1.131 + for (i = 1; i <= nv; i++)
1.132 + link[i] = 0;
1.133 + link[s] = -1, list[1] = s, pos1 = pos2 = 1;
1.134 + /* breadth first search */
1.135 + while (pos1 <= pos2)
1.136 + { /* dequeue node i */
1.137 + i = list[pos1++];
1.138 + /* consider all arcs incident to node i */
1.139 + for (k = ptr[i]; k < ptr[i+1]; k++)
1.140 + { a = arc[k];
1.141 + if (tail[a] == i)
1.142 + { /* a = i->j is a forward arc from s to t */
1.143 + j = head[a];
1.144 + /* if node j has been labelled, skip the arc */
1.145 + if (link[j] != 0) continue;
1.146 + /* if the arc does not allow increasing the flow through
1.147 + it, skip the arc */
1.148 + if (x[a] == cap[a]) continue;
1.149 + }
1.150 + else if (head[a] == i)
1.151 + { /* a = i<-j is a backward arc from s to t */
1.152 + j = tail[a];
1.153 + /* if node j has been labelled, skip the arc */
1.154 + if (link[j] != 0) continue;
1.155 + /* if the arc does not allow decreasing the flow through
1.156 + it, skip the arc */
1.157 + if (x[a] == 0) continue;
1.158 + }
1.159 + else
1.160 + xassert(a != a);
1.161 + /* label node j and enqueue it */
1.162 + link[j] = a, list[++pos2] = j;
1.163 + /* check for breakthrough */
1.164 + if (j == t) goto brkt;
1.165 + }
1.166 + }
1.167 + /* NONBREAKTHROUGH */
1.168 + /* no augmenting path exists; current flow is maximal */
1.169 + /* store minimal cut information, if necessary */
1.170 + if (cut != NULL)
1.171 + { for (i = 1; i <= nv; i++)
1.172 + cut[i] = (char)(link[i] != 0);
1.173 + }
1.174 + goto done;
1.175 +brkt: /* BREAKTHROUGH */
1.176 + /* walk through arcs of the augmenting path (s, ..., t) found in
1.177 + the reverse order and determine maximal change of the flow */
1.178 + delta = 0;
1.179 + for (j = t; j != s; j = i)
1.180 + { /* arc a immediately precedes node j in the path */
1.181 + a = link[j];
1.182 + if (head[a] == j)
1.183 + { /* a = i->j is a forward arc of the cycle */
1.184 + i = tail[a];
1.185 + /* x[a] may be increased until its upper bound */
1.186 + temp = cap[a] - x[a];
1.187 + }
1.188 + else if (tail[a] == j)
1.189 + { /* a = i<-j is a backward arc of the cycle */
1.190 + i = head[a];
1.191 + /* x[a] may be decreased until its lower bound */
1.192 + temp = x[a];
1.193 + }
1.194 + else
1.195 + xassert(a != a);
1.196 + if (delta == 0 || delta > temp) delta = temp;
1.197 + }
1.198 + xassert(delta > 0);
1.199 + /* increase the flow along the path */
1.200 + for (j = t; j != s; j = i)
1.201 + { /* arc a immediately precedes node j in the path */
1.202 + a = link[j];
1.203 + if (head[a] == j)
1.204 + { /* a = i->j is a forward arc of the cycle */
1.205 + i = tail[a];
1.206 + x[a] += delta;
1.207 + }
1.208 + else if (tail[a] == j)
1.209 + { /* a = i<-j is a backward arc of the cycle */
1.210 + i = head[a];
1.211 + x[a] -= delta;
1.212 + }
1.213 + else
1.214 + xassert(a != a);
1.215 + }
1.216 + goto loop;
1.217 +done: /* free working arrays */
1.218 + xfree(ptr);
1.219 + xfree(arc);
1.220 + xfree(link);
1.221 + xfree(list);
1.222 + return;
1.223 +}
1.224 +
1.225 +/* eof */