src/glpios03.c
author Alpar Juttner <alpar@cs.elte.hu>
Mon, 06 Dec 2010 13:09:21 +0100
changeset 1 c445c931472f
permissions -rw-r--r--
Import glpk-4.45

- Generated files and doc/notes are removed
alpar@1
     1
/* glpios03.c (branch-and-cut driver) */
alpar@1
     2
alpar@1
     3
/***********************************************************************
alpar@1
     4
*  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1
     5
*
alpar@1
     6
*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1
     7
*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1
     8
*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1
     9
*  E-mail: <mao@gnu.org>.
alpar@1
    10
*
alpar@1
    11
*  GLPK is free software: you can redistribute it and/or modify it
alpar@1
    12
*  under the terms of the GNU General Public License as published by
alpar@1
    13
*  the Free Software Foundation, either version 3 of the License, or
alpar@1
    14
*  (at your option) any later version.
alpar@1
    15
*
alpar@1
    16
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1
    17
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1
    18
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1
    19
*  License for more details.
alpar@1
    20
*
alpar@1
    21
*  You should have received a copy of the GNU General Public License
alpar@1
    22
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1
    23
***********************************************************************/
alpar@1
    24
alpar@1
    25
#include "glpios.h"
alpar@1
    26
alpar@1
    27
/***********************************************************************
alpar@1
    28
*  show_progress - display current progress of the search
alpar@1
    29
*
alpar@1
    30
*  This routine displays some information about current progress of the
alpar@1
    31
*  search.
alpar@1
    32
*
alpar@1
    33
*  The information includes:
alpar@1
    34
*
alpar@1
    35
*  the current number of iterations performed by the simplex solver;
alpar@1
    36
*
alpar@1
    37
*  the objective value for the best known integer feasible solution,
alpar@1
    38
*  which is upper (minimization) or lower (maximization) global bound
alpar@1
    39
*  for optimal solution of the original mip problem;
alpar@1
    40
*
alpar@1
    41
*  the best local bound for active nodes, which is lower (minimization)
alpar@1
    42
*  or upper (maximization) global bound for optimal solution of the
alpar@1
    43
*  original mip problem;
alpar@1
    44
*
alpar@1
    45
*  the relative mip gap, in percents;
alpar@1
    46
*
alpar@1
    47
*  the number of open (active) subproblems;
alpar@1
    48
*
alpar@1
    49
*  the number of completely explored subproblems, i.e. whose nodes have
alpar@1
    50
*  been removed from the tree. */
alpar@1
    51
alpar@1
    52
static void show_progress(glp_tree *T, int bingo)
alpar@1
    53
{     int p;
alpar@1
    54
      double temp;
alpar@1
    55
      char best_mip[50], best_bound[50], *rho, rel_gap[50];
alpar@1
    56
      /* format the best known integer feasible solution */
alpar@1
    57
      if (T->mip->mip_stat == GLP_FEAS)
alpar@1
    58
         sprintf(best_mip, "%17.9e", T->mip->mip_obj);
alpar@1
    59
      else
alpar@1
    60
         sprintf(best_mip, "%17s", "not found yet");
alpar@1
    61
      /* determine reference number of an active subproblem whose local
alpar@1
    62
         bound is best */
alpar@1
    63
      p = ios_best_node(T);
alpar@1
    64
      /* format the best bound */
alpar@1
    65
      if (p == 0)
alpar@1
    66
         sprintf(best_bound, "%17s", "tree is empty");
alpar@1
    67
      else
alpar@1
    68
      {  temp = T->slot[p].node->bound;
alpar@1
    69
         if (temp == -DBL_MAX)
alpar@1
    70
            sprintf(best_bound, "%17s", "-inf");
alpar@1
    71
         else if (temp == +DBL_MAX)
alpar@1
    72
            sprintf(best_bound, "%17s", "+inf");
alpar@1
    73
         else
alpar@1
    74
            sprintf(best_bound, "%17.9e", temp);
alpar@1
    75
      }
alpar@1
    76
      /* choose the relation sign between global bounds */
alpar@1
    77
      if (T->mip->dir == GLP_MIN)
alpar@1
    78
         rho = ">=";
alpar@1
    79
      else if (T->mip->dir == GLP_MAX)
alpar@1
    80
         rho = "<=";
alpar@1
    81
      else
alpar@1
    82
         xassert(T != T);
alpar@1
    83
      /* format the relative mip gap */
alpar@1
    84
      temp = ios_relative_gap(T);
alpar@1
    85
      if (temp == 0.0)
alpar@1
    86
         sprintf(rel_gap, "  0.0%%");
alpar@1
    87
      else if (temp < 0.001)
alpar@1
    88
         sprintf(rel_gap, "< 0.1%%");
alpar@1
    89
      else if (temp <= 9.999)
alpar@1
    90
         sprintf(rel_gap, "%5.1f%%", 100.0 * temp);
alpar@1
    91
      else
alpar@1
    92
         sprintf(rel_gap, "%6s", "");
alpar@1
    93
      /* display progress of the search */
alpar@1
    94
      xprintf("+%6d: %s %s %s %s %s (%d; %d)\n",
alpar@1
    95
         T->mip->it_cnt, bingo ? ">>>>>" : "mip =", best_mip, rho,
alpar@1
    96
         best_bound, rel_gap, T->a_cnt, T->t_cnt - T->n_cnt);
alpar@1
    97
      T->tm_lag = xtime();
alpar@1
    98
      return;
alpar@1
    99
}
alpar@1
   100
alpar@1
   101
/***********************************************************************
alpar@1
   102
*  is_branch_hopeful - check if specified branch is hopeful
alpar@1
   103
*
alpar@1
   104
*  This routine checks if the specified subproblem can have an integer
alpar@1
   105
*  optimal solution which is better than the best known one.
alpar@1
   106
*
alpar@1
   107
*  The check is based on comparison of the local objective bound stored
alpar@1
   108
*  in the subproblem descriptor and the incumbent objective value which
alpar@1
   109
*  is the global objective bound.
alpar@1
   110
*
alpar@1
   111
*  If there is a chance that the specified subproblem can have a better
alpar@1
   112
*  integer optimal solution, the routine returns non-zero. Otherwise, if
alpar@1
   113
*  the corresponding branch can pruned, zero is returned. */
alpar@1
   114
alpar@1
   115
static int is_branch_hopeful(glp_tree *T, int p)
alpar@1
   116
{     xassert(1 <= p && p <= T->nslots);
alpar@1
   117
      xassert(T->slot[p].node != NULL);
alpar@1
   118
      return ios_is_hopeful(T, T->slot[p].node->bound);
alpar@1
   119
}
alpar@1
   120
alpar@1
   121
/***********************************************************************
alpar@1
   122
*  check_integrality - check integrality of basic solution
alpar@1
   123
*
alpar@1
   124
*  This routine checks if the basic solution of LP relaxation of the
alpar@1
   125
*  current subproblem satisfies to integrality conditions, i.e. that all
alpar@1
   126
*  variables of integer kind have integral primal values. (The solution
alpar@1
   127
*  is assumed to be optimal.)
alpar@1
   128
*
alpar@1
   129
*  For each variable of integer kind the routine computes the following
alpar@1
   130
*  quantity:
alpar@1
   131
*
alpar@1
   132
*     ii(x[j]) = min(x[j] - floor(x[j]), ceil(x[j]) - x[j]),         (1)
alpar@1
   133
*
alpar@1
   134
*  which is a measure of the integer infeasibility (non-integrality) of
alpar@1
   135
*  x[j] (for example, ii(2.1) = 0.1, ii(3.7) = 0.3, ii(5.0) = 0). It is
alpar@1
   136
*  understood that 0 <= ii(x[j]) <= 0.5, and variable x[j] is integer
alpar@1
   137
*  feasible if ii(x[j]) = 0. However, due to floating-point arithmetic
alpar@1
   138
*  the routine checks less restrictive condition:
alpar@1
   139
*
alpar@1
   140
*     ii(x[j]) <= tol_int,                                           (2)
alpar@1
   141
*
alpar@1
   142
*  where tol_int is a given tolerance (small positive number) and marks
alpar@1
   143
*  each variable which does not satisfy to (2) as integer infeasible by
alpar@1
   144
*  setting its fractionality flag.
alpar@1
   145
*
alpar@1
   146
*  In order to characterize integer infeasibility of the basic solution
alpar@1
   147
*  in the whole the routine computes two parameters: ii_cnt, which is
alpar@1
   148
*  the number of variables with the fractionality flag set, and ii_sum,
alpar@1
   149
*  which is the sum of integer infeasibilities (1). */
alpar@1
   150
