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

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <mao@mai2.rcnet.ru>.
*
*  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 "glplib.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 */