/* glpnet06.c (out-of-kilter algorithm) */

 /***********************************************************************

 *  This code is part of GLPK (GNU Linear Programming Kit).

 *

 *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,

 *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,

 *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.

 *  E-mail: <mao@gnu.org>.

 *

 *  GLPK is free software: you can redistribute it and/or modify it

 *  under the terms of the GNU General Public License as published by

 *  the Free Software Foundation, either version 3 of the License, or

 *  (at your option) any later version.

 *

 *  GLPK is distributed in the hope that it will be useful, but WITHOUT

 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

 *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

 *  License for more details.

 *

 *  You should have received a copy of the GNU General Public License

 *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.

 ***********************************************************************/

 #include "glpenv.h"

 #include "glpnet.h"

 /***********************************************************************

 *  NAME

 *

 *  okalg - out-of-kilter algorithm

 *

 *  SYNOPSIS

 *

 *  #include "glpnet.h"

 *  int okalg(int nv, int na, const int tail[], const int head[],

 *     const int low[], const int cap[], const int cost[], int x[],

 *     int pi[]);

 *

 *  DESCRIPTION

 *

 *  The routine okalg implements the out-of-kilter algorithm to find a

 *  minimal-cost circulation in the specified flow network.

 *

 *  INPUT PARAMETERS

 *

 *  nv is the number of nodes, nv >= 0.

 *

 *  na is the number of arcs, na >= 0.

 *

 *  tail[a], a = 1,...,na, is the index of tail node of arc a.

 *

 *  head[a], a = 1,...,na, is the index of head node of arc a.

 *

 *  low[a], a = 1,...,na, is an lower bound to the flow through arc a.

 *

 *  cap[a], a = 1,...,na, is an upper bound to the flow through arc a,

 *  which is the capacity of the arc.

 *

 *  cost[a], a = 1,...,na, is a per-unit cost of the flow through arc a.

 *

 *  NOTES

 *

 *  1. Multiple arcs are allowed, but self-loops are not allowed.

 *

 *  2. It is required that 0 <= low[a] <= cap[a] for all arcs.

 *

 *  3. Arc costs may have any sign.

 *

 *  OUTPUT PARAMETERS

 *

 *  x[a], a = 1,...,na, is optimal value of the flow through arc a.

 *

 *  pi[i], i = 1,...,nv, is Lagrange multiplier for flow conservation

 *  equality constraint corresponding to node i (the node potential).

 *

 *  RETURNS

 *

 *  0  optimal circulation found;

 *

 *  1  there is no feasible circulation;

 *

 *  2  integer overflow occured;

 *

 *  3  optimality test failed (logic error).

 *

 *  REFERENCES

 *

 *  L.R.Ford, Jr., and D.R.Fulkerson, "Flows in Networks," The RAND

 *  Corp., Report R-375-PR (August 1962), Chap. III "Minimal Cost Flow

 *  Problems," pp.113-26. */

 static int overflow(int u, int v)

 {     /* check for integer overflow on computing u + v */

       if (u > 0 && v > 0 && u + v < 0) return 1;

       if (u < 0 && v < 0 && u + v > 0) return 1;

       return 0;

 }

 int okalg(int nv, int na, const int tail[], const int head[],

       const int low[], const int cap[], const int cost[], int x[],

       int pi[])

 {     int a, aok, delta, i, j, k, lambda, pos1, pos2, s, t, temp, ret,

          *ptr, *arc, *link, *list;

       /* sanity checks */

       xassert(nv >= 0);

       xassert(na >= 0);

       for (a = 1; a <= na; a++)

       {  i = tail[a], j = head[a];

          xassert(1 <= i && i <= nv);

          xassert(1 <= j && j <= nv);

          xassert(i != j);

          xassert(0 <= low[a] && low[a] <= cap[a]);

       }

       /* allocate working arrays */

       ptr = xcalloc(1+nv+1, sizeof(int));

       arc = xcalloc(1+na+na, sizeof(int));

       link = xcalloc(1+nv, sizeof(int));

       list = xcalloc(1+nv, sizeof(int));

       /* ptr[i] := (degree of node i) */

       for (i = 1; i <= nv; i++)

          ptr[i] = 0;

       for (a = 1; a <= na; a++)

       {  ptr[tail[a]]++;

          ptr[head[a]]++;

       }

       /* initialize arc pointers */

       ptr[1]++;

       for (i = 1; i < nv; i++)

          ptr[i+1] += ptr[i];

       ptr[nv+1] = ptr[nv];

       /* build arc lists */

       for (a = 1; a <= na; a++)

       {  arc[--ptr[tail[a]]] = a;

          arc[--ptr[head[a]]] = a;

       }

       xassert(ptr[1] == 1);

       xassert(ptr[nv+1] == na+na+1);

       /* now the indices of arcs incident to node i are stored in

          locations arc[ptr[i]], arc[ptr[i]+1], ..., arc[ptr[i+1]-1] */

       /* initialize arc flows and node potentials */

       for (a = 1; a <= na; a++)

          x[a] = 0;

       for (i = 1; i <= nv; i++)

          pi[i] = 0;

 loop: /* main loop starts here */

       /* find out-of-kilter arc */

       aok = 0;

       for (a = 1; a <= na; a++)

       {  i = tail[a], j = head[a];

          if (overflow(cost[a], pi[i] - pi[j]))

          {  ret = 2;

             goto done;

          }

          lambda = cost[a] + (pi[i] - pi[j]);

          if (x[a] < low[a] || lambda < 0 && x[a] < cap[a])

          {  /* arc a = i->j is out of kilter, and we need to increase

                the flow through this arc */

             aok = a, s = j, t = i;

             break;

          }

          if (x[a] > cap[a] || lambda > 0 && x[a] > low[a])

          {  /* arc a = i->j is out of kilter, and we need to decrease

                the flow through this arc */

             aok = a, s = i, t = j;

             break;

          }

       }

       if (aok == 0)

       {  /* all arcs are in kilter */

          /* check for feasibility */

          for (a = 1; a <= na; a++)

          {  if (!(low[a] <= x[a] && x[a] <= cap[a]))

             {  ret = 3;

                goto done;

             }

          }

          for (i = 1; i <= nv; i++)

          {  temp = 0;

             for (k = ptr[i]; k < ptr[i+1]; k++)

             {  a = arc[k];

                if (tail[a] == i)

                {  /* a is outgoing arc */

                   temp += x[a];

                }

                else if (head[a] == i)

                {  /* a is incoming arc */

                   temp -= x[a];

                }

                else

                   xassert(a != a);

             }

             if (temp != 0)

             {  ret = 3;

                goto done;

             }

          }

          /* check for optimality */

          for (a = 1; a <= na; a++)

          {  i = tail[a], j = head[a];

             lambda = cost[a] + (pi[i] - pi[j]);

             if (lambda > 0 && x[a] != low[a] ||

                 lambda < 0 && x[a] != cap[a])

             {  ret = 3;

                goto done;

             }

          }

          /* current circulation is optimal */

          ret = 0;

          goto done;

       }

       /* now we need to find a cycle (t, a, s, ..., t), which allows

          increasing the flow along it, where a is the out-of-kilter arc

          just found */

       /* link[i] = 0 means that node i is not labelled yet;

          link[i] = a means that arc a immediately precedes node i */

       /* initially only node s is labelled */

       for (i = 1; i <= nv; i++)

          link[i] = 0;

       link[s] = aok, list[1] = s, pos1 = pos2 = 1;

       /* breadth first search */

       while (pos1 <= pos2)

       {  /* dequeue node i */

          i = list[pos1++];

          /* consider all arcs incident to node i */

          for (k = ptr[i]; k < ptr[i+1]; k++)

          {  a = arc[k];

             if (tail[a] == i)

             {  /* a = i->j is a forward arc from s to t */

                j = head[a];

                /* if node j has been labelled, skip the arc */

                if (link[j] != 0) continue;

                /* if the arc does not allow increasing the flow through

                   it, skip the arc */

                if (x[a] >= cap[a]) continue;

                if (overflow(cost[a], pi[i] - pi[j]))

                {  ret = 2;

                   goto done;

                }

                lambda = cost[a] + (pi[i] - pi[j]);

                if (lambda > 0 && x[a] >= low[a]) continue;

             }

             else if (head[a] == i)

             {  /* a = i<-j is a backward arc from s to t */

                j = tail[a];

                /* if node j has been labelled, skip the arc */

                if (link[j] != 0) continue;

                /* if the arc does not allow decreasing the flow through

                   it, skip the arc */

                if (x[a] <= low[a]) continue;

                if (overflow(cost[a], pi[j] - pi[i]))

                {  ret = 2;

                   goto done;

                }

                lambda = cost[a] + (pi[j] - pi[i]);

                if (lambda < 0 && x[a] <= cap[a]) continue;

             }

             else

                xassert(a != a);

             /* label node j and enqueue it */

             link[j] = a, list[++pos2] = j;

             /* check for breakthrough */

             if (j == t) goto brkt;

          }

       }

       /* NONBREAKTHROUGH */

       /* consider all arcs, whose one endpoint is labelled and other is

          not, and determine maximal change of node potentials */

       delta = 0;

       for (a = 1; a <= na; a++)

       {  i = tail[a], j = head[a];

          if (link[i] != 0 && link[j] == 0)

          {  /* a = i->j, where node i is labelled, node j is not */

             if (overflow(cost[a], pi[i] - pi[j]))

             {  ret = 2;

                goto done;

             }

             lambda = cost[a] + (pi[i] - pi[j]);

             if (x[a] <= cap[a] && lambda > 0)

                if (delta == 0 || delta > + lambda) delta = + lambda;

          }

          else if (link[i] == 0 && link[j] != 0)

          {  /* a = j<-i, where node j is labelled, node i is not */

             if (overflow(cost[a], pi[i] - pi[j]))

             {  ret = 2;

                goto done;

             }

             lambda = cost[a] + (pi[i] - pi[j]);

             if (x[a] >= low[a] && lambda < 0)

                if (delta == 0 || delta > - lambda) delta = - lambda;

          }

       }

       if (delta == 0)

       {  /* there is no feasible circulation */

          ret = 1;

          goto done;

       }

       /* increase potentials of all unlabelled nodes */

       for (i = 1; i <= nv; i++)

       {  if (link[i] == 0)

          {  if (overflow(pi[i], delta))

             {  ret = 2;

                goto done;

             }

             pi[i] += delta;

          }

       }

       goto loop;

 brkt: /* BREAKTHROUGH */

       /* walk through arcs of the cycle (t, a, s, ..., t) found in the

          reverse order and determine maximal change of the flow */

       delta = 0;

       for (j = t;; j = i)

       {  /* arc a immediately precedes node j in the cycle */

          a = link[j];

          if (head[a] == j)

          {  /* a = i->j is a forward arc of the cycle */

             i = tail[a];

             lambda = cost[a] + (pi[i] - pi[j]);

             if (lambda > 0 && x[a] < low[a])

             {  /* x[a] may be increased until its lower bound */

                temp = low[a] - x[a];

             }

             else if (lambda <= 0 && x[a] < cap[a])

             {  /* x[a] may be increased until its upper bound */

                temp = cap[a] - x[a];

             }

             else

                xassert(a != a);

          }

          else if (tail[a] == j)

          {  /* a = i<-j is a backward arc of the cycle */

             i = head[a];

             lambda = cost[a] + (pi[j] - pi[i]);

             if (lambda < 0 && x[a] > cap[a])

             {  /* x[a] may be decreased until its upper bound */

                temp = x[a] - cap[a];

             }

             else if (lambda >= 0 && x[a] > low[a])

             {  /* x[a] may be decreased until its lower bound */

                temp = x[a] - low[a];

             }

             else

                xassert(a != a);

          }

          else

             xassert(a != a);

          if (delta == 0 || delta > temp) delta = temp;

          /* check for end of the cycle */

          if (i == t) break;

       }

       xassert(delta > 0);

       /* increase the flow along the cycle */

       for (j = t;; j = i)

       {  /* arc a immediately precedes node j in the cycle */

          a = link[j];

          if (head[a] == j)

          {  /* a = i->j is a forward arc of the cycle */

             i = tail[a];

             /* overflow cannot occur */

             x[a] += delta;

          }

          else if (tail[a] == j)

          {  /* a = i<-j is a backward arc of the cycle */

             i = head[a];

             /* overflow cannot occur */

             x[a] -= delta;

          }

          else

             xassert(a != a);

          /* check for end of the cycle */

          if (i == t) break;

       }

       goto loop;

 done: /* free working arrays */

       xfree(ptr);

       xfree(arc);

       xfree(link);

       xfree(list);

       return ret;

 }

 /* eof */