mxf.cpp

/* SAS modified this file. */
/* (C) Copyright 2003 Jens Lysgaard. All rights reserved. */
/* OSI Certified Open Source Software */
/* This software is licensed under the Common Public License Version 1.0 */

#include <stdlib.h>
#include <stdio.h>
#include <cassert>
#include <float.h>
#include "mxf.h"
#include "memmod.h"

typedef struct
{
  double R; /* Residual capacity of the arc */
  double U; /* Capacity of the arc */
  int Tail;
  int Head;
  int Mate;
  int NextInArc;
  int NextOutArc;
} MXF_ArcRec;
typedef MXF_ArcRec *MXF_ArcPtr;

typedef struct
{
  int DLabel;
   double Excess;
  int FirstInArc;
  int LastInArc;
  int FirstOutArc;
  int LastOutArc;
  int CurrentArc;
  int PrevBNode;    /* Previous Bucket Node */
  int NextBNode;    /* Next -               */
  int PrevDLNode;   /* Previous Distance Label Node */
  int NextDLNode;   /* Next -                       */
} MXF_NodeRec;
typedef MXF_NodeRec *MXF_NodePtr;

typedef struct
{
  MXF_ArcPtr ArcPtr;
  MXF_NodePtr NodePtr;
  int ArcListDim;
  int ArcListSize;
  int NodeListDim;
  int NodeListSize;
  int *FirstInBPtr;   /* First Node in Bucket[i] */
  int *FirstInDLPtr;  /* First Node in DistanceLabel[i] */
  int **ArcMapPtr;
  char *CVWrk1;
  int *IVWrk1;
  int *IVWrk2;
  int *IVWrk3;
  int *IVWrk4;
  int *IVWrk5;
} MXF_Rec;
typedef MXF_Rec *MXF_Ptr;

void MXF_WriteArcList(MaxFlowPtr Ptr)
{
  int i;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  printf("ArcList: (ArcListSize=%d, ArcListDim=%d)\n",
          P->ArcListSize,P->ArcListDim);
  printf("[NodeList: (NodeListSize=%d,NodeListDim=%d]\n",
          P->NodeListSize,P->NodeListDim);
  printf(" Arc#    R    U Tail Head Mate NextOutArc NextInArc\n");
  printf("---------------------------------------------------\n");
  for (i=1; i<=P->ArcListSize; i++)
  printf("%5d%5g%5g%5d%5d%5d%7d%10d\n",
         i,
         P->ArcPtr[i].R,
         P->ArcPtr[i].U,
         P->ArcPtr[i].Tail,
         P->ArcPtr[i].Head,
         P->ArcPtr[i].Mate,
         P->ArcPtr[i].NextOutArc,
         P->ArcPtr[i].NextInArc);
  printf("----------------------------------------------\n");
}

void MXF_WriteNodeList(MaxFlowPtr Ptr)
{
  int i;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  printf("NodeList: (NodeListSize=%d,NodeListDim=%d)\n",
          P->NodeListSize,P->NodeListDim);
  printf(" Node  DL   Exc FOArc LOArc FIArc LIArc PrevB NextB PrevDL NextDL\n");
  printf("-----------------------------------------------------------------\n");
  for (i=1; i<=P->NodeListSize; i++)
  {
    printf("%5d%4d%6g%6d%6d%6d%6d%10d%10d%10d%10d\n",
           i,
           P->NodePtr[i].DLabel,
           P->NodePtr[i].Excess,
           P->NodePtr[i].FirstOutArc,
           P->NodePtr[i].LastOutArc,
           P->NodePtr[i].FirstInArc,
           P->NodePtr[i].LastInArc,
           P->NodePtr[i].PrevBNode,
           P->NodePtr[i].NextBNode,
           P->NodePtr[i].PrevDLNode,
           P->NodePtr[i].NextDLNode);
  }
  printf("---------------------------------------------------------------------\n");
}

