eg_dptighten.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include "dp_config.h"
#include "eg_mempool.h"
#include "eg_dptighten.h"
#define DPT_ENABLE 1
#if DPT_ENABLE
#include "eg_ugraph.h"
#include "eg_bit.h"
#include "eg_util.h"
#include "eg_bbtree.h"
/* ========================================================================= */
/* various size, note that the maximujm sizes must be of the form 2^n */
#define DPT_MAX_DOM 512U
#define DPT_MAX_DOM_LOG 9U
#define DPT_MAX_NODE 131072U
#define DPT_MAX_NODE_LOG 17U
#define DPT_MAX_DEPTH 1U
/* ========================================================================= */
/* data structure to handle information about the graph */
/* ========================================================================= */

/* ========================================================================= */
typedef struct
{
  unsigned int n_id:DPT_MAX_NODE_LOG; 
  unsigned int move_id:5;             
  unsigned int dom;                   
}
DPTMove_t;

/* ========================================================================= */
typedef struct
{
  DPTMove_t move[DPT_MAX_DEPTH];    
  double val;                       
  unsigned int depth;               
}
DPTFullMove_t;

/* ========================================================================= */
typedef struct
{
  unsigned int n_H[2];
  unsigned int n_dominos:DPT_MAX_DOM_LOG;
  unsigned int n_nodes:DPT_MAX_NODE_LOG;  
  unsigned int n_A[DPT_MAX_DOM];          
  unsigned int n_B[DPT_MAX_DOM];          
  unsigned int n_Tc[DPT_MAX_DOM];     
}
DPTGdata_t;

/* ========================================================================= */
typedef struct
{
  unsigned int n_AtoB:DPT_MAX_DOM_LOG;  
  unsigned int n_dT:DPT_MAX_DOM_LOG;  
  unsigned int in_F:1;        
}
DPTEdata_t;

/* ========================================================================= */
struct DPTNdata_t
{
  /* ======================================================================= */
  /* link to the internal information */
  /* ======================================================================= */
  EGbbtreeNode_t *tree_cn;      
  /* ======================================================================= */
  /* cut description, this HAVE to be keeped wright at every iteration */
  /* ======================================================================= */
  EGbitset_t Ia[DPT_MAX_DOM >> __EGBIT_SHIFT__];
  EGbitset_t Ib[DPT_MAX_DOM >> __EGBIT_SHIFT__];
  EGbitset_t added[DPT_MAX_DOM >> __EGBIT_SHIFT__];
  EGbitset_t flipA[DPT_MAX_DOM >> __EGBIT_SHIFT__];
  EGbitset_t flipB[DPT_MAX_DOM >> __EGBIT_SHIFT__];
  unsigned int n_in_A:DPT_MAX_DOM_LOG;
  unsigned int n_in_T:DPT_MAX_DOM_LOG;
  unsigned int in_H:1;          
  unsigned int added_H:1;       
  /* ======================================================================= */
  /* move information, this may be changed later on */
  /* ======================================================================= */
  DPTFullMove_t full_move;      
};
typedef struct DPTNdata_t DPTNdata_t;

/* ========================================================================= */
/* static variables, these variables hold the graphData, node data and edge
 * data, the later two indexed by node->id and edge->id respectivelly. */
static DPTNdata_t *nodeData = 0;
static DPTEdata_t *edgeData = 0;
static DPTGdata_t graphData;
static EGuGraph_t *G = 0;
static EGuGraphNode_t **all_nodes;
static EGuGraphEdge_t **all_edges;
static double const *weight;
static EGbbtree_t *tree;
static EGbitset_t node_update[DPT_MAX_NODE >> __EGBIT_SHIFT__];

/* ========================================================================= */
/* some macros to make access easy */
/* ========================================================================= */
/* the minimum improvement is for the local search, if negative it means that
 * it will allow movements that may worsen the inequality */
#define DPTMinImprovement (-1e-6)

/* constants for the moves */
enum move
{
  DPT_Tc_A,                     /* move a node from A to Tc */
  DPT_Tc_B,                     /* move a node from B to Tc */
  DPT_Hc_H,                     /* move a node from H to Hc */
  DPT_A_B,                      /* move a node from A to B */
  DPT_B_A,                      /* move a node from A to B */
  DPT_B_A_flipH,                /* move a node from A to B */
  DPT_A_B_flipH,                /* move a node from A to B */
  DPT_A_Tc,                     /* move a node from A to Tc */
  DPT_B_Tc,                     /* move a node from B to Tc */
  DPT_H_Hc,                     /* move a node from H to Hc */
  DPT_no_move                   /* no feasible move found */
};

/* string names for the moves */
static const char move_name[11][20] = {[DPT_no_move] = "DPT_no_move",
  [DPT_A_Tc] = "DPT_A_Tc",
  [DPT_B_Tc] = "DPT_B_Tc",
  [DPT_A_B] = "DPT_A_B",
  [DPT_H_Hc] = "DPT_H_Hc",
  [DPT_Tc_A] = "DPT_Tc_A",
  [DPT_Tc_B] = "DPT_Tc_B",
  [DPT_B_A] = "DPT_B_A",
  [DPT_Hc_H] = "DPT_Hc_H",
  [DPT_A_B_flipH] = "DPT_A_B_flipH",
  [DPT_B_A_flipH] = "DPT_B_A_flipH"
};
/* inverse move for every move */
static unsigned char DPT_inv_move[11] = {[DPT_no_move] = DPT_no_move,
  [DPT_A_Tc] = DPT_Tc_A,
  [DPT_B_Tc] = DPT_Tc_B,
  [DPT_A_B] = DPT_B_A,
  [DPT_H_Hc] = DPT_Hc_H,
  [DPT_Tc_A] = DPT_A_Tc,
  [DPT_Tc_B] = DPT_B_Tc,
  [DPT_B_A] = DPT_A_B,
  [DPT_Hc_H] = DPT_H_Hc,
  [DPT_A_B_flipH] = DPT_B_A_flipH,
  [DPT_B_A_flipH] = DPT_A_B_flipH
};

/* ========================================================================= */
/* some macros */
#define DPTgetEdge(__e_it) ((EGuGraphEdge_t*)(__e_it->this))
#define DPTgetOtherEndId(__n_id,__e_it) (__n_id == DPTgetEdge(__e_it)->head->id ? DPTgetEdge(__e_it)->tail->id : DPTgetEdge(__e_it)->head->id)
#define DPTgetEdgeId(__e_it) (DPTgetEdge(__e_it)->id)
#define DPTNdata(__this) ((DPTNdata_t*)(__this))

/* ========================================================================= */
/* local static functions */
/* ========================================================================= */

#if DPT_VRB <= DEBUG || DPT_DBG < DEBUG || DPT_EDBG < DEBUG
/* ========================================================================= */
/* this function display on the screen what move are we performing */
/* ========================================================================= */
static inline void DPdisplaySingleMove (const DPTMove_t * const move)
{
  /* local variables */
  EGlistNode_t *e_it;
  unsigned no_id,
    e_id,
    dom = move->dom;
  /* printing information */
  fprintf (stderr, "\tNode %u:%s:%u", move->n_id, move_name[move->move_id],
           dom);
  switch (move->move_id)
  {
  case DPT_A_Tc:
  case DPT_Tc_A:
    fprintf (stderr, ":%u(%u,%u,%u):%s", dom, graphData.n_A[dom],
             graphData.n_B[dom], graphData.n_Tc[dom],
             EGbitTest (nodeData[move->n_id].Ia, dom) ? "A" :
             EGbitTest (nodeData[move->n_id].Ib, dom) ? "B" : "T^c");
    /* printing adjacent edges information */
    for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
    {
      no_id = DPTgetOtherEndId (move->n_id, e_it);
      e_id = DPTgetEdgeId (e_it);
      fprintf (stderr, "\n\t\te=(%u:%s,%u:%s)[%8.6lf]F[%u]",
               move->n_id, EGbitTest (nodeData[move->n_id].Ia, dom) ? "A" :
               EGbitTest (nodeData[move->n_id].Ib, dom) ? "B" : "T^c",
               no_id, EGbitTest (nodeData[no_id].Ia, dom) ? "A" :
               EGbitTest (nodeData[no_id].Ib, dom) ? "B" : "T^c",
               weight[e_id], edgeData[e_id].in_F);
    }
    break;
  case DPT_B_Tc:
  case DPT_Tc_B:
    fprintf (stderr, ":%u(%u,%u,%u):%s", dom, graphData.n_B[dom],
             graphData.n_A[dom], graphData.n_Tc[dom],
             EGbitTest (nodeData[move->n_id].Ia, dom) ? "A" :
             EGbitTest (nodeData[move->n_id].Ib, dom) ? "B" : "T^c");
    /* printing adjacent edges information */
    for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
    {
      no_id = DPTgetOtherEndId (move->n_id, e_it);
      e_id = DPTgetEdgeId (e_it);
      fprintf (stderr, "\n\t\te=(%u:%s,%u:%s)[%8.6lf]F[%u]",
               move->n_id,
               EGbitTest (nodeData[move->n_id].Ia, dom) ? "A" :
               EGbitTest (nodeData[move->n_id].Ib, dom) ? "B" : "T^c",
               no_id, EGbitTest (nodeData[no_id].Ia, dom) ? "A" :
               EGbitTest (nodeData[no_id].Ib, dom) ? "B" : "T^c",
               weight[e_id], edgeData[e_id].in_F);
    }
    break;
  case DPT_A_B_flipH:
  case DPT_B_A_flipH:
  case DPT_A_B:
  case DPT_B_A:
    fprintf (stderr, ":%u(%u,%u,%u):%s", dom, graphData.n_A[dom],
             graphData.n_B[dom], graphData.n_Tc[dom],
             EGbitTest (nodeData[move->n_id].Ia, dom) ? "A" : "B");
    /* printing adjacent edges information */
    for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
    {
      no_id = DPTgetOtherEndId (move->n_id, e_it);
      e_id = DPTgetEdgeId (e_it);
      fprintf (stderr, "\n\t\te=(%u:%s,%u:%s)[%8.6lf]F[%u]",
               move->n_id, EGbitTest (nodeData[move->n_id].Ia, dom) ? "A" :
               EGbitTest (nodeData[move->n_id].Ib, dom) ? "B" : "T^c",
               no_id, EGbitTest (nodeData[no_id].Ia, dom) ? "A" :
               EGbitTest (nodeData[no_id].Ib, dom) ? "B" : "T^c",
               weight[e_id], edgeData[e_id].in_F);
    }
    break;
  case DPT_Hc_H:
  case DPT_H_Hc:
    fprintf (stderr, "(%u,%u):%s", graphData.n_H[1], graphData.n_H[0],
             (nodeData[move->n_id].in_H & 1U) ? "H" : "H^c");
    /* printing adjacent edges information */
    for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
    {
      no_id = DPTgetOtherEndId (move->n_id, e_it);
      e_id = DPTgetEdgeId (e_it);
      fprintf (stderr, "\n\t\te=(%u,%u)[%8.6lf]F[%u]",
               move->n_id, no_id, weight[e_id], edgeData[e_id].in_F);
    }
    break;
  }
  fprintf (stderr, "\n");
}

