lemon-project-template-glpk
diff deps/glpk/src/glpios01.c @ 9:33de93886c88
Import GLPK 4.47
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 20:59:10 +0100 |
| parents | |
| children |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/deps/glpk/src/glpios01.c Sun Nov 06 20:59:10 2011 +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, 2011 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 */