alpar@1
   151
static void check_integrality(glp_tree *T)
alpar@1
   152
{     glp_prob *mip = T->mip;
alpar@1
   153
      int j, type, ii_cnt = 0;
alpar@1
   154
      double lb, ub, x, temp1, temp2, ii_sum = 0.0;
alpar@1
   155
      /* walk through the set of columns (structural variables) */
alpar@1
   156
      for (j = 1; j <= mip->n; j++)
alpar@1
   157
      {  GLPCOL *col = mip->col[j];
alpar@1
   158
         T->non_int[j] = 0;
alpar@1
   159
         /* if the column is not integer, skip it */
alpar@1
   160
         if (col->kind != GLP_IV) continue;
alpar@1
   161
         /* if the column is non-basic, it is integer feasible */
alpar@1
   162
         if (col->stat != GLP_BS) continue;
alpar@1
   163
         /* obtain the type and bounds of the column */
alpar@1
   164
         type = col->type, lb = col->lb, ub = col->ub;
alpar@1
   165
         /* obtain value of the column in optimal basic solution */
alpar@1
   166
         x = col->prim;
alpar@1
   167
         /* if the column's primal value is close to the lower bound,
alpar@1
   168
            the column is integer feasible within given tolerance */
alpar@1
   169
         if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
alpar@1
   170
         {  temp1 = lb - T->parm->tol_int;
alpar@1
   171
            temp2 = lb + T->parm->tol_int;
alpar@1
   172
            if (temp1 <= x && x <= temp2) continue;
alpar@1
   173
#if 0
alpar@1
   174
            /* the lower bound must not be violated */
alpar@1
   175
            xassert(x >= lb);
alpar@1
   176
#else
alpar@1
   177
            if (x < lb) continue;
alpar@1
   178
#endif
alpar@1
   179
         }
alpar@1
   180
         /* if the column's primal value is close to the upper bound,
alpar@1
   181
            the column is integer feasible within given tolerance */
alpar@1
   182
         if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
alpar@1
   183
         {  temp1 = ub - T->parm->tol_int;
alpar@1
   184
            temp2 = ub + T->parm->tol_int;
alpar@1
   185
            if (temp1 <= x && x <= temp2) continue;
alpar@1
   186
#if 0
alpar@1
   187
            /* the upper bound must not be violated */
alpar@1
   188
            xassert(x <= ub);
alpar@1
   189
#else
alpar@1
   190
            if (x > ub) continue;
alpar@1
   191
#endif
alpar@1
   192
         }
alpar@1
   193
         /* if the column's primal value is close to nearest integer,
alpar@1
   194
            the column is integer feasible within given tolerance */
alpar@1
   195
         temp1 = floor(x + 0.5) - T->parm->tol_int;
alpar@1
   196
         temp2 = floor(x + 0.5) + T->parm->tol_int;
alpar@1
   197
         if (temp1 <= x && x <= temp2) continue;
alpar@1
   198
         /* otherwise the column is integer infeasible */
alpar@1
   199
         T->non_int[j] = 1;
alpar@1
   200
         /* increase the number of fractional-valued columns */
alpar@1
   201
         ii_cnt++;
alpar@1
   202
         /* compute the sum of integer infeasibilities */
alpar@1
   203
         temp1 = x - floor(x);
alpar@1
   204
         temp2 = ceil(x) - x;
alpar@1
   205
         xassert(temp1 > 0.0 && temp2 > 0.0);
alpar@1
   206
         ii_sum += (temp1 <= temp2 ? temp1 : temp2);
alpar@1
   207
      }
alpar@1
   208
      /* store ii_cnt and ii_sum to the current problem descriptor */
alpar@1
   209
      xassert(T->curr != NULL);
alpar@1
   210
      T->curr->ii_cnt = ii_cnt;
alpar@1
   211
      T->curr->ii_sum = ii_sum;
alpar@1
   212
      /* and also display these parameters */
alpar@1
   213
      if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   214
      {  if (ii_cnt == 0)
alpar@1
   215
            xprintf("There are no fractional columns\n");
alpar@1
   216
         else if (ii_cnt == 1)
alpar@1
   217
            xprintf("There is one fractional column, integer infeasibil"
alpar@1
   218
               "ity is %.3e\n", ii_sum);
alpar@1
   219
         else
alpar@1
   220
            xprintf("There are %d fractional columns, integer infeasibi"
alpar@1
   221
               "lity is %.3e\n", ii_cnt, ii_sum);
alpar@1
   222
      }
alpar@1
   223
      return;
alpar@1
   224
}
alpar@1
   225
alpar@1
   226
/***********************************************************************
alpar@1
   227
*  record_solution - record better integer feasible solution
alpar@1
   228
*
alpar@1
   229
*  This routine records optimal basic solution of LP relaxation of the
alpar@1
   230
*  current subproblem, which being integer feasible is better than the
alpar@1
   231
*  best known integer feasible solution. */
alpar@1
   232
alpar@1
   233
static void record_solution(glp_tree *T)
alpar@1
   234
{     glp_prob *mip = T->mip;
alpar@1
   235
      int i, j;
alpar@1
   236
      mip->mip_stat = GLP_FEAS;
alpar@1
   237
      mip->mip_obj = mip->obj_val;
alpar@1
   238
      for (i = 1; i <= mip->m; i++)
alpar@1
   239
      {  GLPROW *row = mip->row[i];
alpar@1
   240
         row->mipx = row->prim;
alpar@1
   241
      }
alpar@1
   242
      for (j = 1; j <= mip->n; j++)
alpar@1
   243
      {  GLPCOL *col = mip->col[j];
alpar@1
   244
         if (col->kind == GLP_CV)
alpar@1
   245
            col->mipx = col->prim;
alpar@1
   246
         else if (col->kind == GLP_IV)
alpar@1
   247
         {  /* value of the integer column must be integral */
alpar@1
   248
            col->mipx = floor(col->prim + 0.5);
alpar@1
   249
         }
alpar@1
   250
         else
alpar@1
   251
            xassert(col != col);
alpar@1
   252
      }
alpar@1
   253
      T->sol_cnt++;
alpar@1
   254
      return;
alpar@1
   255
}
alpar@1
   256
alpar@1
   257
/***********************************************************************
alpar@1
   258
*  fix_by_red_cost - fix non-basic integer columns by reduced costs
alpar@1
   259
*
alpar@1
   260
*  This routine fixes some non-basic integer columns if their reduced
alpar@1
   261
*  costs indicate that increasing (decreasing) the column at least by
alpar@1
   262
*  one involves the objective value becoming worse than the incumbent
alpar@1
   263
*  objective value. */
alpar@1
   264
alpar@1
   265
