eg_pdp.c

Go to the documentation of this file.

#include "eg_pdp.h"
/* offsets for EGmengerEdgeData_t */
static const size_t os[EG_MENGER_NUMOS] = {
  [EG_MENGER_PI] = offsetof (EGmengerNodeData_t, pi),
  [EG_MENGER_DIST] = offsetof (EGmengerNodeData_t, dist),
  [EG_MENGER_ORIG_DIST] = offsetof (EGmengerNodeData_t, orig_dist),
  [EG_MENGER_NDIST] = offsetof (EGmengerNodeData_t, ndist),
  [EG_MENGER_ORIG_NDIST] = offsetof (EGmengerNodeData_t, orig_ndist),
  [EG_MENGER_MARK] = offsetof (EGmengerNodeData_t, marker),
  [EG_MENGER_ORIG_MARK] = offsetof (EGmengerNodeData_t, orig_marker),
  [EG_MENGER_FATHER] = offsetof (EGmengerNodeData_t, father),
  [EG_MENGER_ORIG_FATHER] = offsetof (EGmengerNodeData_t, orig_father),
  [EG_MENGER_HEAP_CONNECTOR] = offsetof (EGmengerNodeData_t, hc),
  [EG_MENGER_ECOST] = offsetof (EGmengerEdgeData_t, cost),
  [EG_MENGER_REDUCED_ECOST] = offsetof (EGmengerEdgeData_t, reduced_cost),
  [EG_MENGER_IS_IN_SOL] = offsetof (EGmengerEdgeData_t, is_in_solution),
  [EG_MENGER_EDATA] = offsetof (EGmengerEdgeData_t, data)
};
/* this macro return the 'edge' ponter asociated to the dual arc 'dedge' */
#define EGd2dGetEdge(dedge) ((graphNodeP)(\
    ((EGmengerEdgeData_t*)(((EGdGraphEdge_t*)(dedge))->data))->data ))

/* this macro return the dijksra cost of a dual edge */
#define EGd2dGetEdgeCost(dedge) (\
    ((EGmengerEdgeData_t*)(((EGdGraphEdge_t*)(dedge))->data))->cost)
/* this macro return the 'edge' ponter asociated to the 'b'-th arc in path 'a'
 * of the domino 'dom'. It takes care of different types of dominos in an easy
 * form */
#define EGddGetEdge(DOM,a,b) ((graphNodeP)(((DOM)->DDtype == DOM_CC_ONE) ? \
    (DOM)->path[a][b]:((EGmengerEdgeData_t*)(((EGdGraphEdge_t*)\
    ((DOM)->path[a][b]))->data))->data))

/* given a ddom->path[i] array, a domino, and an index, return the edge pointer
 * asociated with that arc (taking care of the type of the domino ) */
#define EGddPathGetEdge(dom,path,b) ((graphNodeP)((dom->DDtype == \
    DOM_DUAL_NORM) ? ((EGmengerEdgeData_t*)(((EGdGraphEdge_t*)\
    (path[b]))->data))->data : path[b]))

/* return the first edge in the embeding list of a boyer's graph for node
 * 'curroot' */
#define GETROOTEDGE(curroot, Graph) (Graph->G + Graph->G[curroot].link[1])

/* return the next edge (cyclic) after 'curedge' in the embeding of a boyer's
 * graph for node 'nroot' */
#define GETNEXTEDGE(curedge,curroot, Graph) ((curedge->link[1] < Graph->N) ? GETROOTEDGE(curroot,Graph) : (Graph->G + curedge->link[1]))

/* return the other end-pint of an edge (in boyer's graph format */
#define GETOTHEREND(cedge,G) (G->G[gp_GetTwinArc(G,cedge - G->G)].v)

/* static functions and types */
static inline int getNextA (unsigned int *nextA,
                            unsigned int curA,
                            unsigned int frontier,
                            unsigned int bounds[3][2])
{
  /* if the borders haven't yet been set, we set them */
  if (bounds[curA][0] == frontier)
  {
    *nextA = curA ? curA - 1 : 2;
    return 0;
  }
  if (bounds[curA][1] == frontier)
  {
    *nextA = curA == 2 ? 0 : curA + 1;
    return 0;
  }
  if (!bounds[curA][0])
  {
    bounds[curA][0] = frontier;
    bounds[curA ? curA - 1 : 2][1] = frontier;
    *nextA = curA ? curA - 1 : 2;
    return 0;
  }
  if (!bounds[curA][1])
  {
    bounds[curA][1] = frontier;
    bounds[curA == 2 ? 0 : curA + 1][0] = frontier;
    *nextA = curA == 2 ? 0 : curA + 1;
    return 0;
  }
  MESSAGE (0, "Found crossed paths");
  return 1;

}

