glpios03.c

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

            else
               app++;
         }
      }
      if (gmi + mir + cov + clq + app > 0)
      {  if (gmi > 0) xprintf(" gmi:%3d", gmi);
         if (mir > 0) xprintf(" mir:%3d", mir);
         if (cov > 0) xprintf(" cov:%3d", cov);
         if (clq > 0) xprintf(" clq:%3d", clq);
         if (app > 0) xprintf(" app:%3d", app);
         xprintf("\n");
      }
      return;
}

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

/***********************************************************************
*  NAME
*
*  mip_driver - branch-and-bound driver
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  int ios_driver(glp_tree *tree);
*
*  DESCRIPTION
*
*  The routine mip_driver is a branch-and-bound driver that manages the
*  process of solving MIP problem instance.
*
*  RETURNS
*
*  0  The MIP problem instance has been successfully solved. This code
*     does not necessarily mean that the solver has found optimal
*     solution. It only means that the solution process was successful.
*
*  GLP_EFAIL
*     The search was prematurely terminated due to the solver failure.
*
*  GLP_EMIPGAP
*     The search was prematurely terminated, because the relative mip
*     gap tolerance has been reached.
*
*  GLP_ETMLIM
*     The search was prematurely terminated, because the time limit has
*     been exceeded.
*
*  GLP_ESTOP
*     The search was prematurely terminated by application. */

static void generate_cuts(glp_tree *tree);
static void separation(glp_tree *tree);