void MXF_InitMem(MaxFlowPtr *Ptr,
                 int TotalNodes,
                 int TotalArcs)
{
  MXF_Ptr P;
  int i;

  (*Ptr) = (MaxFlowPtr) MemGet(sizeof(MXF_Rec));

  P = (MXF_Ptr) (*Ptr);

  P->ArcPtr  = NULL;
  P->NodePtr = NULL;

  P->ArcMapPtr = NULL;

  P->ArcListDim   = 0;
  P->ArcListSize  = 0;
  P->NodeListDim  = 0;
  P->NodeListSize = 0;

  P->ArcPtr = (MXF_ArcPtr) MemGet(sizeof(MXF_ArcRec) * (TotalArcs+1));
  P->ArcListDim = TotalArcs; /* Space for 1,...,TotalArcs (excl. #0) */
  P->ArcListSize = 0; /* Currently no arcs in the network */

  P->NodePtr = (MXF_NodePtr) MemGet(sizeof(MXF_NodeRec) * (TotalNodes+1));
  P->NodeListDim = TotalNodes;
  P->NodeListSize = 0; /* Currently no nodes in the network */

  P->FirstInBPtr  = MemGetIV(TotalNodes+1);
  P->FirstInDLPtr = MemGetIV(TotalNodes+1);

  for (i=0; i<=TotalNodes; i++)
  {
    P->FirstInBPtr[i]  = 0;
    P->FirstInDLPtr[i] = 0;
  }

  P->CVWrk1 = MemGetCV(TotalNodes+1);
  P->IVWrk1 = MemGetIV(TotalNodes+1);
  P->IVWrk2 = MemGetIV(TotalNodes+1);
  P->IVWrk3 = MemGetIV(TotalNodes+1);
  P->IVWrk4 = MemGetIV(TotalNodes+1);
  P->IVWrk5 = MemGetIV(TotalNodes+1);
}

void MXF_FreeMem(MaxFlowPtr Ptr)
{
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  MemFree(P->CVWrk1);
  MemFree(P->IVWrk1);
  MemFree(P->IVWrk2);
  MemFree(P->IVWrk3);
  MemFree(P->IVWrk4);
  MemFree(P->IVWrk5);

  MemFree(P->FirstInBPtr);
  MemFree(P->FirstInDLPtr);

  MemFree(P->NodePtr);
  MemFree(P->ArcPtr);

  if (P->ArcMapPtr != NULL)
  MemFreeIM(P->ArcMapPtr,P->NodeListSize+1);

  MemFree(P);
}


void MXF_ClearNodeList(MaxFlowPtr Ptr)
{
  int i;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  for (i=0; i<=P->NodeListDim; i++)
  {
    P->NodePtr[i].DLabel      = 0;
    P->NodePtr[i].Excess      = 0;
    P->NodePtr[i].FirstInArc  = 0;
    P->NodePtr[i].LastInArc   = 0;
    P->NodePtr[i].FirstOutArc = 0;
    P->NodePtr[i].LastOutArc  = 0;
  }

  P->NodeListSize = 0;
}

void MXF_ClearArcList(MaxFlowPtr Ptr)
{
  int i;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  for (i=0; i<=P->ArcListDim; i++)
  {
    P->ArcPtr[i].R          = 0;
    P->ArcPtr[i].U          = 0;
    P->ArcPtr[i].Tail       = 0;
    P->ArcPtr[i].Head       = 0;
    P->ArcPtr[i].Mate       = 0;
    P->ArcPtr[i].NextInArc  = 0;
    P->ArcPtr[i].NextOutArc = 0;
  }

  P->ArcListSize = 0;
}

void MXF_SetNodeListSize(MaxFlowPtr Ptr,
                         int TotalNodes)
{
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  if (TotalNodes > P->NodeListDim)
  {
    /* Allocate more memory. */
    printf("Insufficient memory allocated:\n");
    printf("MXF_SetNodeListSize: TotalNodes > NodeListDim (%d > %d)\n",
            TotalNodes,P->NodeListDim);
    assert(0);
    return;
    //exit(0);
  }

  P->NodeListSize = TotalNodes;
}