/* fix the value of the dual dominos */
int EGfixDualDominos (EGlist_t * domino_list, /* dual dominos to be fixed */

                      int nnodes, /* number of nodes in primal graph */

                      int nedges, /* number of edges in primal graph */

                      double *weight, /* weight of all non small edges */

                      int *edges, /* non small edges in cook's format */

                      int *elim_edges,  /* index of eliminated edges */

                      int nelim_edges,  /* size of the array elim_edges */

                      graphP G, /* pointer to the dual graph in 
                                 * boyer's format */

                      EGmemPool_t * mem)  /* memory pool */
{
  /* local variables */
  EGlistNode_t *dIt;
  unsigned char *primalMarks,
   *nodeMarks,
   *edgeMarks,
    pth;
  EGlist_t *Adfs[3];
  graphNodeP curEdge;
  int curNode;
  EGddomino_t *ddom;
  register unsigned int i,
    j;
  unsigned int curA = 0,
    bounds[3][2],
    nextA = 0;
  int rval;

  /* looking for memory */
  nodeMarks = EGmemPoolSMalloc (mem, unsigned char,
                                nnodes);
  edgeMarks = EGmemPoolSMalloc (mem, unsigned char,
                                nedges);
  primalMarks = EGmemPoolSMalloc (mem, unsigned char,
                                  nnodes);
  Adfs[0] = EGnewList (mem);
  Adfs[1] = EGnewList (mem);
  Adfs[2] = EGnewList (mem);

  /* loop through all dominos */
  for (dIt = domino_list->begin; dIt; dIt = dIt->next)
  {
    ddom = (EGddomino_t *) dIt->this;
    memset (edgeMarks, 0, sizeof (unsigned char) * nedges);
    memset (nodeMarks, 0, sizeof (unsigned char) * nnodes);
    memset (primalMarks, 0, sizeof (unsigned char) * nnodes);
    memset (bounds, 0, sizeof (bounds));

    /* mark all edges in the paths to powers of two */
    for (pth = 0; pth < 3; pth++)
    {
      for (i = ddom->npath[pth]; i--;)
      {
        edgeMarks[(EGddGetEdge (ddom, pth, i)->id)] = 1 << pth;
      }                         /* end for edges in path[pth] */
    }                           /* end loop through paths */

    /* up to this point all edges in the paths are marked 1,2 and 4
     * respectivelly, all other edges are zero */
    /* compute A_1, A_2, A_3 */
    curNode = EGddGetEdge (ddom, 0, 0)->v;
    EGlistPushHead (Adfs[curA], (void *) curNode);
    MESSAGE (EGFDOM_LVL, "\nProcessing Domino %u", j);
    MESSAGE (EGFDOM_LVL, "Puting Node %d in Set %u", curNode, curA);
    while (Adfs[0]->size || Adfs[1]->size || Adfs[2]->size)
    {
      if (Adfs[0]->size)
        curA = 0;
      else if (Adfs[1]->size)
        curA = 1;
      else
        curA = 2;
      curNode = (int) (Adfs[curA]->begin->this);
      EGlistErase (Adfs[curA], Adfs[curA]->begin, nullFree);
      if (nodeMarks[curNode])
        continue;
      MESSAGE (EGFDOM_LVL, "Checking Node %d in Set %u", curNode, curA);
      nodeMarks[curNode] = 1 << (curA);
      primalMarks[curNode] = 1 << (curA);
      curEdge = GETROOTEDGE (curNode, G);
      /* loop through the neighbours of the node */
      do
      {
        /* if edgemarks >=8 then is a traversed edge in the DFS and thus we can
         * discard it and keep going. */
        if (edgeMarks[curEdge->id] < 8)
        {
          /* the edge has no marks, then the node go to the current DFSlist, if
           * the edge has values 1,2 or 4, this implies that we are doing a
           * change of set */
          switch (edgeMarks[curEdge->id])
          {
          case 0:
            nextA = curA;
            break;
          case 1:
          case 2:
          case 4:
            rval = getNextA (&nextA, curA, edgeMarks[curEdge->id], bounds);
            CHECKRVAL (rval);
            break;
          default:
            EXIT (1, "This should not happen");
          }
          /* acording to the end we put the other end in the corespondig DFS
           * list */
          EGlistPushHead (Adfs[nextA], (void *) curEdge->v);
          MESSAGE (EGFDOM_LVL, "Puting Node %d in Set %u", curEdge->v, nextA);
          /* set the edge as traversed */
          edgeMarks[curEdge->id] = 8;
        }
        curEdge = GETNEXTEDGE (curEdge, curNode, G);
      } while (curEdge != GETROOTEDGE (curNode, G));
    }                           /* end DFS */
    TEST (Adfs[0]->size, "Adfs[0] is not empty");
    TEST (Adfs[1]->size, "Adfs[1] is not empty");
    TEST (Adfs[2]->size, "Adfs[2] is not empty");

    /* now we start computing the corrected value of the domino */
    ddom->primalValue = ddom->value;
    for (i = nelim_edges; i--;)
    {
      TESTL (1, edges[elim_edges[i] << 1] >= nnodes,
             "Accessing memory outside bounds");
      TESTL (1, edges[(elim_edges[i] << 1) | 1] >= nnodes,
             "Accessing memory outside bounds");
      if ((primalMarks[edges[elim_edges[i] << 1]] !=
           primalMarks[edges[(elim_edges[i] << 1) | 1]]))
      {
        ddom->primalValue =
          EGdijkstraCostAdd (ddom->primalValue,
                             EGdijkstraToCost (weight[elim_edges[i]]));
      }
    }                           /* end for i */
    /* now we check that the price of the 2-Pie is not bellow the threshold, if
     * so we set it to the reset value. This is to avoid problems solving the
     * master problem */
    if (EGdijkstraCostIsLess (ddom->primalValue, EGdijkstraToCost (1e-9)))
    {
      if (EGdijkstraCostIsLess
          (ddom->primalValue, EGdijkstraToCost (DP_RESET_VALUE)))
      {
        fprintf (stderr, "FIX_DDOM: WARNING! SEC constraints not satisfied.\n");
        fprintf (stderr, "FIX_DDOM: Found a 1-domino with value %lf.\n",
                 EGdijkstraCostToLf (ddom->primalValue));
        fprintf (stderr,
                 "FIX_DDOM: This means there exists at least one SEC violated by at least %lf.\n",
                 (-2.0 * EGdijkstraCostToLf (ddom->primalValue)) / 3.0);
        fprintf (stderr,
                 "FIX_DDOM: Setting this ddom to 0 (will introduce some errors).\n");
        fprintf (stderr,
                 "FIX_DDOM: Error written to file \"ERROR_SEC.txt\".\n");
        FILE *efile;
        efile = fopen ("ERROR_SEC.txt", "a");
        fprintf (efile, "FIX_DDOM: WARNING! SEC constraints not satisfied.\n");
        fprintf (efile, "FIX_DDOM: Found a 1-domino with value %lf.\n",
                 EGdijkstraCostToLf (ddom->primalValue));
        fprintf (efile,
                 "FIX_DDOM: This means there exists at least one SEC violated by at least %lf.\n",
                 (-2.0 * EGdijkstraCostToLf (ddom->primalValue)) / 3.0);
        fclose (efile);
      }
      ddom->primalValue = EGdijkstraToCost (1e-9);
    }
  }                             /* end for ddominos */

  /* ending */
  EGfreeList (Adfs[0]);
  EGfreeList (Adfs[1]);
  EGfreeList (Adfs[2]);
  EGmemPoolSFree (primalMarks, unsigned char,
                  nnodes,
                  mem);
  EGmemPoolSFree (nodeMarks, unsigned char,
                  nnodes,
                  mem);
  EGmemPoolSFree (edgeMarks, unsigned char,
                  nedges,
                  mem);
  return 0;
}

