glpios01.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

                  best = node;
            break;
         default:
            xassert(tree != tree);
      }
      return best == NULL ? 0 : best->p;
}

/***********************************************************************
*  NAME
*
*  ios_relative_gap - compute relative mip gap
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  double ios_relative_gap(glp_tree *tree);
*
*  DESCRIPTION
*
*  The routine ios_relative_gap computes the relative mip gap using the
*  formula:
*
*     gap = |best_mip - best_bnd| / (|best_mip| + DBL_EPSILON),
*
*  where best_mip is the best integer feasible solution found so far,
*  best_bnd is the best (global) bound. If no integer feasible solution
*  has been found yet, rel_gap is set to DBL_MAX.
*
*  RETURNS
*
*  The routine ios_relative_gap returns the relative mip gap. */

double ios_relative_gap(glp_tree *tree)
{     glp_prob *mip = tree->mip;
      int p;
      double best_mip, best_bnd, gap;
      if (mip->mip_stat == GLP_FEAS)
      {  best_mip = mip->mip_obj;
         p = ios_best_node(tree);
         if (p == 0)
         {  /* the tree is empty */
            gap = 0.0;
         }
         else
         {  best_bnd = tree->slot[p].node->bound;
            gap = fabs(best_mip - best_bnd) / (fabs(best_mip) +
               DBL_EPSILON);
         }
      }
      else
      {  /* no integer feasible solution has been found yet */
         gap = DBL_MAX;
      }
      return gap;
}

/***********************************************************************
*  NAME
*
*  ios_solve_node - solve LP relaxation of current subproblem
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  int ios_solve_node(glp_tree *tree);
*
*  DESCRIPTION
*
*  The routine ios_solve_node re-optimizes LP relaxation of the current
*  subproblem using the dual simplex method.
*
*  RETURNS
*
*  The routine returns the code which is reported by glp_simplex. */

int ios_solve_node(glp_tree *tree)
{     glp_prob *mip = tree->mip;
      glp_smcp parm;
      int ret;
      /* the current subproblem must exist */
      xassert(tree->curr != NULL);
      /* set some control parameters */
      glp_init_smcp(&parm);
      switch (tree->parm->msg_lev)
      {  case GLP_MSG_OFF:
            parm.msg_lev = GLP_MSG_OFF; break;
         case GLP_MSG_ERR:
            parm.msg_lev = GLP_MSG_ERR; break;
         case GLP_MSG_ON:
         case GLP_MSG_ALL:
            parm.msg_lev = GLP_MSG_ON; break;
         case GLP_MSG_DBG:
            parm.msg_lev = GLP_MSG_ALL; break;
         default:
            xassert(tree != tree);
      }
      parm.meth = GLP_DUALP;
      if (tree->parm->msg_lev < GLP_MSG_DBG)
         parm.out_dly = tree->parm->out_dly;
      else
         parm.out_dly = 0;
      /* if the incumbent objective value is already known, use it to
         prematurely terminate the dual simplex search */
      if (mip->mip_stat == GLP_FEAS)
      {  switch (tree->mip->dir)
         {  case GLP_MIN:
               parm.obj_ul = mip->mip_obj;
               break;
            case GLP_MAX:
               parm.obj_ll = mip->mip_obj;
               break;
            default:
               xassert(mip != mip);
         }
      }
      /* try to solve/re-optimize the LP relaxation */
      ret = glp_simplex(mip, &parm);
#if 0
      xprintf("ret = %d; status = %d; pbs = %d; dbs = %d; some = %d\n",
         ret, glp_get_status(mip), mip->pbs_stat, mip->dbs_stat,
         mip->some);
      lpx_print_sol(mip, "sol");
#endif
      return ret;
}

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

IOSPOOL *ios_create_pool(glp_tree *tree)
{     /* create cut pool */
      IOSPOOL *pool;
#if 0
      pool = dmp_get_atom(tree->pool, sizeof(IOSPOOL));
#else
      xassert(tree == tree);
      pool = xmalloc(sizeof(IOSPOOL));
#endif
      pool->size = 0;
      pool->head = pool->tail = NULL;
      pool->ord = 0, pool->curr = NULL;
      return pool;
}