void LMXF_AddArc(MaxFlowPtr Ptr,
                 int Tail,
                 int Head,
                 double Capacity,
                 int *Index)
{
  int i,j;
  MXF_Ptr P;
  int64_t TotMem;

  P = (MXF_Ptr) Ptr;

  if ((Tail<=0) || (Tail>P->NodeListSize) ||
      (Head<=0) || (Head>P->NodeListSize) ||
      (Capacity < 0))
  {
    printf("Error in input to MXF_AddArc(NodeListSize=%d)\n",
           P->NodeListSize);
    //exit(0);
    assert(0);
    return;
  }

  i = P->ArcListSize + 1;
  if (i > P->ArcListDim)
  {
    /* Allocate more memory. */
    P->ArcListDim = i + 100;
    TotMem = sizeof(MXF_ArcRec) * (P->ArcListDim + 1);
    P->ArcPtr = (MXF_ArcPtr) MemReGet(P->ArcPtr,TotMem);
  }

  (P->ArcListSize)++;

  P->ArcPtr[i].Tail       = Tail;
  P->ArcPtr[i].Head       = Head;
  P->ArcPtr[i].U          = Capacity;
  P->ArcPtr[i].R          = Capacity;
  P->ArcPtr[i].NextInArc  = 0;
  P->ArcPtr[i].NextOutArc = 0;

  /* Update references. */

  if (P->NodePtr[Tail].FirstOutArc == 0)
  {
    P->NodePtr[Tail].FirstOutArc = i;
    P->NodePtr[Tail].LastOutArc  = i;
  }
  else
  {
    j = P->NodePtr[Tail].LastOutArc;
    P->ArcPtr[j].NextOutArc = i;

    P->NodePtr[Tail].LastOutArc = i;
  }

  if (P->NodePtr[Head].FirstInArc == 0)
  {
    P->NodePtr[Head].FirstInArc = i;
    P->NodePtr[Head].LastInArc  = i;
  }
  else
  {
    j = P->NodePtr[Head].LastInArc;
    P->ArcPtr[j].NextInArc = i;

    P->NodePtr[Head].LastInArc = i;
  }

  *Index = i;

}

void MXF_AddArc(MaxFlowPtr Ptr,
                int Tail,
                int Head,
                double Capacity)
{
  int i;
  LMXF_AddArc(Ptr,Tail,Head,Capacity,&i);
}

void MXF_ChgArcCap(MaxFlowPtr Ptr,
                   int Tail,
                   int Head,
                   double Capacity)
{
   int i,ArcNr;
   double Delta;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  if (P->ArcMapPtr == NULL)
  {
    ArcNr = 0;
    for (i=1; i<=P->ArcListSize; i++)
    {
      if ((P->ArcPtr[i].Tail == Tail) && (P->ArcPtr[i].Head == Head))
      {
        ArcNr = i;
        break;
      }
    }
  }
  else
  {
    ArcNr = P->ArcMapPtr[Tail][Head];
  }

  if (ArcNr > 0)
  {
    Delta = Capacity - P->ArcPtr[ArcNr].U;
    P->ArcPtr[ArcNr].R += Delta;

    P->ArcPtr[ArcNr].U = Capacity;
  }
  else
  {
    printf("MXF_ChgArcCap: Arc (%d,%d) not found => stop.\n",
            Tail,Head);
    //exit(0);
    assert(0);
    return;
  }
}

void MXF_UpdateDLList(MaxFlowPtr Ptr)
{
  MXF_Ptr P;
  int i,n,DL,FirstNode;

  P = (MXF_Ptr) Ptr;

  n = P->NodeListSize;

  for (i=1; i<=n; i++)
  {
    P->NodePtr[i].PrevDLNode = 0;
    P->NodePtr[i].NextDLNode = 0;
  }

  for (i=0; i<=n; i++)
  {
    P->FirstInDLPtr[i] = 0;
  }

  for (i=1; i<=n; i++)
  { /* Insert node i */
    DL = P->NodePtr[i].DLabel;

    FirstNode = P->FirstInDLPtr[DL];

    P->NodePtr[i].PrevDLNode = 0;
    P->NodePtr[i].NextDLNode = FirstNode;

    if (FirstNode > 0)
    P->NodePtr[FirstNode].PrevDLNode = i;

    P->FirstInDLPtr[DL] = i;
  }
}

