1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpios01.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,1611 @@
1.4 +/* glpios01.c */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpios.h"
1.29 +
1.30 +/***********************************************************************
1.31 +* NAME
1.32 +*
1.33 +* ios_create_tree - create branch-and-bound tree
1.34 +*
1.35 +* SYNOPSIS
1.36 +*
1.37 +* #include "glpios.h"
1.38 +* glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm);
1.39 +*
1.40 +* DESCRIPTION
1.41 +*
1.42 +* The routine ios_create_tree creates the branch-and-bound tree.
1.43 +*
1.44 +* Being created the tree consists of the only root subproblem whose
1.45 +* reference number is 1. Note that initially the root subproblem is in
1.46 +* frozen state and therefore needs to be revived.
1.47 +*
1.48 +* RETURNS
1.49 +*
1.50 +* The routine returns a pointer to the tree created. */
1.51 +
1.52 +static IOSNPD *new_node(glp_tree *tree, IOSNPD *parent);
1.53 +
1.54 +glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm)
1.55 +{ int m = mip->m;
1.56 + int n = mip->n;
1.57 + glp_tree *tree;
1.58 + int i, j;
1.59 + xassert(mip->tree == NULL);
1.60 + mip->tree = tree = xmalloc(sizeof(glp_tree));
1.61 + tree->pool = dmp_create_pool();
1.62 + tree->n = n;
1.63 + /* save original problem components */
1.64 + tree->orig_m = m;
1.65 + tree->orig_type = xcalloc(1+m+n, sizeof(char));
1.66 + tree->orig_lb = xcalloc(1+m+n, sizeof(double));
1.67 + tree->orig_ub = xcalloc(1+m+n, sizeof(double));
1.68 + tree->orig_stat = xcalloc(1+m+n, sizeof(char));
1.69 + tree->orig_prim = xcalloc(1+m+n, sizeof(double));
1.70 + tree->orig_dual = xcalloc(1+m+n, sizeof(double));
1.71 + for (i = 1; i <= m; i++)
1.72 + { GLPROW *row = mip->row[i];
1.73 + tree->orig_type[i] = (char)row->type;
1.74 + tree->orig_lb[i] = row->lb;
1.75 + tree->orig_ub[i] = row->ub;
1.76 + tree->orig_stat[i] = (char)row->stat;
1.77 + tree->orig_prim[i] = row->prim;
1.78 + tree->orig_dual[i] = row->dual;
1.79 + }
1.80 + for (j = 1; j <= n; j++)
1.81 + { GLPCOL *col = mip->col[j];
1.82 + tree->orig_type[m+j] = (char)col->type;
1.83 + tree->orig_lb[m+j] = col->lb;
1.84 + tree->orig_ub[m+j] = col->ub;
1.85 + tree->orig_stat[m+j] = (char)col->stat;
1.86 + tree->orig_prim[m+j] = col->prim;
1.87 + tree->orig_dual[m+j] = col->dual;
1.88 + }
1.89 + tree->orig_obj = mip->obj_val;
1.90 + /* initialize the branch-and-bound tree */
1.91 + tree->nslots = 0;
1.92 + tree->avail = 0;
1.93 + tree->slot = NULL;
1.94 + tree->head = tree->tail = NULL;
1.95 + tree->a_cnt = tree->n_cnt = tree->t_cnt = 0;
1.96 + /* the root subproblem is not solved yet, so its final components
1.97 + are unknown so far */
1.98 + tree->root_m = 0;
1.99 + tree->root_type = NULL;
1.100 + tree->root_lb = tree->root_ub = NULL;
1.101 + tree->root_stat = NULL;
1.102 + /* the current subproblem does not exist yet */
1.103 + tree->curr = NULL;
1.104 + tree->mip = mip;
1.105 + /*tree->solved = 0;*/
1.106 + tree->non_int = xcalloc(1+n, sizeof(char));
1.107 + memset(&tree->non_int[1], 0, n);
1.108 + /* arrays to save parent subproblem components will be allocated
1.109 + later */
1.110 + tree->pred_m = tree->pred_max = 0;
1.111 + tree->pred_type = NULL;
1.112 + tree->pred_lb = tree->pred_ub = NULL;
1.113 + tree->pred_stat = NULL;
1.114 + /* cut generator */
1.115 + tree->local = ios_create_pool(tree);
1.116 + /*tree->first_attempt = 1;*/
1.117 + /*tree->max_added_cuts = 0;*/
1.118 + /*tree->min_eff = 0.0;*/
1.119 + /*tree->miss = 0;*/
1.120 + /*tree->just_selected = 0;*/
1.121 + tree->mir_gen = NULL;
1.122 + tree->clq_gen = NULL;
1.123 + /*tree->round = 0;*/
1.124 +#if 0
1.125 + /* create the conflict graph */
1.126 + tree->n_ref = xcalloc(1+n, sizeof(int));
1.127 + memset(&tree->n_ref[1], 0, n * sizeof(int));
1.128 + tree->c_ref = xcalloc(1+n, sizeof(int));
1.129 + memset(&tree->c_ref[1], 0, n * sizeof(int));
1.130 + tree->g = scg_create_graph(0);
1.131 + tree->j_ref = xcalloc(1+tree->g->n_max, sizeof(int));
1.132 +#endif
1.133 + /* pseudocost branching */
1.134 + tree->pcost = NULL;
1.135 + tree->iwrk = xcalloc(1+n, sizeof(int));
1.136 + tree->dwrk = xcalloc(1+n, sizeof(double));
1.137 + /* initialize control parameters */
1.138 + tree->parm = parm;
1.139 + tree->tm_beg = xtime();
1.140 + tree->tm_lag = xlset(0);
1.141 + tree->sol_cnt = 0;
1.142 + /* initialize advanced solver interface */
1.143 + tree->reason = 0;
1.144 + tree->reopt = 0;
1.145 + tree->reinv = 0;
1.146 + tree->br_var = 0;
1.147 + tree->br_sel = 0;
1.148 + tree->child = 0;
1.149 + tree->next_p = 0;
1.150 + /*tree->btrack = NULL;*/
1.151 + tree->stop = 0;
1.152 + /* create the root subproblem, which initially is identical to
1.153 + the original MIP */
1.154 + new_node(tree, NULL);
1.155 + return tree;
1.156 +}
1.157 +
1.158 +/***********************************************************************
1.159 +* NAME
1.160 +*
1.161 +* ios_revive_node - revive specified subproblem
1.162 +*
1.163 +* SYNOPSIS
1.164 +*
1.165 +* #include "glpios.h"
1.166 +* void ios_revive_node(glp_tree *tree, int p);
1.167 +*
1.168 +* DESCRIPTION
1.169 +*
1.170 +* The routine ios_revive_node revives the specified subproblem, whose
1.171 +* reference number is p, and thereby makes it the current subproblem.
1.172 +* Note that the specified subproblem must be active. Besides, if the
1.173 +* current subproblem already exists, it must be frozen before reviving
1.174 +* another subproblem. */
1.175 +
1.176 +void ios_revive_node(glp_tree *tree, int p)
1.177 +{ glp_prob *mip = tree->mip;
1.178 + IOSNPD *node, *root;
1.179 + /* obtain pointer to the specified subproblem */
1.180 + xassert(1 <= p && p <= tree->nslots);
1.181 + node = tree->slot[p].node;
1.182 + xassert(node != NULL);
1.183 + /* the specified subproblem must be active */
1.184 + xassert(node->count == 0);
1.185 + /* the current subproblem must not exist */
1.186 + xassert(tree->curr == NULL);
1.187 + /* the specified subproblem becomes current */
1.188 + tree->curr = node;
1.189 + /*tree->solved = 0;*/
1.190 + /* obtain pointer to the root subproblem */
1.191 + root = tree->slot[1].node;
1.192 + xassert(root != NULL);
1.193 + /* at this point problem object components correspond to the root
1.194 + subproblem, so if the root subproblem should be revived, there
1.195 + is nothing more to do */
1.196 + if (node == root) goto done;
1.197 + xassert(mip->m == tree->root_m);
1.198 + /* build path from the root to the current node */
1.199 + node->temp = NULL;
1.200 + for (node = node; node != NULL; node = node->up)
1.201 + { if (node->up == NULL)
1.202 + xassert(node == root);
1.203 + else
1.204 + node->up->temp = node;
1.205 + }
1.206 + /* go down from the root to the current node and make necessary
1.207 + changes to restore components of the current subproblem */
1.208 + for (node = root; node != NULL; node = node->temp)
1.209 + { int m = mip->m;
1.210 + int n = mip->n;
1.211 + /* if the current node is reached, the problem object at this
1.212 + point corresponds to its parent, so save attributes of rows
1.213 + and columns for the parent subproblem */
1.214 + if (node->temp == NULL)
1.215 + { int i, j;
1.216 + tree->pred_m = m;
1.217 + /* allocate/reallocate arrays, if necessary */
1.218 + if (tree->pred_max < m + n)
1.219 + { int new_size = m + n + 100;
1.220 + if (tree->pred_type != NULL) xfree(tree->pred_type);
1.221 + if (tree->pred_lb != NULL) xfree(tree->pred_lb);
1.222 + if (tree->pred_ub != NULL) xfree(tree->pred_ub);
1.223 + if (tree->pred_stat != NULL) xfree(tree->pred_stat);
1.224 + tree->pred_max = new_size;
1.225 + tree->pred_type = xcalloc(1+new_size, sizeof(char));
1.226 + tree->pred_lb = xcalloc(1+new_size, sizeof(double));
1.227 + tree->pred_ub = xcalloc(1+new_size, sizeof(double));
1.228 + tree->pred_stat = xcalloc(1+new_size, sizeof(char));
1.229 + }
1.230 + /* save row attributes */
1.231 + for (i = 1; i <= m; i++)
1.232 + { GLPROW *row = mip->row[i];
1.233 + tree->pred_type[i] = (char)row->type;
1.234 + tree->pred_lb[i] = row->lb;
1.235 + tree->pred_ub[i] = row->ub;
1.236 + tree->pred_stat[i] = (char)row->stat;
1.237 + }
1.238 + /* save column attributes */
1.239 + for (j = 1; j <= n; j++)
1.240 + { GLPCOL *col = mip->col[j];
1.241 + tree->pred_type[mip->m+j] = (char)col->type;
1.242 + tree->pred_lb[mip->m+j] = col->lb;
1.243 + tree->pred_ub[mip->m+j] = col->ub;
1.244 + tree->pred_stat[mip->m+j] = (char)col->stat;
1.245 + }
1.246 + }
1.247 + /* change bounds of rows and columns */
1.248 + { IOSBND *b;
1.249 + for (b = node->b_ptr; b != NULL; b = b->next)
1.250 + { if (b->k <= m)
1.251 + glp_set_row_bnds(mip, b->k, b->type, b->lb, b->ub);
1.252 + else
1.253 + glp_set_col_bnds(mip, b->k-m, b->type, b->lb, b->ub);
1.254 + }
1.255 + }
1.256 + /* change statuses of rows and columns */
1.257 + { IOSTAT *s;
1.258 + for (s = node->s_ptr; s != NULL; s = s->next)
1.259 + { if (s->k <= m)
1.260 + glp_set_row_stat(mip, s->k, s->stat);
1.261 + else
1.262 + glp_set_col_stat(mip, s->k-m, s->stat);
1.263 + }
1.264 + }
1.265 + /* add new rows */
1.266 + if (node->r_ptr != NULL)
1.267 + { IOSROW *r;
1.268 + IOSAIJ *a;
1.269 + int i, len, *ind;
1.270 + double *val;
1.271 + ind = xcalloc(1+n, sizeof(int));
1.272 + val = xcalloc(1+n, sizeof(double));
1.273 + for (r = node->r_ptr; r != NULL; r = r->next)
1.274 + { i = glp_add_rows(mip, 1);
1.275 + glp_set_row_name(mip, i, r->name);
1.276 +#if 1 /* 20/IX-2008 */
1.277 + xassert(mip->row[i]->level == 0);
1.278 + mip->row[i]->level = node->level;
1.279 + mip->row[i]->origin = r->origin;
1.280 + mip->row[i]->klass = r->klass;
1.281 +#endif
1.282 + glp_set_row_bnds(mip, i, r->type, r->lb, r->ub);
1.283 + len = 0;
1.284 + for (a = r->ptr; a != NULL; a = a->next)
1.285 + len++, ind[len] = a->j, val[len] = a->val;
1.286 + glp_set_mat_row(mip, i, len, ind, val);
1.287 + glp_set_rii(mip, i, r->rii);
1.288 + glp_set_row_stat(mip, i, r->stat);
1.289 + }
1.290 + xfree(ind);
1.291 + xfree(val);
1.292 + }
1.293 +#if 0
1.294 + /* add new edges to the conflict graph */
1.295 + /* add new cliques to the conflict graph */
1.296 + /* (not implemented yet) */
1.297 + xassert(node->own_nn == 0);
1.298 + xassert(node->own_nc == 0);
1.299 + xassert(node->e_ptr == NULL);
1.300 +#endif
1.301 + }
1.302 + /* the specified subproblem has been revived */
1.303 + node = tree->curr;
1.304 + /* delete its bound change list */
1.305 + while (node->b_ptr != NULL)
1.306 + { IOSBND *b;
1.307 + b = node->b_ptr;
1.308 + node->b_ptr = b->next;
1.309 + dmp_free_atom(tree->pool, b, sizeof(IOSBND));
1.310 + }
1.311 + /* delete its status change list */
1.312 + while (node->s_ptr != NULL)
1.313 + { IOSTAT *s;
1.314 + s = node->s_ptr;
1.315 + node->s_ptr = s->next;
1.316 + dmp_free_atom(tree->pool, s, sizeof(IOSTAT));
1.317 + }
1.318 +#if 1 /* 20/XI-2009 */
1.319 + /* delete its row addition list (additional rows may appear, for
1.320 + example, due to branching on GUB constraints */
1.321 + while (node->r_ptr != NULL)
1.322 + { IOSROW *r;
1.323 + r = node->r_ptr;
1.324 + node->r_ptr = r->next;
1.325 + xassert(r->name == NULL);
1.326 + while (r->ptr != NULL)
1.327 + { IOSAIJ *a;
1.328 + a = r->ptr;
1.329 + r->ptr = a->next;
1.330 + dmp_free_atom(tree->pool, a, sizeof(IOSAIJ));
1.331 + }
1.332 + dmp_free_atom(tree->pool, r, sizeof(IOSROW));
1.333 + }
1.334 +#endif
1.335 +done: return;
1.336 +}
1.337 +
1.338 +/***********************************************************************
1.339 +* NAME
1.340 +*
1.341 +* ios_freeze_node - freeze current subproblem
1.342 +*
1.343 +* SYNOPSIS
1.344 +*
1.345 +* #include "glpios.h"
1.346 +* void ios_freeze_node(glp_tree *tree);
1.347 +*
1.348 +* DESCRIPTION
1.349 +*
1.350 +* The routine ios_freeze_node freezes the current subproblem. */
1.351 +
1.352 +void ios_freeze_node(glp_tree *tree)
1.353 +{ glp_prob *mip = tree->mip;
1.354 + int m = mip->m;
1.355 + int n = mip->n;
1.356 + IOSNPD *node;
1.357 + /* obtain pointer to the current subproblem */
1.358 + node = tree->curr;
1.359 + xassert(node != NULL);
1.360 + if (node->up == NULL)
1.361 + { /* freeze the root subproblem */
1.362 + int k;
1.363 + xassert(node->p == 1);
1.364 + xassert(tree->root_m == 0);
1.365 + xassert(tree->root_type == NULL);
1.366 + xassert(tree->root_lb == NULL);
1.367 + xassert(tree->root_ub == NULL);
1.368 + xassert(tree->root_stat == NULL);
1.369 + tree->root_m = m;
1.370 + tree->root_type = xcalloc(1+m+n, sizeof(char));
1.371 + tree->root_lb = xcalloc(1+m+n, sizeof(double));
1.372 + tree->root_ub = xcalloc(1+m+n, sizeof(double));
1.373 + tree->root_stat = xcalloc(1+m+n, sizeof(char));
1.374 + for (k = 1; k <= m+n; k++)
1.375 + { if (k <= m)
1.376 + { GLPROW *row = mip->row[k];
1.377 + tree->root_type[k] = (char)row->type;
1.378 + tree->root_lb[k] = row->lb;
1.379 + tree->root_ub[k] = row->ub;
1.380 + tree->root_stat[k] = (char)row->stat;
1.381 + }
1.382 + else
1.383 + { GLPCOL *col = mip->col[k-m];
1.384 + tree->root_type[k] = (char)col->type;
1.385 + tree->root_lb[k] = col->lb;
1.386 + tree->root_ub[k] = col->ub;
1.387 + tree->root_stat[k] = (char)col->stat;
1.388 + }
1.389 + }
1.390 + }
1.391 + else
1.392 + { /* freeze non-root subproblem */
1.393 + int root_m = tree->root_m;
1.394 + int pred_m = tree->pred_m;
1.395 + int i, j, k;
1.396 + xassert(pred_m <= m);
1.397 + /* build change lists for rows and columns which exist in the
1.398 + parent subproblem */
1.399 + xassert(node->b_ptr == NULL);
1.400 + xassert(node->s_ptr == NULL);
1.401 + for (k = 1; k <= pred_m + n; k++)
1.402 + { int pred_type, pred_stat, type, stat;
1.403 + double pred_lb, pred_ub, lb, ub;
1.404 + /* determine attributes in the parent subproblem */
1.405 + pred_type = tree->pred_type[k];
1.406 + pred_lb = tree->pred_lb[k];
1.407 + pred_ub = tree->pred_ub[k];
1.408 + pred_stat = tree->pred_stat[k];
1.409 + /* determine attributes in the current subproblem */
1.410 + if (k <= pred_m)
1.411 + { GLPROW *row = mip->row[k];
1.412 + type = row->type;
1.413 + lb = row->lb;
1.414 + ub = row->ub;
1.415 + stat = row->stat;
1.416 + }
1.417 + else
1.418 + { GLPCOL *col = mip->col[k - pred_m];
1.419 + type = col->type;
1.420 + lb = col->lb;
1.421 + ub = col->ub;
1.422 + stat = col->stat;
1.423 + }
1.424 + /* save type and bounds of a row/column, if changed */
1.425 + if (!(pred_type == type && pred_lb == lb && pred_ub == ub))
1.426 + { IOSBND *b;
1.427 + b = dmp_get_atom(tree->pool, sizeof(IOSBND));
1.428 + b->k = k;
1.429 + b->type = (unsigned char)type;
1.430 + b->lb = lb;
1.431 + b->ub = ub;
1.432 + b->next = node->b_ptr;
1.433 + node->b_ptr = b;
1.434 + }
1.435 + /* save status of a row/column, if changed */
1.436 + if (pred_stat != stat)
1.437 + { IOSTAT *s;
1.438 + s = dmp_get_atom(tree->pool, sizeof(IOSTAT));
1.439 + s->k = k;
1.440 + s->stat = (unsigned char)stat;
1.441 + s->next = node->s_ptr;
1.442 + node->s_ptr = s;
1.443 + }
1.444 + }
1.445 + /* save new rows added to the current subproblem */
1.446 + xassert(node->r_ptr == NULL);
1.447 + if (pred_m < m)
1.448 + { int i, len, *ind;
1.449 + double *val;
1.450 + ind = xcalloc(1+n, sizeof(int));
1.451 + val = xcalloc(1+n, sizeof(double));
1.452 + for (i = m; i > pred_m; i--)
1.453 + { GLPROW *row = mip->row[i];
1.454 + IOSROW *r;
1.455 + const char *name;
1.456 + r = dmp_get_atom(tree->pool, sizeof(IOSROW));
1.457 + name = glp_get_row_name(mip, i);
1.458 + if (name == NULL)
1.459 + r->name = NULL;
1.460 + else
1.461 + { r->name = dmp_get_atom(tree->pool, strlen(name)+1);
1.462 + strcpy(r->name, name);
1.463 + }
1.464 +#if 1 /* 20/IX-2008 */
1.465 + r->origin = row->origin;
1.466 + r->klass = row->klass;
1.467 +#endif
1.468 + r->type = (unsigned char)row->type;
1.469 + r->lb = row->lb;
1.470 + r->ub = row->ub;
1.471 + r->ptr = NULL;
1.472 + len = glp_get_mat_row(mip, i, ind, val);
1.473 + for (k = 1; k <= len; k++)
1.474 + { IOSAIJ *a;
1.475 + a = dmp_get_atom(tree->pool, sizeof(IOSAIJ));
1.476 + a->j = ind[k];
1.477 + a->val = val[k];
1.478 + a->next = r->ptr;
1.479 + r->ptr = a;
1.480 + }
1.481 + r->rii = row->rii;
1.482 + r->stat = (unsigned char)row->stat;
1.483 + r->next = node->r_ptr;
1.484 + node->r_ptr = r;
1.485 + }
1.486 + xfree(ind);
1.487 + xfree(val);
1.488 + }
1.489 + /* remove all rows missing in the root subproblem */
1.490 + if (m != root_m)
1.491 + { int nrs, *num;
1.492 + nrs = m - root_m;
1.493 + xassert(nrs > 0);
1.494 + num = xcalloc(1+nrs, sizeof(int));
1.495 + for (i = 1; i <= nrs; i++) num[i] = root_m + i;
1.496 + glp_del_rows(mip, nrs, num);
1.497 + xfree(num);
1.498 + }
1.499 + m = mip->m;
1.500 + /* and restore attributes of all rows and columns for the root
1.501 + subproblem */
1.502 + xassert(m == root_m);
1.503 + for (i = 1; i <= m; i++)
1.504 + { glp_set_row_bnds(mip, i, tree->root_type[i],
1.505 + tree->root_lb[i], tree->root_ub[i]);
1.506 + glp_set_row_stat(mip, i, tree->root_stat[i]);
1.507 + }
1.508 + for (j = 1; j <= n; j++)
1.509 + { glp_set_col_bnds(mip, j, tree->root_type[m+j],
1.510 + tree->root_lb[m+j], tree->root_ub[m+j]);
1.511 + glp_set_col_stat(mip, j, tree->root_stat[m+j]);
1.512 + }
1.513 +#if 1
1.514 + /* remove all edges and cliques missing in the conflict graph
1.515 + for the root subproblem */
1.516 + /* (not implemented yet) */
1.517 +#endif
1.518 + }
1.519 + /* the current subproblem has been frozen */
1.520 + tree->curr = NULL;
1.521 + return;
1.522 +}
1.523 +
1.524 +/***********************************************************************
1.525 +* NAME
1.526 +*
1.527 +* ios_clone_node - clone specified subproblem
1.528 +*
1.529 +* SYNOPSIS
1.530 +*
1.531 +* #include "glpios.h"
1.532 +* void ios_clone_node(glp_tree *tree, int p, int nnn, int ref[]);
1.533 +*
1.534 +* DESCRIPTION
1.535 +*
1.536 +* The routine ios_clone_node clones the specified subproblem, whose
1.537 +* reference number is p, creating its nnn exact copies. Note that the
1.538 +* specified subproblem must be active and must be in the frozen state
1.539 +* (i.e. it must not be the current subproblem).
1.540 +*
1.541 +* Each clone, an exact copy of the specified subproblem, becomes a new
1.542 +* active subproblem added to the end of the active list. After cloning
1.543 +* the specified subproblem becomes inactive.
1.544 +*
1.545 +* The reference numbers of clone subproblems are stored to locations
1.546 +* ref[1], ..., ref[nnn]. */
1.547 +
1.548 +static int get_slot(glp_tree *tree)
1.549 +{ int p;
1.550 + /* if no free slots are available, increase the room */
1.551 + if (tree->avail == 0)
1.552 + { int nslots = tree->nslots;
1.553 + IOSLOT *save = tree->slot;
1.554 + if (nslots == 0)
1.555 + tree->nslots = 20;
1.556 + else
1.557 + { tree->nslots = nslots + nslots;
1.558 + xassert(tree->nslots > nslots);
1.559 + }
1.560 + tree->slot = xcalloc(1+tree->nslots, sizeof(IOSLOT));
1.561 + if (save != NULL)
1.562 + { memcpy(&tree->slot[1], &save[1], nslots * sizeof(IOSLOT));
1.563 + xfree(save);
1.564 + }
1.565 + /* push more free slots into the stack */
1.566 + for (p = tree->nslots; p > nslots; p--)
1.567 + { tree->slot[p].node = NULL;
1.568 + tree->slot[p].next = tree->avail;
1.569 + tree->avail = p;
1.570 + }
1.571 + }
1.572 + /* pull a free slot from the stack */
1.573 + p = tree->avail;
1.574 + tree->avail = tree->slot[p].next;
1.575 + xassert(tree->slot[p].node == NULL);
1.576 + tree->slot[p].next = 0;
1.577 + return p;
1.578 +}
1.579 +
1.580 +static IOSNPD *new_node(glp_tree *tree, IOSNPD *parent)
1.581 +{ IOSNPD *node;
1.582 + int p;
1.583 + /* pull a free slot for the new node */
1.584 + p = get_slot(tree);
1.585 + /* create descriptor of the new subproblem */
1.586 + node = dmp_get_atom(tree->pool, sizeof(IOSNPD));
1.587 + tree->slot[p].node = node;
1.588 + node->p = p;
1.589 + node->up = parent;
1.590 + node->level = (parent == NULL ? 0 : parent->level + 1);
1.591 + node->count = 0;
1.592 + node->b_ptr = NULL;
1.593 + node->s_ptr = NULL;
1.594 + node->r_ptr = NULL;
1.595 + node->solved = 0;
1.596 +#if 0
1.597 + node->own_nn = node->own_nc = 0;
1.598 + node->e_ptr = NULL;
1.599 +#endif
1.600 +#if 1 /* 04/X-2008 */
1.601 + node->lp_obj = (parent == NULL ? (tree->mip->dir == GLP_MIN ?
1.602 + -DBL_MAX : +DBL_MAX) : parent->lp_obj);
1.603 +#endif
1.604 + node->bound = (parent == NULL ? (tree->mip->dir == GLP_MIN ?
1.605 + -DBL_MAX : +DBL_MAX) : parent->bound);
1.606 + node->br_var = 0;
1.607 + node->br_val = 0.0;
1.608 + node->ii_cnt = 0;
1.609 + node->ii_sum = 0.0;
1.610 +#if 1 /* 30/XI-2009 */
1.611 + node->changed = 0;
1.612 +#endif
1.613 + if (tree->parm->cb_size == 0)
1.614 + node->data = NULL;
1.615 + else
1.616 + { node->data = dmp_get_atom(tree->pool, tree->parm->cb_size);
1.617 + memset(node->data, 0, tree->parm->cb_size);
1.618 + }
1.619 + node->temp = NULL;
1.620 + node->prev = tree->tail;
1.621 + node->next = NULL;
1.622 + /* add the new subproblem to the end of the active list */
1.623 + if (tree->head == NULL)
1.624 + tree->head = node;
1.625 + else
1.626 + tree->tail->next = node;
1.627 + tree->tail = node;
1.628 + tree->a_cnt++;
1.629 + tree->n_cnt++;
1.630 + tree->t_cnt++;
1.631 + /* increase the number of child subproblems */
1.632 + if (parent == NULL)
1.633 + xassert(p == 1);
1.634 + else
1.635 + parent->count++;
1.636 + return node;
1.637 +}
1.638 +
1.639 +void ios_clone_node(glp_tree *tree, int p, int nnn, int ref[])
1.640 +{ IOSNPD *node;
1.641 + int k;
1.642 + /* obtain pointer to the subproblem to be cloned */
1.643 + xassert(1 <= p && p <= tree->nslots);
1.644 + node = tree->slot[p].node;
1.645 + xassert(node != NULL);
1.646 + /* the specified subproblem must be active */
1.647 + xassert(node->count == 0);
1.648 + /* and must be in the frozen state */
1.649 + xassert(tree->curr != node);
1.650 + /* remove the specified subproblem from the active list, because
1.651 + it becomes inactive */
1.652 + if (node->prev == NULL)
1.653 + tree->head = node->next;
1.654 + else
1.655 + node->prev->next = node->next;
1.656 + if (node->next == NULL)
1.657 + tree->tail = node->prev;
1.658 + else
1.659 + node->next->prev = node->prev;
1.660 + node->prev = node->next = NULL;
1.661 + tree->a_cnt--;
1.662 + /* create clone subproblems */
1.663 + xassert(nnn > 0);
1.664 + for (k = 1; k <= nnn; k++)
1.665 + ref[k] = new_node(tree, node)->p;
1.666 + return;
1.667 +}
1.668 +
1.669 +/***********************************************************************
1.670 +* NAME
1.671 +*
1.672 +* ios_delete_node - delete specified subproblem
1.673 +*
1.674 +* SYNOPSIS
1.675 +*
1.676 +* #include "glpios.h"
1.677 +* void ios_delete_node(glp_tree *tree, int p);
1.678 +*
1.679 +* DESCRIPTION
1.680 +*
1.681 +* The routine ios_delete_node deletes the specified subproblem, whose
1.682 +* reference number is p. The subproblem must be active and must be in
1.683 +* the frozen state (i.e. it must not be the current subproblem).
1.684 +*
1.685 +* Note that deletion is performed recursively, i.e. if a subproblem to
1.686 +* be deleted is the only child of its parent, the parent subproblem is
1.687 +* also deleted, etc. */
1.688 +
1.689 +void ios_delete_node(glp_tree *tree, int p)
1.690 +{ IOSNPD *node, *temp;
1.691 + /* obtain pointer to the subproblem to be deleted */
1.692 + xassert(1 <= p && p <= tree->nslots);
1.693 + node = tree->slot[p].node;
1.694 + xassert(node != NULL);
1.695 + /* the specified subproblem must be active */
1.696 + xassert(node->count == 0);
1.697 + /* and must be in the frozen state */
1.698 + xassert(tree->curr != node);
1.699 + /* remove the specified subproblem from the active list, because
1.700 + it is gone from the tree */
1.701 + if (node->prev == NULL)
1.702 + tree->head = node->next;
1.703 + else
1.704 + node->prev->next = node->next;
1.705 + if (node->next == NULL)
1.706 + tree->tail = node->prev;
1.707 + else
1.708 + node->next->prev = node->prev;
1.709 + node->prev = node->next = NULL;
1.710 + tree->a_cnt--;
1.711 +loop: /* recursive deletion starts here */
1.712 + /* delete the bound change list */
1.713 + { IOSBND *b;
1.714 + while (node->b_ptr != NULL)
1.715 + { b = node->b_ptr;
1.716 + node->b_ptr = b->next;
1.717 + dmp_free_atom(tree->pool, b, sizeof(IOSBND));
1.718 + }
1.719 + }
1.720 + /* delete the status change list */
1.721 + { IOSTAT *s;
1.722 + while (node->s_ptr != NULL)
1.723 + { s = node->s_ptr;
1.724 + node->s_ptr = s->next;
1.725 + dmp_free_atom(tree->pool, s, sizeof(IOSTAT));
1.726 + }
1.727 + }
1.728 + /* delete the row addition list */
1.729 + while (node->r_ptr != NULL)
1.730 + { IOSROW *r;
1.731 + r = node->r_ptr;
1.732 + if (r->name != NULL)
1.733 + dmp_free_atom(tree->pool, r->name, strlen(r->name)+1);
1.734 + while (r->ptr != NULL)
1.735 + { IOSAIJ *a;
1.736 + a = r->ptr;
1.737 + r->ptr = a->next;
1.738 + dmp_free_atom(tree->pool, a, sizeof(IOSAIJ));
1.739 + }
1.740 + node->r_ptr = r->next;
1.741 + dmp_free_atom(tree->pool, r, sizeof(IOSROW));
1.742 + }
1.743 +#if 0
1.744 + /* delete the edge addition list */
1.745 + /* delete the clique addition list */
1.746 + /* (not implemented yet) */
1.747 + xassert(node->own_nn == 0);
1.748 + xassert(node->own_nc == 0);
1.749 + xassert(node->e_ptr == NULL);
1.750 +#endif
1.751 + /* free application-specific data */
1.752 + if (tree->parm->cb_size == 0)
1.753 + xassert(node->data == NULL);
1.754 + else
1.755 + dmp_free_atom(tree->pool, node->data, tree->parm->cb_size);
1.756 + /* free the corresponding node slot */
1.757 + p = node->p;
1.758 + xassert(tree->slot[p].node == node);
1.759 + tree->slot[p].node = NULL;
1.760 + tree->slot[p].next = tree->avail;
1.761 + tree->avail = p;
1.762 + /* save pointer to the parent subproblem */
1.763 + temp = node->up;
1.764 + /* delete the subproblem descriptor */
1.765 + dmp_free_atom(tree->pool, node, sizeof(IOSNPD));
1.766 + tree->n_cnt--;
1.767 + /* take pointer to the parent subproblem */
1.768 + node = temp;
1.769 + if (node != NULL)
1.770 + { /* the parent subproblem exists; decrease the number of its
1.771 + child subproblems */
1.772 + xassert(node->count > 0);
1.773 + node->count--;
1.774 + /* if now the parent subproblem has no childs, it also must be
1.775 + deleted */
1.776 + if (node->count == 0) goto loop;
1.777 + }
1.778 + return;
1.779 +}
1.780 +
1.781 +/***********************************************************************
1.782 +* NAME
1.783 +*
1.784 +* ios_delete_tree - delete branch-and-bound tree
1.785 +*
1.786 +* SYNOPSIS
1.787 +*
1.788 +* #include "glpios.h"
1.789 +* void ios_delete_tree(glp_tree *tree);
1.790 +*
1.791 +* DESCRIPTION
1.792 +*
1.793 +* The routine ios_delete_tree deletes the branch-and-bound tree, which
1.794 +* the parameter tree points to, and frees all the memory allocated to
1.795 +* this program object.
1.796 +*
1.797 +* On exit components of the problem object are restored to correspond
1.798 +* to the original MIP passed to the routine ios_create_tree. */
1.799 +
1.800 +void ios_delete_tree(glp_tree *tree)
1.801 +{ glp_prob *mip = tree->mip;
1.802 + int i, j;
1.803 + int m = mip->m;
1.804 + int n = mip->n;
1.805 + xassert(mip->tree == tree);
1.806 + /* remove all additional rows */
1.807 + if (m != tree->orig_m)
1.808 + { int nrs, *num;
1.809 + nrs = m - tree->orig_m;
1.810 + xassert(nrs > 0);
1.811 + num = xcalloc(1+nrs, sizeof(int));
1.812 + for (i = 1; i <= nrs; i++) num[i] = tree->orig_m + i;
1.813 + glp_del_rows(mip, nrs, num);
1.814 + xfree(num);
1.815 + }
1.816 + m = tree->orig_m;
1.817 + /* restore original attributes of rows and columns */
1.818 + xassert(m == tree->orig_m);
1.819 + xassert(n == tree->n);
1.820 + for (i = 1; i <= m; i++)
1.821 + { glp_set_row_bnds(mip, i, tree->orig_type[i],
1.822 + tree->orig_lb[i], tree->orig_ub[i]);
1.823 + glp_set_row_stat(mip, i, tree->orig_stat[i]);
1.824 + mip->row[i]->prim = tree->orig_prim[i];
1.825 + mip->row[i]->dual = tree->orig_dual[i];
1.826 + }
1.827 + for (j = 1; j <= n; j++)
1.828 + { glp_set_col_bnds(mip, j, tree->orig_type[m+j],
1.829 + tree->orig_lb[m+j], tree->orig_ub[m+j]);
1.830 + glp_set_col_stat(mip, j, tree->orig_stat[m+j]);
1.831 + mip->col[j]->prim = tree->orig_prim[m+j];
1.832 + mip->col[j]->dual = tree->orig_dual[m+j];
1.833 + }
1.834 + mip->pbs_stat = mip->dbs_stat = GLP_FEAS;
1.835 + mip->obj_val = tree->orig_obj;
1.836 + /* delete the branch-and-bound tree */
1.837 + xassert(tree->local != NULL);
1.838 + ios_delete_pool(tree, tree->local);
1.839 + dmp_delete_pool(tree->pool);
1.840 + xfree(tree->orig_type);
1.841 + xfree(tree->orig_lb);
1.842 + xfree(tree->orig_ub);
1.843 + xfree(tree->orig_stat);
1.844 + xfree(tree->orig_prim);
1.845 + xfree(tree->orig_dual);
1.846 + xfree(tree->slot);
1.847 + if (tree->root_type != NULL) xfree(tree->root_type);
1.848 + if (tree->root_lb != NULL) xfree(tree->root_lb);
1.849 + if (tree->root_ub != NULL) xfree(tree->root_ub);
1.850 + if (tree->root_stat != NULL) xfree(tree->root_stat);
1.851 + xfree(tree->non_int);
1.852 +#if 0
1.853 + xfree(tree->n_ref);
1.854 + xfree(tree->c_ref);
1.855 + xfree(tree->j_ref);
1.856 +#endif
1.857 + if (tree->pcost != NULL) ios_pcost_free(tree);
1.858 + xfree(tree->iwrk);
1.859 + xfree(tree->dwrk);
1.860 +#if 0
1.861 + scg_delete_graph(tree->g);
1.862 +#endif
1.863 + if (tree->pred_type != NULL) xfree(tree->pred_type);
1.864 + if (tree->pred_lb != NULL) xfree(tree->pred_lb);
1.865 + if (tree->pred_ub != NULL) xfree(tree->pred_ub);
1.866 + if (tree->pred_stat != NULL) xfree(tree->pred_stat);
1.867 +#if 0
1.868 + xassert(tree->cut_gen == NULL);
1.869 +#endif
1.870 + xassert(tree->mir_gen == NULL);
1.871 + xassert(tree->clq_gen == NULL);
1.872 + xfree(tree);
1.873 + mip->tree = NULL;
1.874 + return;
1.875 +}
1.876 +
1.877 +/***********************************************************************
1.878 +* NAME
1.879 +*
1.880 +* ios_eval_degrad - estimate obj. degrad. for down- and up-branches
1.881 +*
1.882 +* SYNOPSIS
1.883 +*
1.884 +* #include "glpios.h"
1.885 +* void ios_eval_degrad(glp_tree *tree, int j, double *dn, double *up);
1.886 +*
1.887 +* DESCRIPTION
1.888 +*
1.889 +* Given optimal basis to LP relaxation of the current subproblem the
1.890 +* routine ios_eval_degrad performs the dual ratio test to compute the
1.891 +* objective values in the adjacent basis for down- and up-branches,
1.892 +* which are stored in locations *dn and *up, assuming that x[j] is a
1.893 +* variable chosen to branch upon. */
1.894 +
1.895 +void ios_eval_degrad(glp_tree *tree, int j, double *dn, double *up)
1.896 +{ glp_prob *mip = tree->mip;
1.897 + int m = mip->m, n = mip->n;
1.898 + int len, kase, k, t, stat;
1.899 + double alfa, beta, gamma, delta, dz;
1.900 + int *ind = tree->iwrk;
1.901 + double *val = tree->dwrk;
1.902 + /* current basis must be optimal */
1.903 + xassert(glp_get_status(mip) == GLP_OPT);
1.904 + /* basis factorization must exist */
1.905 + xassert(glp_bf_exists(mip));
1.906 + /* obtain (fractional) value of x[j] in optimal basic solution
1.907 + to LP relaxation of the current subproblem */
1.908 + xassert(1 <= j && j <= n);
1.909 + beta = mip->col[j]->prim;
1.910 + /* since the value of x[j] is fractional, it is basic; compute
1.911 + corresponding row of the simplex table */
1.912 + len = lpx_eval_tab_row(mip, m+j, ind, val);
1.913 + /* kase < 0 means down-branch; kase > 0 means up-branch */
1.914 + for (kase = -1; kase <= +1; kase += 2)
1.915 + { /* for down-branch we introduce new upper bound floor(beta)
1.916 + for x[j]; similarly, for up-branch we introduce new lower
1.917 + bound ceil(beta) for x[j]; in the current basis this new
1.918 + upper/lower bound is violated, so in the adjacent basis
1.919 + x[j] will leave the basis and go to its new upper/lower
1.920 + bound; we need to know which non-basic variable x[k] should
1.921 + enter the basis to keep dual feasibility */
1.922 +#if 0 /* 23/XI-2009 */
1.923 + k = lpx_dual_ratio_test(mip, len, ind, val, kase, 1e-7);
1.924 +#else
1.925 + k = lpx_dual_ratio_test(mip, len, ind, val, kase, 1e-9);
1.926 +#endif
1.927 + /* if no variable has been chosen, current basis being primal
1.928 + infeasible due to the new upper/lower bound of x[j] is dual
1.929 + unbounded, therefore, LP relaxation to corresponding branch
1.930 + has no primal feasible solution */
1.931 + if (k == 0)
1.932 + { if (mip->dir == GLP_MIN)
1.933 + { if (kase < 0)
1.934 + *dn = +DBL_MAX;
1.935 + else
1.936 + *up = +DBL_MAX;
1.937 + }
1.938 + else if (mip->dir == GLP_MAX)
1.939 + { if (kase < 0)
1.940 + *dn = -DBL_MAX;
1.941 + else
1.942 + *up = -DBL_MAX;
1.943 + }
1.944 + else
1.945 + xassert(mip != mip);
1.946 + continue;
1.947 + }
1.948 + xassert(1 <= k && k <= m+n);
1.949 + /* row of the simplex table corresponding to specified basic
1.950 + variable x[j] is the following:
1.951 + x[j] = ... + alfa * x[k] + ... ;
1.952 + we need to know influence coefficient, alfa, at non-basic
1.953 + variable x[k] chosen with the dual ratio test */
1.954 + for (t = 1; t <= len; t++)
1.955 + if (ind[t] == k) break;
1.956 + xassert(1 <= t && t <= len);
1.957 + alfa = val[t];
1.958 + /* determine status and reduced cost of variable x[k] */
1.959 + if (k <= m)
1.960 + { stat = mip->row[k]->stat;
1.961 + gamma = mip->row[k]->dual;
1.962 + }
1.963 + else
1.964 + { stat = mip->col[k-m]->stat;
1.965 + gamma = mip->col[k-m]->dual;
1.966 + }
1.967 + /* x[k] cannot be basic or fixed non-basic */
1.968 + xassert(stat == GLP_NL || stat == GLP_NU || stat == GLP_NF);
1.969 + /* if the current basis is dual degenerative, some reduced
1.970 + costs, which are close to zero, may have wrong sign due to
1.971 + round-off errors, so correct the sign of gamma */
1.972 + if (mip->dir == GLP_MIN)
1.973 + { if (stat == GLP_NL && gamma < 0.0 ||
1.974 + stat == GLP_NU && gamma > 0.0 ||
1.975 + stat == GLP_NF) gamma = 0.0;
1.976 + }
1.977 + else if (mip->dir == GLP_MAX)
1.978 + { if (stat == GLP_NL && gamma > 0.0 ||
1.979 + stat == GLP_NU && gamma < 0.0 ||
1.980 + stat == GLP_NF) gamma = 0.0;
1.981 + }
1.982 + else
1.983 + xassert(mip != mip);
1.984 + /* determine the change of x[j] in the adjacent basis:
1.985 + delta x[j] = new x[j] - old x[j] */
1.986 + delta = (kase < 0 ? floor(beta) : ceil(beta)) - beta;
1.987 + /* compute the change of x[k] in the adjacent basis:
1.988 + delta x[k] = new x[k] - old x[k] = delta x[j] / alfa */
1.989 + delta /= alfa;
1.990 + /* compute the change of the objective in the adjacent basis:
1.991 + delta z = new z - old z = gamma * delta x[k] */
1.992 + dz = gamma * delta;
1.993 + if (mip->dir == GLP_MIN)
1.994 + xassert(dz >= 0.0);
1.995 + else if (mip->dir == GLP_MAX)
1.996 + xassert(dz <= 0.0);
1.997 + else
1.998 + xassert(mip != mip);
1.999 + /* compute the new objective value in the adjacent basis:
1.1000 + new z = old z + delta z */
1.1001 + if (kase < 0)
1.1002 + *dn = mip->obj_val + dz;
1.1003 + else
1.1004 + *up = mip->obj_val + dz;
1.1005 + }
1.1006 + /*xprintf("obj = %g; dn = %g; up = %g\n",
1.1007 + mip->obj_val, *dn, *up);*/
1.1008 + return;
1.1009 +}
1.1010 +
1.1011 +/***********************************************************************
1.1012 +* NAME
1.1013 +*
1.1014 +* ios_round_bound - improve local bound by rounding
1.1015 +*
1.1016 +* SYNOPSIS
1.1017 +*
1.1018 +* #include "glpios.h"
1.1019 +* double ios_round_bound(glp_tree *tree, double bound);
1.1020 +*
1.1021 +* RETURNS
1.1022 +*
1.1023 +* For the given local bound for any integer feasible solution to the
1.1024 +* current subproblem the routine ios_round_bound returns an improved
1.1025 +* local bound for the same integer feasible solution.
1.1026 +*
1.1027 +* BACKGROUND
1.1028 +*
1.1029 +* Let the current subproblem has the following objective function:
1.1030 +*
1.1031 +* z = sum c[j] * x[j] + s >= b, (1)
1.1032 +* j in J
1.1033 +*
1.1034 +* where J = {j: c[j] is non-zero and integer, x[j] is integer}, s is
1.1035 +* the sum of terms corresponding to fixed variables, b is an initial
1.1036 +* local bound (minimization).
1.1037 +*
1.1038 +* From (1) it follows that:
1.1039 +*
1.1040 +* d * sum (c[j] / d) * x[j] + s >= b, (2)
1.1041 +* j in J
1.1042 +*
1.1043 +* or, equivalently,
1.1044 +*
1.1045 +* sum (c[j] / d) * x[j] >= (b - s) / d = h, (3)
1.1046 +* j in J
1.1047 +*
1.1048 +* where d = gcd(c[j]). Since the left-hand side of (3) is integer,
1.1049 +* h = (b - s) / d can be rounded up to the nearest integer:
1.1050 +*
1.1051 +* h' = ceil(h) = (b' - s) / d, (4)
1.1052 +*
1.1053 +* that gives an rounded, improved local bound:
1.1054 +*
1.1055 +* b' = d * h' + s. (5)
1.1056 +*
1.1057 +* In case of maximization '>=' in (1) should be replaced by '<=' that
1.1058 +* leads to the following formula:
1.1059 +*
1.1060 +* h' = floor(h) = (b' - s) / d, (6)
1.1061 +*
1.1062 +* which should used in the same way as (4).
1.1063 +*
1.1064 +* NOTE: If b is a valid local bound for a child of the current
1.1065 +* subproblem, b' is also valid for that child subproblem. */
1.1066 +
1.1067 +double ios_round_bound(glp_tree *tree, double bound)
1.1068 +{ glp_prob *mip = tree->mip;
1.1069 + int n = mip->n;
1.1070 + int d, j, nn, *c = tree->iwrk;
1.1071 + double s, h;
1.1072 + /* determine c[j] and compute s */
1.1073 + nn = 0, s = mip->c0, d = 0;
1.1074 + for (j = 1; j <= n; j++)
1.1075 + { GLPCOL *col = mip->col[j];
1.1076 + if (col->coef == 0.0) continue;
1.1077 + if (col->type == GLP_FX)
1.1078 + { /* fixed variable */
1.1079 + s += col->coef * col->prim;
1.1080 + }
1.1081 + else
1.1082 + { /* non-fixed variable */
1.1083 + if (col->kind != GLP_IV) goto skip;
1.1084 + if (col->coef != floor(col->coef)) goto skip;
1.1085 + if (fabs(col->coef) <= (double)INT_MAX)
1.1086 + c[++nn] = (int)fabs(col->coef);
1.1087 + else
1.1088 + d = 1;
1.1089 + }
1.1090 + }
1.1091 + /* compute d = gcd(c[1],...c[nn]) */
1.1092 + if (d == 0)
1.1093 + { if (nn == 0) goto skip;
1.1094 + d = gcdn(nn, c);
1.1095 + }
1.1096 + xassert(d > 0);
1.1097 + /* compute new local bound */
1.1098 + if (mip->dir == GLP_MIN)
1.1099 + { if (bound != +DBL_MAX)
1.1100 + { h = (bound - s) / (double)d;
1.1101 + if (h >= floor(h) + 0.001)
1.1102 + { /* round up */
1.1103 + h = ceil(h);
1.1104 + /*xprintf("d = %d; old = %g; ", d, bound);*/
1.1105 + bound = (double)d * h + s;
1.1106 + /*xprintf("new = %g\n", bound);*/
1.1107 + }
1.1108 + }
1.1109 + }
1.1110 + else if (mip->dir == GLP_MAX)
1.1111 + { if (bound != -DBL_MAX)
1.1112 + { h = (bound - s) / (double)d;
1.1113 + if (h <= ceil(h) - 0.001)
1.1114 + { /* round down */
1.1115 + h = floor(h);
1.1116 + bound = (double)d * h + s;
1.1117 + }
1.1118 + }
1.1119 + }
1.1120 + else
1.1121 + xassert(mip != mip);
1.1122 +skip: return bound;
1.1123 +}
1.1124 +
1.1125 +/***********************************************************************
1.1126 +* NAME
1.1127 +*
1.1128 +* ios_is_hopeful - check if subproblem is hopeful
1.1129 +*
1.1130 +* SYNOPSIS
1.1131 +*
1.1132 +* #include "glpios.h"
1.1133 +* int ios_is_hopeful(glp_tree *tree, double bound);
1.1134 +*
1.1135 +* DESCRIPTION
1.1136 +*
1.1137 +* Given the local bound of a subproblem the routine ios_is_hopeful
1.1138 +* checks if the subproblem can have an integer optimal solution which
1.1139 +* is better than the best one currently known.
1.1140 +*
1.1141 +* RETURNS
1.1142 +*
1.1143 +* If the subproblem can have a better integer optimal solution, the
1.1144 +* routine returns non-zero; otherwise, if the corresponding branch can
1.1145 +* be pruned, the routine returns zero. */
1.1146 +
1.1147 +int ios_is_hopeful(glp_tree *tree, double bound)
1.1148 +{ glp_prob *mip = tree->mip;
1.1149 + int ret = 1;
1.1150 + double eps;
1.1151 + if (mip->mip_stat == GLP_FEAS)
1.1152 + { eps = tree->parm->tol_obj * (1.0 + fabs(mip->mip_obj));
1.1153 + switch (mip->dir)
1.1154 + { case GLP_MIN:
1.1155 + if (bound >= mip->mip_obj - eps) ret = 0;
1.1156 + break;
1.1157 + case GLP_MAX:
1.1158 + if (bound <= mip->mip_obj + eps) ret = 0;
1.1159 + break;
1.1160 + default:
1.1161 + xassert(mip != mip);
1.1162 + }
1.1163 + }
1.1164 + else
1.1165 + { switch (mip->dir)
1.1166 + { case GLP_MIN:
1.1167 + if (bound == +DBL_MAX) ret = 0;
1.1168 + break;
1.1169 + case GLP_MAX:
1.1170 + if (bound == -DBL_MAX) ret = 0;
1.1171 + break;
1.1172 + default:
1.1173 + xassert(mip != mip);
1.1174 + }
1.1175 + }
1.1176 + return ret;
1.1177 +}
1.1178 +
1.1179 +/***********************************************************************
1.1180 +* NAME
1.1181 +*
1.1182 +* ios_best_node - find active node with best local bound
1.1183 +*
1.1184 +* SYNOPSIS
1.1185 +*
1.1186 +* #include "glpios.h"
1.1187 +* int ios_best_node(glp_tree *tree);
1.1188 +*
1.1189 +* DESCRIPTION
1.1190 +*
1.1191 +* The routine ios_best_node finds an active node whose local bound is
1.1192 +* best among other active nodes.
1.1193 +*
1.1194 +* It is understood that the integer optimal solution of the original
1.1195 +* mip problem cannot be better than the best bound, so the best bound
1.1196 +* is an lower (minimization) or upper (maximization) global bound for
1.1197 +* the original problem.
1.1198 +*
1.1199 +* RETURNS
1.1200 +*
1.1201 +* The routine ios_best_node returns the subproblem reference number
1.1202 +* for the best node. However, if the tree is empty, it returns zero. */
1.1203 +
1.1204 +int ios_best_node(glp_tree *tree)
1.1205 +{ IOSNPD *node, *best = NULL;
1.1206 + switch (tree->mip->dir)
1.1207 + { case GLP_MIN:
1.1208 + /* minimization */
1.1209 + for (node = tree->head; node != NULL; node = node->next)
1.1210 + if (best == NULL || best->bound > node->bound)
1.1211 + best = node;
1.1212 + break;
1.1213 + case GLP_MAX:
1.1214 + /* maximization */
1.1215 + for (node = tree->head; node != NULL; node = node->next)
1.1216 + if (best == NULL || best->bound < node->bound)
1.1217 + best = node;
1.1218 + break;
1.1219 + default:
1.1220 + xassert(tree != tree);
1.1221 + }
1.1222 + return best == NULL ? 0 : best->p;
1.1223 +}
1.1224 +
1.1225 +/***********************************************************************
1.1226 +* NAME
1.1227 +*
1.1228 +* ios_relative_gap - compute relative mip gap
1.1229 +*
1.1230 +* SYNOPSIS
1.1231 +*
1.1232 +* #include "glpios.h"
1.1233 +* double ios_relative_gap(glp_tree *tree);
1.1234 +*
1.1235 +* DESCRIPTION
1.1236 +*
1.1237 +* The routine ios_relative_gap computes the relative mip gap using the
1.1238 +* formula:
1.1239 +*
1.1240 +* gap = |best_mip - best_bnd| / (|best_mip| + DBL_EPSILON),
1.1241 +*
1.1242 +* where best_mip is the best integer feasible solution found so far,
1.1243 +* best_bnd is the best (global) bound. If no integer feasible solution
1.1244 +* has been found yet, rel_gap is set to DBL_MAX.
1.1245 +*
1.1246 +* RETURNS
1.1247 +*
1.1248 +* The routine ios_relative_gap returns the relative mip gap. */
1.1249 +
1.1250 +double ios_relative_gap(glp_tree *tree)
1.1251 +{ glp_prob *mip = tree->mip;
1.1252 + int p;
1.1253 + double best_mip, best_bnd, gap;
1.1254 + if (mip->mip_stat == GLP_FEAS)
1.1255 + { best_mip = mip->mip_obj;
1.1256 + p = ios_best_node(tree);
1.1257 + if (p == 0)
1.1258 + { /* the tree is empty */
1.1259 + gap = 0.0;
1.1260 + }
1.1261 + else
1.1262 + { best_bnd = tree->slot[p].node->bound;
1.1263 + gap = fabs(best_mip - best_bnd) / (fabs(best_mip) +
1.1264 + DBL_EPSILON);
1.1265 + }
1.1266 + }
1.1267 + else
1.1268 + { /* no integer feasible solution has been found yet */
1.1269 + gap = DBL_MAX;
1.1270 + }
1.1271 + return gap;
1.1272 +}
1.1273 +
1.1274 +/***********************************************************************
1.1275 +* NAME
1.1276 +*
1.1277 +* ios_solve_node - solve LP relaxation of current subproblem
1.1278 +*
1.1279 +* SYNOPSIS
1.1280 +*
1.1281 +* #include "glpios.h"
1.1282 +* int ios_solve_node(glp_tree *tree);
1.1283 +*
1.1284 +* DESCRIPTION
1.1285 +*
1.1286 +* The routine ios_solve_node re-optimizes LP relaxation of the current
1.1287 +* subproblem using the dual simplex method.
1.1288 +*
1.1289 +* RETURNS
1.1290 +*
1.1291 +* The routine returns the code which is reported by glp_simplex. */
1.1292 +
1.1293 +int ios_solve_node(glp_tree *tree)
1.1294 +{ glp_prob *mip = tree->mip;
1.1295 + glp_smcp parm;
1.1296 + int ret;
1.1297 + /* the current subproblem must exist */
1.1298 + xassert(tree->curr != NULL);
1.1299 + /* set some control parameters */
1.1300 + glp_init_smcp(&parm);
1.1301 + switch (tree->parm->msg_lev)
1.1302 + { case GLP_MSG_OFF:
1.1303 + parm.msg_lev = GLP_MSG_OFF; break;
1.1304 + case GLP_MSG_ERR:
1.1305 + parm.msg_lev = GLP_MSG_ERR; break;
1.1306 + case GLP_MSG_ON:
1.1307 + case GLP_MSG_ALL:
1.1308 + parm.msg_lev = GLP_MSG_ON; break;
1.1309 + case GLP_MSG_DBG:
1.1310 + parm.msg_lev = GLP_MSG_ALL; break;
1.1311 + default:
1.1312 + xassert(tree != tree);
1.1313 + }
1.1314 + parm.meth = GLP_DUALP;
1.1315 + if (tree->parm->msg_lev < GLP_MSG_DBG)
1.1316 + parm.out_dly = tree->parm->out_dly;
1.1317 + else
1.1318 + parm.out_dly = 0;
1.1319 + /* if the incumbent objective value is already known, use it to
1.1320 + prematurely terminate the dual simplex search */
1.1321 + if (mip->mip_stat == GLP_FEAS)
1.1322 + { switch (tree->mip->dir)
1.1323 + { case GLP_MIN:
1.1324 + parm.obj_ul = mip->mip_obj;
1.1325 + break;
1.1326 + case GLP_MAX:
1.1327 + parm.obj_ll = mip->mip_obj;
1.1328 + break;
1.1329 + default:
1.1330 + xassert(mip != mip);
1.1331 + }
1.1332 + }
1.1333 + /* try to solve/re-optimize the LP relaxation */
1.1334 + ret = glp_simplex(mip, &parm);
1.1335 + tree->curr->solved++;
1.1336 +#if 0
1.1337 + xprintf("ret = %d; status = %d; pbs = %d; dbs = %d; some = %d\n",
1.1338 + ret, glp_get_status(mip), mip->pbs_stat, mip->dbs_stat,
1.1339 + mip->some);
1.1340 + lpx_print_sol(mip, "sol");
1.1341 +#endif
1.1342 + return ret;
1.1343 +}
1.1344 +
1.1345 +/**********************************************************************/
1.1346 +
1.1347 +IOSPOOL *ios_create_pool(glp_tree *tree)
1.1348 +{ /* create cut pool */
1.1349 + IOSPOOL *pool;
1.1350 +#if 0
1.1351 + pool = dmp_get_atom(tree->pool, sizeof(IOSPOOL));
1.1352 +#else
1.1353 + xassert(tree == tree);
1.1354 + pool = xmalloc(sizeof(IOSPOOL));
1.1355 +#endif
1.1356 + pool->size = 0;
1.1357 + pool->head = pool->tail = NULL;
1.1358 + pool->ord = 0, pool->curr = NULL;
1.1359 + return pool;
1.1360 +}
1.1361 +
1.1362 +int ios_add_row(glp_tree *tree, IOSPOOL *pool,
1.1363 + const char *name, int klass, int flags, int len, const int ind[],
1.1364 + const double val[], int type, double rhs)
1.1365 +{ /* add row (constraint) to the cut pool */
1.1366 + IOSCUT *cut;
1.1367 + IOSAIJ *aij;
1.1368 + int k;
1.1369 + xassert(pool != NULL);
1.1370 + cut = dmp_get_atom(tree->pool, sizeof(IOSCUT));
1.1371 + if (name == NULL || name[0] == '\0')
1.1372 + cut->name = NULL;
1.1373 + else
1.1374 + { for (k = 0; name[k] != '\0'; k++)
1.1375 + { if (k == 256)
1.1376 + xerror("glp_ios_add_row: cut name too long\n");
1.1377 + if (iscntrl((unsigned char)name[k]))
1.1378 + xerror("glp_ios_add_row: cut name contains invalid chara"
1.1379 + "cter(s)\n");
1.1380 + }
1.1381 + cut->name = dmp_get_atom(tree->pool, strlen(name)+1);
1.1382 + strcpy(cut->name, name);
1.1383 + }
1.1384 + if (!(0 <= klass && klass <= 255))
1.1385 + xerror("glp_ios_add_row: klass = %d; invalid cut class\n",
1.1386 + klass);
1.1387 + cut->klass = (unsigned char)klass;
1.1388 + if (flags != 0)
1.1389 + xerror("glp_ios_add_row: flags = %d; invalid cut flags\n",
1.1390 + flags);
1.1391 + cut->ptr = NULL;
1.1392 + if (!(0 <= len && len <= tree->n))
1.1393 + xerror("glp_ios_add_row: len = %d; invalid cut length\n",
1.1394 + len);
1.1395 + for (k = 1; k <= len; k++)
1.1396 + { aij = dmp_get_atom(tree->pool, sizeof(IOSAIJ));
1.1397 + if (!(1 <= ind[k] && ind[k] <= tree->n))
1.1398 + xerror("glp_ios_add_row: ind[%d] = %d; column index out of "
1.1399 + "range\n", k, ind[k]);
1.1400 + aij->j = ind[k];
1.1401 + aij->val = val[k];
1.1402 + aij->next = cut->ptr;
1.1403 + cut->ptr = aij;
1.1404 + }
1.1405 + if (!(type == GLP_LO || type == GLP_UP || type == GLP_FX))
1.1406 + xerror("glp_ios_add_row: type = %d; invalid cut type\n",
1.1407 + type);
1.1408 + cut->type = (unsigned char)type;
1.1409 + cut->rhs = rhs;
1.1410 + cut->prev = pool->tail;
1.1411 + cut->next = NULL;
1.1412 + if (cut->prev == NULL)
1.1413 + pool->head = cut;
1.1414 + else
1.1415 + cut->prev->next = cut;
1.1416 + pool->tail = cut;
1.1417 + pool->size++;
1.1418 + return pool->size;
1.1419 +}
1.1420 +
1.1421 +IOSCUT *ios_find_row(IOSPOOL *pool, int i)
1.1422 +{ /* find row (constraint) in the cut pool */
1.1423 + /* (smart linear search) */
1.1424 + xassert(pool != NULL);
1.1425 + xassert(1 <= i && i <= pool->size);
1.1426 + if (pool->ord == 0)
1.1427 + { xassert(pool->curr == NULL);
1.1428 + pool->ord = 1;
1.1429 + pool->curr = pool->head;
1.1430 + }
1.1431 + xassert(pool->curr != NULL);
1.1432 + if (i < pool->ord)
1.1433 + { if (i < pool->ord - i)
1.1434 + { pool->ord = 1;
1.1435 + pool->curr = pool->head;
1.1436 + while (pool->ord != i)
1.1437 + { pool->ord++;
1.1438 + xassert(pool->curr != NULL);
1.1439 + pool->curr = pool->curr->next;
1.1440 + }
1.1441 + }
1.1442 + else
1.1443 + { while (pool->ord != i)
1.1444 + { pool->ord--;
1.1445 + xassert(pool->curr != NULL);
1.1446 + pool->curr = pool->curr->prev;
1.1447 + }
1.1448 + }
1.1449 + }
1.1450 + else if (i > pool->ord)
1.1451 + { if (i - pool->ord < pool->size - i)
1.1452 + { while (pool->ord != i)
1.1453 + { pool->ord++;
1.1454 + xassert(pool->curr != NULL);
1.1455 + pool->curr = pool->curr->next;
1.1456 + }
1.1457 + }
1.1458 + else
1.1459 + { pool->ord = pool->size;
1.1460 + pool->curr = pool->tail;
1.1461 + while (pool->ord != i)
1.1462 + { pool->ord--;
1.1463 + xassert(pool->curr != NULL);
1.1464 + pool->curr = pool->curr->prev;
1.1465 + }
1.1466 + }
1.1467 + }
1.1468 + xassert(pool->ord == i);
1.1469 + xassert(pool->curr != NULL);
1.1470 + return pool->curr;
1.1471 +}
1.1472 +
1.1473 +void ios_del_row(glp_tree *tree, IOSPOOL *pool, int i)
1.1474 +{ /* remove row (constraint) from the cut pool */
1.1475 + IOSCUT *cut;
1.1476 + IOSAIJ *aij;
1.1477 + xassert(pool != NULL);
1.1478 + if (!(1 <= i && i <= pool->size))
1.1479 + xerror("glp_ios_del_row: i = %d; cut number out of range\n",
1.1480 + i);
1.1481 + cut = ios_find_row(pool, i);
1.1482 + xassert(pool->curr == cut);
1.1483 + if (cut->next != NULL)
1.1484 + pool->curr = cut->next;
1.1485 + else if (cut->prev != NULL)
1.1486 + pool->ord--, pool->curr = cut->prev;
1.1487 + else
1.1488 + pool->ord = 0, pool->curr = NULL;
1.1489 + if (cut->name != NULL)
1.1490 + dmp_free_atom(tree->pool, cut->name, strlen(cut->name)+1);
1.1491 + if (cut->prev == NULL)
1.1492 + { xassert(pool->head == cut);
1.1493 + pool->head = cut->next;
1.1494 + }
1.1495 + else
1.1496 + { xassert(cut->prev->next == cut);
1.1497 + cut->prev->next = cut->next;
1.1498 + }
1.1499 + if (cut->next == NULL)
1.1500 + { xassert(pool->tail == cut);
1.1501 + pool->tail = cut->prev;
1.1502 + }
1.1503 + else
1.1504 + { xassert(cut->next->prev == cut);
1.1505 + cut->next->prev = cut->prev;
1.1506 + }
1.1507 + while (cut->ptr != NULL)
1.1508 + { aij = cut->ptr;
1.1509 + cut->ptr = aij->next;
1.1510 + dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ));
1.1511 + }
1.1512 + dmp_free_atom(tree->pool, cut, sizeof(IOSCUT));
1.1513 + pool->size--;
1.1514 + return;
1.1515 +}
1.1516 +
1.1517 +void ios_clear_pool(glp_tree *tree, IOSPOOL *pool)
1.1518 +{ /* remove all rows (constraints) from the cut pool */
1.1519 + xassert(pool != NULL);
1.1520 + while (pool->head != NULL)
1.1521 + { IOSCUT *cut = pool->head;
1.1522 + pool->head = cut->next;
1.1523 + if (cut->name != NULL)
1.1524 + dmp_free_atom(tree->pool, cut->name, strlen(cut->name)+1);
1.1525 + while (cut->ptr != NULL)
1.1526 + { IOSAIJ *aij = cut->ptr;
1.1527 + cut->ptr = aij->next;
1.1528 + dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ));
1.1529 + }
1.1530 + dmp_free_atom(tree->pool, cut, sizeof(IOSCUT));
1.1531 + }
1.1532 + pool->size = 0;
1.1533 + pool->head = pool->tail = NULL;
1.1534 + pool->ord = 0, pool->curr = NULL;
1.1535 + return;
1.1536 +}
1.1537 +
1.1538 +void ios_delete_pool(glp_tree *tree, IOSPOOL *pool)
1.1539 +{ /* delete cut pool */
1.1540 + xassert(pool != NULL);
1.1541 + ios_clear_pool(tree, pool);
1.1542 + xfree(pool);
1.1543 + return;
1.1544 +}
1.1545 +
1.1546 +/**********************************************************************/
1.1547 +
1.1548 +#if 0
1.1549 +static int refer_to_node(glp_tree *tree, int j)
1.1550 +{ /* determine node number corresponding to binary variable x[j] or
1.1551 + its complement */
1.1552 + glp_prob *mip = tree->mip;
1.1553 + int n = mip->n;
1.1554 + int *ref;
1.1555 + if (j > 0)
1.1556 + ref = tree->n_ref;
1.1557 + else
1.1558 + ref = tree->c_ref, j = - j;
1.1559 + xassert(1 <= j && j <= n);
1.1560 + if (ref[j] == 0)
1.1561 + { /* new node is needed */
1.1562 + SCG *g = tree->g;
1.1563 + int n_max = g->n_max;
1.1564 + ref[j] = scg_add_nodes(g, 1);
1.1565 + if (g->n_max > n_max)
1.1566 + { int *save = tree->j_ref;
1.1567 + tree->j_ref = xcalloc(1+g->n_max, sizeof(int));
1.1568 + memcpy(&tree->j_ref[1], &save[1], g->n * sizeof(int));
1.1569 + xfree(save);
1.1570 + }
1.1571 + xassert(ref[j] == g->n);
1.1572 + tree->j_ref[ref[j]] = j;
1.1573 + xassert(tree->curr != NULL);
1.1574 + if (tree->curr->level > 0) tree->curr->own_nn++;
1.1575 + }
1.1576 + return ref[j];
1.1577 +}
1.1578 +#endif
1.1579 +
1.1580 +#if 0
1.1581 +void ios_add_edge(glp_tree *tree, int j1, int j2)
1.1582 +{ /* add new edge to the conflict graph */
1.1583 + glp_prob *mip = tree->mip;
1.1584 + int n = mip->n;
1.1585 + SCGRIB *e;
1.1586 + int first, i1, i2;
1.1587 + xassert(-n <= j1 && j1 <= +n && j1 != 0);
1.1588 + xassert(-n <= j2 && j2 <= +n && j2 != 0);
1.1589 + xassert(j1 != j2);
1.1590 + /* determine number of the first node, which was added for the
1.1591 + current subproblem */
1.1592 + xassert(tree->curr != NULL);
1.1593 + first = tree->g->n - tree->curr->own_nn + 1;
1.1594 + /* determine node numbers for both endpoints */
1.1595 + i1 = refer_to_node(tree, j1);
1.1596 + i2 = refer_to_node(tree, j2);
1.1597 + /* add edge (i1,i2) to the conflict graph */
1.1598 + e = scg_add_edge(tree->g, i1, i2);
1.1599 + /* if the current subproblem is not the root and both endpoints
1.1600 + were created on some previous levels, save the edge */
1.1601 + if (tree->curr->level > 0 && i1 < first && i2 < first)
1.1602 + { IOSRIB *rib;
1.1603 + rib = dmp_get_atom(tree->pool, sizeof(IOSRIB));
1.1604 + rib->j1 = j1;
1.1605 + rib->j2 = j2;
1.1606 + rib->e = e;
1.1607 + rib->next = tree->curr->e_ptr;
1.1608 + tree->curr->e_ptr = rib;
1.1609 + }
1.1610 + return;
1.1611 +}
1.1612 +#endif
1.1613 +
1.1614 +/* eof */