int ios_add_row(glp_tree *tree, IOSPOOL *pool,
      const char *name, int klass, int flags, int len, const int ind[],
      const double val[], int type, double rhs)
{     /* add row (constraint) to the cut pool */
      IOSCUT *cut;
      IOSAIJ *aij;
      int k;
      xassert(pool != NULL);
      cut = dmp_get_atom(tree->pool, sizeof(IOSCUT));
      if (name == NULL || name[0] == '\0')
         cut->name = NULL;
      else
      {  for (k = 0; name[k] != '\0'; k++)
         {  if (k == 256)
               xerror("glp_ios_add_row: cut name too long\n");
            if (iscntrl((unsigned char)name[k]))
               xerror("glp_ios_add_row: cut name contains invalid chara"
                  "cter(s)\n");
         }
         cut->name = dmp_get_atom(tree->pool, strlen(name)+1);
         strcpy(cut->name, name);
      }
      if (!(0 <= klass && klass <= 255))
         xerror("glp_ios_add_row: klass = %d; invalid cut class\n",
            klass);
      cut->klass = (unsigned char)klass;
      if (flags != 0)
         xerror("glp_ios_add_row: flags = %d; invalid cut flags\n",
            flags);
      cut->ptr = NULL;
      if (!(0 <= len && len <= tree->n))
         xerror("glp_ios_add_row: len = %d; invalid cut length\n",
            len);
      for (k = 1; k <= len; k++)
      {  aij = dmp_get_atom(tree->pool, sizeof(IOSAIJ));
         if (!(1 <= ind[k] && ind[k] <= tree->n))
            xerror("glp_ios_add_row: ind[%d] = %d; column index out of "
               "range\n", k, ind[k]);
         aij->j = ind[k];
         aij->val = val[k];
         aij->next = cut->ptr;
         cut->ptr = aij;
      }
      if (!(type == GLP_LO || type == GLP_UP || type == GLP_FX))
         xerror("glp_ios_add_row: type = %d; invalid cut type\n",
            type);
      cut->type = type;
      cut->rhs = rhs;
      cut->prev = pool->tail;
      cut->next = NULL;
      if (cut->prev == NULL)
         pool->head = cut;
      else
         cut->prev->next = cut;
      pool->tail = cut;
      pool->size++;
      return pool->size;
}

IOSCUT *ios_find_row(IOSPOOL *pool, int i)
{     /* find row (constraint) in the cut pool */
      /* (smart linear search) */
      xassert(pool != NULL);
      xassert(1 <= i && i <= pool->size);
      if (pool->ord == 0)
      {  xassert(pool->curr == NULL);
         pool->ord = 1;
         pool->curr = pool->head;
      }
      xassert(pool->curr != NULL);
      if (i < pool->ord)
      {  if (i < pool->ord - i)
         {  pool->ord = 1;
            pool->curr = pool->head;
            while (pool->ord != i)
            {  pool->ord++;
               xassert(pool->curr != NULL);
               pool->curr = pool->curr->next;
            }
         }
         else
         {  while (pool->ord != i)
            {  pool->ord--;
               xassert(pool->curr != NULL);
               pool->curr = pool->curr->prev;
            }
         }
      }
      else if (i > pool->ord)
      {  if (i - pool->ord < pool->size - i)
         {  while (pool->ord != i)
            {  pool->ord++;
               xassert(pool->curr != NULL);
               pool->curr = pool->curr->next;
            }
         }
         else
         {  pool->ord = pool->size;
            pool->curr = pool->tail;
            while (pool->ord != i)
            {  pool->ord--;
               xassert(pool->curr != NULL);
               pool->curr = pool->curr->prev;
            }
         }
      }
      xassert(pool->ord == i);
      xassert(pool->curr != NULL);
      return pool->curr;
}