void MXF_ComputeDLabels(MaxFlowPtr Ptr,
                        int Source,
                        int Sink)
{
  int i,n,Index,Tail,Head,Arc;
  int CurrentLabel,LabeledNodes;
  char *Labeled;
  int *Node;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  if ((Source<=0) || (Source>P->NodeListSize) ||
      (Sink<=0) || (Sink>P->NodeListSize))
  {
    printf("Error in input to MXF_ComputeDLabels\n");
    //exit(0);
    assert(0);
    return;
  }

  n = P->NodeListSize;

  Labeled = P->CVWrk1;
  Node    = P->IVWrk1;

  for (i=1; i<=n; i++)
  {
    Node[i]              = 0;
    Labeled[i]           = 0;
    P->NodePtr[i].DLabel = n;
  }

  P->NodePtr[Sink].DLabel = 0;
  Labeled[Sink] = 1;

  Node[1] = Sink;
  LabeledNodes = 1;

  Index = 1;

  do
  {
    Head         = Node[Index];
    CurrentLabel = P->NodePtr[Head].DLabel + 1;

    Arc  = P->NodePtr[Head].FirstInArc;
    while (Arc > 0)
    {
      Tail = P->ArcPtr[Arc].Tail;

      if ((Labeled[Tail]==0) && (P->ArcPtr[Arc].R > 0))
      {
        P->NodePtr[Tail].DLabel = CurrentLabel;
        Labeled[Tail] = 1;
        Node[++LabeledNodes] = Tail;
      }

      Arc = P->ArcPtr[Arc].NextInArc;
    }

    Index++;
  } while (Index <= LabeledNodes);

}


void MXF_CreateMates(MaxFlowPtr Ptr)
{
  int i,k,n;
  int Arc, ReverseArc;
  int *InArc, *OutArc, *InArcIndex, *OutArcIndex;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;
  n = P->NodeListSize;

  InArc       = P->IVWrk1;
  OutArc      = P->IVWrk2;
  InArcIndex  = P->IVWrk3;
  OutArcIndex = P->IVWrk4;

  for (i=1; i<=P->ArcListSize; i++) { P->ArcPtr[i].Mate = 0; }
  for (i=1; i<=n; i++) { InArc[i] = OutArc[i] = 0; }

  for (k=1; k<=n; k++)
  { /* Create mates for node k. */

    Arc = P->NodePtr[k].FirstOutArc;
    while (Arc > 0)
    {
      i = P->ArcPtr[Arc].Head;

      OutArc[i] = k;
      OutArcIndex[i] = Arc;
      Arc = P->ArcPtr[Arc].NextOutArc;
    }

    Arc = P->NodePtr[k].FirstInArc;
    while (Arc > 0)
    {
      i = P->ArcPtr[Arc].Tail;

      InArc[i] = k;
      InArcIndex[i] = Arc;
      Arc = P->ArcPtr[Arc].NextInArc;
    }

    /* Scan the lists to find those indices where not both
       inarc and outarc are present. */

    Arc = P->NodePtr[k].FirstOutArc;
    while (Arc > 0)
    {
      i = P->ArcPtr[Arc].Head;

      if (InArc[i] != k)
      {
        LMXF_AddArc(Ptr,i,k,0,&ReverseArc);
      }
      else
      ReverseArc = InArcIndex[i];

      P->ArcPtr[Arc].Mate = ReverseArc;
      P->ArcPtr[ReverseArc].Mate = Arc;

      Arc = P->ArcPtr[Arc].NextOutArc;
    }

    Arc = P->NodePtr[k].FirstInArc;
    while (Arc > 0)
    {
      if (P->ArcPtr[Arc].Mate == 0)
      {
        i = P->ArcPtr[Arc].Tail;
        LMXF_AddArc(Ptr,k,i,0,&ReverseArc);

        P->ArcPtr[Arc].Mate = ReverseArc;
        P->ArcPtr[ReverseArc].Mate = Arc;
      }

      Arc = P->ArcPtr[Arc].NextInArc;
    }
  }

}