static void fix_by_red_cost(glp_tree *T)
alpar@1
   266
{     glp_prob *mip = T->mip;
alpar@1
   267
      int j, stat, fixed = 0;
alpar@1
   268
      double obj, lb, ub, dj;
alpar@1
   269
      /* the global bound must exist */
alpar@1
   270
      xassert(T->mip->mip_stat == GLP_FEAS);
alpar@1
   271
      /* basic solution of LP relaxation must be optimal */
alpar@1
   272
      xassert(mip->pbs_stat == GLP_FEAS && mip->dbs_stat == GLP_FEAS);
alpar@1
   273
      /* determine the objective function value */
alpar@1
   274
      obj = mip->obj_val;
alpar@1
   275
      /* walk through the column list */
alpar@1
   276
      for (j = 1; j <= mip->n; j++)
alpar@1
   277
      {  GLPCOL *col = mip->col[j];
alpar@1
   278
         /* if the column is not integer, skip it */
alpar@1
   279
         if (col->kind != GLP_IV) continue;
alpar@1
   280
         /* obtain bounds of j-th column */
alpar@1
   281
         lb = col->lb, ub = col->ub;
alpar@1
   282
         /* and determine its status and reduced cost */
alpar@1
   283
         stat = col->stat, dj = col->dual;
alpar@1
   284
         /* analyze the reduced cost */
alpar@1
   285
         switch (mip->dir)
alpar@1
   286
         {  case GLP_MIN:
alpar@1
   287
               /* minimization */
alpar@1
   288
               if (stat == GLP_NL)
alpar@1
   289
               {  /* j-th column is non-basic on its lower bound */
alpar@1
   290
                  if (dj < 0.0) dj = 0.0;
alpar@1
   291
                  if (obj + dj >= mip->mip_obj)
alpar@1
   292
                     glp_set_col_bnds(mip, j, GLP_FX, lb, lb), fixed++;
alpar@1
   293
               }
alpar@1
   294
               else if (stat == GLP_NU)
alpar@1
   295
               {  /* j-th column is non-basic on its upper bound */
alpar@1
   296
                  if (dj > 0.0) dj = 0.0;
alpar@1
   297
                  if (obj - dj >= mip->mip_obj)
alpar@1
   298
                     glp_set_col_bnds(mip, j, GLP_FX, ub, ub), fixed++;
alpar@1
   299
               }
alpar@1
   300
               break;
alpar@1
   301
            case GLP_MAX:
alpar@1
   302
               /* maximization */
alpar@1
   303
               if (stat == GLP_NL)
alpar@1
   304
               {  /* j-th column is non-basic on its lower bound */
alpar@1
   305
                  if (dj > 0.0) dj = 0.0;
alpar@1
   306
                  if (obj + dj <= mip->mip_obj)
alpar@1
   307
                     glp_set_col_bnds(mip, j, GLP_FX, lb, lb), fixed++;
alpar@1
   308
               }
alpar@1
   309
               else if (stat == GLP_NU)
alpar@1
   310
               {  /* j-th column is non-basic on its upper bound */
alpar@1
   311
                  if (dj < 0.0) dj = 0.0;
alpar@1
   312
                  if (obj - dj <= mip->mip_obj)
alpar@1
   313
                     glp_set_col_bnds(mip, j, GLP_FX, ub, ub), fixed++;
alpar@1
   314
               }
alpar@1
   315
               break;
alpar@1
   316
            default:
alpar@1
   317
               xassert(T != T);
alpar@1
   318
         }
alpar@1
   319
      }
alpar@1
   320
      if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   321
      {  if (fixed == 0)
alpar@1
   322
            /* nothing to say */;
alpar@1
   323
         else if (fixed == 1)
alpar@1
   324
            xprintf("One column has been fixed by reduced cost\n");
alpar@1
   325
         else
alpar@1
   326
            xprintf("%d columns have been fixed by reduced costs\n",
alpar@1
   327
               fixed);
alpar@1
   328
      }
alpar@1
   329
      /* fixing non-basic columns on their current bounds does not
alpar@1
   330
         change the basic solution */
alpar@1
   331
      xassert(mip->pbs_stat == GLP_FEAS && mip->dbs_stat == GLP_FEAS);
alpar@1
   332
      return;
alpar@1
   333
}
alpar@1
   334
alpar@1
   335
/***********************************************************************
alpar@1
   336
*  branch_on - perform branching on specified variable
alpar@1
   337
*
alpar@1
   338
*  This routine performs branching on j-th column (structural variable)
alpar@1
   339
*  of the current subproblem. The specified column must be of integer
alpar@1
   340
*  kind and must have a fractional value in optimal basic solution of
alpar@1
   341
*  LP relaxation of the current subproblem (i.e. only columns for which
alpar@1
   342
*  the flag non_int[j] is set are valid candidates to branch on).
alpar@1
   343
*
alpar@1
   344
*  Let x be j-th structural variable, and beta be its primal fractional
alpar@1
   345
*  value in the current basic solution. Branching on j-th variable is
alpar@1
   346
*  dividing the current subproblem into two new subproblems, which are
alpar@1
   347
*  identical to the current subproblem with the following exception: in
alpar@1
   348
*  the first subproblem that begins the down-branch x has a new upper
alpar@1
   349
*  bound x <= floor(beta), and in the second subproblem that begins the
alpar@1
   350
*  up-branch x has a new lower bound x >= ceil(beta).
alpar@1
   351
*
alpar@1
   352
*  Depending on estimation of local bounds for down- and up-branches
alpar@1
   353
*  this routine returns the following:
alpar@1
   354
*
alpar@1
   355
*  0 - both branches have been created;
alpar@1
   356
*  1 - one branch is hopeless and has been pruned, so now the current
alpar@1
   357
*      subproblem is other branch;
alpar@1
   358
*  2 - both branches are hopeless and have been pruned; new subproblem
alpar@1
   359
*      selection is needed to continue the search. */
alpar@1
   360
alpar@1
   361
static int branch_on(glp_tree *T, int j, int next)
alpar@1
   362
{     glp_prob *mip = T->mip;
alpar@1
   363
      IOSNPD *node;
alpar@1
   364
      int m = mip->m;
alpar@1
   365
      int n = mip->n;
alpar@1
   366
      int type, dn_type, up_type, dn_bad, up_bad, p, ret, clone[1+2];
alpar@1
   367
      double lb, ub, beta, new_ub, new_lb, dn_lp, up_lp, dn_bnd, up_bnd;
alpar@1
   368
      /* determine bounds and value of x[j] in optimal solution to LP
alpar@1
   369
         relaxation of the current subproblem */
alpar@1
   370
      xassert(1 <= j && j <= n);
alpar@1
   371
      type = mip->col[j]->type;
alpar@1
   372
      lb = mip->col[j]->lb;
alpar@1
   373
      ub = mip->col[j]->ub;
alpar@1
   374
      beta = mip->col[j]->prim;
alpar@1
   375
      /* determine new bounds of x[j] for down- and up-branches */
alpar@1
   376
      new_ub = floor(beta);
alpar@1
   377
      new_lb = ceil(beta);
alpar@1
   378
      switch (type)
alpar@1
   379
      {  case GLP_FR:
alpar@1
   380
            dn_type = GLP_UP;
alpar@1
   381
            up_type = GLP_LO;
alpar@1
   382
            break;
alpar@1
   383
         case GLP_LO:
alpar@1
   384
            xassert(lb <= new_ub);
alpar@1
   385
            dn_type = (lb == new_ub ? GLP_FX : GLP_DB);
alpar@1
   386
            xassert(lb + 1.0 <= new_lb);
alpar@1
   387
            up_type = GLP_LO;
alpar@1
   388
            break;
alpar@1
   389
         case GLP_UP:
alpar@1
   390
            xassert(new_ub <= ub - 1.0);
alpar@1
   391
            dn_type = GLP_UP;
alpar@1
   392
            xassert(new_lb <= ub);
alpar@1
   393
            up_type = (new_lb == ub ? GLP_FX : GLP_DB);
alpar@1
   394
            break;
alpar@1
   395
         case GLP_DB:
alpar@1
   396
            xassert(lb <= new_ub && new_ub <= ub - 1.0);
alpar@1
   397
            dn_type = (lb == new_ub ? GLP_FX : GLP_DB);
alpar@1
   398
            xassert(lb + 1.0 <= new_lb && new_lb <= ub);
alpar@1
   399
            up_type = (new_lb == ub ? GLP_FX : GLP_DB);
alpar@1
   400
            break;
alpar@1
   401
         default:
alpar@1
   402
            xassert(type != type);
alpar@1
   403
      }
alpar@1
   404
      /* compute local bounds to LP relaxation for both branches */
alpar@1
   405
      ios_eval_degrad(T, j, &dn_lp, &up_lp);
alpar@1
   406
      /* and improve them by rounding */
alpar@1
   407
      dn_bnd = ios_round_bound(T, dn_lp);
alpar@1
   408
      up_bnd = ios_round_bound(T, up_lp);
alpar@1
   409
      /* check local bounds for down- and up-branches */
alpar@1
   410
      dn_bad = !ios_is_hopeful(T, dn_bnd);
alpar@1
   411
      up_bad = !ios_is_hopeful(T, up_bnd);
alpar@1
   412
      if (dn_bad && up_bad)
alpar@1
   413
      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   414
            xprintf("Both down- and up-branches are hopeless\n");
alpar@1
   415
         ret = 2;
alpar@1
   416
         goto done;
alpar@1
   417
      }
alpar@1
   418
      else if (up_bad)
alpar@1
   419
      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   420
            xprintf("Up-branch is hopeless\n");
alpar@1
   421
         glp_set_col_bnds(mip, j, dn_type, lb, new_ub);
alpar@1
   422
         T->curr->lp_obj = dn_lp;
alpar@1
   423
         if (mip->dir == GLP_MIN)
alpar@1
   424
         {  if (T->curr->bound < dn_bnd)
alpar@1
   425
                T->curr->bound = dn_bnd;
alpar@1
   426
         }
