alpar@1: /* glpnet06.c (out-of-kilter algorithm) */
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1: *
alpar@1: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1: *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1: *  E-mail: <mao@gnu.org>.
alpar@1: *
alpar@1: *  GLPK is free software: you can redistribute it and/or modify it
alpar@1: *  under the terms of the GNU General Public License as published by
alpar@1: *  the Free Software Foundation, either version 3 of the License, or
alpar@1: *  (at your option) any later version.
alpar@1: *
alpar@1: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1: *  License for more details.
alpar@1: *
alpar@1: *  You should have received a copy of the GNU General Public License
alpar@1: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1: ***********************************************************************/
alpar@1: 
alpar@1: #include "glpenv.h"
alpar@1: #include "glpnet.h"
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  okalg - out-of-kilter algorithm
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpnet.h"
alpar@1: *  int okalg(int nv, int na, const int tail[], const int head[],
alpar@1: *     const int low[], const int cap[], const int cost[], int x[],
alpar@1: *     int pi[]);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine okalg implements the out-of-kilter algorithm to find a
alpar@1: *  minimal-cost circulation in the specified flow network.
alpar@1: *
alpar@1: *  INPUT PARAMETERS
alpar@1: *
alpar@1: *  nv is the number of nodes, nv >= 0.
alpar@1: *
alpar@1: *  na is the number of arcs, na >= 0.
alpar@1: *
alpar@1: *  tail[a], a = 1,...,na, is the index of tail node of arc a.
alpar@1: *
alpar@1: *  head[a], a = 1,...,na, is the index of head node of arc a.
alpar@1: *
alpar@1: *  low[a], a = 1,...,na, is an lower bound to the flow through arc a.
alpar@1: *
alpar@1: *  cap[a], a = 1,...,na, is an upper bound to the flow through arc a,
alpar@1: *  which is the capacity of the arc.
alpar@1: *
alpar@1: *  cost[a], a = 1,...,na, is a per-unit cost of the flow through arc a.
alpar@1: *
alpar@1: *  NOTES
alpar@1: *
alpar@1: *  1. Multiple arcs are allowed, but self-loops are not allowed.
alpar@1: *
alpar@1: *  2. It is required that 0 <= low[a] <= cap[a] for all arcs.
alpar@1: *
alpar@1: *  3. Arc costs may have any sign.
alpar@1: *
alpar@1: *  OUTPUT PARAMETERS
alpar@1: *
alpar@1: *  x[a], a = 1,...,na, is optimal value of the flow through arc a.
alpar@1: *
alpar@1: *  pi[i], i = 1,...,nv, is Lagrange multiplier for flow conservation
alpar@1: *  equality constraint corresponding to node i (the node potential).
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  0  optimal circulation found;
alpar@1: *
alpar@1: *  1  there is no feasible circulation;
alpar@1: *
alpar@1: *  2  integer overflow occured;
alpar@1: *
alpar@1: *  3  optimality test failed (logic error).
alpar@1: *
alpar@1: *  REFERENCES
alpar@1: *
alpar@1: *  L.R.Ford, Jr., and D.R.Fulkerson, "Flows in Networks," The RAND
alpar@1: *  Corp., Report R-375-PR (August 1962), Chap. III "Minimal Cost Flow
alpar@1: *  Problems," pp.113-26. */
alpar@1: 
alpar@1: static int overflow(int u, int v)
alpar@1: {     /* check for integer overflow on computing u + v */
alpar@1:       if (u > 0 && v > 0 && u + v < 0) return 1;
alpar@1:       if (u < 0 && v < 0 && u + v > 0) return 1;
alpar@1:       return 0;
alpar@1: }
alpar@1: 
alpar@1: int okalg(int nv, int na, const int tail[], const int head[],
alpar@1:       const int low[], const int cap[], const int cost[], int x[],
alpar@1:       int pi[])
alpar@1: {     int a, aok, delta, i, j, k, lambda, pos1, pos2, s, t, temp, ret,
alpar@1:          *ptr, *arc, *link, *list;
alpar@1:       /* sanity checks */
alpar@1:       xassert(nv >= 0);
alpar@1:       xassert(na >= 0);
alpar@1:       for (a = 1; a <= na; a++)
alpar@1:       {  i = tail[a], j = head[a];
alpar@1:          xassert(1 <= i && i <= nv);
alpar@1:          xassert(1 <= j && j <= nv);
alpar@1:          xassert(i != j);
alpar@1:          xassert(0 <= low[a] && low[a] <= cap[a]);
alpar@1:       }
alpar@1:       /* allocate working arrays */
alpar@1:       ptr = xcalloc(1+nv+1, sizeof(int));
alpar@1:       arc = xcalloc(1+na+na, sizeof(int));
alpar@1:       link = xcalloc(1+nv, sizeof(int));
alpar@1:       list = xcalloc(1+nv, sizeof(int));
alpar@1:       /* ptr[i] := (degree of node i) */
alpar@1:       for (i = 1; i <= nv; i++)
alpar@1:          ptr[i] = 0;
alpar@1:       for (a = 1; a <= na; a++)
alpar@1:       {  ptr[tail[a]]++;
alpar@1:          ptr[head[a]]++;
alpar@1:       }
alpar@1:       /* initialize arc pointers */
alpar@1:       ptr[1]++;
alpar@1:       for (i = 1; i < nv; i++)
alpar@1:          ptr[i+1] += ptr[i];
alpar@1:       ptr[nv+1] = ptr[nv];
alpar@1:       /* build arc lists */
alpar@1:       for (a = 1; a <= na; a++)
alpar@1:       {  arc[--ptr[tail[a]]] = a;
alpar@1:          arc[--ptr[head[a]]] = a;
alpar@1:       }
alpar@1:       xassert(ptr[1] == 1);
alpar@1:       xassert(ptr[nv+1] == na+na+1);
alpar@1:       /* now the indices of arcs incident to node i are stored in
alpar@1:          locations arc[ptr[i]], arc[ptr[i]+1], ..., arc[ptr[i+1]-1] */
alpar@1:       /* initialize arc flows and node potentials */
alpar@1:       for (a = 1; a <= na; a++)
alpar@1:          x[a] = 0;
alpar@1:       for (i = 1; i <= nv; i++)
alpar@1:          pi[i] = 0;
alpar@1: loop: /* main loop starts here */
alpar@1:       /* find out-of-kilter arc */
alpar@1:       aok = 0;
alpar@1:       for (a = 1; a <= na; a++)
alpar@1:       {  i = tail[a], j = head[a];
alpar@1:          if (overflow(cost[a], pi[i] - pi[j]))
alpar@1:          {  ret = 2;
alpar@1:             goto done;
alpar@1:          }
alpar@1:          lambda = cost[a] + (pi[i] - pi[j]);
alpar@1:          if (x[a] < low[a] || lambda < 0 && x[a] < cap[a])
alpar@1:          {  /* arc a = i->j is out of kilter, and we need to increase
alpar@1:                the flow through this arc */
alpar@1:             aok = a, s = j, t = i;
alpar@1:             break;
alpar@1:          }
alpar@1:          if (x[a] > cap[a] || lambda > 0 && x[a] > low[a])
alpar@1:          {  /* arc a = i->j is out of kilter, and we need to decrease
alpar@1:                the flow through this arc */
alpar@1:             aok = a, s = i, t = j;
alpar@1:             break;
alpar@1:          }
alpar@1:       }
alpar@1:       if (aok == 0)
alpar@1:       {  /* all arcs are in kilter */
alpar@1:          /* check for feasibility */
alpar@1:          for (a = 1; a <= na; a++)
alpar@1:          {  if (!(low[a] <= x[a] && x[a] <= cap[a]))
alpar@1:             {  ret = 3;
alpar@1:                goto done;
alpar@1:             }
alpar@1:          }
alpar@1:          for (i = 1; i <= nv; i++)
alpar@1:          {  temp = 0;
alpar@1:             for (k = ptr[i]; k < ptr[i+1]; k++)
alpar@1:             {  a = arc[k];
alpar@1:                if (tail[a] == i)
alpar@1:                {  /* a is outgoing arc */
alpar@1:                   temp += x[a];
alpar@1:                }
alpar@1:                else if (head[a] == i)
alpar@1:                {  /* a is incoming arc */
alpar@1:                   temp -= x[a];
alpar@1:                }
alpar@1:                else
alpar@1:                   xassert(a != a);
alpar@1:             }
alpar@1:             if (temp != 0)
alpar@1:             {  ret = 3;
alpar@1:                goto done;
alpar@1:             }
alpar@1:          }
alpar@1:          /* check for optimality */
alpar@1:          for (a = 1; a <= na; a++)
alpar@1:          {  i = tail[a], j = head[a];
alpar@1:             lambda = cost[a] + (pi[i] - pi[j]);
alpar@1:             if (lambda > 0 && x[a] != low[a] ||
alpar@1:                 lambda < 0 && x[a] != cap[a])
alpar@1:             {  ret = 3;
alpar@1:                goto done;
alpar@1:             }
alpar@1:          }
alpar@1:          /* current circulation is optimal */
alpar@1:          ret = 0;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* now we need to find a cycle (t, a, s, ..., t), which allows
alpar@1:          increasing the flow along it, where a is the out-of-kilter arc
alpar@1:          just found */
alpar@1:       /* link[i] = 0 means that node i is not labelled yet;
alpar@1:          link[i] = a means that arc a immediately precedes node i */
alpar@1:       /* initially only node s is labelled */
alpar@1:       for (i = 1; i <= nv; i++)
alpar@1:          link[i] = 0;
alpar@1:       link[s] = aok, list[1] = s, pos1 = pos2 = 1;
alpar@1:       /* breadth first search */
alpar@1:       while (pos1 <= pos2)
alpar@1:       {  /* dequeue node i */
alpar@1:          i = list[pos1++];
alpar@1:          /* consider all arcs incident to node i */
alpar@1:          for (k = ptr[i]; k < ptr[i+1]; k++)
alpar@1:          {  a = arc[k];
alpar@1:             if (tail[a] == i)
alpar@1:             {  /* a = i->j is a forward arc from s to t */
alpar@1:                j = head[a];
alpar@1:                /* if node j has been labelled, skip the arc */
alpar@1:                if (link[j] != 0) continue;
alpar@1:                /* if the arc does not allow increasing the flow through
alpar@1:                   it, skip the arc */
alpar@1:                if (x[a] >= cap[a]) continue;
alpar@1:                if (overflow(cost[a], pi[i] - pi[j]))
alpar@1:                {  ret = 2;
alpar@1:                   goto done;
alpar@1:                }
alpar@1:                lambda = cost[a] + (pi[i] - pi[j]);
alpar@1:                if (lambda > 0 && x[a] >= low[a]) continue;
alpar@1:             }
alpar@1:             else if (head[a] == i)
alpar@1:             {  /* a = i<-j is a backward arc from s to t */
alpar@1:                j = tail[a];
alpar@1:                /* if node j has been labelled, skip the arc */
alpar@1:                if (link[j] != 0) continue;
alpar@1:                /* if the arc does not allow decreasing the flow through
alpar@1:                   it, skip the arc */
alpar@1:                if (x[a] <= low[a]) continue;
alpar@1:                if (overflow(cost[a], pi[j] - pi[i]))
alpar@1:                {  ret = 2;
alpar@1:                   goto done;
alpar@1:                }
alpar@1:                lambda = cost[a] + (pi[j] - pi[i]);
alpar@1:                if (lambda < 0 && x[a] <= cap[a]) continue;
alpar@1:             }
alpar@1:             else
alpar@1:                xassert(a != a);
alpar@1:             /* label node j and enqueue it */
alpar@1:             link[j] = a, list[++pos2] = j;
alpar@1:             /* check for breakthrough */
alpar@1:             if (j == t) goto brkt;
alpar@1:          }
alpar@1:       }
alpar@1:       /* NONBREAKTHROUGH */
alpar@1:       /* consider all arcs, whose one endpoint is labelled and other is
alpar@1:          not, and determine maximal change of node potentials */
alpar@1:       delta = 0;
alpar@1:       for (a = 1; a <= na; a++)
alpar@1:       {  i = tail[a], j = head[a];
alpar@1:          if (link[i] != 0 && link[j] == 0)
alpar@1:          {  /* a = i->j, where node i is labelled, node j is not */
alpar@1:             if (overflow(cost[a], pi[i] - pi[j]))
alpar@1:             {  ret = 2;
alpar@1:                goto done;
alpar@1:             }
alpar@1:             lambda = cost[a] + (pi[i] - pi[j]);
alpar@1:             if (x[a] <= cap[a] && lambda > 0)
alpar@1:                if (delta == 0 || delta > + lambda) delta = + lambda;
alpar@1:          }
alpar@1:          else if (link[i] == 0 && link[j] != 0)
alpar@1:          {  /* a = j<-i, where node j is labelled, node i is not */
alpar@1:             if (overflow(cost[a], pi[i] - pi[j]))
alpar@1:             {  ret = 2;
alpar@1:                goto done;
alpar@1:             }
alpar@1:             lambda = cost[a] + (pi[i] - pi[j]);
alpar@1:             if (x[a] >= low[a] && lambda < 0)
alpar@1:                if (delta == 0 || delta > - lambda) delta = - lambda;
alpar@1:          }
alpar@1:       }
alpar@1:       if (delta == 0)
alpar@1:       {  /* there is no feasible circulation */
alpar@1:          ret = 1;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* increase potentials of all unlabelled nodes */
alpar@1:       for (i = 1; i <= nv; i++)
alpar@1:       {  if (link[i] == 0)
alpar@1:          {  if (overflow(pi[i], delta))
alpar@1:             {  ret = 2;
alpar@1:                goto done;
alpar@1:             }
alpar@1:             pi[i] += delta;
alpar@1:          }
alpar@1:       }
alpar@1:       goto loop;
alpar@1: brkt: /* BREAKTHROUGH */
alpar@1:       /* walk through arcs of the cycle (t, a, s, ..., t) found in the
alpar@1:          reverse order and determine maximal change of the flow */
alpar@1:       delta = 0;
alpar@1:       for (j = t;; j = i)
alpar@1:       {  /* arc a immediately precedes node j in the cycle */
alpar@1:          a = link[j];
alpar@1:          if (head[a] == j)
alpar@1:          {  /* a = i->j is a forward arc of the cycle */
alpar@1:             i = tail[a];
alpar@1:             lambda = cost[a] + (pi[i] - pi[j]);
alpar@1:             if (lambda > 0 && x[a] < low[a])
alpar@1:             {  /* x[a] may be increased until its lower bound */
alpar@1:                temp = low[a] - x[a];
alpar@1:             }
alpar@1:             else if (lambda <= 0 && x[a] < cap[a])
alpar@1:             {  /* x[a] may be increased until its upper bound */
alpar@1:                temp = cap[a] - x[a];
alpar@1:             }
alpar@1:             else
alpar@1:                xassert(a != a);
alpar@1:          }
alpar@1:          else if (tail[a] == j)
alpar@1:          {  /* a = i<-j is a backward arc of the cycle */
alpar@1:             i = head[a];
alpar@1:             lambda = cost[a] + (pi[j] - pi[i]);
alpar@1:             if (lambda < 0 && x[a] > cap[a])
alpar@1:             {  /* x[a] may be decreased until its upper bound */
alpar@1:                temp = x[a] - cap[a];
alpar@1:             }
alpar@1:             else if (lambda >= 0 && x[a] > low[a])
alpar@1:             {  /* x[a] may be decreased until its lower bound */
alpar@1:                temp = x[a] - low[a];
alpar@1:             }
alpar@1:             else
alpar@1:                xassert(a != a);
alpar@1:          }
alpar@1:          else
alpar@1:             xassert(a != a);
alpar@1:          if (delta == 0 || delta > temp) delta = temp;
alpar@1:          /* check for end of the cycle */
alpar@1:          if (i == t) break;
alpar@1:       }
alpar@1:       xassert(delta > 0);
alpar@1:       /* increase the flow along the cycle */
alpar@1:       for (j = t;; j = i)
alpar@1:       {  /* arc a immediately precedes node j in the cycle */
alpar@1:          a = link[j];
alpar@1:          if (head[a] == j)
alpar@1:          {  /* a = i->j is a forward arc of the cycle */
alpar@1:             i = tail[a];
alpar@1:             /* overflow cannot occur */
alpar@1:             x[a] += delta;
alpar@1:          }
alpar@1:          else if (tail[a] == j)
alpar@1:          {  /* a = i<-j is a backward arc of the cycle */
alpar@1:             i = head[a];
alpar@1:             /* overflow cannot occur */
alpar@1:             x[a] -= delta;
alpar@1:          }
alpar@1:          else
alpar@1:             xassert(a != a);
alpar@1:          /* check for end of the cycle */
alpar@1:          if (i == t) break;
alpar@1:       }
alpar@1:       goto loop;
alpar@1: done: /* free working arrays */
alpar@1:       xfree(ptr);
alpar@1:       xfree(arc);
alpar@1:       xfree(link);
alpar@1:       xfree(list);
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /* eof */