void MXF_CreateArcMap(MaxFlowPtr Ptr)
{
  int i,j,Size,Tail,Head;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  Size = P->NodeListSize;

  P->ArcMapPtr = MemGetIM(Size+1,Size+1);

  for (i=1; i<=Size; i++)
  for (j=1; j<=Size; j++)
  {
    P->ArcMapPtr[i][j] = 0;
  }

  for (i=1; i<=P->ArcListSize; i++)
  {
    Tail = P->ArcPtr[i].Tail;
    Head = P->ArcPtr[i].Head;

    P->ArcMapPtr[Tail][Head] = i;
  }
}

void LMXF_GetCurrentArc(MXF_Ptr P,
                        int Tail,
                        int *Arc)
{
  /* Returns *Arc = 0 if no admissible arc exists out of Tail. */
  int Head,Label;

  Label = P->NodePtr[Tail].DLabel - 1;
  *Arc  = P->NodePtr[Tail].CurrentArc;

  do
  {
    Head = P->ArcPtr[*Arc].Head;
    if ((P->ArcPtr[*Arc].R > 0) && (P->NodePtr[Head].DLabel == Label))
    {
      P->NodePtr[Tail].CurrentArc = *Arc;
      return;
    }
    (*Arc) = P->ArcPtr[*Arc].NextOutArc;
  } while ((*Arc) > 0);

  P->NodePtr[Tail].CurrentArc = 0;
}

void LMXF_Push(MXF_Ptr P,
               int Arc)
{
   double Delta;
   int Tail,Head,Mate;

  Tail = P->ArcPtr[Arc].Tail;
  Head = P->ArcPtr[Arc].Head;
  Mate = P->ArcPtr[Arc].Mate;

  if (P->NodePtr[Tail].Excess < P->ArcPtr[Arc].R)
  Delta = P->NodePtr[Tail].Excess;
  else
  Delta = P->ArcPtr[Arc].R;

  P->NodePtr[Tail].Excess -= Delta;
  P->ArcPtr[Arc].R -= Delta;

  if(P->NodePtr[Head].Excess < DBL_MAX)
     P->NodePtr[Head].Excess += Delta;
  if(P->ArcPtr[Mate].R < DBL_MAX)
     P->ArcPtr[Mate].R += Delta;
}

void LMXF_ClearBucket(MXF_Ptr P)
{
  int i,n;

  n = P->NodeListSize;

  for (i=1; i<=n; i++)
  {
    P->NodePtr[i].PrevBNode  = 0;
    P->NodePtr[i].NextBNode  = 0;
  }

  for (i=0; i<=n; i++)
  {
    P->FirstInBPtr[i]  = 0;
  }
}

void LMXF_AddToBucket(MXF_Ptr P,
                      int Index,
                      int Node)
{
  int FirstNode;

  FirstNode = P->FirstInBPtr[Index];

  P->NodePtr[Node].PrevBNode = 0;
  P->NodePtr[Node].NextBNode = FirstNode;

  if (FirstNode > 0)
  P->NodePtr[FirstNode].PrevBNode = Node;

  P->FirstInBPtr[Index] = Node;
}

void LMXF_AddToDLList(MXF_Ptr P,
                      int Index,
                      int Node)
{
  int FirstNode;

  FirstNode = P->FirstInDLPtr[Index];

  P->NodePtr[Node].PrevDLNode = 0;
  P->NodePtr[Node].NextDLNode = FirstNode;

  if (FirstNode > 0)
  P->NodePtr[FirstNode].PrevDLNode = Node;

  P->FirstInDLPtr[Index] = Node;
}