alpar@1
   427
         else if (mip->dir == GLP_MAX)
alpar@1
   428
         {  if (T->curr->bound > dn_bnd)
alpar@1
   429
                T->curr->bound = dn_bnd;
alpar@1
   430
         }
alpar@1
   431
         else
alpar@1
   432
            xassert(mip != mip);
alpar@1
   433
         ret = 1;
alpar@1
   434
         goto done;
alpar@1
   435
      }
alpar@1
   436
      else if (dn_bad)
alpar@1
   437
      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   438
            xprintf("Down-branch is hopeless\n");
alpar@1
   439
         glp_set_col_bnds(mip, j, up_type, new_lb, ub);
alpar@1
   440
         T->curr->lp_obj = up_lp;
alpar@1
   441
         if (mip->dir == GLP_MIN)
alpar@1
   442
         {  if (T->curr->bound < up_bnd)
alpar@1
   443
                T->curr->bound = up_bnd;
alpar@1
   444
         }
alpar@1
   445
         else if (mip->dir == GLP_MAX)
alpar@1
   446
         {  if (T->curr->bound > up_bnd)
alpar@1
   447
                T->curr->bound = up_bnd;
alpar@1
   448
         }
alpar@1
   449
         else
alpar@1
   450
            xassert(mip != mip);
alpar@1
   451
         ret = 1;
alpar@1
   452
         goto done;
alpar@1
   453
      }
alpar@1
   454
      /* both down- and up-branches seem to be hopeful */
alpar@1
   455
      if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   456
         xprintf("Branching on column %d, primal value is %.9e\n",
alpar@1
   457
            j, beta);
alpar@1
   458
      /* determine the reference number of the current subproblem */
alpar@1
   459
      xassert(T->curr != NULL);
alpar@1
   460
      p = T->curr->p;
alpar@1
   461
      T->curr->br_var = j;
alpar@1
   462
      T->curr->br_val = beta;
alpar@1
   463
      /* freeze the current subproblem */
alpar@1
   464
      ios_freeze_node(T);
alpar@1
   465
      /* create two clones of the current subproblem; the first clone
alpar@1
   466
         begins the down-branch, the second one begins the up-branch */
alpar@1
   467
      ios_clone_node(T, p, 2, clone);
alpar@1
   468
      if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   469
         xprintf("Node %d begins down branch, node %d begins up branch "
alpar@1
   470
            "\n", clone[1], clone[2]);
alpar@1
   471
      /* set new upper bound of j-th column in the down-branch */
alpar@1
   472
      node = T->slot[clone[1]].node;
alpar@1
   473
      xassert(node != NULL);
alpar@1
   474
      xassert(node->up != NULL);
alpar@1
   475
      xassert(node->b_ptr == NULL);
alpar@1
   476
      node->b_ptr = dmp_get_atom(T->pool, sizeof(IOSBND));
alpar@1
   477
      node->b_ptr->k = m + j;
alpar@1
   478
      node->b_ptr->type = (unsigned char)dn_type;
alpar@1
   479
      node->b_ptr->lb = lb;
alpar@1
   480
      node->b_ptr->ub = new_ub;
alpar@1
   481
      node->b_ptr->next = NULL;
alpar@1
   482
      node->lp_obj = dn_lp;
alpar@1
   483
      if (mip->dir == GLP_MIN)
alpar@1
   484
      {  if (node->bound < dn_bnd)
alpar@1
   485
             node->bound = dn_bnd;
alpar@1
   486
      }
alpar@1
   487
      else if (mip->dir == GLP_MAX)
alpar@1
   488
      {  if (node->bound > dn_bnd)
alpar@1
   489
             node->bound = dn_bnd;
alpar@1
   490
      }
alpar@1
   491
      else
alpar@1
   492
         xassert(mip != mip);
alpar@1
   493
      /* set new lower bound of j-th column in the up-branch */
alpar@1
   494
      node = T->slot[clone[2]].node;
alpar@1
   495
      xassert(node != NULL);
alpar@1
   496
      xassert(node->up != NULL);
alpar@1
   497
      xassert(node->b_ptr == NULL);
alpar@1
   498
      node->b_ptr = dmp_get_atom(T->pool, sizeof(IOSBND));
alpar@1
   499
      node->b_ptr->k = m + j;
alpar@1
   500
      node->b_ptr->type = (unsigned char)up_type;
alpar@1
   501
      node->b_ptr->lb = new_lb;
alpar@1
   502
      node->b_ptr->ub = ub;
alpar@1
   503
      node->b_ptr->next = NULL;
alpar@1
   504
      node->lp_obj = up_lp;
alpar@1
   505
      if (mip->dir == GLP_MIN)
alpar@1
   506
      {  if (node->bound < up_bnd)
alpar@1
   507
             node->bound = up_bnd;
alpar@1
   508
      }
alpar@1
   509
      else if (mip->dir == GLP_MAX)
alpar@1
   510
      {  if (node->bound > up_bnd)
alpar@1
   511
             node->bound = up_bnd;
alpar@1
   512
      }
alpar@1
   513
      else
alpar@1
   514
         xassert(mip != mip);
alpar@1
   515
      /* suggest the subproblem to be solved next */
alpar@1
   516
      xassert(T->child == 0);
alpar@1
   517
      if (next == GLP_NO_BRNCH)
alpar@1
   518
         T->child = 0;
alpar@1
   519
      else if (next == GLP_DN_BRNCH)
alpar@1
   520
         T->child = clone[1];
alpar@1
   521
      else if (next == GLP_UP_BRNCH)
alpar@1
   522
         T->child = clone[2];
alpar@1
   523
      else
alpar@1
   524
         xassert(next != next);
alpar@1
   525
      ret = 0;
alpar@1
   526
done: return ret;
alpar@1
   527
}
alpar@1
   528
alpar@1
   529
/***********************************************************************
alpar@1
   530
*  cleanup_the_tree - prune hopeless branches from the tree
alpar@1
   531
*
alpar@1
   532
*  This routine walks through the active list and checks the local
alpar@1
   533
*  bound for every active subproblem. If the local bound indicates that
alpar@1
   534
*  the subproblem cannot have integer optimal solution better than the
alpar@1
   535
*  incumbent objective value, the routine deletes such subproblem that,
alpar@1
   536
*  in turn, involves pruning the corresponding branch of the tree. */
alpar@1
   537
alpar@1
   538
static void cleanup_the_tree(glp_tree *T)
alpar@1
   539
{     IOSNPD *node, *next_node;
alpar@1
   540
      int count = 0;
alpar@1
   541
      /* the global bound must exist */
alpar@1
   542
      xassert(T->mip->mip_stat == GLP_FEAS);
alpar@1
   543
      /* walk through the list of active subproblems */
alpar@1
   544
      for (node = T->head; node != NULL; node = next_node)
alpar@1
   545
      {  /* deleting some active problem node may involve deleting its
alpar@1
   546
            parents recursively; however, all its parents being created
alpar@1
   547
            *before* it are always *precede* it in the node list, so
alpar@1
   548
            the next problem node is never affected by such deletion */
alpar@1
   549
         next_node = node->next;
alpar@1
   550
         /* if the branch is hopeless, prune it */
alpar@1
   551
         if (!is_branch_hopeful(T, node->p))
alpar@1
   552
            ios_delete_node(T, node->p), count++;
alpar@1
   553
      }
alpar@1
   554
      if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   555
      {  if (count == 1)
alpar@1
   556
            xprintf("One hopeless branch has been pruned\n");
alpar@1
   557
         else if (count > 1)
alpar@1
   558
            xprintf("%d hopeless branches have been pruned\n", count);
alpar@1
   559
      }
alpar@1
   560
      return;
alpar@1
   561
}
alpar@1
   562
alpar@1
   563
/**********************************************************************/
alpar@1
   564
alpar@1
   565