#define DPU_GETID(edge,os) (EGmengerGetEdgeData(edge,os)->id)
/* given a list of dominos of type DOM_DUAL_NORM and the dual (planar) graph we
 * untangle all domino-paths and set path[1] as the longest one. mem is the
 * memory pool from where the dominos where allocated. cook-format edges of the
 * primal graph. */
int EGuntangleAllDomino (EGlist_t * ddom_list,
                         EGdGraph_t * dgraph,
                         EGlist_t ** embed,
                         EGmemPool_t * mem)
{
  /* local variables */
  int rval = 0;
  unsigned int i;
  unsigned int register j;
  EGlistNode_t *ddom_it;
  EGdGraphEdge_t *cur_edge;
  EGddomino_t *cur_ddom;
  unsigned char *node_marks = 0;
  unsigned char *edge_marks = 0;
  unsigned char is_tangled;
  void **path[3] = { 0, 0, 0 };
  unsigned int npath[4];
  unsigned int curpos[3];
  int nroot;
  unsigned int ecount;
  EGlistNode_t *eroot,
   *enexti,
   *enextii;
  EGlist_t *tangled_node = EGnewList (mem);
  int map[3],
    map2[3];

#warning ASSUMING DOMINO PATH ARE DIRECTED PATHS FROM S->T
  MESSAGE (0, "Untangle: untangling %u dual dominos", ddom_list->size);
  edge_marks = EGsMalloc (unsigned char,
                            (dgraph->edge_id_max) >> 1);
  node_marks = EGsMalloc (unsigned char,
                          dgraph->node_id_max);
  /* first we check if there is any intersection */
  for (ddom_it = ddom_list->begin; ddom_it; ddom_it = ddom_it->next)
  {
    MESSAGE (EGPDP_DBGLVL, "Checking domino %p", ddom_it->this);
    memset (node_marks, 0, sizeof (unsigned char) * dgraph->node_id_max);
    memset (edge_marks, 0,
            sizeof (unsigned char) * ((dgraph->edge_id_max) >> 1));
    cur_ddom = (EGddomino_t *) ddom_it->this;
    if ((cur_ddom->DDtype) != (DOM_DUAL_NORM))
      continue;
    /* count edges used in all nodes */
    MESSAGE (EGPDP_DBGLVL, "domino %p s=%u t=%u", ddom_it->this,
             cur_ddom->s->id, cur_ddom->t->id);
    for (i = 3; i--;)
    {
      for (j = 0; j < cur_ddom->npath[i]; j++)
      {
        cur_edge = (EGdGraphEdge_t *) cur_ddom->path[i][j];
        node_marks[cur_edge->head->id]++;
        node_marks[cur_edge->tail->id]++;
        edge_marks[DPU_GETID (cur_edge, os)] |= 1U << i;
        MESSAGE (EGPDP_DBGLVL, "path %u edge %u(%u,%u) id %u(%u) mark %03o",
                 i, j, cur_edge->tail->id, cur_edge->head->id,
                 DPU_GETID (cur_edge, os), cur_edge->id,
                 edge_marks[DPU_GETID (cur_edge, os)]);
      }
    }
    /* find a node with more than 2 edges in the domino touching it that is not
     * the ending points of the domino */
    is_tangled = 0;

    for (j = dgraph->node_id_max; j--;)
    {
      /* we skip the two endpoints of the domino */
      if (j == cur_ddom->s->id)
        continue;
      if (j == cur_ddom->t->id)
        continue;
      if (node_marks[j] > 2)
      {
        TESTG ((rval = (node_marks[j] & 1U)), CLEANUP, "found an node with"
               " odd number of path touching it, this can't be");
        is_tangled = 1;
        break;
      }
    }
    /* if the dual domino is not tangled we continue */
    if (!is_tangled)
      goto SET_MIDDLE;
    MESSAGE (EGPDP_DBGLVL, "domino %p is tangled, untanglin it", ddom_it->this);

    /* now we have a tangled domino in our hands and we need to untangle it,
     * we do so by constructing a new path structure */
    npath[3] = cur_ddom->npath[0] + cur_ddom->npath[1] + cur_ddom->npath[2];
    path[0] = EGmemPoolSMalloc (mem, void *,
                                npath[3]);
    path[1] = EGmemPoolSMalloc (mem, void *,
                                npath[3]);
    path[2] = EGmemPoolSMalloc (mem, void *,
                                npath[3]);
    npath[0] = npath[1] = npath[2] = curpos[0] = curpos[1] = curpos[2] = 0;

    /* we populate the list and create the beginning of the edges */
    EGlistClear (tangled_node, nullFree);
    for (i = 3; i--;)
    {
      for (j = 0; j < cur_ddom->npath[i]; j++)
      {
        cur_edge = (EGdGraphEdge_t *) cur_ddom->path[i][curpos[i]];
        MESSAGE (EGPDP_DBGLVL, "id %3u original mark %03o",
                 DPU_GETID (cur_edge, os),
                 edge_marks[DPU_GETID (cur_edge, os)]);
        edge_marks[DPU_GETID (cur_edge, os)] |= 1 << (3 + i);
        edge_marks[DPU_GETID (cur_edge, os)] &= ~(1U << i);
        MESSAGE (EGPDP_DBGLVL, "id %3u final mark    %03o",
                 DPU_GETID (cur_edge, os),
                 edge_marks[DPU_GETID (cur_edge, os)]);
        path[i][npath[i]++] = cur_edge;
        curpos[i]++;
        if (node_marks[cur_edge->head->id] > 2 && cur_edge->head != cur_ddom->t)
        {
          TESTG ((rval = (node_marks[cur_edge->head->id] & 1U)), CLEANUP,
                 "found an node with odd number of path touching it, this "
                 "can't be");
          EGlistPushBack (tangled_node, (void *) (cur_edge->head->id));
          break;
        }
      }
    }

    /* now we just find the next balanced node in the list of tangled node,
     * compute the local maping, and continue filling out the new paths */
    while ((curpos[0] != cur_ddom->npath[0]) ||
           (curpos[1] != cur_ddom->npath[1]) ||
           (curpos[2] != cur_ddom->npath[2]))
    {
      /* here we find the next map for the next node in the tangled_node 
       * list */
    NEXT_NODE:
      MESSAGE (EGPDP_DBGLVL, "tangled nodes remaining: %u", tangled_node->size);
      TESTG ((rval = !(tangled_node->size)), CLEANUP, "no more tangled nodes "
             "in list but haven't finish covering the paths of the domino");
      nroot = (size_t) tangled_node->begin->this;
      EGlistErase (tangled_node, tangled_node->begin, nullFree);
      /* now we check that half the active edges have been already setted up */
      ecount = node_marks[nroot];
      eroot = EGclistRoot (embed[nroot]);
      enexti = eroot;
      MESSAGE (EGPDP_DBGLVL, "checking if tangled node %d is balanced", nroot);
      do
      {
        if (edge_marks[(size_t) enexti->this] > 7)
          ecount--;
        enexti = EGclistGetNextIt (embed[nroot], enexti);
        MESSAGE (EGPDP_DBGLVL, "checking edge: %u mark %03o",
                 (unsigned) enexti->this, edge_marks[(size_t) enexti->this]);
      } while (enexti != eroot);
      MESSAGE (EGPDP_DBGLVL, "paths:%u edges:%u", node_marks[nroot] - ecount,
               node_marks[nroot]);
      /* if the node is not balanced, we discard it and go to the next node */
      if (ecount << 1 != node_marks[nroot])
        goto NEXT_NODE;
      /* check if ecount is <= 3 */
      TESTG (ecount > 3, CLEANUP, "more than three paths touch this node!");
      /* find eroot such that all next ecount edges (in the domino) are in the 
       * new path (I think that this is true.... but don't have an actual proof
       * for it) */
      MESSAGE (EGPDP_DBGLVL, "Looking for eroot");
      eroot = EGclistRoot (embed[nroot]);
      enexti = eroot;
      do
      {
        ecount = (node_marks[nroot] >> 1);
        enextii = enexti;
        if (edge_marks[(size_t) enexti->this])
        {
          do
          {
            MESSAGE (EGPDP_DBGLVL, "looking edges %u %u %u edge_marks %u",
                     (unsigned) eroot->this, (unsigned) enexti->this,
                     (unsigned) enextii->this,
                     edge_marks[(size_t) enextii->this]);
            if (edge_marks[(size_t) enextii->this])
            {
              if (edge_marks[(size_t) enextii->this] < 8)
                break;
              ecount--;
            }
            if (!ecount)
              break;
            enextii = EGclistGetNextIt (embed[nroot], enextii);
          } while (enextii != enexti);
          if (!ecount)
            break;
        }
        enexti = EGclistGetNextIt (embed[nroot], enexti);
      } while (enexti != eroot);
      TESTG (ecount, CLEANUP, "this should not happen (or may?)");
      eroot = enexti;
      MESSAGE (EGPDP_DBGLVL, "eroot is edge %u", (unsigned) eroot->this);
      /* reset the maping */
      map[0] = map[1] = map[2] = 3;
      map2[0] = map2[1] = map2[2] = 3;
      /* now we just need to fill the map2 and map */
      MESSAGE (EGPDP_DBGLVL, "looking for local map");
      ecount = node_marks[nroot] >> 1;
      for (enexti = eroot, i = 0;
           i < ecount; enexti = EGclistGetNextIt (embed[nroot], enexti))
      {
        MESSAGE (EGPDP_DBGLVL, "edge %3u mark %03o", (unsigned) enexti->this,
                 edge_marks[(size_t) enexti->this]);
        if (edge_marks[(size_t) enexti->this])
        {
          TESTG ((rval = edge_marks[(size_t) enexti->this] < 8), CLEANUP,
                 "this should not happen");
          map2[i++] = edge_marks[(size_t) enexti->this] >> 4;
        }
      }
      MESSAGE (EGPDP_DBGLVL, "map2[0]=%d map2[1]=%d map2[2]=%d i=%u", map2[0],
               map2[1], map2[2], i);
      for (; i; enexti = EGclistGetNextIt (embed[nroot], enexti))
      {
        MESSAGE (EGPDP_DBGLVL, "edge %3u mark %03o", (unsigned) enexti->this,
                 edge_marks[(size_t) enexti->this]);
        if (edge_marks[(size_t) enexti->this])
        {
          TESTG ((rval = edge_marks[(size_t) enexti->this] > 7), CLEANUP,
                 "this should not happen");
          map[map2[--i]] = edge_marks[(size_t) enexti->this] >> 1;
        }
      }
      MESSAGE (EGPDP_DBGLVL, "map[0]=%d map[1]=%d map[2]=%d i=%u", map[0],
               map[1], map[2], i);

      /* now, given the map we keep filling-up the new path */
      for (i = 3; i--;)
      {
        /* the local maping may have only two branches to follow, and we must
         * skip the other, we know wich map to skip because the value is bigger
         * than 2 */
        if (map[i] > 2)
          continue;
        for (j = curpos[map[i]]; j < cur_ddom->npath[map[i]]; j++)
        {
          cur_edge = (EGdGraphEdge_t *) cur_ddom->path[map[i]][curpos[map[i]]];
          edge_marks[DPU_GETID (cur_edge, os)] |= 1 << (3 + i);
          edge_marks[DPU_GETID (cur_edge, os)] &= ~(1U << map[i]);
          path[i][npath[i]++] = cur_edge;
          curpos[map[i]]++;
          if (node_marks[cur_edge->head->id] > 2
              && cur_edge->head != cur_ddom->t)
          {
            TESTG ((rval = (node_marks[cur_edge->head->id] & 1U) &&
                    (cur_edge->head != cur_ddom->t)), CLEANUP,
                   "found an node with odd number of path touching it, this"
                   " can't be");
            EGlistPushBack (tangled_node, (void *) (cur_edge->head->id));
            break;
          }
        }
      }
    }                           /* end computing the new untangled paths */

    /* now we save the untangled paths */
    MESSAGE (EGPDP_DBGLVL, "saving new paths");
    EGmemPoolSFree (cur_ddom->path[0], void *,
                    cur_ddom->npath[0],
                    mem);
    EGmemPoolSFree (cur_ddom->path[1], void *,
                    cur_ddom->npath[1],
                    mem);
    EGmemPoolSFree (cur_ddom->path[2], void *,
                    cur_ddom->npath[2],
                    mem);
    cur_ddom->path[0] = EGmemPoolSMalloc (mem, void *,
                                          npath[0]);
    cur_ddom->path[1] = EGmemPoolSMalloc (mem, void *,
                                          npath[1]);
    cur_ddom->path[2] = EGmemPoolSMalloc (mem, void *,
                                          npath[2]);
    memcpy (cur_ddom->path[0], path[0], sizeof (void *) * npath[0]);
    memcpy (cur_ddom->path[1], path[1], sizeof (void *) * npath[1]);
    memcpy (cur_ddom->path[2], path[2], sizeof (void *) * npath[2]);
    EGmemPoolSFree (path[0], void *,
                    npath[3],
                    mem);
    EGmemPoolSFree (path[1], void *,
                    npath[3],
                    mem);
    EGmemPoolSFree (path[2], void *,
                    npath[3],
                    mem);
    cur_ddom->npath[0] = npath[0];
    cur_ddom->npath[1] = npath[1];
    cur_ddom->npath[2] = npath[2];

  SET_MIDDLE:
    /* here we select wich path is going to be the middle path */
    map[0] = 0;
    map[1] = 1;
    map[2] = 2;
    MESSAGE (EGPDP_DBGLVL, "setting middle path");
#if DP_CHOOSE_MIDDLE == DP_LONGEST
    if (cur_ddom->npath[0] > cur_ddom->npath[1])
    {
      npath[1] = cur_ddom->npath[1];
      cur_ddom->npath[1] = cur_ddom->npath[0];
      cur_ddom->npath[0] = npath[1];
      path[1] = cur_ddom->path[1];
      cur_ddom->path[1] = cur_ddom->path[0];
      cur_ddom->path[0] = path[1];
    }
    if (cur_ddom->npath[2] > cur_ddom->npath[1])
    {
      npath[1] = cur_ddom->npath[1];
      cur_ddom->npath[1] = cur_ddom->npath[2];
      cur_ddom->npath[2] = npath[1];
      path[1] = cur_ddom->path[1];
      cur_ddom->path[1] = cur_ddom->path[2];
      cur_ddom->path[2] = path[1];
    }
#elif DP_CHOOSE_MIDDLE == DP_SHORTEST
    if (cur_ddom->npath[0] < cur_ddom->npath[1])
    {
      npath[1] = cur_ddom->npath[1];
      cur_ddom->npath[1] = cur_ddom->npath[0];
      cur_ddom->npath[0] = npath[1];
      path[1] = cur_ddom->path[1];
      cur_ddom->path[1] = cur_ddom->path[0];
      cur_ddom->path[0] = path[1];
    }
    if (cur_ddom->npath[2] < cur_ddom->npath[1])
    {
      npath[1] = cur_ddom->npath[1];
      cur_ddom->npath[1] = cur_ddom->npath[2];
      cur_ddom->npath[2] = npath[1];
      path[1] = cur_ddom->path[1];
      cur_ddom->path[1] = cur_ddom->path[2];
      cur_ddom->path[2] = path[1];
    }
#elif DP_CHOOSE_MIDDLE == DP_LIGHTEST
    {
      EGdijkstraCost_t vpath[3];
      vpath[0] = EGdijkstraToCost (0.0);
      vpath[1] = EGdijkstraToCost (0.0);
      vpath[2] = EGdijkstraToCost (0.0);
      for (j = cur_ddom->npath[0]; j--;)
        vpath[0] =
          EGdijkstraCostAdd (vpath[0],
                             EGmengerGetEdgeCost (cur_ddom->path[0][j]));
      for (j = cur_ddom->npath[1]; j--;)
        vpath[1] =
          EGdijkstraCostAdd (vpath[1],
                             EGmengerGetEdgeCost (cur_ddom->path[1][j]));
      for (j = cur_ddom->npath[2]; j--;)
        vpath[2] =
          EGdijkstraCostAdd (vpath[2],
                             EGmengerGetEdgeCost (cur_ddom->path[2][j]));
      if (vpath[0] < vpath[1])
      {
        npath[1] = cur_ddom->npath[1];
        cur_ddom->npath[1] = cur_ddom->npath[0];
        cur_ddom->npath[0] = npath[1];
        path[1] = cur_ddom->path[1];
        cur_ddom->path[1] = cur_ddom->path[0];
        cur_ddom->path[0] = path[1];
      }
      if (vpath[2] < vpath[1])
      {
        npath[1] = cur_ddom->npath[1];
        cur_ddom->npath[1] = cur_ddom->npath[2];
        cur_ddom->npath[2] = npath[1];
        path[1] = cur_ddom->path[1];
        cur_ddom->path[1] = cur_ddom->path[2];
        cur_ddom->path[2] = path[1];
      }
    }
#elif DP_CHOOSE_MIDDLE == DP_HEVIEST
    {
      EGdijkstraCost_t vpath[3];
      vpath[0] = EGdijkstraToCost (0.0);
      vpath[1] = EGdijkstraToCost (0.0);
      vpath[2] = EGdijkstraToCost (0.0);
      for (j = cur_ddom->npath[0]; j--;)
        vpath[0] =
          EGdijkstraCostAdd (vpath[0],
                             EGmengerGetEdgeCost (cur_ddom->path[0][j]));
      for (j = cur_ddom->npath[1]; j--;)
        vpath[1] =
          EGdijkstraCostAdd (vpath[1],
                             EGmengerGetEdgeCost (cur_ddom->path[1][j]));
      for (j = cur_ddom->npath[2]; j--;)
        vpath[2] =
          EGdijkstraCostAdd (vpath[2],
                             EGmengerGetEdgeCost (cur_ddom->path[2][j]));
      if (vpath[0] > vpath[1])
      {
        npath[1] = cur_ddom->npath[1];
        cur_ddom->npath[1] = cur_ddom->npath[0];
        cur_ddom->npath[0] = npath[1];
        path[1] = cur_ddom->path[1];
        cur_ddom->path[1] = cur_ddom->path[0];
        cur_ddom->path[0] = path[1];
      }
      if (vpath[2] > vpath[1])
      {
        npath[1] = cur_ddom->npath[1];
        cur_ddom->npath[1] = cur_ddom->npath[2];
        cur_ddom->npath[2] = npath[1];
        path[1] = cur_ddom->path[1];
        cur_ddom->path[1] = cur_ddom->path[2];
        cur_ddom->path[2] = path[1];
      }
    }
#else
#error DP_CHOOSE_MIDDLE Mmust be either DP_LONGEST DP_SHORTEST DP_LIGHTEST or DP_HEVIEST
#endif
  }                             /* end loop through all dual dominos */

  /* ending and cleaning-up */
CLEANUP:
  if (edge_marks)
    EGfree (edge_marks);
  if (node_marks)
    EGfree (node_marks);
  if (tangled_node)
    EGfreeList (tangled_node);
  MESSAGE (0, "Untangle: done");
  return rval;

}