void LMXF_RemoveFromBucket(MXF_Ptr P,
                           int Index,
                           int Node)
{
  int Prev,Next;

  Prev  = P->NodePtr[Node].PrevBNode;
  Next  = P->NodePtr[Node].NextBNode;

  P->NodePtr[Node].PrevBNode = 0;
  P->NodePtr[Node].NextBNode = 0;

  if (Prev > 0) P->NodePtr[Prev].NextBNode = Next;
  if (Next > 0) P->NodePtr[Next].PrevBNode = Prev;

  if (Prev == 0) P->FirstInBPtr[Index] = Next;
}

void LMXF_RemoveFromDLList(MXF_Ptr P,
                           int Index,
                           int Node)
{
  int Prev,Next;

  Prev  = P->NodePtr[Node].PrevDLNode;
  Next  = P->NodePtr[Node].NextDLNode;

  P->NodePtr[Node].PrevDLNode = 0;
  P->NodePtr[Node].NextDLNode = 0;

  if (Prev > 0) P->NodePtr[Prev].NextDLNode = Next;
  if (Next > 0) P->NodePtr[Next].PrevDLNode = Prev;

  if (Prev == 0) P->FirstInDLPtr[Index] = Next;
}


void LMXF_BucketGetMaxLevel(MXF_Ptr P,
                            int InitLevel,
                            int *Level)
{
  *Level = InitLevel;
  while (*Level > 0)
  {
    if (P->FirstInBPtr[*Level] > 0) return;
    (*Level)--;
  }
}