static void generate_cuts(glp_tree *T)
alpar@1
   566
{     /* generate generic cuts with built-in generators */
alpar@1
   567
      if (!(T->parm->mir_cuts == GLP_ON ||
alpar@1
   568
            T->parm->gmi_cuts == GLP_ON ||
alpar@1
   569
            T->parm->cov_cuts == GLP_ON ||
alpar@1
   570
            T->parm->clq_cuts == GLP_ON)) goto done;
alpar@1
   571
#if 1 /* 20/IX-2008 */
alpar@1
   572
      {  int i, max_cuts, added_cuts;
alpar@1
   573
         max_cuts = T->n;
alpar@1
   574
         if (max_cuts < 1000) max_cuts = 1000;
alpar@1
   575
         added_cuts = 0;
alpar@1
   576
         for (i = T->orig_m+1; i <= T->mip->m; i++)
alpar@1
   577
         {  if (T->mip->row[i]->origin == GLP_RF_CUT)
alpar@1
   578
               added_cuts++;
alpar@1
   579
         }
alpar@1
   580
         /* xprintf("added_cuts = %d\n", added_cuts); */
alpar@1
   581
         if (added_cuts >= max_cuts) goto done;
alpar@1
   582
      }
alpar@1
   583
#endif
alpar@1
   584
      /* generate and add to POOL all cuts violated by x* */
alpar@1
   585
      if (T->parm->gmi_cuts == GLP_ON)
alpar@1
   586
      {  if (T->curr->changed < 5)
alpar@1
   587
            ios_gmi_gen(T);
alpar@1
   588
      }
alpar@1
   589
      if (T->parm->mir_cuts == GLP_ON)
alpar@1
   590
      {  xassert(T->mir_gen != NULL);
alpar@1
   591
         ios_mir_gen(T, T->mir_gen);
alpar@1
   592
      }
alpar@1
   593
      if (T->parm->cov_cuts == GLP_ON)
alpar@1
   594
      {  /* cover cuts works well along with mir cuts */
alpar@1
   595
         /*if (T->round <= 5)*/
alpar@1
   596
            ios_cov_gen(T);
alpar@1
   597
      }
alpar@1
   598
      if (T->parm->clq_cuts == GLP_ON)
alpar@1
   599
      {  if (T->clq_gen != NULL)
alpar@1
   600
         {  if (T->curr->level == 0 && T->curr->changed < 50 ||
alpar@1
   601
                T->curr->level >  0 && T->curr->changed < 5)
alpar@1
   602
               ios_clq_gen(T, T->clq_gen);
alpar@1
   603
         }
alpar@1
   604
      }
alpar@1
   605
done: return;
alpar@1
   606
}
alpar@1
   607
alpar@1
   608
/**********************************************************************/
alpar@1
   609
alpar@1
   610
static void remove_cuts(glp_tree *T)
alpar@1
   611
{     /* remove inactive cuts (some valueable globally valid cut might
alpar@1
   612
         be saved in the global cut pool) */
alpar@1
   613
      int i, cnt = 0, *num = NULL;
alpar@1
   614
      xassert(T->curr != NULL);
alpar@1
   615
      for (i = T->orig_m+1; i <= T->mip->m; i++)
alpar@1
   616
      {  if (T->mip->row[i]->origin == GLP_RF_CUT &&
alpar@1
   617
             T->mip->row[i]->level == T->curr->level &&
alpar@1
   618
             T->mip->row[i]->stat == GLP_BS)
alpar@1
   619
         {  if (num == NULL)
alpar@1
   620
               num = xcalloc(1+T->mip->m, sizeof(int));
alpar@1
   621
            num[++cnt] = i;
alpar@1
   622
         }
alpar@1
   623
      }
alpar@1
   624
      if (cnt > 0)
alpar@1
   625
      {  glp_del_rows(T->mip, cnt, num);
alpar@1
   626
#if 0
alpar@1
   627
         xprintf("%d inactive cut(s) removed\n", cnt);
alpar@1
   628
#endif
alpar@1
   629
         xfree(num);
alpar@1
   630
         xassert(glp_factorize(T->mip) == 0);
alpar@1
   631
      }
alpar@1
   632
      return;
alpar@1
   633
}
alpar@1
   634
alpar@1
   635
/**********************************************************************/
alpar@1
   636
alpar@1
   637
static void display_cut_info(glp_tree *T)
alpar@1
   638
{     glp_prob *mip = T->mip;
alpar@1
   639
      int i, gmi = 0, mir = 0, cov = 0, clq = 0, app = 0;
alpar@1
   640
      for (i = mip->m; i > 0; i--)
alpar@1
   641
      {  GLPROW *row;
alpar@1
   642
         row = mip->row[i];
alpar@1
   643
         /* if (row->level < T->curr->level) break; */
alpar@1
   644
         if (row->origin == GLP_RF_CUT)
alpar@1
   645
         {  if (row->klass == GLP_RF_GMI)
alpar@1
   646
               gmi++;
alpar@1
   647
            else if (row->klass == GLP_RF_MIR)
alpar@1
   648
               mir++;
alpar@1
   649
            else if (row->klass == GLP_RF_COV)
alpar@1
   650
               cov++;
alpar@1
   651
            else if (row->klass == GLP_RF_CLQ)
alpar@1
   652
               clq++;
alpar@1
   653
            else
alpar@1
   654
               app++;
alpar@1
   655
         }
alpar@1
   656
      }
alpar@1
   657
      xassert(T->curr != NULL);
alpar@1
   658
      if (gmi + mir + cov + clq + app > 0)
alpar@1
   659
      {  xprintf("Cuts on level %d:", T->curr->level);
alpar@1
   660
         if (gmi > 0) xprintf(" gmi = %d;", gmi);
alpar@1
   661
         if (mir > 0) xprintf(" mir = %d;", mir);
alpar@1
   662
         if (cov > 0) xprintf(" cov = %d;", cov);
alpar@1
   663
         if (clq > 0) xprintf(" clq = %d;", clq);
alpar@1
   664
         if (app > 0) xprintf(" app = %d;", app);
alpar@1
   665
         xprintf("\n");
alpar@1
   666
      }
alpar@1
   667
      return;
alpar@1
   668
}
alpar@1
   669
alpar@1
   670
/***********************************************************************
alpar@1
   671
*  NAME
alpar@1
   672
*
alpar@1
   673
*  ios_driver - branch-and-cut driver
alpar@1
   674
*
alpar@1
   675
*  SYNOPSIS
alpar@1
   676
*
alpar@1
   677
*  #include "glpios.h"
alpar@1
   678
*  int ios_driver(glp_tree *T);
alpar@1
   679
*
alpar@1
   680
*  DESCRIPTION
alpar@1
   681
*
alpar@1
   682
*  The routine ios_driver is a branch-and-cut driver. It controls the
alpar@1
   683
*  MIP solution process.
alpar@1
   684
*
alpar@1
   685
*  RETURNS
alpar@1
   686
*
alpar@1
   687
*  0  The MIP problem instance has been successfully solved. This code
alpar@1
   688
*     does not necessarily mean that the solver has found optimal
alpar@1
   689
*     solution. It only means that the solution process was successful.
alpar@1
   690
*
alpar@1
   691
*  GLP_EFAIL
alpar@1
   692
*     The search was prematurely terminated due to the solver failure.
alpar@1
   693
*
alpar@1
   694
*  GLP_EMIPGAP
alpar@1
   695
*     The search was prematurely terminated, because the relative mip
alpar@1
   696
*     gap tolerance has been reached.
alpar@1
   697
*
alpar@1
   698
*  GLP_ETMLIM
alpar@1
   699
*     The search was prematurely terminated, because the time limit has
alpar@1
   700
*     been exceeded.
alpar@1
   701
*
alpar@1
   702
*  GLP_ESTOP
alpar@1
   703
*     The search was prematurely terminated by application. */
alpar@1
   704
alpar@1
   705
int ios_driver(glp_tree *T)
alpar@1
   706
{     int p, curr_p, p_stat, d_stat, ret;
alpar@1
   707
#if 1 /* carry out to glp_tree */
alpar@1
   708
      int pred_p = 0;
alpar@1
   709
      /* if the current subproblem has been just created due to
alpar@1
   710
         branching, pred_p is the reference number of its parent
alpar@1
   711
         subproblem, otherwise pred_p is zero */
alpar@1
   712
#endif
alpar@1
   713
      glp_long ttt = T->tm_beg;
alpar@1
   714
#if 0
alpar@1
   715
      ((glp_iocp *)T->parm)->msg_lev = GLP_MSG_DBG;
alpar@1
   716
#endif
alpar@1
   717
      /* on entry to the B&B driver it is assumed that the active list
alpar@1
   718
         contains the only active (i.e. root) subproblem, which is the
alpar@1
   719
         original MIP problem to be solved */
alpar@1
   720
loop: /* main loop starts here */
alpar@1
   721
      /* at this point the current subproblem does not exist */
alpar@1
   722
      xassert(T->curr == NULL);
alpar@1
   723
      /* if the active list is empty, the search is finished */
alpar@1
   724
      if (T->head == NULL)
alpar@1
   725
      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   726
            xprintf("Active list is empty!\n");
alpar@1
   727
         xassert(dmp_in_use(T->pool).lo == 0);
alpar@1
   728
         ret = 0;
alpar@1
   729
         goto done;
alpar@1
   730
      }