int EGgetPrimalDP (EGddpConstraint_t * ddp, /* dual DP constraint */

                   int nnodes,  /* number of nodes in primal graph */

                   int nedges,  /* number of edges in primal graph */

                   int *const edges,  /* non small edges in cook's format */

                   int *const elim_edges, /* index of eliminated edges */

                   int nelim_edges, /* size of the array elim_edges */

                   double *const weight,  /* weight of all non small edges */

                   int *const ndom, /* store number of dominos */

                   int **const nAset, /* store the size of A each sets */

                   int **const nBset, /* store the size of B each sets */

                   int *const nHandle,  /* store the size of the handle */

                   int ***const Aset, /* store each A set */

                   int ***const Bset, /* store each B set */

                   int **const Handle,  /* store the handle */

                   double *const violation, /* primal violation of the cut */

                   graphP G,    /* pointer to the dual graph in 
                                 * boyer's format */

                   EGmemPool_t * const mem) /* memory pool */
{

  /* local variables */
  unsigned char *primalMarks,
   *nodeMarks,
   *edgeMarks,
   *nodeOld,
   *edgeOld,
    mA,
    mB,
    pth;
  const int topMark = INT_MAX | 1;
  int rval,
    curNode;
  register int j,
    k;
  unsigned int *elimMarks = 0,
    bounds[3][2],
    curA = 0,
    nextA = 0,
    nABedges,
    npt;
  void **ABedges = 0;
  graphNodeP curEdge,
    ep1,
    ep2,
    ep3;
  EGddomino_t *ddom;
  EGlistNode_t *lIt;
  EGlist_t *Hlist[2],
   *H,
   *Alist[3],
   *Adfs[3],
   *HdfsList[2],
   *fA,
   *fB;
  EGdijkstraCost_t p_cost[3] = { EGdijkstraToCost (0), EGdijkstraToCost (0),
    EGdijkstraToCost (0)
  };

  /* some non-null test */
  TEST (!ddp, "This shouldn't happen");
  TEST (!edges, "This shouldn't happen");
  TEST (!weight, "This shouldn't happen");
  TEST (!ndom, "This shouldn't happen");
  TEST (!nAset, "This shouldn't happen");
  TEST (!nBset, "This shouldn't happen");
  TEST (!nHandle, "This shouldn't happen");
  TEST (!violation, "This shouldn't happen");
  TEST (!Aset, "This shouldn't happen");
  TEST (!Bset, "This shouldn't happen");
  TEST (!Handle, "This shouldn't happen");

  /* looking for memory */
  HdfsList[1] = EGnewList (mem);
  HdfsList[0] = EGnewList (mem);
  Hlist[1] = EGnewList (mem);
  Hlist[0] = EGnewList (mem);
  Alist[0] = EGnewList (mem);
  Alist[1] = EGnewList (mem);
  Alist[2] = EGnewList (mem);
  Adfs[0] = EGnewList (mem);
  Adfs[1] = EGnewList (mem);
  Adfs[2] = EGnewList (mem);
  nodeMarks = EGmemPoolSMalloc (mem, unsigned char,
                                nnodes);
  edgeMarks = EGmemPoolSMalloc (mem, unsigned char,
                                nedges);
  primalMarks = EGmemPoolSMalloc (mem, unsigned char,
                                  nnodes);
  memset (primalMarks, 0, sizeof (unsigned char) * nnodes);
  if (nelim_edges)
  {
    elimMarks = EGmemPoolSMalloc (mem, unsigned,
                                  nelim_edges);
    memset (elimMarks, 0, sizeof (unsigned int) * nelim_edges);
  }

  /* alloc memory for the arrays */
  TEST (ddp->ndom < 1, "no dominos");
  *ndom = ddp->ndom;
  *Bset = EGsMalloc (int *,
                     ddp->ndom);
  *Aset = EGsMalloc (int *,
                     ddp->ndom);
  *nAset = EGsMalloc (int,
                      ddp->ndom);
  *nBset = EGsMalloc (int,
                      ddp->ndom);

  /* reset violation to zero */
  *violation = 0;

  /* we use the field 'old' in graph to maintain the parity information to
   * build H */
  nodeOld = EGmemPoolSMalloc (mem, unsigned char,
                              nnodes);
  edgeOld = EGmemPoolSMalloc (mem, unsigned char,
                              nedges);
  memset (nodeOld, 0, sizeof (unsigned char) * nnodes);
  memset (edgeOld, 0, sizeof (unsigned char) * nedges);

  /* now we update the parity of the edges in F */
  for (k = ddp->nF; k--;)
  {
    TESTL (1, ((graphNodeP) (((EGmengerEdgeData_t *)
                              (ddp->F[k]->data))->data))->id >= nedges,
           "Accessing memory outside bounds");
    edgeOld[((graphNodeP) (((EGmengerEdgeData_t *)
                            (ddp->F[k]->data))->data))->id] += 1;
    *violation +=
      EGdijkstraCostToLf (((EGmengerEdgeData_t *) (ddp->F[k]->data))->cost);
  }

  /* here we compute the primal dominoes, and loop through all of them */
  for (j = ddp->ndom; j--;)
  {
    ddom = ddp->ddom[j];
    *violation += EGdijkstraCostToLf (ddom->primalValue);

    /* reset all marks to zero */
    memset (nodeMarks, 0, sizeof (unsigned char) * nnodes);
    memset (edgeMarks, 0, sizeof (unsigned char) * nedges);
    memset (bounds, 0, sizeof (bounds));
    p_cost[0] = EGdijkstraToCost (0);
    p_cost[1] = EGdijkstraToCost (0);
    p_cost[2] = EGdijkstraToCost (0);

    /* mark all paths with 2**i */
    for (pth = 3; pth--;)
    {
      for (k = ddom->npath[pth]; k--;)
      {
        edgeMarks[EGddGetEdge (ddom, pth, k)->id] = 1 << pth;
        p_cost[pth] = EGdijkstraCostAdd (p_cost[pth],
                                         EGd2dGetEdgeCost (ddom->path[pth][k]));
        //EGmengerGetEdgeCost(EGddGetEdge(ddom,pth,k),os));
      }                         /* end for k (edges in path) */
    }                           /* end for pth */

    /* compute A_1, A_2, A_3 */
    curNode = EGddGetEdge (ddom, 0, 0)->v;
    curA = 0;
    EGlistPushHead (Adfs[curA], (void *) curNode);
    MESSAGE (EGFDOM_LVL, "\nProcessing Domino %u", j);
    MESSAGE (EGFDOM_LVL, "Puting Node %d in Set %u", curNode, curA);
    while (Adfs[0]->size || Adfs[1]->size || Adfs[2]->size)
    {
      if (Adfs[0]->size)
        curA = 0;
      else if (Adfs[1]->size)
        curA = 1;
      else
        curA = 2;
      curNode = (int) (Adfs[curA]->begin->this);
      EGlistErase (Adfs[curA], Adfs[curA]->begin, nullFree);
      if (nodeMarks[curNode])
        continue;
      MESSAGE (EGFDOM_LVL, "Checking Node %d in Set %u", curNode, curA);
      EGlistPushBack (Alist[curA], (void *) curNode);
      nodeMarks[curNode] = 1 << (curA);
      primalMarks[curNode] = 1 << (curA);
      curEdge = GETROOTEDGE (curNode, G);
      /* loop through the neighbours of the node */
      do
      {
        /* if edgemarks >=8 then is a traversed edge in the DFS and thus we can
         * discard it and keep going. */
        if (edgeMarks[curEdge->id] < 8)
        {
          /* the edge has no marks, then the node go to the current DFSlist, if
           * the edge has values 1,2 or 4, this implies that we are doing a
           * change of set */
          switch (edgeMarks[curEdge->id])
          {
          case 0:
            nextA = curA;
            break;
          case 1:
          case 2:
          case 4:
            rval = getNextA (&nextA, curA, edgeMarks[curEdge->id], bounds);
            CHECKRVAL (rval);
            break;
          default:
            EXIT (1, "This should not happen");
          }
          /* acording to the end we put the other end in the corespondig DFS
           * list */
          EGlistPushHead (Adfs[nextA], (void *) curEdge->v);
          MESSAGE (EGFDOM_LVL, "Puting Node %d in Set %u", curEdge->v, nextA);
          /* set the edge as traversed */
          edgeMarks[curEdge->id] = 8;
        }
        curEdge = GETNEXTEDGE (curEdge, curNode, G);
      } while (curEdge != GETROOTEDGE (curNode, G));
    }                           /* end DFS */
    TEST (Adfs[0]->size, "Adfs[0] is not empty");
    TEST (Adfs[1]->size, "Adfs[1] is not empty");
    TEST (Adfs[2]->size, "Adfs[2] is not empty");
    TEST (!Alist[0]->size, "Aset[0] is empty");
    TEST (!Alist[1]->size, "Aset[1] is empty");
    TEST (!Alist[2]->size, "Aset[2] is empty");

    /* compute the first edge of Pi */
    ep1 = EGddGetEdge (ddom, 0, 0);
    ep2 = EGddGetEdge (ddom, 1, 0);
    ep3 = EGddGetEdge (ddom, 2, 0);

    /* this sets the first try to find a 'good' domino path */
    /* up to now we have A, B and C, but we have to choose a domino path that
     * does not eliminate all 'odd' edges (i.e. that leaves a cut for H) */
    /* now we compute fA and fB as the two smallest sets */
    if (EGdijkstraCostIsLess (p_cost[0], p_cost[1]) &&
        EGdijkstraCostIsLess (p_cost[0], p_cost[2]))
    {
      mA = nodeMarks[ep1->v];
      mB = nodeMarks[GETOTHEREND (ep1, G)];
      ABedges = ddom->path[0];
      nABedges = ddom->npath[0];
    }
    else if (EGdijkstraCostIsLess (p_cost[1], p_cost[0]) &&
             EGdijkstraCostIsLess (p_cost[1], p_cost[2]))
    {
      mA = nodeMarks[ep2->v];
      mB = nodeMarks[GETOTHEREND (ep2, G)];
      ABedges = ddom->path[1];
      nABedges = ddom->npath[1];
    }
    else
    {
      mA = nodeMarks[ep3->v];
      mB = nodeMarks[GETOTHEREND (ep3, G)];
      ABedges = ddom->path[2];
      nABedges = ddom->npath[2];
    }
    fA = Alist[mA / 2];
    fB = Alist[mB / 2];
    TEST (nABedges == 0, "one domino path has zero length");

    /* now we update the 'odd' status of the eliminated edges */
    for (k = nelim_edges; k--;)
    {
      TESTL (1, edges[elim_edges[k] << 1] >= nnodes,
             "Accessing memory outside bounds");
      TESTL (1, edges[(elim_edges[k] << 1) | 1] >= nnodes,
             "Accessing memory outside bounds");
      if (primalMarks[edges[elim_edges[k] << 1]] +
          primalMarks[edges[(elim_edges[k] << 1) | 1]] == mA + mB)
        elimMarks[k]++;
    }                           /* end for k (elimEdges) */

    /* now we set-up nAset, nBset, Aset, Bset */
    (*nAset)[j] = fA->size;
    (*nBset)[j] = fB->size;
    (*Aset)[j] = EGsMalloc (int,
                            fA->size);
    (*Bset)[j] = EGsMalloc (int,
                            fB->size);
    /* copy set A */
    for (lIt = fA->begin, k = 0; lIt; lIt = lIt->next, k++)
      (*Aset)[j][k] = (int) (lIt->this);
    /* copy set B */
    for (lIt = fB->begin, k = 0; lIt; lIt = lIt->next, k++)
      (*Bset)[j][k] = (int) (lIt->this);

    /* update parity of the ABedges */
    for (k = nABedges; k--;)
    {
      edgeOld[EGddPathGetEdge (ddom, ABedges, k)->id] += 1;
    }

    /* clear all A, B, C, DFSlist is already empty */
    EGlistClear (Alist[0], nullFree);
    EGlistClear (Alist[1], nullFree);
    EGlistClear (Alist[2], nullFree);
  }                             /* end for j (ddominos) */

  /* now we compute the handle, to do that we compute handle and co-handle,
   * note that those sets might be non-contiguous */
  npt = 0;
  EGlistPushHead (HdfsList[npt & 1], (void *) 0);
  while (HdfsList[0]->size || HdfsList[1]->size)
  {
    /* this allow to change current partition */
    if (!(HdfsList[npt & 1]->size))
      npt++;
    curNode = ((int) (HdfsList[npt & 1]->begin->this));
    EGlistErase (HdfsList[npt & 1], HdfsList[npt & 1]->begin, nullFree);
    if (nodeOld[curNode])
      continue;
    EGlistPushBack (Hlist[npt & 1], (void *) curNode);
    nodeOld[curNode] = npt + 1;
    TESTL (1, curNode >= nnodes, "Accessing memory outside bounds");
    primalMarks[curNode] = npt + 1;
    curEdge = GETROOTEDGE (curNode, G);
    do
    {
      /* if the edge is odd or it has an odd mark, we don't look at it */
      TESTL (1, curEdge->id >= nedges, "Accessing memory outside bounds");
      if (!(edgeOld[curEdge->id] & 1))
      {

        EGlistPushHead (HdfsList[npt & 1], (void *) curEdge->v);
        TESTL (1, curEdge->id >= nedges, "Accessing memory outside bounds");
        edgeOld[curEdge->id] = topMark;
      }                         /* end current partition */
      /* if it is odd, (but not topMarked) we add the other end to the other
       * HdfsList */
      else if (edgeOld[curEdge->id] < topMark)
      {
        EGlistPushHead (HdfsList[(npt + 1) & 1], (void *) curEdge->v);
      }                         /* end complement */
      curEdge = GETNEXTEDGE (curEdge, curNode, G);
    } while (curEdge != GETROOTEDGE (curNode, G));
  }                             /* end computing H and H^c */
  EXIT (Hlist[0]->size + Hlist[1]->size != (unsigned) nnodes,
        "The primal graph is not connected.");

  /* now we update the 'odd' status of the eliminated edges */
  for (k = nelim_edges; k--;)
  {
    TESTL (1, edges[elim_edges[k] << 1] >= nnodes,
           "Accessing memory outside bounds");
    TESTL (1, edges[(elim_edges[k] << 1) | 1] >= nnodes,
           "Accessing memory outside bounds");
    if (((primalMarks[edges[elim_edges[k] << 1]]) & 1) !=
        ((primalMarks[edges[(elim_edges[k] << 1) | 1]]) & 1))
      elimMarks[k]++;
    if (elimMarks[k] & 1)
      *violation += weight[elim_edges[k]];
  }                             /* end for k (elimEdges) */

  /* final update of the violation */
  *violation -= 1;
  //*violation /= 2;

  /* now we define H as the smallest Hlist */
  H = (Hlist[0]->size < Hlist[1]->size) ? Hlist[0] : Hlist[1];
  *nHandle = H->size;
  if (H->size)
    *Handle = EGsMalloc (int,
                         H->size);
  for (lIt = H->begin, k = 0; lIt; lIt = lIt->next, k++)
    (*Handle)[k] = (int) (lIt->this);

  /* ending */
  EGmemPoolSFree (primalMarks, unsigned char,
                  nnodes,
                  mem);
  EGmemPoolSFree (nodeOld, unsigned char,
                  nnodes,
                  mem);
  EGmemPoolSFree (nodeMarks, unsigned char,
                  nnodes,
                  mem);
  EGmemPoolSFree (edgeOld, unsigned char,
                  nedges,
                  mem);
  EGmemPoolSFree (edgeMarks, unsigned char,
                  nedges,
                  mem);
  if (nelim_edges)
    EGmemPoolSFree (elimMarks, unsigned int,
                    nelim_edges,
                    mem);
  EGfreeList (Alist[0]);
  EGfreeList (Alist[1]);
  EGfreeList (Alist[2]);
  EGfreeList (HdfsList[0]);
  EGfreeList (HdfsList[1]);
  EGfreeList (Hlist[0]);
  EGfreeList (Hlist[1]);
  EGfreeList (Adfs[0]);
  EGfreeList (Adfs[1]);
  EGfreeList (Adfs[2]);
  return 0;
}

---