void MXF_SolveMaxFlow(MaxFlowPtr Ptr,
                      char InitByZeroFlow,
                      int Source,
                      int Sink,
                      double *CutValue,
                      char GetSinkSide,
                      int *NodeListSize,
                      int *NodeList)
{
  /* CreateMates must be called before calling this routine */
  char InsertHead;
  int i,j,n,Node,Tail,Head,Mate,NextInDList,Size;
  double Delta;
  int Arc;
  int DLabel,MinLabel,Level,CurrentLevel,ReLabels,ReLabelsLimit;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;
  n = P->NodeListSize;

  for (i=1; i<=n; i++)
  P->NodePtr[i].CurrentArc = P->NodePtr[i].FirstOutArc;

  if (InitByZeroFlow)
  {
    for (i=1; i<=P->ArcListSize; i++)
    P->ArcPtr[i].R = P->ArcPtr[i].U; /* Zero flow; Flow equals U-R */

    for (i=1; i<=n; i++) P->NodePtr[i].Excess = 0;
  }

  if (InitByZeroFlow)
  {
    MXF_ComputeDLabels(Ptr,Source,Sink);
    P->NodePtr[Source].DLabel = n;

    MXF_UpdateDLList(Ptr);
    LMXF_ClearBucket(P);

    Level = 0;

    Tail = Source;
    Arc = P->NodePtr[Tail].FirstOutArc;
    while (Arc > 0)
    {
      Head = P->ArcPtr[Arc].Head;
      Mate = P->ArcPtr[Arc].Mate;
      P->ArcPtr[Arc].R = 0;
      P->ArcPtr[Mate].R += P->ArcPtr[Arc].U;
      P->NodePtr[Head].Excess = P->ArcPtr[Arc].U;

      if (P->NodePtr[Head].Excess > 0)
      {
        DLabel = P->NodePtr[Head].DLabel;
        LMXF_AddToBucket(P,DLabel,Head);
        if (DLabel > Level) Level = DLabel;
      }

      Arc = P->ArcPtr[Arc].NextOutArc;
    }

  } /* InitByZeroFlow */
  else
  { /* Start from advanced flow */

    /* saturate all arcs out of the source */
    Tail = Source;
    Arc = P->NodePtr[Tail].FirstOutArc;
    while (Arc > 0)
    {
      Head = P->ArcPtr[Arc].Head;
      Mate = P->ArcPtr[Arc].Mate;

      Delta = P->ArcPtr[Arc].R;

      P->NodePtr[Head].Excess += Delta;
      P->ArcPtr[Arc].R = 0;
      P->ArcPtr[Mate].R += Delta;

      Arc = P->ArcPtr[Arc].NextOutArc;
    }

    /* all arcs out of the source are now saturated */

    /* Perform a global relabeling */

    MXF_ComputeDLabels(Ptr,Source,Sink);
    P->NodePtr[Source].DLabel = n;

    MXF_UpdateDLList(Ptr);
    LMXF_ClearBucket(P);

    Level = 0;

    for (Node=1; Node<=n; Node++)
    {
      if (Node == Sink) continue;

      if (P->NodePtr[Node].Excess > 0)
      {
        DLabel = P->NodePtr[Node].DLabel;

        if (DLabel < n)
        {
          LMXF_AddToBucket(P,DLabel,Node);
          if (DLabel > Level) Level = DLabel;
        }
      }

      P->NodePtr[Node].CurrentArc = P->NodePtr[Node].FirstOutArc;
    }
  }

  ReLabels = 0;
  ReLabelsLimit = n;

  BeginMainLoop:

  while (Level > 0)
  {
    Tail = P->FirstInBPtr[Level];

    if (Tail == 0)
    {
      /* No active nodes on current level. */
      Level--;
      continue;
    }

    LMXF_GetCurrentArc(P,Tail,&Arc);
    if (Arc > 0)
    {
      Head = P->ArcPtr[Arc].Head;
      InsertHead = (P->NodePtr[Head].Excess == 0);

      LMXF_Push(P,Arc);
      if (P->NodePtr[Tail].Excess == 0)
      {
        /* Remove Tail from list of active nodes */
        LMXF_RemoveFromBucket(P,Level,Tail);
      }

      if (InsertHead)
      LMXF_AddToBucket(P,P->NodePtr[Head].DLabel,Head);
    }
    else
    { /* Relabel */
      ReLabels++;
      if (ReLabels == ReLabelsLimit)
      { /* Global Relabeling */

        MXF_ComputeDLabels(Ptr,Source,Sink);
        P->NodePtr[Source].DLabel = n;

        MXF_UpdateDLList(Ptr);
        LMXF_ClearBucket(P);

        Level = 0;

        for (Node=1; Node<=n; Node++)
        {
          if (Node == Sink) continue;

          if (P->NodePtr[Node].Excess > 0)
          {
            DLabel = P->NodePtr[Node].DLabel;

            if (DLabel < n)
            {
              LMXF_AddToBucket(P,DLabel,Node);
              if (DLabel > Level) Level = DLabel;
            }
          }

          P->NodePtr[Node].CurrentArc = P->NodePtr[Node].FirstOutArc;
        }

        ReLabels = 0;

        goto BeginMainLoop;
      }


      MinLabel = n;
      Arc = P->NodePtr[Tail].FirstOutArc;

      while (Arc > 0)
      {
        if (P->ArcPtr[Arc].R > 0)
        {
          Head = P->ArcPtr[Arc].Head;
          if (P->NodePtr[Head].DLabel < MinLabel)
          MinLabel = P->NodePtr[Head].DLabel;
        }
        Arc = P->ArcPtr[Arc].NextOutArc;
      }

      P->NodePtr[Tail].CurrentArc = P->NodePtr[Tail].FirstOutArc;

      if (MinLabel < n) MinLabel++;

      LMXF_RemoveFromDLList(P,Level,Tail);
      LMXF_AddToDLList(P,MinLabel,Tail);
      P->NodePtr[Tail].DLabel = MinLabel;

      LMXF_RemoveFromBucket(P,Level,Tail);
      if (MinLabel < n) LMXF_AddToBucket(P,MinLabel,Tail);

      /* Check for Gap Relabeling condition: */
      if (P->FirstInDLPtr[Level] == 0)
      { /* Gap Relabeling */

        for (i=Level+1; i<n; i++)
        {
          j = P->FirstInDLPtr[i];
          while (j > 0)
          {
            NextInDList = P->NodePtr[j].NextDLNode;

            LMXF_RemoveFromDLList(P,i,j);
            LMXF_AddToDLList(P,n,j);

            if (P->NodePtr[j].Excess > 0)
            {
              LMXF_RemoveFromBucket(P,i,j);
            }

            P->NodePtr[j].DLabel = n;

            j = NextInDList;
          }
        }

        CurrentLevel = Level;
        LMXF_BucketGetMaxLevel(P,CurrentLevel,&Level);
      }
      else
      if (MinLabel < n)
      {
        Level = MinLabel;
      }
      else
      {
        CurrentLevel = Level;
        LMXF_BucketGetMaxLevel(P,CurrentLevel,&Level);
      }
    }

  }

  (*CutValue) = P->NodePtr[Sink].Excess;

  MXF_ComputeDLabels(Ptr,Source,Sink);

  Size=0;
  if (GetSinkSide)
  {
    for (Node=1; Node<=n; Node++)
    {
      DLabel = P->NodePtr[Node].DLabel;
      if (DLabel < n) NodeList[++Size] = Node;
    }
  }
  else
  {
    for (Node=1; Node<=n; Node++)
    {
      DLabel = P->NodePtr[Node].DLabel;
      if (DLabel == n) NodeList[++Size] = Node;
    }
  }

  (*NodeListSize) = Size;
}