alpar@1
   731
      /* select some active subproblem to continue the search */
alpar@1
   732
      xassert(T->next_p == 0);
alpar@1
   733
      /* let the application program select subproblem */
alpar@1
   734
      if (T->parm->cb_func != NULL)
alpar@1
   735
      {  xassert(T->reason == 0);
alpar@1
   736
         T->reason = GLP_ISELECT;
alpar@1
   737
         T->parm->cb_func(T, T->parm->cb_info);
alpar@1
   738
         T->reason = 0;
alpar@1
   739
         if (T->stop)
alpar@1
   740
         {  ret = GLP_ESTOP;
alpar@1
   741
            goto done;
alpar@1
   742
         }
alpar@1
   743
      }
alpar@1
   744
      if (T->next_p != 0)
alpar@1
   745
      {  /* the application program has selected something */
alpar@1
   746
         ;
alpar@1
   747
      }
alpar@1
   748
      else if (T->a_cnt == 1)
alpar@1
   749
      {  /* the only active subproblem exists, so select it */
alpar@1
   750
         xassert(T->head->next == NULL);
alpar@1
   751
         T->next_p = T->head->p;
alpar@1
   752
      }
alpar@1
   753
      else if (T->child != 0)
alpar@1
   754
      {  /* select one of branching childs suggested by the branching
alpar@1
   755
            heuristic */
alpar@1
   756
         T->next_p = T->child;
alpar@1
   757
      }
alpar@1
   758
      else
alpar@1
   759
      {  /* select active subproblem as specified by the backtracking
alpar@1
   760
            technique option */
alpar@1
   761
         T->next_p = ios_choose_node(T);
alpar@1
   762
      }
alpar@1
   763
      /* the active subproblem just selected becomes current */
alpar@1
   764
      ios_revive_node(T, T->next_p);
alpar@1
   765
      T->next_p = T->child = 0;
alpar@1
   766
      /* invalidate pred_p, if it is not the reference number of the
alpar@1
   767
         parent of the current subproblem */
alpar@1
   768
      if (T->curr->up != NULL && T->curr->up->p != pred_p) pred_p = 0;
alpar@1
   769
      /* determine the reference number of the current subproblem */
alpar@1
   770
      p = T->curr->p;
alpar@1
   771
      if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   772
      {  xprintf("-----------------------------------------------------"
alpar@1
   773
            "-------------------\n");
alpar@1
   774
         xprintf("Processing node %d at level %d\n", p, T->curr->level);
alpar@1
   775
      }
alpar@1
   776
      /* if it is the root subproblem, initialize cut generators */
alpar@1
   777
      if (p == 1)
alpar@1
   778
      {  if (T->parm->gmi_cuts == GLP_ON)
alpar@1
   779
         {  if (T->parm->msg_lev >= GLP_MSG_ALL)
alpar@1
   780
               xprintf("Gomory's cuts enabled\n");
alpar@1
   781
         }
alpar@1
   782
         if (T->parm->mir_cuts == GLP_ON)
alpar@1
   783
         {  if (T->parm->msg_lev >= GLP_MSG_ALL)
alpar@1
   784
               xprintf("MIR cuts enabled\n");
alpar@1
   785
            xassert(T->mir_gen == NULL);
alpar@1
   786
            T->mir_gen = ios_mir_init(T);
alpar@1
   787
         }
alpar@1
   788
         if (T->parm->cov_cuts == GLP_ON)
alpar@1
   789
         {  if (T->parm->msg_lev >= GLP_MSG_ALL)
alpar@1
   790
               xprintf("Cover cuts enabled\n");
alpar@1
   791
         }
alpar@1
   792
         if (T->parm->clq_cuts == GLP_ON)
alpar@1
   793
         {  xassert(T->clq_gen == NULL);
alpar@1
   794
            if (T->parm->msg_lev >= GLP_MSG_ALL)
alpar@1
   795
               xprintf("Clique cuts enabled\n");
alpar@1
   796
            T->clq_gen = ios_clq_init(T);
alpar@1
   797
         }
alpar@1
   798
      }
alpar@1
   799
more: /* minor loop starts here */
alpar@1
   800
      /* at this point the current subproblem needs either to be solved
alpar@1
   801
         for the first time or re-optimized due to reformulation */
alpar@1
   802
      /* display current progress of the search */
alpar@1
   803
      if (T->parm->msg_lev >= GLP_MSG_DBG ||
alpar@1
   804
          T->parm->msg_lev >= GLP_MSG_ON &&
alpar@1
   805
        (double)(T->parm->out_frq - 1) <=
alpar@1
   806
            1000.0 * xdifftime(xtime(), T->tm_lag))
alpar@1
   807
         show_progress(T, 0);
alpar@1
   808
      if (T->parm->msg_lev >= GLP_MSG_ALL &&
alpar@1
   809
            xdifftime(xtime(), ttt) >= 60.0)
alpar@1
   810
      {  glp_long total;
alpar@1
   811
         glp_mem_usage(NULL, NULL, &total, NULL);
alpar@1
   812
         xprintf("Time used: %.1f secs.  Memory used: %.1f Mb.\n",
alpar@1
   813
            xdifftime(xtime(), T->tm_beg), xltod(total) / 1048576.0);
alpar@1
   814
         ttt = xtime();
alpar@1
   815
      }
alpar@1
   816
      /* check the mip gap */
alpar@1
   817
      if (T->parm->mip_gap > 0.0 &&
alpar@1
   818
          ios_relative_gap(T) <= T->parm->mip_gap)
alpar@1
   819
      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   820
            xprintf("Relative gap tolerance reached; search terminated "
alpar@1
   821
               "\n");
alpar@1
   822
         ret = GLP_EMIPGAP;
alpar@1
   823
         goto done;
alpar@1
   824
      }
alpar@1
   825
      /* check if the time limit has been exhausted */
alpar@1
   826
      if (T->parm->tm_lim < INT_MAX &&
alpar@1
   827
         (double)(T->parm->tm_lim - 1) <=
alpar@1
   828
         1000.0 * xdifftime(xtime(), T->tm_beg))
alpar@1
   829
      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   830
            xprintf("Time limit exhausted; search terminated\n");
alpar@1
   831
         ret = GLP_ETMLIM;
alpar@1
   832
         goto done;
alpar@1
   833
      }
alpar@1
   834
      /* let the application program preprocess the subproblem */
alpar@1
   835
      if (T->parm->cb_func != NULL)
alpar@1
   836
      {  xassert(T->reason == 0);
alpar@1
   837
         T->reason = GLP_IPREPRO;
alpar@1
   838
         T->parm->cb_func(T, T->parm->cb_info);
alpar@1
   839
         T->reason = 0;
alpar@1
   840
         if (T->stop)
alpar@1
   841
         {  ret = GLP_ESTOP;
alpar@1
   842
            goto done;
alpar@1
   843
         }
alpar@1
   844
      }
alpar@1
   845
      /* perform basic preprocessing */
alpar@1
   846
      if (T->parm->pp_tech == GLP_PP_NONE)
alpar@1
   847
         ;
alpar@1
   848
      else if (T->parm->pp_tech == GLP_PP_ROOT)
alpar@1
   849
      {  if (T->curr->level == 0)
alpar@1
   850
         {  if (ios_preprocess_node(T, 100))
alpar@1
   851
               goto fath;
alpar@1
   852
         }
alpar@1
   853
      }
alpar@1
   854
      else if (T->parm->pp_tech == GLP_PP_ALL)
alpar@1
   855
      {  if (ios_preprocess_node(T, T->curr->level == 0 ? 100 : 10))
alpar@1
   856
            goto fath;
alpar@1
   857
      }
alpar@1
   858
      else
alpar@1
   859
         xassert(T != T);
alpar@1
   860
      /* preprocessing may improve the global bound */
alpar@1
   861
      if (!is_branch_hopeful(T, p))
alpar@1
   862
      {  xprintf("*** not tested yet ***\n");
alpar@1
   863
         goto fath;
alpar@1
   864
      }
alpar@1
   865
      /* solve LP relaxation of the current subproblem */
alpar@1
   866
      if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   867
         xprintf("Solving LP relaxation...\n");