/* ========================================================================= */
/* this function display on the screen what move are we performing */
/* ========================================================================= */
static inline int DPdisplayMove (unsigned int const nc_id)
{
  unsigned depth = nodeData[nc_id].full_move.depth;
  fprintf (stderr, "Move Stored for node %u:%u\n", nc_id,
           nodeData[nc_id].full_move.depth);
  while (depth--)
    DPdisplaySingleMove (nodeData[nc_id].full_move.move + depth);
  fprintf (stderr, "\tValue %lf\n", nodeData[nc_id].full_move.val);
  fprintf (stderr, "\t");
  __EG_PRINTLOC2__;
  return 1;
}
#endif

/* ========================================================================= */
/* this function compute the violation of the curently stored constraint */
/* ========================================================================= */
static inline void DPTpriceConstraint (double *l_violation)
{
  register unsigned int i;
  /* finally we compute the in_F field for all edges */
  for (i = G->nedges; i--;)
  {
    (*l_violation) += (edgeData[i].in_F + edgeData[i].n_AtoB +
                       edgeData[i].n_dT) * weight[i];
  }
  return;
}

#if (DPT_VRB <= DEBUG-1) || (DPT_EDBG <= DEBUG)
/* ========================================================================= */
/* ========================================================================= */
static inline int DPTDisplayEdges (void)
{
  register unsigned int i;
  register unsigned j;
  unsigned h_id,
    t_id;
  /* finally we compute the in_F field for all edges */
  for (i = G->nedges; i--;)
  {
    h_id = all_edges[i]->head->id;
    t_id = all_edges[i]->tail->id;
    if (weight[i] > 1e-3)
    {
      fprintf (stderr, "e(%u)[%u,%u][%8.6lf] bellongs to: ", i,
               h_id, t_id, weight[i]);
      for (j = graphData.n_dominos; j--;)
      {
        if (EGbitTest (nodeData[h_id].Ia, j) !=
            EGbitTest (nodeData[t_id].Ia, j))
          fprintf (stderr, "delta(A_%u) ", j);
        if (EGbitTest (nodeData[h_id].Ib, j) !=
            EGbitTest (nodeData[t_id].Ib, j))
          fprintf (stderr, "delta(B_%u) ", j);
      }
      if (nodeData[h_id].in_H != nodeData[t_id].in_H)
        fprintf (stderr, "delta(H) ");
      if (edgeData[i].in_F)
        fprintf (stderr, "F ");
      fprintf (stderr, "\n");
    }
  }
  return 1;
}
#endif