void ios_del_row(glp_tree *tree, IOSPOOL *pool, int i)
{     /* remove row (constraint) from the cut pool */
      IOSCUT *cut;
      IOSAIJ *aij;
      xassert(pool != NULL);
      if (!(1 <= i && i <= pool->size))
         xerror("glp_ios_del_row: i = %d; cut number out of range\n",
            i);
      cut = ios_find_row(pool, i);
      xassert(pool->curr == cut);
      if (cut->next != NULL)
         pool->curr = cut->next;
      else if (cut->prev != NULL)
         pool->ord--, pool->curr = cut->prev;
      else
         pool->ord = 0, pool->curr = NULL;
      if (cut->name != NULL)
         dmp_free_atom(tree->pool, cut->name, strlen(cut->name)+1);
      if (cut->prev == NULL)
      {  xassert(pool->head == cut);
         pool->head = cut->next;
      }
      else
      {  xassert(cut->prev->next == cut);
         cut->prev->next = cut->next;
      }
      if (cut->next == NULL)
      {  xassert(pool->tail == cut);
         pool->tail = cut->prev;
      }
      else
      {  xassert(cut->next->prev == cut);
         cut->next->prev = cut->prev;
      }
      while (cut->ptr != NULL)
      {  aij = cut->ptr;
         cut->ptr = aij->next;
         dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ));
      }
      dmp_free_atom(tree->pool, cut, sizeof(IOSCUT));
      pool->size--;
      return;
}

void ios_clear_pool(glp_tree *tree, IOSPOOL *pool)
{     /* remove all rows (constraints) from the cut pool */
      xassert(pool != NULL);
      while (pool->head != NULL)
      {  IOSCUT *cut = pool->head;
         pool->head = cut->next;
         if (cut->name != NULL)
            dmp_free_atom(tree->pool, cut->name, strlen(cut->name)+1);
         while (cut->ptr != NULL)
         {  IOSAIJ *aij = cut->ptr;
            cut->ptr = aij->next;
            dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ));
         }
         dmp_free_atom(tree->pool, cut, sizeof(IOSCUT));
      }
      pool->size = 0;
      pool->head = pool->tail = NULL;
      pool->ord = 0, pool->curr = NULL;
      return;
}

void ios_delete_pool(glp_tree *tree, IOSPOOL *pool)
{     /* delete cut pool */
      xassert(pool != NULL);
      ios_clear_pool(tree, pool);
      xfree(pool);
      return;
}

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

static int refer_to_node(glp_tree *tree, int j)
{     /* determine node number corresponding to binary variable x[j] or
         its complement */
      glp_prob *mip = tree->mip;
      int n = mip->n;
      int *ref;
      if (j > 0)
         ref = tree->n_ref;
      else
         ref = tree->c_ref, j = - j;
      xassert(1 <= j && j <= n);
      if (ref[j] == 0)
      {  /* new node is needed */
         SCG *g = tree->g;
         int n_max = g->n_max;
         ref[j] = scg_add_nodes(g, 1);
         if (g->n_max > n_max)
         {  int *save = tree->j_ref;
            tree->j_ref = xcalloc(1+g->n_max, sizeof(int));
            memcpy(&tree->j_ref[1], &save[1], g->n * sizeof(int));
            xfree(save);
         }
         xassert(ref[j] == g->n);
         tree->j_ref[ref[j]] = j;
         xassert(tree->curr != NULL);
         if (tree->curr->level > 0) tree->curr->own_nn++;
      }
      return ref[j];
}

void ios_add_edge(glp_tree *tree, int j1, int j2)
{     /* add new edge to the conflict graph */
      glp_prob *mip = tree->mip;
      int n = mip->n;
      SCGRIB *e;
      int first, i1, i2;
      xassert(-n <= j1 && j1 <= +n && j1 != 0);
      xassert(-n <= j2 && j2 <= +n && j2 != 0);
      xassert(j1 != j2);
      /* determine number of the first node, which was added for the
         current subproblem */
      xassert(tree->curr != NULL);
      first = tree->g->n - tree->curr->own_nn + 1;
      /* determine node numbers for both endpoints */
      i1 = refer_to_node(tree, j1);
      i2 = refer_to_node(tree, j2);
      /* add edge (i1,i2) to the conflict graph */
      e = scg_add_edge(tree->g, i1, i2);
      /* if the current subproblem is not the root and both endpoints
         were created on some previous levels, save the edge */
      if (tree->curr->level > 0 && i1 < first && i2 < first)
      {  IOSRIB *rib;
         rib = dmp_get_atom(tree->pool, sizeof(IOSRIB));
         rib->j1 = j1;
         rib->j2 = j2;
         rib->e = e;
         rib->next = tree->curr->e_ptr;
         tree->curr->e_ptr = rib;
      }
      return;
}

/* eof */