alpar@1
   868
      ret = ios_solve_node(T);
alpar@1
   869
      if (!(ret == 0 || ret == GLP_EOBJLL || ret == GLP_EOBJUL))
alpar@1
   870
      {  if (T->parm->msg_lev >= GLP_MSG_ERR)
alpar@1
   871
            xprintf("ios_driver: unable to solve current LP relaxation;"
alpar@1
   872
               " glp_simplex returned %d\n", ret);
alpar@1
   873
         ret = GLP_EFAIL;
alpar@1
   874
         goto done;
alpar@1
   875
      }
alpar@1
   876
      /* analyze status of the basic solution to LP relaxation found */
alpar@1
   877
      p_stat = T->mip->pbs_stat;
alpar@1
   878
      d_stat = T->mip->dbs_stat;
alpar@1
   879
      if (p_stat == GLP_FEAS && d_stat == GLP_FEAS)
alpar@1
   880
      {  /* LP relaxation has optimal solution */
alpar@1
   881
         if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   882
            xprintf("Found optimal solution to LP relaxation\n");
alpar@1
   883
      }
alpar@1
   884
      else if (d_stat == GLP_NOFEAS)
alpar@1
   885
      {  /* LP relaxation has no dual feasible solution */
alpar@1
   886
         /* since the current subproblem cannot have a larger feasible
alpar@1
   887
            region than its parent, there is something wrong */
alpar@1
   888
         if (T->parm->msg_lev >= GLP_MSG_ERR)
alpar@1
   889
            xprintf("ios_driver: current LP relaxation has no dual feas"
alpar@1
   890
               "ible solution\n");
alpar@1
   891
         ret = GLP_EFAIL;
alpar@1
   892
         goto done;
alpar@1
   893
      }
alpar@1
   894
      else if (p_stat == GLP_INFEAS && d_stat == GLP_FEAS)
alpar@1
   895
      {  /* LP relaxation has no primal solution which is better than
alpar@1
   896
            the incumbent objective value */
alpar@1
   897
         xassert(T->mip->mip_stat == GLP_FEAS);
alpar@1
   898
         if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   899
            xprintf("LP relaxation has no solution better than incumben"
alpar@1
   900
               "t objective value\n");
alpar@1
   901
         /* prune the branch */
alpar@1
   902
         goto fath;
alpar@1
   903
      }
alpar@1
   904
      else if (p_stat == GLP_NOFEAS)
alpar@1
   905
      {  /* LP relaxation has no primal feasible solution */
alpar@1
   906
         if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   907
            xprintf("LP relaxation has no feasible solution\n");
alpar@1
   908
         /* prune the branch */
alpar@1
   909
         goto fath;
alpar@1
   910
      }
alpar@1
   911
      else
alpar@1
   912
      {  /* other cases cannot appear */
alpar@1
   913
         xassert(T->mip != T->mip);
alpar@1
   914
      }
alpar@1
   915
      /* at this point basic solution to LP relaxation of the current
alpar@1
   916
         subproblem is optimal */
alpar@1
   917
      xassert(p_stat == GLP_FEAS && d_stat == GLP_FEAS);
alpar@1
   918
      xassert(T->curr != NULL);
alpar@1
   919
      T->curr->lp_obj = T->mip->obj_val;
alpar@1
   920
      /* thus, it defines a local bound to integer optimal solution of
alpar@1
   921
         the current subproblem */
alpar@1
   922
      {  double bound = T->mip->obj_val;
alpar@1
   923
         /* some local bound to the current subproblem could be already
alpar@1
   924
            set before, so we should only improve it */
alpar@1
   925
         bound = ios_round_bound(T, bound);
alpar@1
   926
         if (T->mip->dir == GLP_MIN)
alpar@1
   927
         {  if (T->curr->bound < bound)
alpar@1
   928
               T->curr->bound = bound;
alpar@1
   929
         }
alpar@1
   930
         else if (T->mip->dir == GLP_MAX)
alpar@1
   931
         {  if (T->curr->bound > bound)
alpar@1
   932
               T->curr->bound = bound;
alpar@1
   933
         }
alpar@1
   934
         else
alpar@1
   935
            xassert(T->mip != T->mip);
alpar@1
   936
         if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   937
            xprintf("Local bound is %.9e\n", bound);
alpar@1
   938
      }
alpar@1
   939
      /* if the local bound indicates that integer optimal solution of
alpar@1
   940
         the current subproblem cannot be better than the global bound,
alpar@1
   941
         prune the branch */
alpar@1
   942
      if (!is_branch_hopeful(T, p))
alpar@1
   943
      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   944
            xprintf("Current branch is hopeless and can be pruned\n");
alpar@1
   945
         goto fath;
alpar@1
   946
      }
alpar@1
   947
      /* let the application program generate additional rows ("lazy"
alpar@1
   948
         constraints) */
alpar@1
   949
      xassert(T->reopt == 0);
alpar@1
   950
      xassert(T->reinv == 0);
alpar@1
   951
      if (T->parm->cb_func != NULL)
alpar@1
   952
      {  xassert(T->reason == 0);
alpar@1
   953
         T->reason = GLP_IROWGEN;
alpar@1
   954
         T->parm->cb_func(T, T->parm->cb_info);
alpar@1
   955
         T->reason = 0;
alpar@1
   956
         if (T->stop)
alpar@1
   957
         {  ret = GLP_ESTOP;
alpar@1
   958
            goto done;
alpar@1
   959
         }
alpar@1
   960
         if (T->reopt)
alpar@1
   961
         {  /* some rows were added; re-optimization is needed */
alpar@1
   962
            T->reopt = T->reinv = 0;
alpar@1
   963
            goto more;
alpar@1
   964
         }
alpar@1
   965
         if (T->reinv)
alpar@1
   966
         {  /* no rows were added, however, some inactive rows were
alpar@1
   967
               removed */
alpar@1
   968
            T->reinv = 0;
alpar@1
   969
            xassert(glp_factorize(T->mip) == 0);
alpar@1
   970
         }
alpar@1
   971
      }
alpar@1
   972
      /* check if the basic solution is integer feasible */
alpar@1
   973
      check_integrality(T);
alpar@1
   974
      /* if the basic solution satisfies to all integrality conditions,
alpar@1
   975
         it is a new, better integer feasible solution */
alpar@1
   976
      if (T->curr->ii_cnt == 0)
alpar@1
   977
      {  if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
   978
            xprintf("New integer feasible solution found\n");
alpar@1
   979
         if (T->parm->msg_lev >= GLP_MSG_ALL)
alpar@1
   980
            display_cut_info(T);
alpar@1
   981
         record_solution(T);
alpar@1
   982
         if (T->parm->msg_lev >= GLP_MSG_ON)
alpar@1
   983
            show_progress(T, 1);
alpar@1
   984
         /* make the application program happy */
alpar@1
   985
         if (T->parm->cb_func != NULL)
alpar@1
   986
         {  xassert(T->reason == 0);
alpar@1
   987
            T->reason = GLP_IBINGO;
alpar@1
   988
            T->parm->cb_func(T, T->parm->cb_info);
alpar@1
   989
            T->reason = 0;
alpar@1
   990
            if (T->stop)
alpar@1
   991
            {  ret = GLP_ESTOP;
alpar@1
   992
               goto done;
alpar@1
   993
            }
alpar@1
   994
         }
alpar@1
   995
         /* since the current subproblem has been fathomed, prune its
alpar@1
   996
            branch */
alpar@1
   997
         goto fath;
alpar@1
   998
      }
alpar@1
   999
      /* at this point basic solution to LP relaxation of the current
alpar@1
  1000
         subproblem is optimal, but integer infeasible */
alpar@1
  1001
      /* try to fix some non-basic structural variables of integer kind
alpar@1
  1002
         on their current bounds due to reduced costs */
alpar@1
  1003
      if (T->mip->mip_stat == GLP_FEAS)
alpar@1
  1004
         fix_by_red_cost(T);
alpar@1
  1005
      /* let the application program try to find some solution to the
alpar@1
  1006
         original MIP with a primal heuristic */
alpar@1
  1007
      if (T->parm->cb_func != NULL)
