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