/* ========================================================================= */
/* ========================================================================= */
static int EGdpTightNdataCompare (const void *N1,
                                  const void *N2)
{
  /* we first check if the values are near-zero, if so, we sort according to 
   * the type of the movement */
  if ((fabs (((const DPTNdata_t *) N1)->full_move.val) <
       fabs (DPTMinImprovement))
      && (fabs (((const DPTNdata_t *) N2)->full_move.val) <
          fabs (DPTMinImprovement)))
  {
    /* Sort by first move order */
    if (((const DPTNdata_t *) N1)->full_move.move[0].move_id <
        ((const DPTNdata_t *) N2)->full_move.move[0].move_id)
      return -1;
    if (((const DPTNdata_t *) N1)->full_move.move[0].move_id >
        ((const DPTNdata_t *) N2)->full_move.move[0].move_id)
      return 1;
    /* if they have the same type. choose by simplicity */
    if (((const DPTNdata_t *) N1)->full_move.depth <
        ((const DPTNdata_t *) N2)->full_move.depth)
      return -1;
    if (((const DPTNdata_t *) N1)->full_move.depth >
        ((const DPTNdata_t *) N2)->full_move.depth)
      return 1;
    if( N1 < N2) return -1;
    if( N1 > N2) return 1;
    return 0;
  }
  /* otherwise we sort by value */
  if (((const DPTNdata_t *) N1)->full_move.val <
      ((const DPTNdata_t *) N2)->full_move.val)
    return -1;
  if (((const DPTNdata_t *) N1)->full_move.val >
      ((const DPTNdata_t *) N2)->full_move.val)
    return 1;
  if( N1 < N2) return -1;
  if( N1 > N2) return 1;
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
#if DPT_EDBG < DEBUG
#define DPTupdateMove(__cmove,__bmove) __DPTupdateMove(__cmove,__bmove,__FILE__,__LINE__,__func__)
#else
#define DPTupdateMove(__cmove,__bmove) __DPTupdateMove(__cmove,__bmove)
#endif
static inline int __DPTupdateMove (const DPTFullMove_t * const cur_move,
                                   DPTFullMove_t * const best_move
#if DPT_EDBG < DEBUG
                                   ,
                                   const char *file,
                                   const int line,
                                   const char *function
#endif
  )
{
  if ((cur_move->val + fabs (DPTMinImprovement) < best_move->val) ||
      (cur_move->move[0].move_id < best_move->move[0].move_id &&
       (fabs (cur_move->val - best_move->val) < fabs (DPTMinImprovement))))
  {
#if DPT_VRB+1900<DEBUG
    fprintf (stderr, "Storing Move for node %u\n", cur_move->move[0].n_id);
    {
      unsigned l_depth = cur_move->depth + 1;
      while (l_depth--)
        DPdisplaySingleMove (cur_move->move + l_depth);
      fprintf (stderr, "\tValue %lf\n", cur_move->val);
    }
#endif
#if DPT_EDBG < DEBUG
    {
      unsigned int l_depth = cur_move->depth;
      while (l_depth--)
        EXIT (cur_move->move[l_depth].move_id == DPT_no_move,
              "Storing DPT_no_move for depth %u, called from %s at %s:%d",
              l_depth, function, file, line);
    }
#endif
    memcpy (best_move, cur_move, sizeof (DPTFullMove_t));
  }
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTpriceTcA (const DPTMove_t * const cur_move,
                               double *const move_val)
{
  /* local variables */
  EGlistNode_t *e_it = all_nodes[cur_move->n_id]->edges->begin;
  unsigned e_id,
    no_id,
    pos;

  /* we have to compute every move */
  for (; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (cur_move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, cur_move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, cur_move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case the only change is that we have to pay for this edge in the
       * Delta(T) */
      (*move_val) += weight[e_id];
      break;
      /* in this case we save one Delta(T), pay one E(A:B) and flip in_F because
       * we are changing parity of \sum E(A:B) */
    case 1U:
      if (edgeData[e_id].in_F)
        (*move_val) -= weight[e_id];
      else
        (*move_val) += weight[e_id];
      break;
      /* in this case we save one Delta(T) */
    case 2U:
      (*move_val) -= weight[e_id];
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }                             /* end loop through all incident edges */
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTpriceTcB (const DPTMove_t * const cur_move,
                               double *const move_val)
{
  /* local variables */
  EGlistNode_t *e_it = all_nodes[cur_move->n_id]->edges->begin;
  unsigned e_id,
    no_id,
    pos;

  /* we have to compute every move */
  for (; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (cur_move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, cur_move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, cur_move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case the only change is that we have to pay for this edge in the
       * Delta(T) */
      (*move_val) += weight[e_id];
      break;
      /* in this case we save one Delta(T), pay one E(A:B) and flip in_F because
       * we are changing parity of \sum E(A:B) */
    case 2U:
      if (edgeData[e_id].in_F)
        (*move_val) -= weight[e_id];
      else
        (*move_val) += weight[e_id];
      break;
      /* in this case we save one Delta(T) */
    case 1U:
      (*move_val) -= weight[e_id];
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }                             /* end loop through all incident edges */
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTpriceAB (const DPTMove_t * const cur_move,
                              double *const move_val)
{
  /* local variables */
  EGlistNode_t *e_it = all_nodes[cur_move->n_id]->edges->begin;
  unsigned e_id,
    no_id,
    pos;

  /* we have to compute every move */
  for (; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (cur_move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, cur_move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, cur_move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case nathing change */
      break;
      /* in this case we save one E(A:B) but have to flip in_F */
    case 1U:
      if (edgeData[e_id].in_F)
        (*move_val) -= 2 * weight[e_id];
      break;
      /* in this case we pay one E(A:B) but have to flip in_F */
    case 2U:
      if (!edgeData[e_id].in_F)
        (*move_val) += 2 * weight[e_id];
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }                             /* end loop through all incident edges */
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTpriceABflipH (DPTMove_t const *const cur_move,
                                   double *const move_val)
{
  /* local variables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* we have to compute every move */
  for (e_it = all_nodes[cur_move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (cur_move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, cur_move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, cur_move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case nathing change */
      if (edgeData[e_id].in_F)
        (*move_val) -= weight[e_id];
      else
        (*move_val) += weight[e_id];
      break;
      /* in this case we save one E(A:B) but have to double flip in_F */
    case 1U:
      (*move_val) -= weight[e_id];
      break;
      /* in this case we pay one E(A:B) but have to double flip in_F */
    case 2U:
      (*move_val) += weight[e_id];
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }                             /* end loop through all incident edges */
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTpriceBA (DPTMove_t const *const cur_move,
                              double *const move_val)
{
  /* local variables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* we have to compute every move */
  for (e_it = all_nodes[cur_move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (cur_move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, cur_move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, cur_move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case nathing change */
      break;
      /* in this case we save one E(A:B) but have to flip in_F */
    case 2U:
      if (edgeData[e_id].in_F)
        (*move_val) -= 2 * weight[e_id];
      break;
      /* in this case we pay one E(A:B) but have to flip in_F */
    case 1U:
      if (!edgeData[e_id].in_F)
        (*move_val) += 2 * weight[e_id];
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }                             /* end loop through all incident edges */
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTpriceBAflipH (DPTMove_t const *const cur_move,
                                   double *const move_val)
{
  /* local variables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* we have to compute every move */
  for (e_it = all_nodes[cur_move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (cur_move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, cur_move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, cur_move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case nathing change */
      if (edgeData[e_id].in_F)
        (*move_val) -= weight[e_id];
      else
        (*move_val) += weight[e_id];
      break;
      /* in this case we save one E(A:B) but have to double flip in_F */
    case 2U:
      (*move_val) -= weight[e_id];
      break;
      /* in this case we pay one E(A:B) but have to double flip in_F */
    case 1U:
      (*move_val) += weight[e_id];
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }                             /* end loop through all incident edges */
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTpriceHHc (DPTMove_t const *const cur_move,
                               double *const move_val)
{
  /* local variables */
  EGlistNode_t *e_it;
  unsigned e_id;

  /* we have to compute every move */
  for (e_it = all_nodes[cur_move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    if (edgeData[e_id].in_F)
      (*move_val) -= weight[e_id];
    else
      (*move_val) += weight[e_id];
  }                             /* end loop through all incident edges */
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
#define DPTpriceHcH(__id,__val) DPTpriceHHc(__id,__val)

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTpriceBTc (DPTMove_t const *const cur_move,
                               double *const move_val)
{
  /* local variables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* we have to compute every move */
  for (e_it = all_nodes[cur_move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (cur_move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, cur_move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, cur_move->dom) ? 1U : 0U);
    switch (pos)
    {
      /* in this case we save one Delta(T) */
    case 0U:
      (*move_val) -= weight[e_id];
      break;
      /* in this case we save one E(A:B), pay one Delta(T) and flip the F */
    case 2U:
      if (edgeData[e_id].in_F)
        (*move_val) -= weight[e_id];
      else
        (*move_val) += weight[e_id];
      break;
      /* in this case we pay one Delta(T) */
    case 1U:
      (*move_val) += weight[e_id];
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }                             /* end loop through all incident edges */
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTpriceATc (DPTMove_t const *const cur_move,
                               double *const move_val)
{
  /* local variables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* we have to compute every move */
  for (e_it = all_nodes[cur_move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (cur_move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, cur_move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, cur_move->dom) ? 1U : 0U);
    switch (pos)
    {
      /* in this case we save one Delta(T) */
    case 0U:
      (*move_val) -= weight[e_id];
      break;
      /* in this case we save one E(A:B), pay one Delta(T) and flip the F */
    case 1U:
      if (edgeData[e_id].in_F)
        (*move_val) -= weight[e_id];
      else
        (*move_val) += weight[e_id];
      break;
      /* in this case we pay one Delta(T) */
    case 2U:
      (*move_val) += weight[e_id];
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }                             /* end loop through all incident edges */
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTmakeMoveTcA (DPTMove_t const *const move,
                                  const unsigned int update_flags)
{
  /* local vaariables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* we have to change all related edges */
  for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case the only change is that we have to pay for this edge in the
       * Delta(T) */
      edgeData[e_id].n_dT++;
      break;
      /* in this case we save one Delta(T), pay one E(A:B) and flip in_F because
       * we are changing parity of \sum E(A:B) */
    case 1U:
      edgeData[e_id].n_dT--;
      edgeData[e_id].n_AtoB++;
      edgeData[e_id].in_F++;
      break;
      /* in this case we save one Delta(T) */
    case 2U:
      edgeData[e_id].n_dT--;
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }

  /* move the node */
  graphData.n_Tc[move->dom]--;
  graphData.n_A[move->dom]++;
  EGbitSet (nodeData[move->n_id].Ia, move->dom);
  nodeData[move->n_id].n_in_A++;
  nodeData[move->n_id].n_in_T++;
  if (update_flags)
    EGbitSet (nodeData[move->n_id].added, move->dom);
  /* ending */
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTmakeMoveTcB (DPTMove_t const *const move,
                                  const unsigned int update_flags)
{
  /* local vaariables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* move the node */
  graphData.n_Tc[move->dom]--;
  graphData.n_B[move->dom]++;
  EGbitSet (nodeData[move->n_id].Ib, move->dom);
  nodeData[move->n_id].n_in_T++;
  if (update_flags)
    EGbitSet (nodeData[move->n_id].added, move->dom);

  /* we have to change all related edges */
  for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case the only change is that we have to pay for this edge in the
       * Delta(T) */
      edgeData[e_id].n_dT++;
      break;
      /* in this case we save one Delta(T), pay one E(A:B) and flip in_F because
       * we are changing parity of \sum E(A:B) */
    case 2U:
      edgeData[e_id].n_dT--;
      edgeData[e_id].n_AtoB++;
      edgeData[e_id].in_F++;
      break;
      /* in this case we save one Delta(T) */
    case 1U:
      edgeData[e_id].n_dT--;
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTmakeMoveAB (DPTMove_t const *const move,
                                 const unsigned int update_flags)
{
  /* local vaariables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* move the node */
  graphData.n_A[move->dom]--;
  graphData.n_B[move->dom]++;
  EGbitUnset (nodeData[move->n_id].Ia, move->dom);
  EGbitSet (nodeData[move->n_id].Ib, move->dom);
  nodeData[move->n_id].n_in_A--;
  if (update_flags)
    EGbitSet (nodeData[move->n_id].flipA, move->dom);

  /* we have to change all related edges */
  for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case nathing change */
      break;
      /* in this case we save one E(A:B) but have to flip in_F */
    case 1U:
      edgeData[e_id].in_F++;
      edgeData[e_id].n_AtoB--;
      break;
      /* in this case we pay one E(A:B) but have to flip in_F */
    case 2U:
      edgeData[e_id].in_F++;
      edgeData[e_id].n_AtoB++;
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTmakeMoveBA (DPTMove_t const *const move,
                                 const unsigned int update_flags)
{
  /* local vaariables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* move the node */
  graphData.n_B[move->dom]--;
  graphData.n_A[move->dom]++;
  EGbitUnset (nodeData[move->n_id].Ib, move->dom);
  EGbitSet (nodeData[move->n_id].Ia, move->dom);
  nodeData[move->n_id].n_in_A++;
  if (update_flags)
    EGbitSet (nodeData[move->n_id].flipB, move->dom);

  /* we have to change all related edges */
  for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case nathing change */
      break;
      /* in this case we save one E(A:B) but have to flip in_F */
    case 2U:
      edgeData[e_id].in_F++;
      edgeData[e_id].n_AtoB--;
      break;
      /* in this case we pay one E(A:B) but have to flip in_F */
    case 1U:
      edgeData[e_id].in_F++;
      edgeData[e_id].n_AtoB++;
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTmakeMoveBAflipH (DPTMove_t const *const move,
                                      const unsigned int update_flags)
{
  /* local vaariables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* we have to change all related edges */
  for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case nathing change */
      edgeData[e_id].in_F++;
      break;
      /* in this case we save one E(A:B) but have to double flip in_F */
    case 2U:
      edgeData[e_id].n_AtoB--;
      break;
      /* in this case we pay one E(A:B) but have to double flip in_F */
    case 1U:
      edgeData[e_id].n_AtoB++;
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }

  /* update flags if necesary */
  if (update_flags)
  {
    EGbitSet (nodeData[move->n_id].flipB, move->dom);
  }
  /* flip the node in the handle */
  graphData.n_H[nodeData[move->n_id].in_H]--;
  nodeData[move->n_id].in_H++;
  graphData.n_H[nodeData[move->n_id].in_H]++;
  /* change the node from B to A */
  graphData.n_B[move->dom]--;
  graphData.n_A[move->dom]++;
  EGbitUnset (nodeData[move->n_id].Ib, move->dom);
  EGbitSet (nodeData[move->n_id].Ia, move->dom);
  nodeData[move->n_id].n_in_A++;

  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTmakeMoveABflipH (DPTMove_t const *const move,
                                      const unsigned int update_flags)
{
  /* local vaariables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* we have to change all related edges */
  for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, move->dom) ? 1U : 0U);
    switch (pos)
    {
    case 0U:
      /* in this case nathing change */
      edgeData[e_id].in_F++;
      break;
      /* in this case we save one E(A:B) but have to double flip in_F */
    case 1U:
      edgeData[e_id].n_AtoB--;
      break;
      /* in this case we pay one E(A:B) but have to double flip in_F */
    case 2U:
      edgeData[e_id].n_AtoB++;
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }
  /* update flags if necesary */
  if (update_flags)
  {
    EGbitSet (nodeData[move->n_id].flipA, move->dom);
  }
  /* flip the node in the handle */
  graphData.n_H[nodeData[move->n_id].in_H]--;
  nodeData[move->n_id].in_H++;
  graphData.n_H[nodeData[move->n_id].in_H]++;
  /* change the node from B to A */
  graphData.n_A[move->dom]--;
  graphData.n_B[move->dom]++;
  EGbitUnset (nodeData[move->n_id].Ia, move->dom);
  EGbitSet (nodeData[move->n_id].Ib, move->dom);
  nodeData[move->n_id].n_in_A--;

  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTmakeMoveHHc (DPTMove_t const *const move,
                                  const unsigned int update_flags)
{
  /* local vaariables */
  EGlistNode_t *e_it;

  /* update flags if necesary */
  if (update_flags)
    if (nodeData[move->n_id].in_H)
      nodeData[move->n_id].added_H = 1U;
  /* flip the node in the handle */
  graphData.n_H[nodeData[move->n_id].in_H]--;
  nodeData[move->n_id].in_H++;
  graphData.n_H[nodeData[move->n_id].in_H]++;

  /* we have to change all related edges */
  for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
    edgeData[DPTgetEdgeId (e_it)].in_F++;
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
#define DPTmakeMoveHcH(__move,__flags) DPTmakeMoveHHc(__move,__flags)

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTmakeMoveBTc (DPTMove_t const *const move)
{
  /* local vaariables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* change the node from B to Tc */
  graphData.n_B[move->dom]--;
  graphData.n_Tc[move->dom]++;
  EGbitUnset (nodeData[move->n_id].Ib, move->dom);
  nodeData[move->n_id].n_in_T--;

  /* we have to change all related edges */
  for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, move->dom) ? 1U : 0U);
    switch (pos)
    {
      /* in this case we save one Delta(T) */
    case 0U:
      edgeData[e_id].n_dT--;
      break;
      /* in this case we save one E(A:B), pay one Delta(T) and flip the F */
    case 2U:
      edgeData[e_id].in_F++;
      edgeData[e_id].n_AtoB--;
      edgeData[e_id].n_dT++;
      break;
      /* in this case we pay one Delta(T) */
    case 1U:
      edgeData[e_id].n_dT++;
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTTestEdges (void)
{
  /* local variables */
  register unsigned int i = G->nedges,
    j;
  unsigned int n_AtoB,
    n_dT,
    N1_id,
    N2_id,
    N1_inA,
    N1_inB,
    N1_inTc,
    N2_inA,
    N2_inB,
    N2_inTc,
    in_F,
    pos;
  while (i--)
  {
    N1_id = all_edges[i]->head->id;
    N2_id = all_edges[i]->tail->id;
    n_AtoB = n_dT = in_F = 0;
    for (j = graphData.n_dominos; j--;)
    {
      N1_inA = EGbitTest (nodeData[N1_id].Ia, j);
      N1_inB = EGbitTest (nodeData[N1_id].Ib, j);
      N2_inA = EGbitTest (nodeData[N2_id].Ia, j);
      N2_inB = EGbitTest (nodeData[N2_id].Ib, j);
      N1_inTc = !(N1_inA || N1_inB);
      N2_inTc = !(N2_inA || N2_inB);
      if ((!N1_inTc && N2_inTc) || (N1_inTc && !N2_inTc))
        n_dT++;
      if ((N1_inA && N2_inB) || (N1_inB && N2_inA))
        n_AtoB++;
    }
    pos =
      1U & ((nodeData[N1_id].in_H ^ nodeData[N2_id].in_H) ^ edgeData[i].n_AtoB);
    if (pos)
      in_F = 1U;
    TEST (in_F != edgeData[i].in_F || n_dT != edgeData[i].n_dT ||
          n_AtoB != edgeData[i].n_AtoB, "in_F(%u,%u) n_dT(%u,%u) n_AtoB(%u,%u)"
          " don't match for edge %u", in_F, edgeData[i].in_F, n_dT,
          edgeData[i].n_dT, n_AtoB, edgeData[i].n_AtoB, i);
  }
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTmakeMoveATc (DPTMove_t const *const move)
{
  /* local vaariables */
  EGlistNode_t *e_it;
  unsigned e_id,
    no_id,
    pos;

  /* change the node from B to Tc */
  graphData.n_A[move->dom]--;
  graphData.n_Tc[move->dom]++;
  EGbitUnset (nodeData[move->n_id].Ia, move->dom);
  nodeData[move->n_id].n_in_T--;
  nodeData[move->n_id].n_in_A--;

  /* we have to change all related edges */
  for (e_it = all_nodes[move->n_id]->edges->begin; e_it; e_it = e_it->next)
  {
    e_id = DPTgetEdgeId (e_it);
    no_id = DPTgetOtherEndId (move->n_id, e_it);
    pos = (EGbitTest (nodeData[no_id].Ia, move->dom) ? 2U : 0U) |
      (EGbitTest (nodeData[no_id].Ib, move->dom) ? 1U : 0U);
    switch (pos)
    {
      /* in this case we save one Delta(T) */
    case 0U:
      edgeData[e_id].n_dT--;
      break;
      /* in this case we save one E(A:B), pay one Delta(T) and flip the F */
    case 1U:
      edgeData[e_id].in_F++;
      edgeData[e_id].n_AtoB--;
      edgeData[e_id].n_dT++;
      break;
      /* in this case we pay one Delta(T) */
    case 2U:
      edgeData[e_id].n_dT++;
      break;
    default:
      EXIT (1, "This should not happen");
    }                           /* end switch */
  }
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
#define DPTmakeMove(__move,__flags) __DPTmakeMove(__move,__flags,__LINE__,__FILE__,__func__)
static int __DPTmakeMove (DPTMove_t const *const move,
                          const unsigned int update_flags,
                          const int line,
                          const char *file,
                          const char *function)
{
  /* we just switch the correct call */
  MESSAGE (DPT_VRB + 19, "doing move %s for node %u from %s at %s:%d",
           move_name[move->move_id], move->n_id, function, file, line);
  switch (move->move_id)
  {
  case DPT_A_Tc:
    DPTmakeMoveATc (move);
    break;
  case DPT_B_Tc:
    DPTmakeMoveBTc (move);
    break;
  case DPT_H_Hc:
    DPTmakeMoveHHc (move, update_flags);
    break;
  case DPT_Hc_H:
    DPTmakeMoveHcH (move, update_flags);
    break;
  case DPT_A_B:
    DPTmakeMoveAB (move, update_flags);
    break;
  case DPT_B_A:
    DPTmakeMoveBA (move, update_flags);
    break;
  case DPT_A_B_flipH:
    DPTmakeMoveABflipH (move, update_flags);
    break;
  case DPT_B_A_flipH:
    DPTmakeMoveBAflipH (move, update_flags);
    break;
  case DPT_Tc_A:
    DPTmakeMoveTcA (move, update_flags);
    break;
  case DPT_Tc_B:
    DPTmakeMoveTcB (move, update_flags);
    break;
  case DPT_no_move:
    break;
  default:
    EXIT (1, "unknown move %d", move->move_id);
  }
  /* ending */
  //#if DPT_EDBG <= DEBUG
  //EXITRVAL(DPTTestEdges());
  //#endif
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static int DPTmakeInvMove (DPTMove_t const *const move,
                           const unsigned int update_flags)
{
  /* we just switch the correct call */
  MESSAGE (DPT_VRB + 19, "doing inv move %s for node %u",
           move_name[move->move_id], move->n_id);
  switch (DPT_inv_move[move->move_id])
  {
  case DPT_A_Tc:
    DPTmakeMoveATc (move);
    break;
  case DPT_B_Tc:
    DPTmakeMoveBTc (move);
    break;
  case DPT_H_Hc:
    DPTmakeMoveHHc (move, update_flags);
    break;
  case DPT_Hc_H:
    DPTmakeMoveHcH (move, update_flags);
    break;
  case DPT_A_B:
    DPTmakeMoveAB (move, update_flags);
    break;
  case DPT_B_A:
    DPTmakeMoveBA (move, update_flags);
    break;
  case DPT_A_B_flipH:
    DPTmakeMoveABflipH (move, update_flags);
    break;
  case DPT_B_A_flipH:
    DPTmakeMoveBAflipH (move, update_flags);
    break;
  case DPT_Tc_A:
    DPTmakeMoveTcA (move, update_flags);
    break;
  case DPT_Tc_B:
    DPTmakeMoveTcB (move, update_flags);
    break;
  case DPT_no_move:
    break;
  default:
    EXIT (1, "unknown move %d", move->move_id);
  }
  /* ending */
  //#if DPT_EDBG <= DEBUG
  //EXITRVAL(DPTTestEdges());
  //#endif
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTmakeFullMove (DPTFullMove_t const *const full_move)
{
  unsigned int depth = 0;
  for (; depth < full_move->depth; depth++)
    DPTmakeMove (full_move->move + depth, 1);
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
#define DPTgoDeeper(__depth) if(__depth+1<DPT_MAX_DEPTH){\
  /* now we call for deeper moves */\
  DPTmakeMove(base_move->move+__depth,0);\
  for( e_it = all_nodes[nc_id]->edges->begin; e_it ; e_it=e_it->next){\
    /*if(DPTgetOtherEndId(nc_id,e_it)<nc_id)*/\
    DPTSetBestMove( DPTgetOtherEndId( nc_id, e_it), best_move, \
                    base_move, __depth + 1);}\
  /* undo the last move */\
  DPTmakeInvMove(base_move->move+__depth,0);\
  }


/* ========================================================================= */
/* ========================================================================= */
static int DPTSetBestMove (const unsigned int nc_id,
                           DPTFullMove_t * const best_move,
                           DPTFullMove_t * const base_move,
                           const unsigned int depth)
{
  /* local variables */
  unsigned nc_inA,
    nc_inB,
    nc_inTc;
  const double move_val = base_move->val;
  EGlistNode_t *e_it;
  MESSAGE (DPT_VRB + 19, "entering depth %u", depth);
  /* check that the depth is OK, if no, we return back */
  if (depth >= DPT_MAX_DEPTH)
    return 0;
  /* basic set-up */
  base_move->move[depth].n_id = nc_id;
  base_move->move[depth].move_id = DPT_no_move;
  base_move->depth = depth + 1;
  /* now we check the adding moves if they can be considered */
  base_move->move[depth].dom = graphData.n_dominos;
  for (; base_move->move[depth].dom--;)
  {
    nc_inA = EGbitTest (nodeData[nc_id].Ia, base_move->move[depth].dom);
    nc_inB = EGbitTest (nodeData[nc_id].Ib, base_move->move[depth].dom);
    nc_inTc = !(nc_inA || nc_inB);
    /* first check for addind moves, but only if they have not been done */
    if (nc_inTc &&
        !EGbitTest (nodeData[nc_id].added, base_move->move[depth].dom) &&
        (graphData.n_Tc[base_move->move[depth].dom] > 1))
    {
      /* ------------------------------------------------------------------- */
      /* Check DPT_Tc_A */
      base_move->val = move_val;
      base_move->move[depth].move_id = DPT_Tc_A;
      DPTpriceTcA (base_move->move + depth, &(base_move->val));
      DPTupdateMove (base_move, best_move);
      /* now we call for deeper moves */
      DPTgoDeeper (depth);
      /* ------------------------------------------------------------------- */
      /* Check DPT_Tc_B */
      base_move->val = move_val;
      base_move->move[depth].move_id = DPT_Tc_B;
      DPTpriceTcB (base_move->move + depth, &(base_move->val));
      DPTupdateMove (base_move, best_move);
      /* now we call for deeper moves */
      DPTgoDeeper (depth);
      /* ------------------------------------------------------------------- */
    }                           /* end checking adding moves */
  }                             /* end loop through all dominos */
  /* now we check if this node can be added to the handle */
  if (!nodeData[nc_id].in_H && !nodeData[nc_id].added_H
      && (graphData.n_H[0] > 1))
  {
    /* --------------------------------------------------------------------- */
    /* Check DPT_Hc_H */
    base_move->val = move_val;
    base_move->move[depth].move_id = DPT_Hc_H;
    DPTpriceHcH (base_move->move + depth, &(base_move->val));
    DPTupdateMove (base_move, best_move);
    /* now we call for deeper moves */
    DPTgoDeeper (depth);
    /* --------------------------------------------------------------------- */
  }                             /* end checking move to add this node into H */
  /* now we check for moves that flip this node from A to B */
  base_move->move[depth].dom = graphData.n_dominos;
  for (; base_move->move[depth].dom--;)
  {
    nc_inA = EGbitTest (nodeData[nc_id].Ia, base_move->move[depth].dom);
    nc_inB = EGbitTest (nodeData[nc_id].Ib, base_move->move[depth].dom);
    nc_inTc = !(nc_inA || nc_inB);
    /* --------------------------------------------------------------------- */
    /* check DPT_A_B */
    if (nc_inA && (graphData.n_A[base_move->move[depth].dom] > 1) &&
        !EGbitTest (nodeData[nc_id].flipA, base_move->move[depth].dom))
    {
      base_move->val = move_val;
      base_move->move[depth].move_id = DPT_A_B;
      DPTpriceAB (base_move->move + depth, &(base_move->val));
      DPTupdateMove (base_move, best_move);
      /* now we call for deeper moves */
      DPTgoDeeper (depth);
      /* --------------------------------------------------------------------- */
      /* check DPT_A_B_flipH */
      if (graphData.n_H[nodeData[nc_id].in_H] > 1)
      {
        base_move->val = move_val;
        base_move->move[depth].move_id = DPT_A_B_flipH;
        DPTpriceABflipH (base_move->move + depth, &(base_move->val));
        DPTupdateMove (base_move, best_move);
        /* now we call for deeper moves */
        DPTgoDeeper (depth);
      }
    }
    /* --------------------------------------------------------------------- */
    /* check DPT_B_A */
    else if (nc_inB && (graphData.n_B[base_move->move[depth].dom] > 1) &&
             !EGbitTest (nodeData[nc_id].flipB, base_move->move[depth].dom))
    {
      base_move->val = move_val;
      base_move->move[depth].move_id = DPT_B_A;
      DPTpriceBA (base_move->move + depth, &(base_move->val));
      DPTupdateMove (base_move, best_move);
      /* now we call for deeper moves */
      DPTgoDeeper (depth);
      /* --------------------------------------------------------------------- */
      /* check DPT_B_A_flipH */
      if (graphData.n_H[nodeData[nc_id].in_H] > 1)
      {
        base_move->val = move_val;
        base_move->move[depth].move_id = DPT_B_A_flipH;
        DPTpriceBAflipH (base_move->move + depth, &(base_move->val));
        DPTupdateMove (base_move, best_move);
        /* now we call for deeper moves */
        DPTgoDeeper (depth);
      }
    }
    /* --------------------------------------------------------------------- */
  }                             /* end checking non-growing moves */
  /* now we check for moves that may get this node out of a teeth */
  base_move->move[depth].dom = graphData.n_dominos;
  for (; base_move->move[depth].dom--;)
  {
    nc_inA = EGbitTest (nodeData[nc_id].Ia, base_move->move[depth].dom);
    nc_inB = EGbitTest (nodeData[nc_id].Ib, base_move->move[depth].dom);
    nc_inTc = !(nc_inA || nc_inB);
    /* --------------------------------------------------------------------- */
    /* check DPT_A_Tc */
    if (nc_inA && (graphData.n_A[base_move->move[depth].dom] > 1))
    {
      base_move->val = move_val;
      base_move->move[depth].move_id = DPT_A_Tc;
      DPTpriceATc (base_move->move + depth, &(base_move->val));
      DPTupdateMove (base_move, best_move);
      /* now we call for deeper moves */
      DPTgoDeeper (depth);
    }
    /* --------------------------------------------------------------------- */
    /* check DPT_B_Tc */
    else if (nc_inB && (graphData.n_B[base_move->move[depth].dom] > 1))
    {
      base_move->val = move_val;
      base_move->move[depth].move_id = DPT_B_Tc;
      DPTpriceBTc (base_move->move + depth, &(base_move->val));
      DPTupdateMove (base_move, best_move);
      /* now we call for deeper moves */
      DPTgoDeeper (depth);
    }
    /* --------------------------------------------------------------------- */
  }                             /*end checking moves to shrink T */
  /* now we check if this node can be shrinked from the handle */
  if (nodeData[nc_id].in_H && (graphData.n_H[1] > 1))
  {
    /* --------------------------------------------------------------------- */
    /* check DPT_H_Hc */
    base_move->val = move_val;
    base_move->move[depth].move_id = DPT_H_Hc;
    DPTpriceHHc (base_move->move + depth, &(base_move->val));
    DPTupdateMove (base_move, best_move);
    /* now we call for deeper moves */
    DPTgoDeeper (depth);
  }                             /* end checking for shrinking moves for H */
  //base_move->move[depth].move_id = DPT_no_move;
  base_move->val = move_val;
  base_move->depth = depth;
  MESSAGE (DPT_VRB + 19, "done");
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTisMoveFeasible (DPTMove_t const *const move,
                                     const unsigned int update_flags)
{
  int rval = 0;
  /* we just switch the correct call */
  switch (move->move_id)
  {
  case DPT_A_Tc:
    rval = (EGbitTest (nodeData[move->n_id].Ia, move->dom) &&
            graphData.n_A[move->dom] > 1);
    break;
  case DPT_B_Tc:
    rval = (EGbitTest (nodeData[move->n_id].Ib, move->dom) &&
            graphData.n_B[move->dom] > 1);
    break;
  case DPT_H_Hc:
    rval = (nodeData[move->n_id].in_H && graphData.n_H[1] > 1);
    break;
  case DPT_Hc_H:
    rval = (!nodeData[move->n_id].in_H && graphData.n_H[0] > 1 &&
            (!update_flags || !nodeData[move->n_id].added_H));
    break;
  case DPT_A_B:
    rval = (EGbitTest (nodeData[move->n_id].Ia, move->dom) &&
            graphData.n_A[move->dom] > 1 &&
            (!update_flags
             || !EGbitTest (nodeData[move->n_id].flipA, move->dom)));
    break;
  case DPT_B_A:
    rval = (EGbitTest (nodeData[move->n_id].Ib, move->dom) &&
            graphData.n_B[move->dom] > 1 &&
            (!update_flags
             || !EGbitTest (nodeData[move->n_id].flipB, move->dom)));
    break;
  case DPT_A_B_flipH:
    rval = (EGbitTest (nodeData[move->n_id].Ia, move->dom) &&
            graphData.n_A[move->dom] > 1 &&
            graphData.n_H[nodeData[move->n_id].in_H] > 1 &&
            (!update_flags
             || (!EGbitTest (nodeData[move->n_id].flipA, move->dom))));
    break;
  case DPT_B_A_flipH:
    rval = (EGbitTest (nodeData[move->n_id].Ib, move->dom) &&
            graphData.n_B[move->dom] > 1 &&
            graphData.n_H[nodeData[move->n_id].in_H] > 1 &&
            (!update_flags
             || (!EGbitTest (nodeData[move->n_id].flipB, move->dom))));
    break;
  case DPT_Tc_A:
    rval = (!EGbitTest (nodeData[move->n_id].Ia, move->dom) &&
            !EGbitTest (nodeData[move->n_id].Ib, move->dom) &&
            graphData.n_Tc[move->dom] > 1 &&
            (!update_flags
             || !EGbitTest (nodeData[move->n_id].added, move->dom)));
    break;
  case DPT_Tc_B:
    rval = (!EGbitTest (nodeData[move->n_id].Ia, move->dom) &&
            !EGbitTest (nodeData[move->n_id].Ib, move->dom) &&
            graphData.n_Tc[move->dom] > 1 &&
            (!update_flags
             || !EGbitTest (nodeData[move->n_id].added, move->dom)));
    break;
  case DPT_no_move:
    rval = 1;
    break;
  default:
    EXIT (1, "unknown move %d", move->move_id);
  }
  /* ending */
  return rval;
}

/* ========================================================================= */
static inline int DPTisFullMoveFeasible (DPTFullMove_t const *const full_move)
{
  int rval = 1;
  unsigned int depth = 0;
  MESSAGE (DPT_VRB + 19, "entering");
  /* check if each move is feasible, make it, and check next move */
  for (; depth + 1 < full_move->depth &&
       (rval = DPTisMoveFeasible (full_move->move + depth, 1)); depth++)
    DPTmakeMove (full_move->move + depth, 0);
  WARNINGL (DPT_VRB + 19, !rval, "move level %u is infeasible", depth);
  if (rval)
    rval = DPTisMoveFeasible (full_move->move + depth, 1);
  WARNINGL (DPT_VRB + 19, !rval, "move level %u is infeasible", depth);
  /* undo the moves */
  while (depth--)
  {
    DPTmakeInvMove (full_move->move + depth, 0);
  }
  MESSAGE (DPT_VRB + 19, "done, return %d", rval);
  return rval;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTresetFlags (void)
{
  /* local variables */
  register unsigned int i = G->nodes->size;
  /* loop through all nodes */
  while (i--)
  {
    memset (nodeData[i].added, 0, DPT_MAX_DOM >> 3);
    memset (nodeData[i].flipA, 0, DPT_MAX_DOM >> 3);
    memset (nodeData[i].flipB, 0, DPT_MAX_DOM >> 3);
    nodeData[i].added_H = 0;
  }
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTStoreFullMove (const unsigned int nc_id,
                                    DPTFullMove_t const *const full_move)
{
  unsigned int l_depth = full_move->depth;
  memcpy (&(nodeData[nc_id].full_move), full_move, sizeof (DPTFullMove_t));
  while (l_depth--)
  {
#if DPT_EDBG < DEBUG
    EXIT (full_move->move[l_depth].move_id == DPT_no_move, "Storing "
          "DPT_no_move for node %u, depth %u", nc_id, l_depth);
#endif
  }
  EGbbtreeRemove (tree, nodeData[nc_id].tree_cn);
  nodeData[nc_id].tree_cn = EGbbtreeAdd (tree, nodeData + nc_id);
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTResetMove (DPTFullMove_t * move)
{
  memset (move, 0, sizeof (DPTFullMove_t));
  return 0;
}

/* ========================================================================= */
/* ========================================================================= */
static inline int DPTupdateNeighMove (DPTFullMove_t * full_move)
{
  unsigned int const max_depth = full_move->depth + DPT_MAX_DEPTH;
  unsigned int no_id[2 * DPT_MAX_DEPTH + 1],
    n_depth = 0;
  EGlistNode_t *e_it[DPT_MAX_DEPTH * 2];
  DPTFullMove_t base_move,
    cur_move;
  memset (&base_move, 0, sizeof (base_move));
  memset (&cur_move, 0, sizeof (base_move));
  DPTResetMove (&cur_move);
  cur_move.val = DBL_MAX;
  base_move.val = 0;
  no_id[0] = full_move->move[0].n_id;
  DPTSetBestMove (no_id[0], &cur_move, &base_move, 0);
  DPTStoreFullMove (no_id[0], &cur_move);
  EGbitSet (node_update, no_id[0]);
#if DPT_VRB+19<DEBUG
  DPdisplayMove (no_id);
#endif
  e_it[0] = all_nodes[no_id[0]]->edges->begin;
  n_depth = 0;
  /* for each depth move, we update the neighbours up to depth
   * DPT_MAX_DEPTH+full_move->depth */
  while (e_it[0])
  {
    no_id[n_depth + 1] = DPTgetOtherEndId (no_id[n_depth], e_it[n_depth]);
    if (!EGbitTest (node_update, no_id[n_depth + 1]))
    {
      DPTResetMove (&cur_move);
      cur_move.val = DBL_MAX;
      base_move.val = 0;
      DPTSetBestMove (no_id[n_depth + 1], &cur_move, &base_move, 0);
      DPTStoreFullMove (no_id[n_depth + 1], &cur_move);
#if DPT_VRB+19<DEBUG
      DPdisplayMove (no_id);
#endif
      EGbitSet (node_update, no_id[n_depth + 1]);
    }
    /* now we look at the next neighbour, if we are not at the end of the
     * recursion, we go down one more level */
    if (n_depth + 1 < max_depth)
    {
      n_depth++;
      e_it[n_depth] = all_nodes[no_id[n_depth]]->edges->begin;
    }
    /* otherwise we move horizontally in the current level */
    else
    {
      e_it[n_depth] = e_it[n_depth]->next;
      /* if the next e_it is null, we have to move up one level, 
       * if posible */
      while (n_depth && !e_it[n_depth])
      {
        n_depth--;
        e_it[n_depth] = e_it[n_depth]->next;
      }
    }                           /* end looking at next neighbour */
  }                             /* end updating for the given simple move */
/* ========================================================================= */
/*
#warning Note that here we update a node too many times, we could \
  decrease this by adding a flag for each node and update only if they \
  haven't been updated so far. of course, this would need to reset the \
  flags after we are done.
*/
/* ========================================================================= */
  memset (node_update, 0, sizeof (node_update));
  return 0;
}

/* ========================================================================= */
/* given a residual graph G^*, and a valid DP inequality, it will try to find a
 * 'more' violated DP-inequality constraint; the function will not touch the
 * data abour the graph nor about the original inequality, but it will alloc
 * memory for all data in the new returned inequality. Note too that the number
 * of dominos of the new inequality is the same as in the original inequality.
 * Finaly, the violation (if it is positive) of the new inequality is returned
 * in (*violation). If an error occurs the function will return 1, zero on
 * success (this doesn't imply that we were able to find a violated constraint
 * from the original onw */
/* ========================================================================= */
int DPtighten (                 /* graph G* data */
                int n_nodes,
                int n_edges,
                int *edges,
                double const *const external_weight,
                /* original constrain to tighten */
                int n_dominos,
                int *n_aset,
                int *n_bset,
                int n_handle,
                int **aset,
                int **bset,
                int *handle,
                /* new generated constraint */
                int **new_n_aset,
                int **new_n_bset,
                int *new_n_handle,
                int ***new_aset,
                int ***new_bset,
                int **new_handle,
                double *violation)
{
  /* --------------------------------------------------------------------- */
  /* local variables */
  EGmemPool_t *mem;
  int rval = 0;
  unsigned pos,
    max_iter = n_nodes * 5,
    in_H = 0,
    N1_inA,
    N1_inB,
    N2_inA,
    N2_inB,
    N1_id,
    N2_id,
    N1_inTc,
    N2_inTc;
  double cost_tmp,
    l_violation,
    o_violation;
  register unsigned int i,
    j;
  DPTNdata_t *nc_data = 0,
   *n_data;
  DPTFullMove_t cur_move,
    base_move;

  /* --------------------------------------------------------------------- */
  /* we test that the necesary input is not NULL */
  MESSAGE (DPT_VRB + 19, "Entering with n_nodes %d n_edges %d n_dominos %d",
           n_nodes, n_edges, n_dominos);
  EXIT ((unsigned) n_dominos > DPT_MAX_DOM,
        "n_dominos (%d) exced max_n_dom (%u)"
        "change the value of DPT_MAX_NODE to an apropiate value", n_dominos,
        DPT_MAX_DOM);
  EXIT ((unsigned) n_nodes > DPT_MAX_NODE,
        "n_nodes (%d) exced max_n_nodes (%u)"
        "change the value of DPT_MAX_NODE to an apropiate value", n_nodes,
        DPT_MAX_NODE);
  TEST (!n_nodes, "graph has no nodes");
  TEST (!n_edges, "graph has no edges");
  TEST (!n_dominos, "inequality has no dominos");
  TEST (!n_handle, "inequality has handle");
  TEST (!edges, "edges is NULL");
  TEST (!external_weight, "weight is NULL");
  TEST (!n_aset, "n_aset is NULL");
  TEST (!n_bset, "n_bset is NULL");
  TEST (!aset, "aset is NULL");
  TEST (!bset, "bset is NULL");
  TEST (!handle, "handle is NULL");
  TEST (!new_n_aset, "new_n_aset is NULL");
  TEST (!new_n_bset, "new_n_bset is NULL");
  TEST (!new_n_handle, "new_n_handle is NULL");
  TEST (!new_aset, "new_aset is NULL");
  TEST (!new_bset, "new_bset is NULL");
  TEST (!new_handle, "new_handle is NULL");
  TEST (!violation, "violation is NULL");
#if LOG_MODE
  /* we save the graph if necesary */
  saveGraph ("graph_tighten.x", n_nodes, n_edges, edges, external_weight);
#endif
  /* basic tests */
  ADVTESTL (DPT_DBG, (n_nodes < 2) || (n_edges < 1) || !(n_dominos & 1) ||
            (n_handle < 1) || (n_nodes - n_handle < 1),
            1, "either n_nodes (%d), n_edges (%d), n_dominos (%d) or n_handle"
            " (%d) has wrong value", n_nodes, n_edges, n_dominos, n_handle);
  for (i = n_dominos; i--;)
    ADVTESTL (DPT_EDBG, (n_aset[i] < 1) || (n_bset[i] < 1) ||
              (n_nodes - n_aset[i] - n_bset[i] < 1), 1, "either "
              "n_aset[%d] (%d) or n_bset[%d] (%d) or n_Tc[%d] (%d) is < 1",
              i, n_aset[i], i, n_bset[i], i, n_nodes - n_aset[i]);

  /* --------------------------------------------------------------------- */
  /* basic set-up */
  weight = external_weight;
  memset (node_update, 0, sizeof (node_update));
  mem = EGnewMemPool (512, EGmemPoolNewSize, 4096);
  tree = EGnewBbtree (mem, EGdpTightNdataCompare);
  nodeData = EGsMalloc (DPTNdata_t, n_nodes);
  edgeData = EGsMalloc (DPTEdata_t, n_edges);
  all_nodes = EGsMalloc (EGuGraphNode_t *, n_nodes);
  all_edges = EGsMalloc (EGuGraphEdge_t *, n_edges);
  memset (all_nodes, 0, sizeof (EGuGraphNode_t *) * n_nodes);
  memset (nodeData, 0, sizeof (DPTNdata_t) * n_nodes);
  memset (edgeData, 0, sizeof (DPTEdata_t) * n_edges);
  memset (all_edges, 0, sizeof (EGuGraphEdge_t *) * n_nodes);
  memset (&cur_move, 0, sizeof (DPTFullMove_t));
  memset (&base_move, 0, sizeof (DPTFullMove_t));
  l_violation = ((3 * n_dominos + 1));
  (*new_n_aset) = 0;
  (*new_n_bset) = 0;
  (*new_n_handle) = 0;
  (*new_aset) = 0;
  (*new_bset) = 0;
  (*new_handle) = 0;

  /* first we build the graph with its internal data */
  memset (&graphData, 0, sizeof (DPTGdata_t));
  graphData.n_dominos = n_dominos;
  graphData.n_nodes = n_nodes;
  G = EGnewUGraph (mem);
  G->data = &graphData;

  /* now we initialize all the data structures */
  MESSAGE (DPT_VRB + 19, "Initializing");
  /* add all nodes */
  for (i = 0; i < (unsigned) n_nodes; i++)
  {
    nodeData[i].full_move.val = DBL_MAX;
    all_nodes[i] = EGuGraphAddNode (G, nodeData + i);
    nodeData[i].tree_cn = EGbbtreeAdd (tree, nodeData + i);
  }
  /* add all edges */
  for (i = 0; i < (unsigned) n_edges; i++)
  {
    all_edges[i] = EGuGraphAddEdge (G, edgeData + i,
                                    all_nodes[edges[i << 1]],
                                    all_nodes[edges[(i << 1) | 1]]);
  }
  /* compute handle stuff */
  for (i = n_handle; i--;)
  {
    nodeData[handle[i]].in_H = 1U;
    graphData.n_H[1]++;
  }
  graphData.n_H[0] = n_nodes - graphData.n_H[1];
  /* now we update information related to the dominoes */
  for (i = n_dominos; i--;)
  {
    graphData.n_A[i] = n_aset[i];
    graphData.n_B[i] = n_bset[i];
    graphData.n_Tc[i] = n_nodes - n_aset[i] - n_bset[i];
    for (j = n_aset[i]; j--;)
    {
      pos = aset[i][j];
      EGbitSet (nodeData[pos].Ia, i);
      nodeData[pos].n_in_A++;
      nodeData[pos].n_in_T++;
    }
    for (j = n_bset[i]; j--;)
    {
      pos = bset[i][j];
      ADVTESTL (DPT_EDBG, EGbitTest (nodeData[pos].Ia, i), 1, "a node (%d)"
                "is in A_%d and B_%d at the same time", pos, i, i);
      EGbitSet (nodeData[pos].Ib, i);
      nodeData[pos].n_in_T++;
    }
  }
  /* now we update the information for all edges */
  for (i = n_edges; i--;)
  {
    N1_id = all_edges[i]->head->id;
    N2_id = all_edges[i]->tail->id;
    for (j = n_dominos; j--;)
    {
      N1_inA = EGbitTest (nodeData[N1_id].Ia, j);
      N1_inB = EGbitTest (nodeData[N1_id].Ib, j);
      N2_inA = EGbitTest (nodeData[N2_id].Ia, j);
      N2_inB = EGbitTest (nodeData[N2_id].Ib, j);
      N1_inTc = !(N1_inA || N1_inB);
      N2_inTc = !(N2_inA || N2_inB);
      if ((!N1_inTc && N2_inTc) || (N1_inTc && !N2_inTc))
        edgeData[i].n_dT++;
      if ((N1_inA && N2_inB) || (N1_inB && N2_inA))
        edgeData[i].n_AtoB++;
    }
    pos =
      1U & ((nodeData[N1_id].in_H ^ nodeData[N2_id].in_H) ^ edgeData[i].n_AtoB);
    if (pos)
      edgeData[i].in_F = 1U;
  }
#if DPT_VRB <= DEBUG-1
  //DPTDisplayEdges();
#endif
  /* --------------------------------------------------------------------- */
  /* now we need to get the best move for each node */
  MESSAGE (DPT_VRB + 19, "Setting up best move for all nodes");
  /* note that we only need to update moves for nodes whose edges has an edge
   * with coefficient of non-zero */
  for (i = n_edges; i--;)
  {
    /* if the edge has a non-zero coefficient, we update it's both ends if they
     * haven't been update yet. */
    if (edgeData[i].n_dT || edgeData[i].n_AtoB || edgeData[i].in_F)
    {
      if (!EGbitTest (node_update, all_edges[i]->tail->id))
      {
        DPTResetMove (&cur_move);
        cur_move.val = DBL_MAX;
        base_move.val = 0;
        DPTSetBestMove (all_edges[i]->tail->id, &cur_move, &base_move, 0);
        DPTStoreFullMove (all_edges[i]->tail->id, &cur_move);
        EGbitSet (node_update, all_edges[i]->tail->id);
#if DPT_VRB+19<DEBUG
        DPdisplayMove (all_edges[i]->tail->id);
#endif
      }
      if (!EGbitTest (node_update, all_edges[i]->head->id))
      {
        DPTResetMove (&cur_move);
        cur_move.val = DBL_MAX;
        base_move.val = 0;
        DPTSetBestMove (all_edges[i]->head->id, &cur_move, &base_move, 0);
        DPTStoreFullMove (all_edges[i]->head->id, &cur_move);
        EGbitSet (node_update, all_edges[i]->head->id);
#if DPT_VRB+19<DEBUG
        DPdisplayMove (all_edges[i]->head->id);
#endif
      }
    }
  }                             /* end checking all edges */
  /* reset the node update flags */
  memset (node_update, 0, sizeof (node_update));
#if 0
  for (i = n_nodes; i--;)
  {
    DPTResetMove (&cur_move);
    cur_move.val = DBL_MAX;
    base_move.val = 0;
    DPTSetBestMove (i, &cur_move, &base_move, 0);
    DPTStoreFullMove (i, &cur_move);
#if DPT_VRB+19<DEBUG
    DPdisplayMove (i);
#endif
  }
#endif
  MESSAGE (DPT_VRB + 19, "Best move set for all nodes");

  /* --------------------------------------------------------------------- */
  /* now we compute the original violation */
  o_violation = -(3 * n_dominos + 1);
  DPTpriceConstraint (&o_violation);
  MESSAGE (DPT_VRB + 19, "Computing Original violation %lg", o_violation);
  l_violation = o_violation;

  /* --------------------------------------------------------------------- */
  /* now we enter our main loop, while there is a negative move, we keep going
   * */
  while (((cost_tmp = (nc_data =
                       DPTNdata (EGbbtreeMin (tree)->this))->full_move.val)
          < -DPTMinImprovement) && max_iter--)
  {
    MESSAGE (DPT_VRB + 19, "Try to find best move (current best %8.6lf)",
             cost_tmp);
    memcpy (&cur_move, &(nc_data->full_move), sizeof (DPTFullMove_t));
#if DPT_VRB+19<DEBUG
    DPdisplayMove (cur_move.move[0].n_id);
#endif
    /* now check if the movement is infeasible, if so, update the best move and
     * move-on */
    if (!DPTisFullMoveFeasible (&cur_move))
    {
      i = cur_move.move[0].n_id;
#if DPT_VRB+0<DEBUG
      DPdisplayMove (i);
#endif
      MESSAGE (DPT_VRB + 1, "Move for node %u dom %u infeasible", i,
               cur_move.move[0].dom);
      DPTResetMove (&cur_move);
      cur_move.val = DBL_MAX;
      base_move.val = 0;
      DPTSetBestMove (i, &cur_move, &base_move, 0);
      DPTStoreFullMove (i, &cur_move);
#if DPT_VRB+0<DEBUG
      DPdisplayMove (i);
#endif
      continue;
    }
    /* now, if the move is good (i.e. non-zero) we reset the flags */
    if (cost_tmp < DPTMinImprovement)
      DPTresetFlags ();
    /* else we check if it is an adding move, if no such move, check for a flip
     * move, and finally, if no such move, check for a shrink move. Note that
     * the ordering of the moves assures that all zero-valued moves are well
     * ordered. */
    /* now we perform the move */
    MESSAGE (DPT_VRB + 19, "Imp: %8.6lf Tot: %8.6lf node %u move %s depth %u",
             cost_tmp, l_violation, cur_move.move[0].n_id,
             move_name[cur_move.move[0].move_id], nc_data->full_move.depth);
    DPTmakeFullMove (&cur_move);
    l_violation += cost_tmp;
    /* this check is to be sure if the violation predicted is OK */
#if DPT_EDBG <= DEBUG
    EXITRVAL (DPTTestEdges ());
    cost_tmp = -(3 * n_dominos + 1);
    DPTpriceConstraint (&cost_tmp);
    EXIT ((fabs (cost_tmp - l_violation) > 1e-3) &&
          DPdisplayMove (cur_move.move[0].n_id), "Error computing violation"
          ", predicted %8.6lf computed %8.6lf diference %8.6lf",
          l_violation, cost_tmp, l_violation - cost_tmp);
#endif
    /* now update best move for the node and all its neighbours */
    /* update moves */
    DPTupdateNeighMove (&cur_move);
  }
  /* --------------------------------------------------------------------- */
  cost_tmp = ((DPTNdata_t *) (EGbbtreeMin (tree)->this))->full_move.val;
  MESSAGE (DPT_VRB + 19, "No good moves left (current best %8.6lf)",
           (cost_tmp));

  /* --------------------------------------------------------------------- */
  if (o_violation > l_violation + 1e-4)
    /* now, if we have a violated and better inequality, we save it and send it
     * back to the calling function. */
  {
    MESSAGE (DPT_VRB - 1, "Original: %8.6lf New: %8.6lf Diff: %lg", o_violation,
             l_violation, fabs (o_violation - l_violation));
    /* look for memory */
    (*new_n_aset) = EGsMalloc (int,
                               n_dominos);
    (*new_n_bset) = EGsMalloc (int,
                               n_dominos);
    (*new_aset) = EGsMalloc (int *,
                             n_dominos);
    (*new_bset) = EGsMalloc (int *,
                             n_dominos);
    (*new_n_handle) = 0;
    in_H = graphData.n_H[1] < graphData.n_H[0] ? 1 : 0;
    TESTL (DPT_EDBG, !(graphData.n_H[in_H]), "handle is empty!");
    (*new_handle) = EGsMalloc (int,
                               graphData.n_H[in_H]);
    for (i = n_dominos; i--;)
    {
      TESTL (DPT_EDBG, !(graphData.n_A[i]) || !(graphData.n_B[i]) ||
             !(graphData.n_Tc[i]),
             "domino(%u) with empty component A=%u, B=%u, " "T^c=%u", i,
             graphData.n_A[i], graphData.n_B[i], graphData.n_Tc[i]);
      (*new_aset)[i] = EGsMalloc (int,
                                  graphData.n_A[i]);
      (*new_bset)[i] = EGsMalloc (int,
                                  graphData.n_B[i]);
      (*new_n_aset)[i] = 0;
      (*new_n_bset)[i] = 0;
    }

    /* we loop through all nodes and put them where they should go */
    for (i = n_nodes; i--;)
    {
      n_data = nodeData + i;
      if ((n_data->in_H & 1U) == in_H)
        (*new_handle)[(*new_n_handle)++] = i;
      /* loop through all dominos to check to which one the node belong */
      for (j = n_dominos; j--;)
      {
        if (EGbitTest (n_data->Ia, j))
          (*new_aset)[j][(*new_n_aset)[j]++] = i;
        if (EGbitTest (n_data->Ib, j))
          (*new_bset)[j][(*new_n_bset)[j]++] = i;
      }                         /* end loop through dominos for each node */
    }                           /* end loop every node */

    TESTL (DPT_EDBG, (in_H ? ((*new_n_handle) != (signed) graphData.n_H[1]) :
                      ((*new_n_handle) != (signed) graphData.n_H[0])),
           "handle size " "doesn't agree with what we have found");
    for (i = n_dominos; i--;)
    {
      TESTL (DPT_EDBG, (*new_n_aset)[i] != (signed) graphData.n_A[i],
             "size of A_[%u] doesn't match", i);
      TESTL (DPT_EDBG, (*new_n_bset)[i] != (signed) graphData.n_B[i],
             "size of B_[%u] doesn't match", i);
    }
  }                             /* end saving newly found constraint */
  /* --------------------------------------------------------------------- */
  /* ending */
  MESSAGE (DPT_VRB + 19, "Clearing graph");
  EGuGraphClearMP (G, nullFreeMP, nullFreeMP, nullFreeMP, mem, mem, mem);
  MESSAGE (DPT_VRB + 19, "Freing graph");
  EGfreeUGraph (G);
  MESSAGE (DPT_VRB + 19, "Freing all_nodes");
  EGfree (all_nodes);
  MESSAGE (DPT_VRB + 19, "Freing all_edges");
  EGfree (all_edges);
  MESSAGE (DPT_VRB + 19, "Freing tree");
  EGfreeBbtree (tree);
  MESSAGE (DPT_VRB + 19, "Freing nodeData");
  EGfree (nodeData);
  MESSAGE (DPT_VRB + 19, "Freing edgeData");
  EGfree (edgeData);
  MESSAGE (DPT_VRB + 19, "Freing mempool");
  EGfreeMemPool (mem);
  *violation = l_violation;
  MESSAGE (DPT_VRB + 19, "Ending with violation %8.6lf", *violation);
  return rval;
}
#else
int DPtighten (                 /* graph G* data */
                int n_nodes,
                int n_edges,
                int *edges,
                double *weight,
                /* original constrain to tighten */
                int n_dominos,
                int *n_aset,
                int *n_bset,
                int n_handle,
                int **aset,
                int **bset,
                int *handle,
                /* new generated constraint */
                int **new_n_aset,
                int **new_n_bset,
                int *new_n_handle,
                int ***new_aset,
                int ***new_bset,
                int **new_handle,
                double *violation)
{
  register int i;
  *new_n_aset = EGsMalloc (int,
                           n_dominos);
  *new_n_bset = EGsMalloc (int,
                           n_dominos);
  *new_aset = EGsMalloc (int *,
                         n_dominos);
  *new_bset = EGsMalloc (int *,
                         n_dominos);
  *new_n_handle = n_handle;
  *new_handle = EGsMalloc (int,
                           n_handle);
  memcpy (*new_handle, handle, sizeof (int) * n_handle);
  memcpy (*new_n_aset, n_aset, sizeof (int) * n_dominos);
  memcpy (*new_n_bset, n_bset, sizeof (int) * n_dominos);
  for (i = n_dominos; i--;)
  {
    (*new_aset)[i] = EGsMalloc (int,
                                n_aset[i]);
    (*new_bset)[i] = EGsMalloc (int,
                                n_bset[i]);
    memcpy (((*new_aset)[i]), (aset[i]), sizeof (int) * n_aset[i]);
    memcpy (((*new_bset)[i]), (bset[i]), sizeof (int) * n_bset[i]);
  }
  *violation = 0;
  return 0;
}
#endif

---