alpar@1
  1008
      {  xassert(T->reason == 0);
alpar@1
  1009
         T->reason = GLP_IHEUR;
alpar@1
  1010
         T->parm->cb_func(T, T->parm->cb_info);
alpar@1
  1011
         T->reason = 0;
alpar@1
  1012
         if (T->stop)
alpar@1
  1013
         {  ret = GLP_ESTOP;
alpar@1
  1014
            goto done;
alpar@1
  1015
         }
alpar@1
  1016
         /* check if the current branch became hopeless */
alpar@1
  1017
         if (!is_branch_hopeful(T, p))
alpar@1
  1018
         {  if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
  1019
               xprintf("Current branch became hopeless and can be prune"
alpar@1
  1020
                  "d\n");
alpar@1
  1021
            goto fath;
alpar@1
  1022
         }
alpar@1
  1023
      }
alpar@1
  1024
      /* try to find solution with the feasibility pump heuristic */
alpar@1
  1025
      if (T->parm->fp_heur)
alpar@1
  1026
      {  xassert(T->reason == 0);
alpar@1
  1027
         T->reason = GLP_IHEUR;
alpar@1
  1028
         ios_feas_pump(T);
alpar@1
  1029
         T->reason = 0;
alpar@1
  1030
         /* check if the current branch became hopeless */
alpar@1
  1031
         if (!is_branch_hopeful(T, p))
alpar@1
  1032
         {  if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
  1033
               xprintf("Current branch became hopeless and can be prune"
alpar@1
  1034
                  "d\n");
alpar@1
  1035
            goto fath;
alpar@1
  1036
         }
alpar@1
  1037
      }
alpar@1
  1038
      /* it's time to generate cutting planes */
alpar@1
  1039
      xassert(T->local != NULL);
alpar@1
  1040
      xassert(T->local->size == 0);
alpar@1
  1041
      /* let the application program generate some cuts; note that it
alpar@1
  1042
         can add cuts either to the local cut pool or directly to the
alpar@1
  1043
         current subproblem */
alpar@1
  1044
      if (T->parm->cb_func != NULL)
alpar@1
  1045
      {  xassert(T->reason == 0);
alpar@1
  1046
         T->reason = GLP_ICUTGEN;
alpar@1
  1047
         T->parm->cb_func(T, T->parm->cb_info);
alpar@1
  1048
         T->reason = 0;
alpar@1
  1049
         if (T->stop)
alpar@1
  1050
         {  ret = GLP_ESTOP;
alpar@1
  1051
            goto done;
alpar@1
  1052
         }
alpar@1
  1053
      }
alpar@1
  1054
      /* try to generate generic cuts with built-in generators
alpar@1
  1055
         (as suggested by Matteo Fischetti et al. the built-in cuts
alpar@1
  1056
         are not generated at each branching node; an intense attempt
alpar@1
  1057
         of generating new cuts is only made at the root node, and then
alpar@1
  1058
         a moderate effort is spent after each backtracking step) */
alpar@1
  1059
      if (T->curr->level == 0 || pred_p == 0)
alpar@1
  1060
      {  xassert(T->reason == 0);
alpar@1
  1061
         T->reason = GLP_ICUTGEN;
alpar@1
  1062
         generate_cuts(T);
alpar@1
  1063
         T->reason = 0;
alpar@1
  1064
      }
alpar@1
  1065
      /* if the local cut pool is not empty, select useful cuts and add
alpar@1
  1066
         them to the current subproblem */
alpar@1
  1067
      if (T->local->size > 0)
alpar@1
  1068
      {  xassert(T->reason == 0);
alpar@1
  1069
         T->reason = GLP_ICUTGEN;
alpar@1
  1070
         ios_process_cuts(T);
alpar@1
  1071
         T->reason = 0;
alpar@1
  1072
      }
alpar@1
  1073
      /* clear the local cut pool */
alpar@1
  1074
      ios_clear_pool(T, T->local);
alpar@1
  1075
      /* perform re-optimization, if necessary */
alpar@1
  1076
      if (T->reopt)
alpar@1
  1077
      {  T->reopt = 0;
alpar@1
  1078
         T->curr->changed++;
alpar@1
  1079
         goto more;
alpar@1
  1080
      }
alpar@1
  1081
      /* no cuts were generated; remove inactive cuts */
alpar@1
  1082
      remove_cuts(T);
alpar@1
  1083
      if (T->parm->msg_lev >= GLP_MSG_ALL && T->curr->level == 0)
alpar@1
  1084
         display_cut_info(T);
alpar@1
  1085
      /* update history information used on pseudocost branching */
alpar@1
  1086
      if (T->pcost != NULL) ios_pcost_update(T);
alpar@1
  1087
      /* it's time to perform branching */
alpar@1
  1088
      xassert(T->br_var == 0);
alpar@1
  1089
      xassert(T->br_sel == 0);
alpar@1
  1090
      /* let the application program choose variable to branch on */
alpar@1
  1091
      if (T->parm->cb_func != NULL)
alpar@1
  1092
      {  xassert(T->reason == 0);
alpar@1
  1093
         xassert(T->br_var == 0);
alpar@1
  1094
         xassert(T->br_sel == 0);
alpar@1
  1095
         T->reason = GLP_IBRANCH;
alpar@1
  1096
         T->parm->cb_func(T, T->parm->cb_info);
alpar@1
  1097
         T->reason = 0;
alpar@1
  1098
         if (T->stop)
alpar@1
  1099
         {  ret = GLP_ESTOP;
alpar@1
  1100
            goto done;
alpar@1
  1101
         }
alpar@1
  1102
      }
alpar@1
  1103
      /* if nothing has been chosen, choose some variable as specified
alpar@1
  1104
         by the branching technique option */
alpar@1
  1105
      if (T->br_var == 0)
alpar@1
  1106
         T->br_var = ios_choose_var(T, &T->br_sel);
alpar@1
  1107
      /* perform actual branching */
alpar@1
  1108
      curr_p = T->curr->p;
alpar@1
  1109
      ret = branch_on(T, T->br_var, T->br_sel);
alpar@1
  1110
      T->br_var = T->br_sel = 0;
alpar@1
  1111
      if (ret == 0)
alpar@1
  1112
      {  /* both branches have been created */
alpar@1
  1113
         pred_p = curr_p;
alpar@1
  1114
         goto loop;
alpar@1
  1115
      }
alpar@1
  1116
      else if (ret == 1)
alpar@1
  1117
      {  /* one branch is hopeless and has been pruned, so now the
alpar@1
  1118
            current subproblem is other branch */
alpar@1
  1119
         /* the current subproblem should be considered as a new one,
alpar@1
  1120
            since one bound of the branching variable was changed */
alpar@1
  1121
         T->curr->solved = T->curr->changed = 0;
alpar@1
  1122
         goto more;
alpar@1
  1123
      }
alpar@1
  1124
      else if (ret == 2)
alpar@1
  1125
      {  /* both branches are hopeless and have been pruned; new
alpar@1
  1126
            subproblem selection is needed to continue the search */
alpar@1
  1127
         goto fath;
alpar@1
  1128
      }
alpar@1
  1129
      else
alpar@1
  1130
         xassert(ret != ret);
alpar@1
  1131
fath: /* the current subproblem has been fathomed */
alpar@1
  1132
      if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@1
  1133
         xprintf("Node %d fathomed\n", p);
alpar@1
  1134
      /* freeze the current subproblem */
alpar@1
  1135
      ios_freeze_node(T);
alpar@1
  1136
      /* and prune the corresponding branch of the tree */
alpar@1
  1137
      ios_delete_node(T, p);
alpar@1
  1138
      /* if a new integer feasible solution has just been found, other
alpar@1
  1139
         branches may become hopeless and therefore must be pruned */
alpar@1
  1140
      if (T->mip->mip_stat == GLP_FEAS) cleanup_the_tree(T);
alpar@1
  1141
      /* new subproblem selection is needed due to backtracking */
alpar@1
  1142
      pred_p = 0;
alpar@1
  1143
      goto loop;
alpar@1
  1144
done: /* display progress of the search on exit from the solver */
alpar@1
  1145
      if (T->parm->msg_lev >= GLP_MSG_ON)
alpar@1
  1146
         show_progress(T, 0);
alpar@1
  1147
      if (T->mir_gen != NULL)
alpar@1
  1148
         ios_mir_term(T->mir_gen), T->mir_gen = NULL;
alpar@1
  1149
      if (T->clq_gen != NULL)
alpar@1
  1150
         ios_clq_term(T->clq_gen), T->clq_gen = NULL;
alpar@1
  1151
      /* return to the calling program */
alpar@1
  1152
      return ret;
alpar@1
  1153
}
alpar@1
  1154
alpar@1
  1155
/* eof */