lemon-project-template-glpk

annotate deps/glpk/src/glpios03.c @ 11:4fc6ad2fb8a6

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