int ios_driver(glp_tree *tree)
{     glp_prob *mip = tree->mip;
      int p, p_stat, d_stat, ret;
      xlong_t ttt = tree->tm_beg;
      /* on entry to the B&B driver it is assumed that the active list
         contains the only active (i.e. root) subproblem, which is the
         original MIP problem to be solved */
back: /* at this point the current subproblem does not exist */
      xassert(tree->curr == NULL);
      /* if the active list is empty, the search is finished */
      if (tree->head == NULL)
      {  if (tree->parm->msg_lev >= GLP_MSG_DBG)
            xprintf("Active list is empty!\n");
         xassert(dmp_in_use(tree->pool).lo == 0);
         ret = 0;
         goto done;
      }
#if 1
      tree->just_selected = 1;
#endif
      /* select some active subproblem to be solved next and make it
         current */
      if (tree->parm->cb_func != NULL)
      {  xassert(tree->reason == 0);
         xassert(tree->btrack == NULL);
         tree->reason = GLP_ISELECT;
         tree->parm->cb_func(tree, tree->parm->cb_info);
         tree->reason = 0;
         if (tree->terminate)
         {  ret = GLP_ESTOP;
            goto done;
         }
         if (tree->btrack != NULL)
         {  /* revive the subproblem selected */
            ios_revive_node(tree, tree->btrack->p);
            tree->btrack = NULL;
            goto loop;
         }
      }
      /****************************************************************/
      if (tree->a_cnt == 1)
      {  xassert(tree->head->next == NULL);
         ios_revive_node(tree, tree->head->p);
      }
      else
      switch (tree->parm->bt_tech)
      {  case GLP_BT_DFS:
            /* depth first search */
            xassert(tree->tail != NULL);
            ios_revive_node(tree, tree->tail->p);
            break;
         case GLP_BT_BFS:
            /* breadth first search */
            xassert(tree->head != NULL);
            ios_revive_node(tree, tree->head->p);
            break;
         case GLP_BT_BLB:
            /* select node with best local bound */
#if 0
            ios_revive_node(tree, ios_best_node(tree));
#else
            btrack_best_node(tree);
#endif
            break;
         case GLP_BT_BPH:
            if (mip->mip_stat == GLP_UNDEF)
            {  /* "most integer feasible" subproblem */
               btrack_most_feas(tree);
            }
            else
            {  /* best projection heuristic */
               btrack_best_proj(tree);
            }
            break;
         default:
            xassert(tree != tree);
      }
#if 0
      /* revive the root subproblem */
      ios_revive_node(tree, 1);
#endif
loop: /* main loop starts here; at this point some subproblem has been
         just selected from the active list and made current */
      xassert(tree->curr != NULL);
      xassert(tree->solved == 0);
#if 1
      tree->round = 0;
#endif
      /* determine the reference number of the current subproblem */
      p = tree->curr->p;
      if (tree->parm->msg_lev >= GLP_MSG_DBG)
      {  int level = tree->slot[p].node->level;
         xprintf("-----------------------------------------------------"
            "-------------------\n");
         xprintf("Processing node %d at level %d\n", p, level);
      }
#if 1
      /* if it is the root subproblem, initialize cut generators */
      if (p == 1)
      {  if (tree->parm->gmi_cuts == GLP_ON)
         {  if (tree->parm->msg_lev >= GLP_MSG_ALL)
               xprintf("Gomory's cuts enabled\n");
         }
         if (tree->parm->mir_cuts == GLP_ON)
         {  if (tree->parm->msg_lev >= GLP_MSG_ALL)
               xprintf("MIR cuts enabled\n");
            xassert(tree->mir_gen == NULL);
            tree->mir_gen = ios_mir_init(tree);
         }
         if (tree->parm->cov_cuts == GLP_ON)
         {  if (tree->parm->msg_lev >= GLP_MSG_ALL)
               xprintf("Cover cuts enabled\n");
         }
         if (tree->parm->clq_cuts == GLP_ON)
         {  xassert(tree->clq_gen == NULL);
            if (tree->parm->msg_lev >= GLP_MSG_ALL)
               xprintf("Clique cuts enabled\n");
            tree->clq_gen = ios_clq_init(tree);
         }
      }
#endif
#if 1 /* 09/X-2008 */
more: /* minor loop starts here; at this point the current subproblem
         is either to be solved for the first time or to be reoptimized
         due to reformulation */
      /* display current progress of the search */
      if (tree->parm->msg_lev >= GLP_MSG_DBG ||
          tree->parm->msg_lev >= GLP_MSG_ON &&
        (double)(tree->parm->out_frq - 1) <=
            1000.0 * xdifftime(xtime(), tree->tm_lag))
         show_progress(tree, 0);
      if (tree->parm->msg_lev >= GLP_MSG_ALL &&
            xdifftime(xtime(), ttt) >= 60.0)
      {  xlong_t total;
         lib_mem_usage(NULL, NULL, &total, NULL);
         xprintf("Time used: %.1f secs.  Memory used: %.1f Mb.\n",
            xdifftime(xtime(), tree->tm_beg), xltod(total) / 1048576.0);
         ttt = xtime();
      }
#endif
#if 1
      /* check mip gap */
      if (tree->parm->mip_gap > 0.0 &&
          ios_relative_gap(tree) <= tree->parm->mip_gap)
      {  if (tree->parm->msg_lev >= GLP_MSG_DBG)
            xprintf("Relative gap tolerance reached; search terminated "
               "\n");
         ret = GLP_EMIPGAP;
         goto done;
      }
#endif
      /* check if the time limit has been exhausted */
      if (tree->parm->tm_lim < INT_MAX &&
         (double)(tree->parm->tm_lim - 1) <=
         1000.0 * xdifftime(xtime(), tree->tm_beg))
      {  if (tree->parm->msg_lev >= GLP_MSG_DBG)
            xprintf("Time limit exhausted; search terminated\n");
         ret = GLP_ETMLIM;
         goto done;
      }
#if 1
      /* perform basic preprocessing */
      switch (tree->parm->pp_tech)
      {  case GLP_PP_NONE:
            break;
         case GLP_PP_ROOT:
            if (tree->curr->level == 0)
            {  if (ios_preprocess_node(tree, 99))
                  goto fath;
            }
            break;
         case GLP_PP_ALL:
            if (ios_preprocess_node(tree,
               tree->curr->level == 0 ? 99 : 10)) goto fath;
            break;
         default:
            xassert(tree != tree);
      }
#endif
      /* call the user-defined callback routine for preprocessing */
      if (tree->parm->cb_func != NULL)
      {  xassert(tree->reason == 0);
         tree->reason = GLP_IPREPRO;
         tree->parm->cb_func(tree, tree->parm->cb_info);
         tree->reason = 0;
         if (tree->terminate)
         {  ret = GLP_ESTOP;
            goto done;
         }
      }
#if 1
      if (!is_curr_node_hopeful(tree))
      {  /*xprintf("fathomed by preprocessing\n");*/
         goto fath;
      }
#endif
#if 0 /* 09/X-2008 */
more: /* minor loop starts here; at this point the current subproblem
         is either to be solved for the first time or to be reoptimized
         due to reformulation */
      /* display current progress of the search */
      if (tree->parm->msg_lev >= GLP_MSG_DBG ||
          tree->parm->msg_lev >= GLP_MSG_ON &&
        (double)(tree->parm->out_frq - 1) <=
            1000.0 * xdifftime(xtime(), tree->tm_lag))
         show_progress(tree, 0);
      if (tree->parm->msg_lev >= GLP_MSG_ALL &&
            xdifftime(xtime(), ttt) >= 60.0)
      {  xlong_t total;
         lib_mem_usage(NULL, NULL, &total, NULL);
         xprintf("Time used: %.1f secs.  Memory used: %.1f Mb.\n",
            xdifftime(xtime(), tree->tm_beg), xltod(total) / 1048576.0);
         ttt = xtime();
      }
#endif
      /* solve LP relaxation of the current subproblem */
      if (tree->parm->msg_lev >= GLP_MSG_DBG)
         xprintf("Solving LP relaxation...\n");
      ret = ios_solve_node(tree);
      tree->solved++;
      switch (ret)
      {  case 0:
         case GLP_EOBJLL:
         case GLP_EOBJUL:
            break;
         default:
            if (tree->parm->msg_lev >= GLP_MSG_ERR)
               xprintf("ios_driver: unable to solve current LP relaxati"
                  "on; glp_simplex returned %d\n", ret);
            ret = GLP_EFAIL;
            goto done;
      }
      /* analyze status of the basic solution */
      p_stat = mip->pbs_stat;
      d_stat = mip->dbs_stat;
      if (p_stat == GLP_FEAS && d_stat == GLP_FEAS)
      {  /* LP relaxation has optimal solution */
         if (tree->parm->msg_lev >= GLP_MSG_DBG)
            xprintf("Found optimal solution to LP relaxation\n");
      }
#if 0 /* 19/VIII-2008 */
      else if (p_stat == GLP_FEAS && d_stat == GLP_NOFEAS)
      {  /* LP relaxation has unbounded solution */
#else
      else if (d_stat == GLP_NOFEAS)
      {  /* LP relaxation has no dual feasible solution */
#endif
         /* since the current subproblem cannot have a larger feasible
            region than its parent, there is something wrong */
         if (tree->parm->msg_lev >= GLP_MSG_ERR)
#if 0 /* 19/VIII-2008 */
            xprintf("ios_driver: current LP relaxation has unbounded so"
               "lution\n");
#else
            xprintf("ios_driver: current LP relaxation has no dual feas"
               "ible solution\n");
#endif
         ret = GLP_EFAIL;
         goto done;
      }
      else if (p_stat == GLP_INFEAS && d_stat == GLP_FEAS)
      {  /* LP relaxation has no primal solution which is better than
            the incumbent objective value */
         xassert(mip->mip_stat == GLP_FEAS);
         if (tree->parm->msg_lev >= GLP_MSG_DBG)
            xprintf("LP relaxation has no solution better than incumben"
               "t objective value\n");
         /* prune the branch */
         goto fath;
      }
      else if (p_stat == GLP_NOFEAS)
      {  /* LP relaxation has no primal feasible solution */
         if (tree->parm->msg_lev >= GLP_MSG_DBG)
            xprintf("LP relaxation has no feasible solution\n");
         /* prune the branch */
         goto fath;
      }
      else
      {  /* other cases cannot appear */
         xassert(mip != mip);
      }
      /* at this point basic solution of LP relaxation of the current
         subproblem is optimal */
      xassert(p_stat == GLP_FEAS && d_stat == GLP_FEAS);
#if 1 /* 04/X-2008 */
      xassert(tree->curr != NULL);
      tree->curr->lp_obj = tree->mip->obj_val;
#endif
      /* thus, it defines a local bound for integer optimal solution of
         the current subproblem */
#if 0 /* 12/X-2008 */
      set_local_bound(tree);
#else
      {  double bound = tree->mip->obj_val;