alpar@9: /* glpios01.c */ alpar@9: alpar@9: /*********************************************************************** alpar@9: * This code is part of GLPK (GNU Linear Programming Kit). alpar@9: * alpar@9: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@9: * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, alpar@9: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@9: * E-mail: . alpar@9: * alpar@9: * GLPK is free software: you can redistribute it and/or modify it alpar@9: * under the terms of the GNU General Public License as published by alpar@9: * the Free Software Foundation, either version 3 of the License, or alpar@9: * (at your option) any later version. alpar@9: * alpar@9: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@9: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@9: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@9: * License for more details. alpar@9: * alpar@9: * You should have received a copy of the GNU General Public License alpar@9: * along with GLPK. If not, see . alpar@9: ***********************************************************************/ alpar@9: alpar@9: #include "glpios.h" alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_create_tree - create branch-and-bound tree alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine ios_create_tree creates the branch-and-bound tree. alpar@9: * alpar@9: * Being created the tree consists of the only root subproblem whose alpar@9: * reference number is 1. Note that initially the root subproblem is in alpar@9: * frozen state and therefore needs to be revived. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine returns a pointer to the tree created. */ alpar@9: alpar@9: static IOSNPD *new_node(glp_tree *tree, IOSNPD *parent); alpar@9: alpar@9: glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm) alpar@9: { int m = mip->m; alpar@9: int n = mip->n; alpar@9: glp_tree *tree; alpar@9: int i, j; alpar@9: xassert(mip->tree == NULL); alpar@9: mip->tree = tree = xmalloc(sizeof(glp_tree)); alpar@9: tree->pool = dmp_create_pool(); alpar@9: tree->n = n; alpar@9: /* save original problem components */ alpar@9: tree->orig_m = m; alpar@9: tree->orig_type = xcalloc(1+m+n, sizeof(char)); alpar@9: tree->orig_lb = xcalloc(1+m+n, sizeof(double)); alpar@9: tree->orig_ub = xcalloc(1+m+n, sizeof(double)); alpar@9: tree->orig_stat = xcalloc(1+m+n, sizeof(char)); alpar@9: tree->orig_prim = xcalloc(1+m+n, sizeof(double)); alpar@9: tree->orig_dual = xcalloc(1+m+n, sizeof(double)); alpar@9: for (i = 1; i <= m; i++) alpar@9: { GLPROW *row = mip->row[i]; alpar@9: tree->orig_type[i] = (char)row->type; alpar@9: tree->orig_lb[i] = row->lb; alpar@9: tree->orig_ub[i] = row->ub; alpar@9: tree->orig_stat[i] = (char)row->stat; alpar@9: tree->orig_prim[i] = row->prim; alpar@9: tree->orig_dual[i] = row->dual; alpar@9: } alpar@9: for (j = 1; j <= n; j++) alpar@9: { GLPCOL *col = mip->col[j]; alpar@9: tree->orig_type[m+j] = (char)col->type; alpar@9: tree->orig_lb[m+j] = col->lb; alpar@9: tree->orig_ub[m+j] = col->ub; alpar@9: tree->orig_stat[m+j] = (char)col->stat; alpar@9: tree->orig_prim[m+j] = col->prim; alpar@9: tree->orig_dual[m+j] = col->dual; alpar@9: } alpar@9: tree->orig_obj = mip->obj_val; alpar@9: /* initialize the branch-and-bound tree */ alpar@9: tree->nslots = 0; alpar@9: tree->avail = 0; alpar@9: tree->slot = NULL; alpar@9: tree->head = tree->tail = NULL; alpar@9: tree->a_cnt = tree->n_cnt = tree->t_cnt = 0; alpar@9: /* the root subproblem is not solved yet, so its final components alpar@9: are unknown so far */ alpar@9: tree->root_m = 0; alpar@9: tree->root_type = NULL; alpar@9: tree->root_lb = tree->root_ub = NULL; alpar@9: tree->root_stat = NULL; alpar@9: /* the current subproblem does not exist yet */ alpar@9: tree->curr = NULL; alpar@9: tree->mip = mip; alpar@9: /*tree->solved = 0;*/ alpar@9: tree->non_int = xcalloc(1+n, sizeof(char)); alpar@9: memset(&tree->non_int[1], 0, n); alpar@9: /* arrays to save parent subproblem components will be allocated alpar@9: later */ alpar@9: tree->pred_m = tree->pred_max = 0; alpar@9: tree->pred_type = NULL; alpar@9: tree->pred_lb = tree->pred_ub = NULL; alpar@9: tree->pred_stat = NULL; alpar@9: /* cut generator */ alpar@9: tree->local = ios_create_pool(tree); alpar@9: /*tree->first_attempt = 1;*/ alpar@9: /*tree->max_added_cuts = 0;*/ alpar@9: /*tree->min_eff = 0.0;*/ alpar@9: /*tree->miss = 0;*/ alpar@9: /*tree->just_selected = 0;*/ alpar@9: tree->mir_gen = NULL; alpar@9: tree->clq_gen = NULL; alpar@9: /*tree->round = 0;*/ alpar@9: #if 0 alpar@9: /* create the conflict graph */ alpar@9: tree->n_ref = xcalloc(1+n, sizeof(int)); alpar@9: memset(&tree->n_ref[1], 0, n * sizeof(int)); alpar@9: tree->c_ref = xcalloc(1+n, sizeof(int)); alpar@9: memset(&tree->c_ref[1], 0, n * sizeof(int)); alpar@9: tree->g = scg_create_graph(0); alpar@9: tree->j_ref = xcalloc(1+tree->g->n_max, sizeof(int)); alpar@9: #endif alpar@9: /* pseudocost branching */ alpar@9: tree->pcost = NULL; alpar@9: tree->iwrk = xcalloc(1+n, sizeof(int)); alpar@9: tree->dwrk = xcalloc(1+n, sizeof(double)); alpar@9: /* initialize control parameters */ alpar@9: tree->parm = parm; alpar@9: tree->tm_beg = xtime(); alpar@9: tree->tm_lag = xlset(0); alpar@9: tree->sol_cnt = 0; alpar@9: /* initialize advanced solver interface */ alpar@9: tree->reason = 0; alpar@9: tree->reopt = 0; alpar@9: tree->reinv = 0; alpar@9: tree->br_var = 0; alpar@9: tree->br_sel = 0; alpar@9: tree->child = 0; alpar@9: tree->next_p = 0; alpar@9: /*tree->btrack = NULL;*/ alpar@9: tree->stop = 0; alpar@9: /* create the root subproblem, which initially is identical to alpar@9: the original MIP */ alpar@9: new_node(tree, NULL); alpar@9: return tree; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_revive_node - revive specified subproblem alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * void ios_revive_node(glp_tree *tree, int p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine ios_revive_node revives the specified subproblem, whose alpar@9: * reference number is p, and thereby makes it the current subproblem. alpar@9: * Note that the specified subproblem must be active. Besides, if the alpar@9: * current subproblem already exists, it must be frozen before reviving alpar@9: * another subproblem. */ alpar@9: alpar@9: void ios_revive_node(glp_tree *tree, int p) alpar@9: { glp_prob *mip = tree->mip; alpar@9: IOSNPD *node, *root; alpar@9: /* obtain pointer to the specified subproblem */ alpar@9: xassert(1 <= p && p <= tree->nslots); alpar@9: node = tree->slot[p].node; alpar@9: xassert(node != NULL); alpar@9: /* the specified subproblem must be active */ alpar@9: xassert(node->count == 0); alpar@9: /* the current subproblem must not exist */ alpar@9: xassert(tree->curr == NULL); alpar@9: /* the specified subproblem becomes current */ alpar@9: tree->curr = node; alpar@9: /*tree->solved = 0;*/ alpar@9: /* obtain pointer to the root subproblem */ alpar@9: root = tree->slot[1].node; alpar@9: xassert(root != NULL); alpar@9: /* at this point problem object components correspond to the root alpar@9: subproblem, so if the root subproblem should be revived, there alpar@9: is nothing more to do */ alpar@9: if (node == root) goto done; alpar@9: xassert(mip->m == tree->root_m); alpar@9: /* build path from the root to the current node */ alpar@9: node->temp = NULL; alpar@9: for (node = node; node != NULL; node = node->up) alpar@9: { if (node->up == NULL) alpar@9: xassert(node == root); alpar@9: else alpar@9: node->up->temp = node; alpar@9: } alpar@9: /* go down from the root to the current node and make necessary alpar@9: changes to restore components of the current subproblem */ alpar@9: for (node = root; node != NULL; node = node->temp) alpar@9: { int m = mip->m; alpar@9: int n = mip->n; alpar@9: /* if the current node is reached, the problem object at this alpar@9: point corresponds to its parent, so save attributes of rows alpar@9: and columns for the parent subproblem */ alpar@9: if (node->temp == NULL) alpar@9: { int i, j; alpar@9: tree->pred_m = m; alpar@9: /* allocate/reallocate arrays, if necessary */ alpar@9: if (tree->pred_max < m + n) alpar@9: { int new_size = m + n + 100; alpar@9: if (tree->pred_type != NULL) xfree(tree->pred_type); alpar@9: if (tree->pred_lb != NULL) xfree(tree->pred_lb); alpar@9: if (tree->pred_ub != NULL) xfree(tree->pred_ub); alpar@9: if (tree->pred_stat != NULL) xfree(tree->pred_stat); alpar@9: tree->pred_max = new_size; alpar@9: tree->pred_type = xcalloc(1+new_size, sizeof(char)); alpar@9: tree->pred_lb = xcalloc(1+new_size, sizeof(double)); alpar@9: tree->pred_ub = xcalloc(1+new_size, sizeof(double)); alpar@9: tree->pred_stat = xcalloc(1+new_size, sizeof(char)); alpar@9: } alpar@9: /* save row attributes */ alpar@9: for (i = 1; i <= m; i++) alpar@9: { GLPROW *row = mip->row[i]; alpar@9: tree->pred_type[i] = (char)row->type; alpar@9: tree->pred_lb[i] = row->lb; alpar@9: tree->pred_ub[i] = row->ub; alpar@9: tree->pred_stat[i] = (char)row->stat; alpar@9: } alpar@9: /* save column attributes */ alpar@9: for (j = 1; j <= n; j++) alpar@9: { GLPCOL *col = mip->col[j]; alpar@9: tree->pred_type[mip->m+j] = (char)col->type; alpar@9: tree->pred_lb[mip->m+j] = col->lb; alpar@9: tree->pred_ub[mip->m+j] = col->ub; alpar@9: tree->pred_stat[mip->m+j] = (char)col->stat; alpar@9: } alpar@9: } alpar@9: /* change bounds of rows and columns */ alpar@9: { IOSBND *b; alpar@9: for (b = node->b_ptr; b != NULL; b = b->next) alpar@9: { if (b->k <= m) alpar@9: glp_set_row_bnds(mip, b->k, b->type, b->lb, b->ub); alpar@9: else alpar@9: glp_set_col_bnds(mip, b->k-m, b->type, b->lb, b->ub); alpar@9: } alpar@9: } alpar@9: /* change statuses of rows and columns */ alpar@9: { IOSTAT *s; alpar@9: for (s = node->s_ptr; s != NULL; s = s->next) alpar@9: { if (s->k <= m) alpar@9: glp_set_row_stat(mip, s->k, s->stat); alpar@9: else alpar@9: glp_set_col_stat(mip, s->k-m, s->stat); alpar@9: } alpar@9: } alpar@9: /* add new rows */ alpar@9: if (node->r_ptr != NULL) alpar@9: { IOSROW *r; alpar@9: IOSAIJ *a; alpar@9: int i, len, *ind; alpar@9: double *val; alpar@9: ind = xcalloc(1+n, sizeof(int)); alpar@9: val = xcalloc(1+n, sizeof(double)); alpar@9: for (r = node->r_ptr; r != NULL; r = r->next) alpar@9: { i = glp_add_rows(mip, 1); alpar@9: glp_set_row_name(mip, i, r->name); alpar@9: #if 1 /* 20/IX-2008 */ alpar@9: xassert(mip->row[i]->level == 0); alpar@9: mip->row[i]->level = node->level; alpar@9: mip->row[i]->origin = r->origin; alpar@9: mip->row[i]->klass = r->klass; alpar@9: #endif alpar@9: glp_set_row_bnds(mip, i, r->type, r->lb, r->ub); alpar@9: len = 0; alpar@9: for (a = r->ptr; a != NULL; a = a->next) alpar@9: len++, ind[len] = a->j, val[len] = a->val; alpar@9: glp_set_mat_row(mip, i, len, ind, val); alpar@9: glp_set_rii(mip, i, r->rii); alpar@9: glp_set_row_stat(mip, i, r->stat); alpar@9: } alpar@9: xfree(ind); alpar@9: xfree(val); alpar@9: } alpar@9: #if 0 alpar@9: /* add new edges to the conflict graph */ alpar@9: /* add new cliques to the conflict graph */ alpar@9: /* (not implemented yet) */ alpar@9: xassert(node->own_nn == 0); alpar@9: xassert(node->own_nc == 0); alpar@9: xassert(node->e_ptr == NULL); alpar@9: #endif alpar@9: } alpar@9: /* the specified subproblem has been revived */ alpar@9: node = tree->curr; alpar@9: /* delete its bound change list */ alpar@9: while (node->b_ptr != NULL) alpar@9: { IOSBND *b; alpar@9: b = node->b_ptr; alpar@9: node->b_ptr = b->next; alpar@9: dmp_free_atom(tree->pool, b, sizeof(IOSBND)); alpar@9: } alpar@9: /* delete its status change list */ alpar@9: while (node->s_ptr != NULL) alpar@9: { IOSTAT *s; alpar@9: s = node->s_ptr; alpar@9: node->s_ptr = s->next; alpar@9: dmp_free_atom(tree->pool, s, sizeof(IOSTAT)); alpar@9: } alpar@9: #if 1 /* 20/XI-2009 */ alpar@9: /* delete its row addition list (additional rows may appear, for alpar@9: example, due to branching on GUB constraints */ alpar@9: while (node->r_ptr != NULL) alpar@9: { IOSROW *r; alpar@9: r = node->r_ptr; alpar@9: node->r_ptr = r->next; alpar@9: xassert(r->name == NULL); alpar@9: while (r->ptr != NULL) alpar@9: { IOSAIJ *a; alpar@9: a = r->ptr; alpar@9: r->ptr = a->next; alpar@9: dmp_free_atom(tree->pool, a, sizeof(IOSAIJ)); alpar@9: } alpar@9: dmp_free_atom(tree->pool, r, sizeof(IOSROW)); alpar@9: } alpar@9: #endif alpar@9: done: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_freeze_node - freeze current subproblem alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * void ios_freeze_node(glp_tree *tree); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine ios_freeze_node freezes the current subproblem. */ alpar@9: alpar@9: void ios_freeze_node(glp_tree *tree) alpar@9: { glp_prob *mip = tree->mip; alpar@9: int m = mip->m; alpar@9: int n = mip->n; alpar@9: IOSNPD *node; alpar@9: /* obtain pointer to the current subproblem */ alpar@9: node = tree->curr; alpar@9: xassert(node != NULL); alpar@9: if (node->up == NULL) alpar@9: { /* freeze the root subproblem */ alpar@9: int k; alpar@9: xassert(node->p == 1); alpar@9: xassert(tree->root_m == 0); alpar@9: xassert(tree->root_type == NULL); alpar@9: xassert(tree->root_lb == NULL); alpar@9: xassert(tree->root_ub == NULL); alpar@9: xassert(tree->root_stat == NULL); alpar@9: tree->root_m = m; alpar@9: tree->root_type = xcalloc(1+m+n, sizeof(char)); alpar@9: tree->root_lb = xcalloc(1+m+n, sizeof(double)); alpar@9: tree->root_ub = xcalloc(1+m+n, sizeof(double)); alpar@9: tree->root_stat = xcalloc(1+m+n, sizeof(char)); alpar@9: for (k = 1; k <= m+n; k++) alpar@9: { if (k <= m) alpar@9: { GLPROW *row = mip->row[k]; alpar@9: tree->root_type[k] = (char)row->type; alpar@9: tree->root_lb[k] = row->lb; alpar@9: tree->root_ub[k] = row->ub; alpar@9: tree->root_stat[k] = (char)row->stat; alpar@9: } alpar@9: else alpar@9: { GLPCOL *col = mip->col[k-m]; alpar@9: tree->root_type[k] = (char)col->type; alpar@9: tree->root_lb[k] = col->lb; alpar@9: tree->root_ub[k] = col->ub; alpar@9: tree->root_stat[k] = (char)col->stat; alpar@9: } alpar@9: } alpar@9: } alpar@9: else alpar@9: { /* freeze non-root subproblem */ alpar@9: int root_m = tree->root_m; alpar@9: int pred_m = tree->pred_m; alpar@9: int i, j, k; alpar@9: xassert(pred_m <= m); alpar@9: /* build change lists for rows and columns which exist in the alpar@9: parent subproblem */ alpar@9: xassert(node->b_ptr == NULL); alpar@9: xassert(node->s_ptr == NULL); alpar@9: for (k = 1; k <= pred_m + n; k++) alpar@9: { int pred_type, pred_stat, type, stat; alpar@9: double pred_lb, pred_ub, lb, ub; alpar@9: /* determine attributes in the parent subproblem */ alpar@9: pred_type = tree->pred_type[k]; alpar@9: pred_lb = tree->pred_lb[k]; alpar@9: pred_ub = tree->pred_ub[k]; alpar@9: pred_stat = tree->pred_stat[k]; alpar@9: /* determine attributes in the current subproblem */ alpar@9: if (k <= pred_m) alpar@9: { GLPROW *row = mip->row[k]; alpar@9: type = row->type; alpar@9: lb = row->lb; alpar@9: ub = row->ub; alpar@9: stat = row->stat; alpar@9: } alpar@9: else alpar@9: { GLPCOL *col = mip->col[k - pred_m]; alpar@9: type = col->type; alpar@9: lb = col->lb; alpar@9: ub = col->ub; alpar@9: stat = col->stat; alpar@9: } alpar@9: /* save type and bounds of a row/column, if changed */ alpar@9: if (!(pred_type == type && pred_lb == lb && pred_ub == ub)) alpar@9: { IOSBND *b; alpar@9: b = dmp_get_atom(tree->pool, sizeof(IOSBND)); alpar@9: b->k = k; alpar@9: b->type = (unsigned char)type; alpar@9: b->lb = lb; alpar@9: b->ub = ub; alpar@9: b->next = node->b_ptr; alpar@9: node->b_ptr = b; alpar@9: } alpar@9: /* save status of a row/column, if changed */ alpar@9: if (pred_stat != stat) alpar@9: { IOSTAT *s; alpar@9: s = dmp_get_atom(tree->pool, sizeof(IOSTAT)); alpar@9: s->k = k; alpar@9: s->stat = (unsigned char)stat; alpar@9: s->next = node->s_ptr; alpar@9: node->s_ptr = s; alpar@9: } alpar@9: } alpar@9: /* save new rows added to the current subproblem */ alpar@9: xassert(node->r_ptr == NULL); alpar@9: if (pred_m < m) alpar@9: { int i, len, *ind; alpar@9: double *val; alpar@9: ind = xcalloc(1+n, sizeof(int)); alpar@9: val = xcalloc(1+n, sizeof(double)); alpar@9: for (i = m; i > pred_m; i--) alpar@9: { GLPROW *row = mip->row[i]; alpar@9: IOSROW *r; alpar@9: const char *name; alpar@9: r = dmp_get_atom(tree->pool, sizeof(IOSROW)); alpar@9: name = glp_get_row_name(mip, i); alpar@9: if (name == NULL) alpar@9: r->name = NULL; alpar@9: else alpar@9: { r->name = dmp_get_atom(tree->pool, strlen(name)+1); alpar@9: strcpy(r->name, name); alpar@9: } alpar@9: #if 1 /* 20/IX-2008 */ alpar@9: r->origin = row->origin; alpar@9: r->klass = row->klass; alpar@9: #endif alpar@9: r->type = (unsigned char)row->type; alpar@9: r->lb = row->lb; alpar@9: r->ub = row->ub; alpar@9: r->ptr = NULL; alpar@9: len = glp_get_mat_row(mip, i, ind, val); alpar@9: for (k = 1; k <= len; k++) alpar@9: { IOSAIJ *a; alpar@9: a = dmp_get_atom(tree->pool, sizeof(IOSAIJ)); alpar@9: a->j = ind[k]; alpar@9: a->val = val[k]; alpar@9: a->next = r->ptr; alpar@9: r->ptr = a; alpar@9: } alpar@9: r->rii = row->rii; alpar@9: r->stat = (unsigned char)row->stat; alpar@9: r->next = node->r_ptr; alpar@9: node->r_ptr = r; alpar@9: } alpar@9: xfree(ind); alpar@9: xfree(val); alpar@9: } alpar@9: /* remove all rows missing in the root subproblem */ alpar@9: if (m != root_m) alpar@9: { int nrs, *num; alpar@9: nrs = m - root_m; alpar@9: xassert(nrs > 0); alpar@9: num = xcalloc(1+nrs, sizeof(int)); alpar@9: for (i = 1; i <= nrs; i++) num[i] = root_m + i; alpar@9: glp_del_rows(mip, nrs, num); alpar@9: xfree(num); alpar@9: } alpar@9: m = mip->m; alpar@9: /* and restore attributes of all rows and columns for the root alpar@9: subproblem */ alpar@9: xassert(m == root_m); alpar@9: for (i = 1; i <= m; i++) alpar@9: { glp_set_row_bnds(mip, i, tree->root_type[i], alpar@9: tree->root_lb[i], tree->root_ub[i]); alpar@9: glp_set_row_stat(mip, i, tree->root_stat[i]); alpar@9: } alpar@9: for (j = 1; j <= n; j++) alpar@9: { glp_set_col_bnds(mip, j, tree->root_type[m+j], alpar@9: tree->root_lb[m+j], tree->root_ub[m+j]); alpar@9: glp_set_col_stat(mip, j, tree->root_stat[m+j]); alpar@9: } alpar@9: #if 1 alpar@9: /* remove all edges and cliques missing in the conflict graph alpar@9: for the root subproblem */ alpar@9: /* (not implemented yet) */ alpar@9: #endif alpar@9: } alpar@9: /* the current subproblem has been frozen */ alpar@9: tree->curr = NULL; alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_clone_node - clone specified subproblem alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * void ios_clone_node(glp_tree *tree, int p, int nnn, int ref[]); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine ios_clone_node clones the specified subproblem, whose alpar@9: * reference number is p, creating its nnn exact copies. Note that the alpar@9: * specified subproblem must be active and must be in the frozen state alpar@9: * (i.e. it must not be the current subproblem). alpar@9: * alpar@9: * Each clone, an exact copy of the specified subproblem, becomes a new alpar@9: * active subproblem added to the end of the active list. After cloning alpar@9: * the specified subproblem becomes inactive. alpar@9: * alpar@9: * The reference numbers of clone subproblems are stored to locations alpar@9: * ref[1], ..., ref[nnn]. */ alpar@9: alpar@9: static int get_slot(glp_tree *tree) alpar@9: { int p; alpar@9: /* if no free slots are available, increase the room */ alpar@9: if (tree->avail == 0) alpar@9: { int nslots = tree->nslots; alpar@9: IOSLOT *save = tree->slot; alpar@9: if (nslots == 0) alpar@9: tree->nslots = 20; alpar@9: else alpar@9: { tree->nslots = nslots + nslots; alpar@9: xassert(tree->nslots > nslots); alpar@9: } alpar@9: tree->slot = xcalloc(1+tree->nslots, sizeof(IOSLOT)); alpar@9: if (save != NULL) alpar@9: { memcpy(&tree->slot[1], &save[1], nslots * sizeof(IOSLOT)); alpar@9: xfree(save); alpar@9: } alpar@9: /* push more free slots into the stack */ alpar@9: for (p = tree->nslots; p > nslots; p--) alpar@9: { tree->slot[p].node = NULL; alpar@9: tree->slot[p].next = tree->avail; alpar@9: tree->avail = p; alpar@9: } alpar@9: } alpar@9: /* pull a free slot from the stack */ alpar@9: p = tree->avail; alpar@9: tree->avail = tree->slot[p].next; alpar@9: xassert(tree->slot[p].node == NULL); alpar@9: tree->slot[p].next = 0; alpar@9: return p; alpar@9: } alpar@9: alpar@9: static IOSNPD *new_node(glp_tree *tree, IOSNPD *parent) alpar@9: { IOSNPD *node; alpar@9: int p; alpar@9: /* pull a free slot for the new node */ alpar@9: p = get_slot(tree); alpar@9: /* create descriptor of the new subproblem */ alpar@9: node = dmp_get_atom(tree->pool, sizeof(IOSNPD)); alpar@9: tree->slot[p].node = node; alpar@9: node->p = p; alpar@9: node->up = parent; alpar@9: node->level = (parent == NULL ? 0 : parent->level + 1); alpar@9: node->count = 0; alpar@9: node->b_ptr = NULL; alpar@9: node->s_ptr = NULL; alpar@9: node->r_ptr = NULL; alpar@9: node->solved = 0; alpar@9: #if 0 alpar@9: node->own_nn = node->own_nc = 0; alpar@9: node->e_ptr = NULL; alpar@9: #endif alpar@9: #if 1 /* 04/X-2008 */ alpar@9: node->lp_obj = (parent == NULL ? (tree->mip->dir == GLP_MIN ? alpar@9: -DBL_MAX : +DBL_MAX) : parent->lp_obj); alpar@9: #endif alpar@9: node->bound = (parent == NULL ? (tree->mip->dir == GLP_MIN ? alpar@9: -DBL_MAX : +DBL_MAX) : parent->bound); alpar@9: node->br_var = 0; alpar@9: node->br_val = 0.0; alpar@9: node->ii_cnt = 0; alpar@9: node->ii_sum = 0.0; alpar@9: #if 1 /* 30/XI-2009 */ alpar@9: node->changed = 0; alpar@9: #endif alpar@9: if (tree->parm->cb_size == 0) alpar@9: node->data = NULL; alpar@9: else alpar@9: { node->data = dmp_get_atom(tree->pool, tree->parm->cb_size); alpar@9: memset(node->data, 0, tree->parm->cb_size); alpar@9: } alpar@9: node->temp = NULL; alpar@9: node->prev = tree->tail; alpar@9: node->next = NULL; alpar@9: /* add the new subproblem to the end of the active list */ alpar@9: if (tree->head == NULL) alpar@9: tree->head = node; alpar@9: else alpar@9: tree->tail->next = node; alpar@9: tree->tail = node; alpar@9: tree->a_cnt++; alpar@9: tree->n_cnt++; alpar@9: tree->t_cnt++; alpar@9: /* increase the number of child subproblems */ alpar@9: if (parent == NULL) alpar@9: xassert(p == 1); alpar@9: else alpar@9: parent->count++; alpar@9: return node; alpar@9: } alpar@9: alpar@9: void ios_clone_node(glp_tree *tree, int p, int nnn, int ref[]) alpar@9: { IOSNPD *node; alpar@9: int k; alpar@9: /* obtain pointer to the subproblem to be cloned */ alpar@9: xassert(1 <= p && p <= tree->nslots); alpar@9: node = tree->slot[p].node; alpar@9: xassert(node != NULL); alpar@9: /* the specified subproblem must be active */ alpar@9: xassert(node->count == 0); alpar@9: /* and must be in the frozen state */ alpar@9: xassert(tree->curr != node); alpar@9: /* remove the specified subproblem from the active list, because alpar@9: it becomes inactive */ alpar@9: if (node->prev == NULL) alpar@9: tree->head = node->next; alpar@9: else alpar@9: node->prev->next = node->next; alpar@9: if (node->next == NULL) alpar@9: tree->tail = node->prev; alpar@9: else alpar@9: node->next->prev = node->prev; alpar@9: node->prev = node->next = NULL; alpar@9: tree->a_cnt--; alpar@9: /* create clone subproblems */ alpar@9: xassert(nnn > 0); alpar@9: for (k = 1; k <= nnn; k++) alpar@9: ref[k] = new_node(tree, node)->p; alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_delete_node - delete specified subproblem alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * void ios_delete_node(glp_tree *tree, int p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine ios_delete_node deletes the specified subproblem, whose alpar@9: * reference number is p. The subproblem must be active and must be in alpar@9: * the frozen state (i.e. it must not be the current subproblem). alpar@9: * alpar@9: * Note that deletion is performed recursively, i.e. if a subproblem to alpar@9: * be deleted is the only child of its parent, the parent subproblem is alpar@9: * also deleted, etc. */ alpar@9: alpar@9: void ios_delete_node(glp_tree *tree, int p) alpar@9: { IOSNPD *node, *temp; alpar@9: /* obtain pointer to the subproblem to be deleted */ alpar@9: xassert(1 <= p && p <= tree->nslots); alpar@9: node = tree->slot[p].node; alpar@9: xassert(node != NULL); alpar@9: /* the specified subproblem must be active */ alpar@9: xassert(node->count == 0); alpar@9: /* and must be in the frozen state */ alpar@9: xassert(tree->curr != node); alpar@9: /* remove the specified subproblem from the active list, because alpar@9: it is gone from the tree */ alpar@9: if (node->prev == NULL) alpar@9: tree->head = node->next; alpar@9: else alpar@9: node->prev->next = node->next; alpar@9: if (node->next == NULL) alpar@9: tree->tail = node->prev; alpar@9: else alpar@9: node->next->prev = node->prev; alpar@9: node->prev = node->next = NULL; alpar@9: tree->a_cnt--; alpar@9: loop: /* recursive deletion starts here */ alpar@9: /* delete the bound change list */ alpar@9: { IOSBND *b; alpar@9: while (node->b_ptr != NULL) alpar@9: { b = node->b_ptr; alpar@9: node->b_ptr = b->next; alpar@9: dmp_free_atom(tree->pool, b, sizeof(IOSBND)); alpar@9: } alpar@9: } alpar@9: /* delete the status change list */ alpar@9: { IOSTAT *s; alpar@9: while (node->s_ptr != NULL) alpar@9: { s = node->s_ptr; alpar@9: node->s_ptr = s->next; alpar@9: dmp_free_atom(tree->pool, s, sizeof(IOSTAT)); alpar@9: } alpar@9: } alpar@9: /* delete the row addition list */ alpar@9: while (node->r_ptr != NULL) alpar@9: { IOSROW *r; alpar@9: r = node->r_ptr; alpar@9: if (r->name != NULL) alpar@9: dmp_free_atom(tree->pool, r->name, strlen(r->name)+1); alpar@9: while (r->ptr != NULL) alpar@9: { IOSAIJ *a; alpar@9: a = r->ptr; alpar@9: r->ptr = a->next; alpar@9: dmp_free_atom(tree->pool, a, sizeof(IOSAIJ)); alpar@9: } alpar@9: node->r_ptr = r->next; alpar@9: dmp_free_atom(tree->pool, r, sizeof(IOSROW)); alpar@9: } alpar@9: #if 0 alpar@9: /* delete the edge addition list */ alpar@9: /* delete the clique addition list */ alpar@9: /* (not implemented yet) */ alpar@9: xassert(node->own_nn == 0); alpar@9: xassert(node->own_nc == 0); alpar@9: xassert(node->e_ptr == NULL); alpar@9: #endif alpar@9: /* free application-specific data */ alpar@9: if (tree->parm->cb_size == 0) alpar@9: xassert(node->data == NULL); alpar@9: else alpar@9: dmp_free_atom(tree->pool, node->data, tree->parm->cb_size); alpar@9: /* free the corresponding node slot */ alpar@9: p = node->p; alpar@9: xassert(tree->slot[p].node == node); alpar@9: tree->slot[p].node = NULL; alpar@9: tree->slot[p].next = tree->avail; alpar@9: tree->avail = p; alpar@9: /* save pointer to the parent subproblem */ alpar@9: temp = node->up; alpar@9: /* delete the subproblem descriptor */ alpar@9: dmp_free_atom(tree->pool, node, sizeof(IOSNPD)); alpar@9: tree->n_cnt--; alpar@9: /* take pointer to the parent subproblem */ alpar@9: node = temp; alpar@9: if (node != NULL) alpar@9: { /* the parent subproblem exists; decrease the number of its alpar@9: child subproblems */ alpar@9: xassert(node->count > 0); alpar@9: node->count--; alpar@9: /* if now the parent subproblem has no childs, it also must be alpar@9: deleted */ alpar@9: if (node->count == 0) goto loop; alpar@9: } alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_delete_tree - delete branch-and-bound tree alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * void ios_delete_tree(glp_tree *tree); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine ios_delete_tree deletes the branch-and-bound tree, which alpar@9: * the parameter tree points to, and frees all the memory allocated to alpar@9: * this program object. alpar@9: * alpar@9: * On exit components of the problem object are restored to correspond alpar@9: * to the original MIP passed to the routine ios_create_tree. */ alpar@9: alpar@9: void ios_delete_tree(glp_tree *tree) alpar@9: { glp_prob *mip = tree->mip; alpar@9: int i, j; alpar@9: int m = mip->m; alpar@9: int n = mip->n; alpar@9: xassert(mip->tree == tree); alpar@9: /* remove all additional rows */ alpar@9: if (m != tree->orig_m) alpar@9: { int nrs, *num; alpar@9: nrs = m - tree->orig_m; alpar@9: xassert(nrs > 0); alpar@9: num = xcalloc(1+nrs, sizeof(int)); alpar@9: for (i = 1; i <= nrs; i++) num[i] = tree->orig_m + i; alpar@9: glp_del_rows(mip, nrs, num); alpar@9: xfree(num); alpar@9: } alpar@9: m = tree->orig_m; alpar@9: /* restore original attributes of rows and columns */ alpar@9: xassert(m == tree->orig_m); alpar@9: xassert(n == tree->n); alpar@9: for (i = 1; i <= m; i++) alpar@9: { glp_set_row_bnds(mip, i, tree->orig_type[i], alpar@9: tree->orig_lb[i], tree->orig_ub[i]); alpar@9: glp_set_row_stat(mip, i, tree->orig_stat[i]); alpar@9: mip->row[i]->prim = tree->orig_prim[i]; alpar@9: mip->row[i]->dual = tree->orig_dual[i]; alpar@9: } alpar@9: for (j = 1; j <= n; j++) alpar@9: { glp_set_col_bnds(mip, j, tree->orig_type[m+j], alpar@9: tree->orig_lb[m+j], tree->orig_ub[m+j]); alpar@9: glp_set_col_stat(mip, j, tree->orig_stat[m+j]); alpar@9: mip->col[j]->prim = tree->orig_prim[m+j]; alpar@9: mip->col[j]->dual = tree->orig_dual[m+j]; alpar@9: } alpar@9: mip->pbs_stat = mip->dbs_stat = GLP_FEAS; alpar@9: mip->obj_val = tree->orig_obj; alpar@9: /* delete the branch-and-bound tree */ alpar@9: xassert(tree->local != NULL); alpar@9: ios_delete_pool(tree, tree->local); alpar@9: dmp_delete_pool(tree->pool); alpar@9: xfree(tree->orig_type); alpar@9: xfree(tree->orig_lb); alpar@9: xfree(tree->orig_ub); alpar@9: xfree(tree->orig_stat); alpar@9: xfree(tree->orig_prim); alpar@9: xfree(tree->orig_dual); alpar@9: xfree(tree->slot); alpar@9: if (tree->root_type != NULL) xfree(tree->root_type); alpar@9: if (tree->root_lb != NULL) xfree(tree->root_lb); alpar@9: if (tree->root_ub != NULL) xfree(tree->root_ub); alpar@9: if (tree->root_stat != NULL) xfree(tree->root_stat); alpar@9: xfree(tree->non_int); alpar@9: #if 0 alpar@9: xfree(tree->n_ref); alpar@9: xfree(tree->c_ref); alpar@9: xfree(tree->j_ref); alpar@9: #endif alpar@9: if (tree->pcost != NULL) ios_pcost_free(tree); alpar@9: xfree(tree->iwrk); alpar@9: xfree(tree->dwrk); alpar@9: #if 0 alpar@9: scg_delete_graph(tree->g); alpar@9: #endif alpar@9: if (tree->pred_type != NULL) xfree(tree->pred_type); alpar@9: if (tree->pred_lb != NULL) xfree(tree->pred_lb); alpar@9: if (tree->pred_ub != NULL) xfree(tree->pred_ub); alpar@9: if (tree->pred_stat != NULL) xfree(tree->pred_stat); alpar@9: #if 0 alpar@9: xassert(tree->cut_gen == NULL); alpar@9: #endif alpar@9: xassert(tree->mir_gen == NULL); alpar@9: xassert(tree->clq_gen == NULL); alpar@9: xfree(tree); alpar@9: mip->tree = NULL; alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_eval_degrad - estimate obj. degrad. for down- and up-branches alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * void ios_eval_degrad(glp_tree *tree, int j, double *dn, double *up); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * Given optimal basis to LP relaxation of the current subproblem the alpar@9: * routine ios_eval_degrad performs the dual ratio test to compute the alpar@9: * objective values in the adjacent basis for down- and up-branches, alpar@9: * which are stored in locations *dn and *up, assuming that x[j] is a alpar@9: * variable chosen to branch upon. */ alpar@9: alpar@9: void ios_eval_degrad(glp_tree *tree, int j, double *dn, double *up) alpar@9: { glp_prob *mip = tree->mip; alpar@9: int m = mip->m, n = mip->n; alpar@9: int len, kase, k, t, stat; alpar@9: double alfa, beta, gamma, delta, dz; alpar@9: int *ind = tree->iwrk; alpar@9: double *val = tree->dwrk; alpar@9: /* current basis must be optimal */ alpar@9: xassert(glp_get_status(mip) == GLP_OPT); alpar@9: /* basis factorization must exist */ alpar@9: xassert(glp_bf_exists(mip)); alpar@9: /* obtain (fractional) value of x[j] in optimal basic solution alpar@9: to LP relaxation of the current subproblem */ alpar@9: xassert(1 <= j && j <= n); alpar@9: beta = mip->col[j]->prim; alpar@9: /* since the value of x[j] is fractional, it is basic; compute alpar@9: corresponding row of the simplex table */ alpar@9: len = lpx_eval_tab_row(mip, m+j, ind, val); alpar@9: /* kase < 0 means down-branch; kase > 0 means up-branch */ alpar@9: for (kase = -1; kase <= +1; kase += 2) alpar@9: { /* for down-branch we introduce new upper bound floor(beta) alpar@9: for x[j]; similarly, for up-branch we introduce new lower alpar@9: bound ceil(beta) for x[j]; in the current basis this new alpar@9: upper/lower bound is violated, so in the adjacent basis alpar@9: x[j] will leave the basis and go to its new upper/lower alpar@9: bound; we need to know which non-basic variable x[k] should alpar@9: enter the basis to keep dual feasibility */ alpar@9: #if 0 /* 23/XI-2009 */ alpar@9: k = lpx_dual_ratio_test(mip, len, ind, val, kase, 1e-7); alpar@9: #else alpar@9: k = lpx_dual_ratio_test(mip, len, ind, val, kase, 1e-9); alpar@9: #endif alpar@9: /* if no variable has been chosen, current basis being primal alpar@9: infeasible due to the new upper/lower bound of x[j] is dual alpar@9: unbounded, therefore, LP relaxation to corresponding branch alpar@9: has no primal feasible solution */ alpar@9: if (k == 0) alpar@9: { if (mip->dir == GLP_MIN) alpar@9: { if (kase < 0) alpar@9: *dn = +DBL_MAX; alpar@9: else alpar@9: *up = +DBL_MAX; alpar@9: } alpar@9: else if (mip->dir == GLP_MAX) alpar@9: { if (kase < 0) alpar@9: *dn = -DBL_MAX; alpar@9: else alpar@9: *up = -DBL_MAX; alpar@9: } alpar@9: else alpar@9: xassert(mip != mip); alpar@9: continue; alpar@9: } alpar@9: xassert(1 <= k && k <= m+n); alpar@9: /* row of the simplex table corresponding to specified basic alpar@9: variable x[j] is the following: alpar@9: x[j] = ... + alfa * x[k] + ... ; alpar@9: we need to know influence coefficient, alfa, at non-basic alpar@9: variable x[k] chosen with the dual ratio test */ alpar@9: for (t = 1; t <= len; t++) alpar@9: if (ind[t] == k) break; alpar@9: xassert(1 <= t && t <= len); alpar@9: alfa = val[t]; alpar@9: /* determine status and reduced cost of variable x[k] */ alpar@9: if (k <= m) alpar@9: { stat = mip->row[k]->stat; alpar@9: gamma = mip->row[k]->dual; alpar@9: } alpar@9: else alpar@9: { stat = mip->col[k-m]->stat; alpar@9: gamma = mip->col[k-m]->dual; alpar@9: } alpar@9: /* x[k] cannot be basic or fixed non-basic */ alpar@9: xassert(stat == GLP_NL || stat == GLP_NU || stat == GLP_NF); alpar@9: /* if the current basis is dual degenerative, some reduced alpar@9: costs, which are close to zero, may have wrong sign due to alpar@9: round-off errors, so correct the sign of gamma */ alpar@9: if (mip->dir == GLP_MIN) alpar@9: { if (stat == GLP_NL && gamma < 0.0 || alpar@9: stat == GLP_NU && gamma > 0.0 || alpar@9: stat == GLP_NF) gamma = 0.0; alpar@9: } alpar@9: else if (mip->dir == GLP_MAX) alpar@9: { if (stat == GLP_NL && gamma > 0.0 || alpar@9: stat == GLP_NU && gamma < 0.0 || alpar@9: stat == GLP_NF) gamma = 0.0; alpar@9: } alpar@9: else alpar@9: xassert(mip != mip); alpar@9: /* determine the change of x[j] in the adjacent basis: alpar@9: delta x[j] = new x[j] - old x[j] */ alpar@9: delta = (kase < 0 ? floor(beta) : ceil(beta)) - beta; alpar@9: /* compute the change of x[k] in the adjacent basis: alpar@9: delta x[k] = new x[k] - old x[k] = delta x[j] / alfa */ alpar@9: delta /= alfa; alpar@9: /* compute the change of the objective in the adjacent basis: alpar@9: delta z = new z - old z = gamma * delta x[k] */ alpar@9: dz = gamma * delta; alpar@9: if (mip->dir == GLP_MIN) alpar@9: xassert(dz >= 0.0); alpar@9: else if (mip->dir == GLP_MAX) alpar@9: xassert(dz <= 0.0); alpar@9: else alpar@9: xassert(mip != mip); alpar@9: /* compute the new objective value in the adjacent basis: alpar@9: new z = old z + delta z */ alpar@9: if (kase < 0) alpar@9: *dn = mip->obj_val + dz; alpar@9: else alpar@9: *up = mip->obj_val + dz; alpar@9: } alpar@9: /*xprintf("obj = %g; dn = %g; up = %g\n", alpar@9: mip->obj_val, *dn, *up);*/ alpar@9: return; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_round_bound - improve local bound by rounding alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * double ios_round_bound(glp_tree *tree, double bound); alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * For the given local bound for any integer feasible solution to the alpar@9: * current subproblem the routine ios_round_bound returns an improved alpar@9: * local bound for the same integer feasible solution. alpar@9: * alpar@9: * BACKGROUND alpar@9: * alpar@9: * Let the current subproblem has the following objective function: alpar@9: * alpar@9: * z = sum c[j] * x[j] + s >= b, (1) alpar@9: * j in J alpar@9: * alpar@9: * where J = {j: c[j] is non-zero and integer, x[j] is integer}, s is alpar@9: * the sum of terms corresponding to fixed variables, b is an initial alpar@9: * local bound (minimization). alpar@9: * alpar@9: * From (1) it follows that: alpar@9: * alpar@9: * d * sum (c[j] / d) * x[j] + s >= b, (2) alpar@9: * j in J alpar@9: * alpar@9: * or, equivalently, alpar@9: * alpar@9: * sum (c[j] / d) * x[j] >= (b - s) / d = h, (3) alpar@9: * j in J alpar@9: * alpar@9: * where d = gcd(c[j]). Since the left-hand side of (3) is integer, alpar@9: * h = (b - s) / d can be rounded up to the nearest integer: alpar@9: * alpar@9: * h' = ceil(h) = (b' - s) / d, (4) alpar@9: * alpar@9: * that gives an rounded, improved local bound: alpar@9: * alpar@9: * b' = d * h' + s. (5) alpar@9: * alpar@9: * In case of maximization '>=' in (1) should be replaced by '<=' that alpar@9: * leads to the following formula: alpar@9: * alpar@9: * h' = floor(h) = (b' - s) / d, (6) alpar@9: * alpar@9: * which should used in the same way as (4). alpar@9: * alpar@9: * NOTE: If b is a valid local bound for a child of the current alpar@9: * subproblem, b' is also valid for that child subproblem. */ alpar@9: alpar@9: double ios_round_bound(glp_tree *tree, double bound) alpar@9: { glp_prob *mip = tree->mip; alpar@9: int n = mip->n; alpar@9: int d, j, nn, *c = tree->iwrk; alpar@9: double s, h; alpar@9: /* determine c[j] and compute s */ alpar@9: nn = 0, s = mip->c0, d = 0; alpar@9: for (j = 1; j <= n; j++) alpar@9: { GLPCOL *col = mip->col[j]; alpar@9: if (col->coef == 0.0) continue; alpar@9: if (col->type == GLP_FX) alpar@9: { /* fixed variable */ alpar@9: s += col->coef * col->prim; alpar@9: } alpar@9: else alpar@9: { /* non-fixed variable */ alpar@9: if (col->kind != GLP_IV) goto skip; alpar@9: if (col->coef != floor(col->coef)) goto skip; alpar@9: if (fabs(col->coef) <= (double)INT_MAX) alpar@9: c[++nn] = (int)fabs(col->coef); alpar@9: else alpar@9: d = 1; alpar@9: } alpar@9: } alpar@9: /* compute d = gcd(c[1],...c[nn]) */ alpar@9: if (d == 0) alpar@9: { if (nn == 0) goto skip; alpar@9: d = gcdn(nn, c); alpar@9: } alpar@9: xassert(d > 0); alpar@9: /* compute new local bound */ alpar@9: if (mip->dir == GLP_MIN) alpar@9: { if (bound != +DBL_MAX) alpar@9: { h = (bound - s) / (double)d; alpar@9: if (h >= floor(h) + 0.001) alpar@9: { /* round up */ alpar@9: h = ceil(h); alpar@9: /*xprintf("d = %d; old = %g; ", d, bound);*/ alpar@9: bound = (double)d * h + s; alpar@9: /*xprintf("new = %g\n", bound);*/ alpar@9: } alpar@9: } alpar@9: } alpar@9: else if (mip->dir == GLP_MAX) alpar@9: { if (bound != -DBL_MAX) alpar@9: { h = (bound - s) / (double)d; alpar@9: if (h <= ceil(h) - 0.001) alpar@9: { /* round down */ alpar@9: h = floor(h); alpar@9: bound = (double)d * h + s; alpar@9: } alpar@9: } alpar@9: } alpar@9: else alpar@9: xassert(mip != mip); alpar@9: skip: return bound; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_is_hopeful - check if subproblem is hopeful alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * int ios_is_hopeful(glp_tree *tree, double bound); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * Given the local bound of a subproblem the routine ios_is_hopeful alpar@9: * checks if the subproblem can have an integer optimal solution which alpar@9: * is better than the best one currently known. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * If the subproblem can have a better integer optimal solution, the alpar@9: * routine returns non-zero; otherwise, if the corresponding branch can alpar@9: * be pruned, the routine returns zero. */ alpar@9: alpar@9: int ios_is_hopeful(glp_tree *tree, double bound) alpar@9: { glp_prob *mip = tree->mip; alpar@9: int ret = 1; alpar@9: double eps; alpar@9: if (mip->mip_stat == GLP_FEAS) alpar@9: { eps = tree->parm->tol_obj * (1.0 + fabs(mip->mip_obj)); alpar@9: switch (mip->dir) alpar@9: { case GLP_MIN: alpar@9: if (bound >= mip->mip_obj - eps) ret = 0; alpar@9: break; alpar@9: case GLP_MAX: alpar@9: if (bound <= mip->mip_obj + eps) ret = 0; alpar@9: break; alpar@9: default: alpar@9: xassert(mip != mip); alpar@9: } alpar@9: } alpar@9: else alpar@9: { switch (mip->dir) alpar@9: { case GLP_MIN: alpar@9: if (bound == +DBL_MAX) ret = 0; alpar@9: break; alpar@9: case GLP_MAX: alpar@9: if (bound == -DBL_MAX) ret = 0; alpar@9: break; alpar@9: default: alpar@9: xassert(mip != mip); alpar@9: } alpar@9: } alpar@9: return ret; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_best_node - find active node with best local bound alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * int ios_best_node(glp_tree *tree); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine ios_best_node finds an active node whose local bound is alpar@9: * best among other active nodes. alpar@9: * alpar@9: * It is understood that the integer optimal solution of the original alpar@9: * mip problem cannot be better than the best bound, so the best bound alpar@9: * is an lower (minimization) or upper (maximization) global bound for alpar@9: * the original problem. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine ios_best_node returns the subproblem reference number alpar@9: * for the best node. However, if the tree is empty, it returns zero. */ alpar@9: alpar@9: int ios_best_node(glp_tree *tree) alpar@9: { IOSNPD *node, *best = NULL; alpar@9: switch (tree->mip->dir) alpar@9: { case GLP_MIN: alpar@9: /* minimization */ alpar@9: for (node = tree->head; node != NULL; node = node->next) alpar@9: if (best == NULL || best->bound > node->bound) alpar@9: best = node; alpar@9: break; alpar@9: case GLP_MAX: alpar@9: /* maximization */ alpar@9: for (node = tree->head; node != NULL; node = node->next) alpar@9: if (best == NULL || best->bound < node->bound) alpar@9: best = node; alpar@9: break; alpar@9: default: alpar@9: xassert(tree != tree); alpar@9: } alpar@9: return best == NULL ? 0 : best->p; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_relative_gap - compute relative mip gap alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * double ios_relative_gap(glp_tree *tree); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine ios_relative_gap computes the relative mip gap using the alpar@9: * formula: alpar@9: * alpar@9: * gap = |best_mip - best_bnd| / (|best_mip| + DBL_EPSILON), alpar@9: * alpar@9: * where best_mip is the best integer feasible solution found so far, alpar@9: * best_bnd is the best (global) bound. If no integer feasible solution alpar@9: * has been found yet, rel_gap is set to DBL_MAX. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine ios_relative_gap returns the relative mip gap. */ alpar@9: alpar@9: double ios_relative_gap(glp_tree *tree) alpar@9: { glp_prob *mip = tree->mip; alpar@9: int p; alpar@9: double best_mip, best_bnd, gap; alpar@9: if (mip->mip_stat == GLP_FEAS) alpar@9: { best_mip = mip->mip_obj; alpar@9: p = ios_best_node(tree); alpar@9: if (p == 0) alpar@9: { /* the tree is empty */ alpar@9: gap = 0.0; alpar@9: } alpar@9: else alpar@9: { best_bnd = tree->slot[p].node->bound; alpar@9: gap = fabs(best_mip - best_bnd) / (fabs(best_mip) + alpar@9: DBL_EPSILON); alpar@9: } alpar@9: } alpar@9: else alpar@9: { /* no integer feasible solution has been found yet */ alpar@9: gap = DBL_MAX; alpar@9: } alpar@9: return gap; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * ios_solve_node - solve LP relaxation of current subproblem alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpios.h" alpar@9: * int ios_solve_node(glp_tree *tree); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine ios_solve_node re-optimizes LP relaxation of the current alpar@9: * subproblem using the dual simplex method. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine returns the code which is reported by glp_simplex. */ alpar@9: alpar@9: int ios_solve_node(glp_tree *tree) alpar@9: { glp_prob *mip = tree->mip; alpar@9: glp_smcp parm; alpar@9: int ret; alpar@9: /* the current subproblem must exist */ alpar@9: xassert(tree->curr != NULL); alpar@9: /* set some control parameters */ alpar@9: glp_init_smcp(&parm); alpar@9: switch (tree->parm->msg_lev) alpar@9: { case GLP_MSG_OFF: alpar@9: parm.msg_lev = GLP_MSG_OFF; break; alpar@9: case GLP_MSG_ERR: alpar@9: parm.msg_lev = GLP_MSG_ERR; break; alpar@9: case GLP_MSG_ON: alpar@9: case GLP_MSG_ALL: alpar@9: parm.msg_lev = GLP_MSG_ON; break; alpar@9: case GLP_MSG_DBG: alpar@9: parm.msg_lev = GLP_MSG_ALL; break; alpar@9: default: alpar@9: xassert(tree != tree); alpar@9: } alpar@9: parm.meth = GLP_DUALP; alpar@9: if (tree->parm->msg_lev < GLP_MSG_DBG) alpar@9: parm.out_dly = tree->parm->out_dly; alpar@9: else alpar@9: parm.out_dly = 0; alpar@9: /* if the incumbent objective value is already known, use it to alpar@9: prematurely terminate the dual simplex search */ alpar@9: if (mip->mip_stat == GLP_FEAS) alpar@9: { switch (tree->mip->dir) alpar@9: { case GLP_MIN: alpar@9: parm.obj_ul = mip->mip_obj; alpar@9: break; alpar@9: case GLP_MAX: alpar@9: parm.obj_ll = mip->mip_obj; alpar@9: break; alpar@9: default: alpar@9: xassert(mip != mip); alpar@9: } alpar@9: } alpar@9: /* try to solve/re-optimize the LP relaxation */ alpar@9: ret = glp_simplex(mip, &parm); alpar@9: tree->curr->solved++; alpar@9: #if 0 alpar@9: xprintf("ret = %d; status = %d; pbs = %d; dbs = %d; some = %d\n", alpar@9: ret, glp_get_status(mip), mip->pbs_stat, mip->dbs_stat, alpar@9: mip->some); alpar@9: lpx_print_sol(mip, "sol"); alpar@9: #endif alpar@9: return ret; alpar@9: } alpar@9: alpar@9: /**********************************************************************/ alpar@9: alpar@9: IOSPOOL *ios_create_pool(glp_tree *tree) alpar@9: { /* create cut pool */ alpar@9: IOSPOOL *pool; alpar@9: #if 0 alpar@9: pool = dmp_get_atom(tree->pool, sizeof(IOSPOOL)); alpar@9: #else alpar@9: xassert(tree == tree); alpar@9: pool = xmalloc(sizeof(IOSPOOL)); alpar@9: #endif alpar@9: pool->size = 0; alpar@9: pool->head = pool->tail = NULL; alpar@9: pool->ord = 0, pool->curr = NULL; alpar@9: return pool; alpar@9: } alpar@9: alpar@9: int ios_add_row(glp_tree *tree, IOSPOOL *pool, alpar@9: const char *name, int klass, int flags, int len, const int ind[], alpar@9: const double val[], int type, double rhs) alpar@9: { /* add row (constraint) to the cut pool */ alpar@9: IOSCUT *cut; alpar@9: IOSAIJ *aij; alpar@9: int k; alpar@9: xassert(pool != NULL); alpar@9: cut = dmp_get_atom(tree->pool, sizeof(IOSCUT)); alpar@9: if (name == NULL || name[0] == '\0') alpar@9: cut->name = NULL; alpar@9: else alpar@9: { for (k = 0; name[k] != '\0'; k++) alpar@9: { if (k == 256) alpar@9: xerror("glp_ios_add_row: cut name too long\n"); alpar@9: if (iscntrl((unsigned char)name[k])) alpar@9: xerror("glp_ios_add_row: cut name contains invalid chara" alpar@9: "cter(s)\n"); alpar@9: } alpar@9: cut->name = dmp_get_atom(tree->pool, strlen(name)+1); alpar@9: strcpy(cut->name, name); alpar@9: } alpar@9: if (!(0 <= klass && klass <= 255)) alpar@9: xerror("glp_ios_add_row: klass = %d; invalid cut class\n", alpar@9: klass); alpar@9: cut->klass = (unsigned char)klass; alpar@9: if (flags != 0) alpar@9: xerror("glp_ios_add_row: flags = %d; invalid cut flags\n", alpar@9: flags); alpar@9: cut->ptr = NULL; alpar@9: if (!(0 <= len && len <= tree->n)) alpar@9: xerror("glp_ios_add_row: len = %d; invalid cut length\n", alpar@9: len); alpar@9: for (k = 1; k <= len; k++) alpar@9: { aij = dmp_get_atom(tree->pool, sizeof(IOSAIJ)); alpar@9: if (!(1 <= ind[k] && ind[k] <= tree->n)) alpar@9: xerror("glp_ios_add_row: ind[%d] = %d; column index out of " alpar@9: "range\n", k, ind[k]); alpar@9: aij->j = ind[k]; alpar@9: aij->val = val[k]; alpar@9: aij->next = cut->ptr; alpar@9: cut->ptr = aij; alpar@9: } alpar@9: if (!(type == GLP_LO || type == GLP_UP || type == GLP_FX)) alpar@9: xerror("glp_ios_add_row: type = %d; invalid cut type\n", alpar@9: type); alpar@9: cut->type = (unsigned char)type; alpar@9: cut->rhs = rhs; alpar@9: cut->prev = pool->tail; alpar@9: cut->next = NULL; alpar@9: if (cut->prev == NULL) alpar@9: pool->head = cut; alpar@9: else alpar@9: cut->prev->next = cut; alpar@9: pool->tail = cut; alpar@9: pool->size++; alpar@9: return pool->size; alpar@9: } alpar@9: alpar@9: IOSCUT *ios_find_row(IOSPOOL *pool, int i) alpar@9: { /* find row (constraint) in the cut pool */ alpar@9: /* (smart linear search) */ alpar@9: xassert(pool != NULL); alpar@9: xassert(1 <= i && i <= pool->size); alpar@9: if (pool->ord == 0) alpar@9: { xassert(pool->curr == NULL); alpar@9: pool->ord = 1; alpar@9: pool->curr = pool->head; alpar@9: } alpar@9: xassert(pool->curr != NULL); alpar@9: if (i < pool->ord) alpar@9: { if (i < pool->ord - i) alpar@9: { pool->ord = 1; alpar@9: pool->curr = pool->head; alpar@9: while (pool->ord != i) alpar@9: { pool->ord++; alpar@9: xassert(pool->curr != NULL); alpar@9: pool->curr = pool->curr->next; alpar@9: } alpar@9: } alpar@9: else alpar@9: { while (pool->ord != i) alpar@9: { pool->ord--; alpar@9: xassert(pool->curr != NULL); alpar@9: pool->curr = pool->curr->prev; alpar@9: } alpar@9: } alpar@9: } alpar@9: else if (i > pool->ord) alpar@9: { if (i - pool->ord < pool->size - i) alpar@9: { while (pool->ord != i) alpar@9: { pool->ord++; alpar@9: xassert(pool->curr != NULL); alpar@9: pool->curr = pool->curr->next; alpar@9: } alpar@9: } alpar@9: else alpar@9: { pool->ord = pool->size; alpar@9: pool->curr = pool->tail; alpar@9: while (pool->ord != i) alpar@9: { pool->ord--; alpar@9: xassert(pool->curr != NULL); alpar@9: pool->curr = pool->curr->prev; alpar@9: } alpar@9: } alpar@9: } alpar@9: xassert(pool->ord == i); alpar@9: xassert(pool->curr != NULL); alpar@9: return pool->curr; alpar@9: } alpar@9: alpar@9: void ios_del_row(glp_tree *tree, IOSPOOL *pool, int i) alpar@9: { /* remove row (constraint) from the cut pool */ alpar@9: IOSCUT *cut; alpar@9: IOSAIJ *aij; alpar@9: xassert(pool != NULL); alpar@9: if (!(1 <= i && i <= pool->size)) alpar@9: xerror("glp_ios_del_row: i = %d; cut number out of range\n", alpar@9: i); alpar@9: cut = ios_find_row(pool, i); alpar@9: xassert(pool->curr == cut); alpar@9: if (cut->next != NULL) alpar@9: pool->curr = cut->next; alpar@9: else if (cut->prev != NULL) alpar@9: pool->ord--, pool->curr = cut->prev; alpar@9: else alpar@9: pool->ord = 0, pool->curr = NULL; alpar@9: if (cut->name != NULL) alpar@9: dmp_free_atom(tree->pool, cut->name, strlen(cut->name)+1); alpar@9: if (cut->prev == NULL) alpar@9: { xassert(pool->head == cut); alpar@9: pool->head = cut->next; alpar@9: } alpar@9: else alpar@9: { xassert(cut->prev->next == cut); alpar@9: cut->prev->next = cut->next; alpar@9: } alpar@9: if (cut->next == NULL) alpar@9: { xassert(pool->tail == cut); alpar@9: pool->tail = cut->prev; alpar@9: } alpar@9: else alpar@9: { xassert(cut->next->prev == cut); alpar@9: cut->next->prev = cut->prev; alpar@9: } alpar@9: while (cut->ptr != NULL) alpar@9: { aij = cut->ptr; alpar@9: cut->ptr = aij->next; alpar@9: dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ)); alpar@9: } alpar@9: dmp_free_atom(tree->pool, cut, sizeof(IOSCUT)); alpar@9: pool->size--; alpar@9: return; alpar@9: } alpar@9: alpar@9: void ios_clear_pool(glp_tree *tree, IOSPOOL *pool) alpar@9: { /* remove all rows (constraints) from the cut pool */ alpar@9: xassert(pool != NULL); alpar@9: while (pool->head != NULL) alpar@9: { IOSCUT *cut = pool->head; alpar@9: pool->head = cut->next; alpar@9: if (cut->name != NULL) alpar@9: dmp_free_atom(tree->pool, cut->name, strlen(cut->name)+1); alpar@9: while (cut->ptr != NULL) alpar@9: { IOSAIJ *aij = cut->ptr; alpar@9: cut->ptr = aij->next; alpar@9: dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ)); alpar@9: } alpar@9: dmp_free_atom(tree->pool, cut, sizeof(IOSCUT)); alpar@9: } alpar@9: pool->size = 0; alpar@9: pool->head = pool->tail = NULL; alpar@9: pool->ord = 0, pool->curr = NULL; alpar@9: return; alpar@9: } alpar@9: alpar@9: void ios_delete_pool(glp_tree *tree, IOSPOOL *pool) alpar@9: { /* delete cut pool */ alpar@9: xassert(pool != NULL); alpar@9: ios_clear_pool(tree, pool); alpar@9: xfree(pool); alpar@9: return; alpar@9: } alpar@9: alpar@9: /**********************************************************************/ alpar@9: alpar@9: #if 0 alpar@9: static int refer_to_node(glp_tree *tree, int j) alpar@9: { /* determine node number corresponding to binary variable x[j] or alpar@9: its complement */ alpar@9: glp_prob *mip = tree->mip; alpar@9: int n = mip->n; alpar@9: int *ref; alpar@9: if (j > 0) alpar@9: ref = tree->n_ref; alpar@9: else alpar@9: ref = tree->c_ref, j = - j; alpar@9: xassert(1 <= j && j <= n); alpar@9: if (ref[j] == 0) alpar@9: { /* new node is needed */ alpar@9: SCG *g = tree->g; alpar@9: int n_max = g->n_max; alpar@9: ref[j] = scg_add_nodes(g, 1); alpar@9: if (g->n_max > n_max) alpar@9: { int *save = tree->j_ref; alpar@9: tree->j_ref = xcalloc(1+g->n_max, sizeof(int)); alpar@9: memcpy(&tree->j_ref[1], &save[1], g->n * sizeof(int)); alpar@9: xfree(save); alpar@9: } alpar@9: xassert(ref[j] == g->n); alpar@9: tree->j_ref[ref[j]] = j; alpar@9: xassert(tree->curr != NULL); alpar@9: if (tree->curr->level > 0) tree->curr->own_nn++; alpar@9: } alpar@9: return ref[j]; alpar@9: } alpar@9: #endif alpar@9: alpar@9: #if 0 alpar@9: void ios_add_edge(glp_tree *tree, int j1, int j2) alpar@9: { /* add new edge to the conflict graph */ alpar@9: glp_prob *mip = tree->mip; alpar@9: int n = mip->n; alpar@9: SCGRIB *e; alpar@9: int first, i1, i2; alpar@9: xassert(-n <= j1 && j1 <= +n && j1 != 0); alpar@9: xassert(-n <= j2 && j2 <= +n && j2 != 0); alpar@9: xassert(j1 != j2); alpar@9: /* determine number of the first node, which was added for the alpar@9: current subproblem */ alpar@9: xassert(tree->curr != NULL); alpar@9: first = tree->g->n - tree->curr->own_nn + 1; alpar@9: /* determine node numbers for both endpoints */ alpar@9: i1 = refer_to_node(tree, j1); alpar@9: i2 = refer_to_node(tree, j2); alpar@9: /* add edge (i1,i2) to the conflict graph */ alpar@9: e = scg_add_edge(tree->g, i1, i2); alpar@9: /* if the current subproblem is not the root and both endpoints alpar@9: were created on some previous levels, save the edge */ alpar@9: if (tree->curr->level > 0 && i1 < first && i2 < first) alpar@9: { IOSRIB *rib; alpar@9: rib = dmp_get_atom(tree->pool, sizeof(IOSRIB)); alpar@9: rib->j1 = j1; alpar@9: rib->j2 = j2; alpar@9: rib->e = e; alpar@9: rib->next = tree->curr->e_ptr; alpar@9: tree->curr->e_ptr = rib; alpar@9: } alpar@9: return; alpar@9: } alpar@9: #endif alpar@9: alpar@9: /* eof */