void MXF_GetNetworkSize(MaxFlowPtr Ptr,
                        int *Nodes,
                        int *Arcs)
{
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;
  *Nodes = P->NodeListSize;
  *Arcs  = P->ArcListSize;
}

void MXF_GetCurrentFlow(MaxFlowPtr Ptr,
                        double *ArcResidualCapacity,
                        double *NodeExcess)
{
  int i;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  for (i=1; i<=P->ArcListSize; i++)
  ArcResidualCapacity[i] = P->ArcPtr[i].R;

  for (i=1; i<=P->NodeListSize; i++)
  NodeExcess[i] = P->NodePtr[i].Excess;
}

void MXF_SetFlow(MaxFlowPtr Ptr,
                 double *ArcResidualCapacity,
                 double *NodeExcess)
{
  int i;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;

  for (i=1; i<=P->ArcListSize; i++)
  P->ArcPtr[i].R = ArcResidualCapacity[i];

  for (i=1; i<=P->NodeListSize; i++)
  P->NodePtr[i].Excess = NodeExcess[i];
}

void MXF_ComputeGHCutTree(MaxFlowPtr Ptr,
                          int CenterNode,
                          double *CutValue,
                          int *NextOnPath)
{
  char Shift;
  int i,j,n,Source,Sink,SourceSideSize,SinkNeighbor;
  double CutVal;
  int *SourceSide;
  MXF_Ptr P;

  P = (MXF_Ptr) Ptr;
  n = P->NodeListSize;

  SourceSide = P->IVWrk5;

  CutValue[CenterNode] = 0;
  for (i=1; i<=n; i++) NextOnPath[i] = CenterNode;

  for (Source=1; Source<=n; Source++)
  if (Source != CenterNode)
  {
    Sink = NextOnPath[Source];

    MXF_SolveMaxFlow(Ptr,1,Source,Sink,
                     &CutVal,0,&SourceSideSize,SourceSide);

    CutValue[Source] = CutVal;

    for (i=1; i<=SourceSideSize; i++)
    {
      j = SourceSide[i];
      if ((j != Source) && (NextOnPath[j] == Sink)) NextOnPath[j] = Source;
    }

    SinkNeighbor = NextOnPath[Sink];
    Shift = 0;

    for (i=1; i<=SourceSideSize; i++)
    {
      j = SourceSide[i];
      if (j == SinkNeighbor)
      {
        Shift = 1;
        break;
      }
    }

    if (Shift)
    {
      NextOnPath[Source] = NextOnPath[Sink];
      NextOnPath[Sink] = Source;
      CutValue[Source] = CutValue[Sink];
      CutValue[Sink] = CutVal;
    }
  }

  /*
  printf("CenterNode = %d\n",CenterNode);
  printf("NextOnPath =");
  for (i=1; i<=n; i++) printf(" %d",NextOnPath[i]);
  printf("\n");
  printf("CutValue =");
  for (i=1; i<=n; i++) printf(" %d",CutValue[i]);
  */
}