src/glpios01.c
author Alpar Juttner <alpar@cs.elte.hu>
Mon, 06 Dec 2010 13:09:21 +0100
changeset 1 c445c931472f
permissions -rw-r--r--
Import glpk-4.45

- Generated files and doc/notes are removed
     1 /* glpios01.c */
     2 
     3 /***********************************************************************
     4 *  This code is part of GLPK (GNU Linear Programming Kit).
     5 *
     6 *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
     7 *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
     8 *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
     9 *  E-mail: <mao@gnu.org>.
    10 *
    11 *  GLPK is free software: you can redistribute it and/or modify it
    12 *  under the terms of the GNU General Public License as published by
    13 *  the Free Software Foundation, either version 3 of the License, or
    14 *  (at your option) any later version.
    15 *
    16 *  GLPK is distributed in the hope that it will be useful, but WITHOUT
    17 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    18 *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    19 *  License for more details.
    20 *
    21 *  You should have received a copy of the GNU General Public License
    22 *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    23 ***********************************************************************/
    24 
    25 #include "glpios.h"
    26 
    27 /***********************************************************************
    28 *  NAME
    29 *
    30 *  ios_create_tree - create branch-and-bound tree
    31 *
    32 *  SYNOPSIS
    33 *
    34 *  #include "glpios.h"
    35 *  glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm);
    36 *
    37 *  DESCRIPTION
    38 *
    39 *  The routine ios_create_tree creates the branch-and-bound tree.
    40 *
    41 *  Being created the tree consists of the only root subproblem whose
    42 *  reference number is 1. Note that initially the root subproblem is in
    43 *  frozen state and therefore needs to be revived.
    44 *
    45 *  RETURNS
    46 *
    47 *  The routine returns a pointer to the tree created. */
    48 
    49 static IOSNPD *new_node(glp_tree *tree, IOSNPD *parent);
    50 
    51 glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm)
    52 {     int m = mip->m;
    53       int n = mip->n;
    54       glp_tree *tree;
    55       int i, j;
    56       xassert(mip->tree == NULL);
    57       mip->tree = tree = xmalloc(sizeof(glp_tree));
    58       tree->pool = dmp_create_pool();
    59       tree->n = n;
    60       /* save original problem components */
    61       tree->orig_m = m;
    62       tree->orig_type = xcalloc(1+m+n, sizeof(char));
    63       tree->orig_lb = xcalloc(1+m+n, sizeof(double));
    64       tree->orig_ub = xcalloc(1+m+n, sizeof(double));
    65       tree->orig_stat = xcalloc(1+m+n, sizeof(char));
    66       tree->orig_prim = xcalloc(1+m+n, sizeof(double));
    67       tree->orig_dual = xcalloc(1+m+n, sizeof(double));
    68       for (i = 1; i <= m; i++)
    69       {  GLPROW *row = mip->row[i];
    70          tree->orig_type[i] = (char)row->type;
    71          tree->orig_lb[i] = row->lb;
    72          tree->orig_ub[i] = row->ub;
    73          tree->orig_stat[i] = (char)row->stat;
    74          tree->orig_prim[i] = row->prim;
    75          tree->orig_dual[i] = row->dual;
    76       }
    77       for (j = 1; j <= n; j++)
    78       {  GLPCOL *col = mip->col[j];
    79          tree->orig_type[m+j] = (char)col->type;
    80          tree->orig_lb[m+j] = col->lb;
    81          tree->orig_ub[m+j] = col->ub;
    82          tree->orig_stat[m+j] = (char)col->stat;
    83          tree->orig_prim[m+j] = col->prim;
    84          tree->orig_dual[m+j] = col->dual;
    85       }
    86       tree->orig_obj = mip->obj_val;
    87       /* initialize the branch-and-bound tree */
    88       tree->nslots = 0;
    89       tree->avail = 0;
    90       tree->slot = NULL;
    91       tree->head = tree->tail = NULL;
    92       tree->a_cnt = tree->n_cnt = tree->t_cnt = 0;
    93       /* the root subproblem is not solved yet, so its final components
    94          are unknown so far */
    95       tree->root_m = 0;
    96       tree->root_type = NULL;
    97       tree->root_lb = tree->root_ub = NULL;
    98       tree->root_stat = NULL;
    99       /* the current subproblem does not exist yet */
   100       tree->curr = NULL;
   101       tree->mip = mip;
   102       /*tree->solved = 0;*/
   103       tree->non_int = xcalloc(1+n, sizeof(char));
   104       memset(&tree->non_int[1], 0, n);
   105       /* arrays to save parent subproblem components will be allocated
   106          later */
   107       tree->pred_m = tree->pred_max = 0;
   108       tree->pred_type = NULL;
   109       tree->pred_lb = tree->pred_ub = NULL;
   110       tree->pred_stat = NULL;
   111       /* cut generator */
   112       tree->local = ios_create_pool(tree);
   113       /*tree->first_attempt = 1;*/
   114       /*tree->max_added_cuts = 0;*/
   115       /*tree->min_eff = 0.0;*/
   116       /*tree->miss = 0;*/
   117       /*tree->just_selected = 0;*/
   118       tree->mir_gen = NULL;
   119       tree->clq_gen = NULL;
   120       /*tree->round = 0;*/
   121 #if 0
   122       /* create the conflict graph */
   123       tree->n_ref = xcalloc(1+n, sizeof(int));
   124       memset(&tree->n_ref[1], 0, n * sizeof(int));
   125       tree->c_ref = xcalloc(1+n, sizeof(int));
   126       memset(&tree->c_ref[1], 0, n * sizeof(int));
   127       tree->g = scg_create_graph(0);
   128       tree->j_ref = xcalloc(1+tree->g->n_max, sizeof(int));
   129 #endif
   130       /* pseudocost branching */
   131       tree->pcost = NULL;
   132       tree->iwrk = xcalloc(1+n, sizeof(int));
   133       tree->dwrk = xcalloc(1+n, sizeof(double));
   134       /* initialize control parameters */
   135       tree->parm = parm;
   136       tree->tm_beg = xtime();
   137       tree->tm_lag = xlset(0);
   138       tree->sol_cnt = 0;
   139       /* initialize advanced solver interface */
   140       tree->reason = 0;
   141       tree->reopt = 0;
   142       tree->reinv = 0;
   143       tree->br_var = 0;
   144       tree->br_sel = 0;
   145       tree->child = 0;
   146       tree->next_p = 0;
   147       /*tree->btrack = NULL;*/
   148       tree->stop = 0;
   149       /* create the root subproblem, which initially is identical to
   150          the original MIP */
   151       new_node(tree, NULL);
   152       return tree;
   153 }
   154 
   155 /***********************************************************************
   156 *  NAME
   157 *
   158 *  ios_revive_node - revive specified subproblem
   159 *
   160 *  SYNOPSIS
   161 *
   162 *  #include "glpios.h"
   163 *  void ios_revive_node(glp_tree *tree, int p);
   164 *
   165 *  DESCRIPTION
   166 *
   167 *  The routine ios_revive_node revives the specified subproblem, whose
   168 *  reference number is p, and thereby makes it the current subproblem.
   169 *  Note that the specified subproblem must be active. Besides, if the
   170 *  current subproblem already exists, it must be frozen before reviving
   171 *  another subproblem. */
   172 
   173 void ios_revive_node(glp_tree *tree, int p)
   174 {     glp_prob *mip = tree->mip;
   175       IOSNPD *node, *root;
   176       /* obtain pointer to the specified subproblem */
   177       xassert(1 <= p && p <= tree->nslots);
   178       node = tree->slot[p].node;
   179       xassert(node != NULL);
   180       /* the specified subproblem must be active */
   181       xassert(node->count == 0);
   182       /* the current subproblem must not exist */
   183       xassert(tree->curr == NULL);
   184       /* the specified subproblem becomes current */
   185       tree->curr = node;
   186       /*tree->solved = 0;*/
   187       /* obtain pointer to the root subproblem */
   188       root = tree->slot[1].node;
   189       xassert(root != NULL);
   190       /* at this point problem object components correspond to the root
   191          subproblem, so if the root subproblem should be revived, there
   192          is nothing more to do */
   193       if (node == root) goto done;
   194       xassert(mip->m == tree->root_m);
   195       /* build path from the root to the current node */
   196       node->temp = NULL;
   197       for (node = node; node != NULL; node = node->up)
   198       {  if (node->up == NULL)
   199             xassert(node == root);
   200          else
   201             node->up->temp = node;
   202       }
   203       /* go down from the root to the current node and make necessary
   204          changes to restore components of the current subproblem */
   205       for (node = root; node != NULL; node = node->temp)
   206       {  int m = mip->m;
   207          int n = mip->n;
   208          /* if the current node is reached, the problem object at this
   209             point corresponds to its parent, so save attributes of rows
   210             and columns for the parent subproblem */
   211          if (node->temp == NULL)
   212          {  int i, j;
   213             tree->pred_m = m;
   214             /* allocate/reallocate arrays, if necessary */
   215             if (tree->pred_max < m + n)
   216             {  int new_size = m + n + 100;
   217                if (tree->pred_type != NULL) xfree(tree->pred_type);
   218                if (tree->pred_lb != NULL) xfree(tree->pred_lb);
   219                if (tree->pred_ub != NULL) xfree(tree->pred_ub);
   220                if (tree->pred_stat != NULL) xfree(tree->pred_stat);
   221                tree->pred_max = new_size;
   222                tree->pred_type = xcalloc(1+new_size, sizeof(char));
   223                tree->pred_lb = xcalloc(1+new_size, sizeof(double));
   224                tree->pred_ub = xcalloc(1+new_size, sizeof(double));
   225                tree->pred_stat = xcalloc(1+new_size, sizeof(char));
   226             }
   227             /* save row attributes */
   228             for (i = 1; i <= m; i++)
   229             {  GLPROW *row = mip->row[i];
   230                tree->pred_type[i] = (char)row->type;
   231                tree->pred_lb[i] = row->lb;
   232                tree->pred_ub[i] = row->ub;
   233                tree->pred_stat[i] = (char)row->stat;
   234             }
   235             /* save column attributes */
   236             for (j = 1; j <= n; j++)
   237             {  GLPCOL *col = mip->col[j];
   238                tree->pred_type[mip->m+j] = (char)col->type;
   239                tree->pred_lb[mip->m+j] = col->lb;
   240                tree->pred_ub[mip->m+j] = col->ub;
   241                tree->pred_stat[mip->m+j] = (char)col->stat;
   242             }
   243          }
   244          /* change bounds of rows and columns */
   245          {  IOSBND *b;
   246             for (b = node->b_ptr; b != NULL; b = b->next)
   247             {  if (b->k <= m)
   248                   glp_set_row_bnds(mip, b->k, b->type, b->lb, b->ub);
   249                else
   250                   glp_set_col_bnds(mip, b->k-m, b->type, b->lb, b->ub);
   251             }
   252          }
   253          /* change statuses of rows and columns */
   254          {  IOSTAT *s;
   255             for (s = node->s_ptr; s != NULL; s = s->next)
   256             {  if (s->k <= m)
   257                   glp_set_row_stat(mip, s->k, s->stat);
   258                else
   259                   glp_set_col_stat(mip, s->k-m, s->stat);
   260             }
   261          }
   262          /* add new rows */
   263          if (node->r_ptr != NULL)
   264          {  IOSROW *r;
   265             IOSAIJ *a;
   266             int i, len, *ind;
   267             double *val;
   268             ind = xcalloc(1+n, sizeof(int));
   269             val = xcalloc(1+n, sizeof(double));
   270             for (r = node->r_ptr; r != NULL; r = r->next)
   271             {  i = glp_add_rows(mip, 1);
   272                glp_set_row_name(mip, i, r->name);
   273 #if 1 /* 20/IX-2008 */
   274                xassert(mip->row[i]->level == 0);
   275                mip->row[i]->level = node->level;
   276                mip->row[i]->origin = r->origin;
   277                mip->row[i]->klass = r->klass;
   278 #endif
   279                glp_set_row_bnds(mip, i, r->type, r->lb, r->ub);
   280                len = 0;
   281                for (a = r->ptr; a != NULL; a = a->next)
   282                   len++, ind[len] = a->j, val[len] = a->val;
   283                glp_set_mat_row(mip, i, len, ind, val);
   284                glp_set_rii(mip, i, r->rii);
   285                glp_set_row_stat(mip, i, r->stat);
   286             }
   287             xfree(ind);
   288             xfree(val);
   289          }
   290 #if 0
   291          /* add new edges to the conflict graph */
   292          /* add new cliques to the conflict graph */
   293          /* (not implemented yet) */
   294          xassert(node->own_nn == 0);
   295          xassert(node->own_nc == 0);
   296          xassert(node->e_ptr == NULL);
   297 #endif
   298       }
   299       /* the specified subproblem has been revived */
   300       node = tree->curr;
   301       /* delete its bound change list */
   302       while (node->b_ptr != NULL)
   303       {  IOSBND *b;
   304          b = node->b_ptr;
   305          node->b_ptr = b->next;
   306          dmp_free_atom(tree->pool, b, sizeof(IOSBND));
   307       }
   308       /* delete its status change list */
   309       while (node->s_ptr != NULL)
   310       {  IOSTAT *s;
   311          s = node->s_ptr;
   312          node->s_ptr = s->next;
   313          dmp_free_atom(tree->pool, s, sizeof(IOSTAT));
   314       }
   315 #if 1 /* 20/XI-2009 */
   316       /* delete its row addition list (additional rows may appear, for
   317          example, due to branching on GUB constraints */
   318       while (node->r_ptr != NULL)
   319       {  IOSROW *r;
   320          r = node->r_ptr;
   321          node->r_ptr = r->next;
   322          xassert(r->name == NULL);
   323          while (r->ptr != NULL)
   324          {  IOSAIJ *a;
   325             a = r->ptr;
   326             r->ptr = a->next;
   327             dmp_free_atom(tree->pool, a, sizeof(IOSAIJ));
   328          }
   329          dmp_free_atom(tree->pool, r, sizeof(IOSROW));
   330       }
   331 #endif
   332 done: return;
   333 }
   334 
   335 /***********************************************************************
   336 *  NAME
   337 *
   338 *  ios_freeze_node - freeze current subproblem
   339 *
   340 *  SYNOPSIS
   341 *
   342 *  #include "glpios.h"
   343 *  void ios_freeze_node(glp_tree *tree);
   344 *
   345 *  DESCRIPTION
   346 *
   347 *  The routine ios_freeze_node freezes the current subproblem. */
   348 
   349 void ios_freeze_node(glp_tree *tree)
   350 {     glp_prob *mip = tree->mip;
   351       int m = mip->m;
   352       int n = mip->n;
   353       IOSNPD *node;
   354       /* obtain pointer to the current subproblem */
   355       node = tree->curr;
   356       xassert(node != NULL);
   357       if (node->up == NULL)
   358       {  /* freeze the root subproblem */
   359          int k;
   360          xassert(node->p == 1);
   361          xassert(tree->root_m == 0);
   362          xassert(tree->root_type == NULL);
   363          xassert(tree->root_lb == NULL);
   364          xassert(tree->root_ub == NULL);
   365          xassert(tree->root_stat == NULL);
   366          tree->root_m = m;
   367          tree->root_type = xcalloc(1+m+n, sizeof(char));
   368          tree->root_lb = xcalloc(1+m+n, sizeof(double));
   369          tree->root_ub = xcalloc(1+m+n, sizeof(double));
   370          tree->root_stat = xcalloc(1+m+n, sizeof(char));
   371          for (k = 1; k <= m+n; k++)
   372          {  if (k <= m)
   373             {  GLPROW *row = mip->row[k];
   374                tree->root_type[k] = (char)row->type;
   375                tree->root_lb[k] = row->lb;
   376                tree->root_ub[k] = row->ub;
   377                tree->root_stat[k] = (char)row->stat;
   378             }
   379             else
   380             {  GLPCOL *col = mip->col[k-m];
   381                tree->root_type[k] = (char)col->type;
   382                tree->root_lb[k] = col->lb;
   383                tree->root_ub[k] = col->ub;
   384                tree->root_stat[k] = (char)col->stat;
   385             }
   386          }
   387       }
   388       else
   389       {  /* freeze non-root subproblem */
   390          int root_m = tree->root_m;
   391          int pred_m = tree->pred_m;
   392          int i, j, k;
   393          xassert(pred_m <= m);
   394          /* build change lists for rows and columns which exist in the
   395             parent subproblem */
   396          xassert(node->b_ptr == NULL);
   397          xassert(node->s_ptr == NULL);
   398          for (k = 1; k <= pred_m + n; k++)
   399          {  int pred_type, pred_stat, type, stat;
   400             double pred_lb, pred_ub, lb, ub;
   401             /* determine attributes in the parent subproblem */
   402             pred_type = tree->pred_type[k];
   403             pred_lb = tree->pred_lb[k];
   404             pred_ub = tree->pred_ub[k];
   405             pred_stat = tree->pred_stat[k];
   406             /* determine attributes in the current subproblem */
   407             if (k <= pred_m)
   408             {  GLPROW *row = mip->row[k];
   409                type = row->type;
   410                lb = row->lb;
   411                ub = row->ub;
   412                stat = row->stat;
   413             }
   414             else
   415             {  GLPCOL *col = mip->col[k - pred_m];
   416                type = col->type;
   417                lb = col->lb;
   418                ub = col->ub;
   419                stat = col->stat;
   420             }
   421             /* save type and bounds of a row/column, if changed */
   422             if (!(pred_type == type && pred_lb == lb && pred_ub == ub))
   423             {  IOSBND *b;
   424                b = dmp_get_atom(tree->pool, sizeof(IOSBND));
   425                b->k = k;
   426                b->type = (unsigned char)type;
   427                b->lb = lb;
   428                b->ub = ub;
   429                b->next = node->b_ptr;
   430                node->b_ptr = b;
   431             }
   432             /* save status of a row/column, if changed */
   433             if (pred_stat != stat)
   434             {  IOSTAT *s;
   435                s = dmp_get_atom(tree->pool, sizeof(IOSTAT));
   436                s->k = k;
   437                s->stat = (unsigned char)stat;
   438                s->next = node->s_ptr;
   439                node->s_ptr = s;
   440             }
   441          }
   442          /* save new rows added to the current subproblem */
   443          xassert(node->r_ptr == NULL);
   444          if (pred_m < m)
   445          {  int i, len, *ind;
   446             double *val;
   447             ind = xcalloc(1+n, sizeof(int));
   448             val = xcalloc(1+n, sizeof(double));
   449             for (i = m; i > pred_m; i--)
   450             {  GLPROW *row = mip->row[i];
   451                IOSROW *r;
   452                const char *name;
   453                r = dmp_get_atom(tree->pool, sizeof(IOSROW));
   454                name = glp_get_row_name(mip, i);
   455                if (name == NULL)
   456                   r->name = NULL;
   457                else
   458                {  r->name = dmp_get_atom(tree->pool, strlen(name)+1);
   459                   strcpy(r->name, name);
   460                }
   461 #if 1 /* 20/IX-2008 */
   462                r->origin = row->origin;
   463                r->klass = row->klass;
   464 #endif
   465                r->type = (unsigned char)row->type;
   466                r->lb = row->lb;
   467                r->ub = row->ub;
   468                r->ptr = NULL;
   469                len = glp_get_mat_row(mip, i, ind, val);
   470                for (k = 1; k <= len; k++)
   471                {  IOSAIJ *a;
   472                   a = dmp_get_atom(tree->pool, sizeof(IOSAIJ));
   473                   a->j = ind[k];
   474                   a->val = val[k];
   475                   a->next = r->ptr;
   476                   r->ptr = a;
   477                }
   478                r->rii = row->rii;
   479                r->stat = (unsigned char)row->stat;
   480                r->next = node->r_ptr;
   481                node->r_ptr = r;
   482             }
   483             xfree(ind);
   484             xfree(val);
   485          }
   486          /* remove all rows missing in the root subproblem */
   487          if (m != root_m)
   488          {  int nrs, *num;
   489             nrs = m - root_m;
   490             xassert(nrs > 0);
   491             num = xcalloc(1+nrs, sizeof(int));
   492             for (i = 1; i <= nrs; i++) num[i] = root_m + i;
   493             glp_del_rows(mip, nrs, num);
   494             xfree(num);
   495          }
   496          m = mip->m;
   497          /* and restore attributes of all rows and columns for the root
   498             subproblem */
   499          xassert(m == root_m);
   500          for (i = 1; i <= m; i++)
   501          {  glp_set_row_bnds(mip, i, tree->root_type[i],
   502                tree->root_lb[i], tree->root_ub[i]);
   503             glp_set_row_stat(mip, i, tree->root_stat[i]);
   504          }
   505          for (j = 1; j <= n; j++)
   506          {  glp_set_col_bnds(mip, j, tree->root_type[m+j],
   507                tree->root_lb[m+j], tree->root_ub[m+j]);
   508             glp_set_col_stat(mip, j, tree->root_stat[m+j]);
   509          }
   510 #if 1
   511          /* remove all edges and cliques missing in the conflict graph
   512             for the root subproblem */
   513          /* (not implemented yet) */
   514 #endif
   515       }
   516       /* the current subproblem has been frozen */
   517       tree->curr = NULL;
   518       return;
   519 }
   520 
   521 /***********************************************************************
   522 *  NAME
   523 *
   524 *  ios_clone_node - clone specified subproblem
   525 *
   526 *  SYNOPSIS
   527 *
   528 *  #include "glpios.h"
   529 *  void ios_clone_node(glp_tree *tree, int p, int nnn, int ref[]);
   530 *
   531 *  DESCRIPTION
   532 *
   533 *  The routine ios_clone_node clones the specified subproblem, whose
   534 *  reference number is p, creating its nnn exact copies. Note that the
   535 *  specified subproblem must be active and must be in the frozen state
   536 *  (i.e. it must not be the current subproblem).
   537 *
   538 *  Each clone, an exact copy of the specified subproblem, becomes a new
   539 *  active subproblem added to the end of the active list. After cloning
   540 *  the specified subproblem becomes inactive.
   541 *
   542 *  The reference numbers of clone subproblems are stored to locations
   543 *  ref[1], ..., ref[nnn]. */
   544 
   545 static int get_slot(glp_tree *tree)
   546 {     int p;
   547       /* if no free slots are available, increase the room */
   548       if (tree->avail == 0)
   549       {  int nslots = tree->nslots;
   550          IOSLOT *save = tree->slot;
   551          if (nslots == 0)
   552             tree->nslots = 20;
   553          else
   554          {  tree->nslots = nslots + nslots;
   555             xassert(tree->nslots > nslots);
   556          }
   557          tree->slot = xcalloc(1+tree->nslots, sizeof(IOSLOT));
   558          if (save != NULL)
   559          {  memcpy(&tree->slot[1], &save[1], nslots * sizeof(IOSLOT));
   560             xfree(save);
   561          }
   562          /* push more free slots into the stack */
   563          for (p = tree->nslots; p > nslots; p--)
   564          {  tree->slot[p].node = NULL;
   565             tree->slot[p].next = tree->avail;
   566             tree->avail = p;
   567          }
   568       }
   569       /* pull a free slot from the stack */
   570       p = tree->avail;
   571       tree->avail = tree->slot[p].next;
   572       xassert(tree->slot[p].node == NULL);
   573       tree->slot[p].next = 0;
   574       return p;
   575 }
   576 
   577 static IOSNPD *new_node(glp_tree *tree, IOSNPD *parent)
   578 {     IOSNPD *node;
   579       int p;
   580       /* pull a free slot for the new node */
   581       p = get_slot(tree);
   582       /* create descriptor of the new subproblem */
   583       node = dmp_get_atom(tree->pool, sizeof(IOSNPD));
   584       tree->slot[p].node = node;
   585       node->p = p;
   586       node->up = parent;
   587       node->level = (parent == NULL ? 0 : parent->level + 1);
   588       node->count = 0;
   589       node->b_ptr = NULL;
   590       node->s_ptr = NULL;
   591       node->r_ptr = NULL;
   592       node->solved = 0;
   593 #if 0
   594       node->own_nn = node->own_nc = 0;
   595       node->e_ptr = NULL;
   596 #endif
   597 #if 1 /* 04/X-2008 */
   598       node->lp_obj = (parent == NULL ? (tree->mip->dir == GLP_MIN ?
   599          -DBL_MAX : +DBL_MAX) : parent->lp_obj);
   600 #endif
   601       node->bound = (parent == NULL ? (tree->mip->dir == GLP_MIN ?
   602          -DBL_MAX : +DBL_MAX) : parent->bound);
   603       node->br_var = 0;
   604       node->br_val = 0.0;
   605       node->ii_cnt = 0;
   606       node->ii_sum = 0.0;
   607 #if 1 /* 30/XI-2009 */
   608       node->changed = 0;
   609 #endif
   610       if (tree->parm->cb_size == 0)
   611          node->data = NULL;
   612       else
   613       {  node->data = dmp_get_atom(tree->pool, tree->parm->cb_size);
   614          memset(node->data, 0, tree->parm->cb_size);
   615       }
   616       node->temp = NULL;
   617       node->prev = tree->tail;
   618       node->next = NULL;
   619       /* add the new subproblem to the end of the active list */
   620       if (tree->head == NULL)
   621          tree->head = node;
   622       else
   623          tree->tail->next = node;
   624       tree->tail = node;
   625       tree->a_cnt++;
   626       tree->n_cnt++;
   627       tree->t_cnt++;
   628       /* increase the number of child subproblems */
   629       if (parent == NULL)
   630          xassert(p == 1);
   631       else
   632          parent->count++;
   633       return node;
   634 }
   635 
   636 void ios_clone_node(glp_tree *tree, int p, int nnn, int ref[])
   637 {     IOSNPD *node;
   638       int k;
   639       /* obtain pointer to the subproblem to be cloned */
   640       xassert(1 <= p && p <= tree->nslots);
   641       node = tree->slot[p].node;
   642       xassert(node != NULL);
   643       /* the specified subproblem must be active */
   644       xassert(node->count == 0);
   645       /* and must be in the frozen state */
   646       xassert(tree->curr != node);
   647       /* remove the specified subproblem from the active list, because
   648          it becomes inactive */
   649       if (node->prev == NULL)
   650          tree->head = node->next;
   651       else
   652          node->prev->next = node->next;
   653       if (node->next == NULL)
   654          tree->tail = node->prev;
   655       else
   656          node->next->prev = node->prev;
   657       node->prev = node->next = NULL;
   658       tree->a_cnt--;
   659       /* create clone subproblems */
   660       xassert(nnn > 0);
   661       for (k = 1; k <= nnn; k++)
   662          ref[k] = new_node(tree, node)->p;
   663       return;
   664 }
   665 
   666 /***********************************************************************
   667 *  NAME
   668 *
   669 *  ios_delete_node - delete specified subproblem
   670 *
   671 *  SYNOPSIS
   672 *
   673 *  #include "glpios.h"
   674 *  void ios_delete_node(glp_tree *tree, int p);
   675 *
   676 *  DESCRIPTION
   677 *
   678 *  The routine ios_delete_node deletes the specified subproblem, whose
   679 *  reference number is p. The subproblem must be active and must be in
   680 *  the frozen state (i.e. it must not be the current subproblem).
   681 *
   682 *  Note that deletion is performed recursively, i.e. if a subproblem to
   683 *  be deleted is the only child of its parent, the parent subproblem is
   684 *  also deleted, etc. */
   685 
   686 void ios_delete_node(glp_tree *tree, int p)
   687 {     IOSNPD *node, *temp;
   688       /* obtain pointer to the subproblem to be deleted */
   689       xassert(1 <= p && p <= tree->nslots);
   690       node = tree->slot[p].node;
   691       xassert(node != NULL);
   692       /* the specified subproblem must be active */
   693       xassert(node->count == 0);
   694       /* and must be in the frozen state */
   695       xassert(tree->curr != node);
   696       /* remove the specified subproblem from the active list, because
   697          it is gone from the tree */
   698       if (node->prev == NULL)
   699          tree->head = node->next;
   700       else
   701          node->prev->next = node->next;
   702       if (node->next == NULL)
   703          tree->tail = node->prev;
   704       else
   705          node->next->prev = node->prev;
   706       node->prev = node->next = NULL;
   707       tree->a_cnt--;
   708 loop: /* recursive deletion starts here */
   709       /* delete the bound change list */
   710       {  IOSBND *b;
   711          while (node->b_ptr != NULL)
   712          {  b = node->b_ptr;
   713             node->b_ptr = b->next;
   714             dmp_free_atom(tree->pool, b, sizeof(IOSBND));
   715          }
   716       }
   717       /* delete the status change list */
   718       {  IOSTAT *s;
   719          while (node->s_ptr != NULL)
   720          {  s = node->s_ptr;
   721             node->s_ptr = s->next;
   722             dmp_free_atom(tree->pool, s, sizeof(IOSTAT));
   723          }
   724       }
   725       /* delete the row addition list */
   726       while (node->r_ptr != NULL)
   727       {  IOSROW *r;
   728          r = node->r_ptr;
   729          if (r->name != NULL)
   730             dmp_free_atom(tree->pool, r->name, strlen(r->name)+1);
   731          while (r->ptr != NULL)
   732          {  IOSAIJ *a;
   733             a = r->ptr;
   734             r->ptr = a->next;
   735             dmp_free_atom(tree->pool, a, sizeof(IOSAIJ));
   736          }
   737          node->r_ptr = r->next;
   738          dmp_free_atom(tree->pool, r, sizeof(IOSROW));
   739       }
   740 #if 0
   741       /* delete the edge addition list */
   742       /* delete the clique addition list */
   743       /* (not implemented yet) */
   744       xassert(node->own_nn == 0);
   745       xassert(node->own_nc == 0);
   746       xassert(node->e_ptr == NULL);
   747 #endif
   748       /* free application-specific data */
   749       if (tree->parm->cb_size == 0)
   750          xassert(node->data == NULL);
   751       else
   752          dmp_free_atom(tree->pool, node->data, tree->parm->cb_size);
   753       /* free the corresponding node slot */
   754       p = node->p;
   755       xassert(tree->slot[p].node == node);
   756       tree->slot[p].node = NULL;
   757       tree->slot[p].next = tree->avail;
   758       tree->avail = p;
   759       /* save pointer to the parent subproblem */
   760       temp = node->up;
   761       /* delete the subproblem descriptor */
   762       dmp_free_atom(tree->pool, node, sizeof(IOSNPD));
   763       tree->n_cnt--;
   764       /* take pointer to the parent subproblem */
   765       node = temp;
   766       if (node != NULL)
   767       {  /* the parent subproblem exists; decrease the number of its
   768             child subproblems */
   769          xassert(node->count > 0);
   770          node->count--;
   771          /* if now the parent subproblem has no childs, it also must be
   772             deleted */
   773          if (node->count == 0) goto loop;
   774       }
   775       return;
   776 }
   777 
   778 /***********************************************************************
   779 *  NAME
   780 *
   781 *  ios_delete_tree - delete branch-and-bound tree
   782 *
   783 *  SYNOPSIS
   784 *
   785 *  #include "glpios.h"
   786 *  void ios_delete_tree(glp_tree *tree);
   787 *
   788 *  DESCRIPTION
   789 *
   790 *  The routine ios_delete_tree deletes the branch-and-bound tree, which
   791 *  the parameter tree points to, and frees all the memory allocated to
   792 *  this program object.
   793 *
   794 *  On exit components of the problem object are restored to correspond
   795 *  to the original MIP passed to the routine ios_create_tree. */
   796 
   797 void ios_delete_tree(glp_tree *tree)
   798 {     glp_prob *mip = tree->mip;
   799       int i, j;
   800       int m = mip->m;
   801       int n = mip->n;
   802       xassert(mip->tree == tree);
   803       /* remove all additional rows */
   804       if (m != tree->orig_m)
   805       {  int nrs, *num;
   806          nrs = m - tree->orig_m;
   807          xassert(nrs > 0);
   808          num = xcalloc(1+nrs, sizeof(int));
   809          for (i = 1; i <= nrs; i++) num[i] = tree->orig_m + i;
   810          glp_del_rows(mip, nrs, num);
   811          xfree(num);
   812       }
   813       m = tree->orig_m;
   814       /* restore original attributes of rows and columns */
   815       xassert(m == tree->orig_m);
   816       xassert(n == tree->n);
   817       for (i = 1; i <= m; i++)
   818       {  glp_set_row_bnds(mip, i, tree->orig_type[i],
   819             tree->orig_lb[i], tree->orig_ub[i]);
   820          glp_set_row_stat(mip, i, tree->orig_stat[i]);
   821          mip->row[i]->prim = tree->orig_prim[i];
   822          mip->row[i]->dual = tree->orig_dual[i];
   823       }
   824       for (j = 1; j <= n; j++)
   825       {  glp_set_col_bnds(mip, j, tree->orig_type[m+j],
   826             tree->orig_lb[m+j], tree->orig_ub[m+j]);
   827          glp_set_col_stat(mip, j, tree->orig_stat[m+j]);
   828          mip->col[j]->prim = tree->orig_prim[m+j];
   829          mip->col[j]->dual = tree->orig_dual[m+j];
   830       }
   831       mip->pbs_stat = mip->dbs_stat = GLP_FEAS;
   832       mip->obj_val = tree->orig_obj;
   833       /* delete the branch-and-bound tree */
   834       xassert(tree->local != NULL);
   835       ios_delete_pool(tree, tree->local);
   836       dmp_delete_pool(tree->pool);
   837       xfree(tree->orig_type);
   838       xfree(tree->orig_lb);
   839       xfree(tree->orig_ub);
   840       xfree(tree->orig_stat);
   841       xfree(tree->orig_prim);
   842       xfree(tree->orig_dual);
   843       xfree(tree->slot);
   844       if (tree->root_type != NULL) xfree(tree->root_type);
   845       if (tree->root_lb != NULL) xfree(tree->root_lb);
   846       if (tree->root_ub != NULL) xfree(tree->root_ub);
   847       if (tree->root_stat != NULL) xfree(tree->root_stat);
   848       xfree(tree->non_int);
   849 #if 0
   850       xfree(tree->n_ref);
   851       xfree(tree->c_ref);
   852       xfree(tree->j_ref);
   853 #endif
   854       if (tree->pcost != NULL) ios_pcost_free(tree);
   855       xfree(tree->iwrk);
   856       xfree(tree->dwrk);
   857 #if 0
   858       scg_delete_graph(tree->g);
   859 #endif
   860       if (tree->pred_type != NULL) xfree(tree->pred_type);
   861       if (tree->pred_lb != NULL) xfree(tree->pred_lb);
   862       if (tree->pred_ub != NULL) xfree(tree->pred_ub);
   863       if (tree->pred_stat != NULL) xfree(tree->pred_stat);
   864 #if 0
   865       xassert(tree->cut_gen == NULL);
   866 #endif
   867       xassert(tree->mir_gen == NULL);
   868       xassert(tree->clq_gen == NULL);
   869       xfree(tree);
   870       mip->tree = NULL;
   871       return;
   872 }
   873 
   874 /***********************************************************************
   875 *  NAME
   876 *
   877 *  ios_eval_degrad - estimate obj. degrad. for down- and up-branches
   878 *
   879 *  SYNOPSIS
   880 *
   881 *  #include "glpios.h"
   882 *  void ios_eval_degrad(glp_tree *tree, int j, double *dn, double *up);
   883 *
   884 *  DESCRIPTION
   885 *
   886 *  Given optimal basis to LP relaxation of the current subproblem the
   887 *  routine ios_eval_degrad performs the dual ratio test to compute the
   888 *  objective values in the adjacent basis for down- and up-branches,
   889 *  which are stored in locations *dn and *up, assuming that x[j] is a
   890 *  variable chosen to branch upon. */
   891 
   892 void ios_eval_degrad(glp_tree *tree, int j, double *dn, double *up)
   893 {     glp_prob *mip = tree->mip;
   894       int m = mip->m, n = mip->n;
   895       int len, kase, k, t, stat;
   896       double alfa, beta, gamma, delta, dz;
   897       int *ind = tree->iwrk;
   898       double *val = tree->dwrk;
   899       /* current basis must be optimal */
   900       xassert(glp_get_status(mip) == GLP_OPT);
   901       /* basis factorization must exist */
   902       xassert(glp_bf_exists(mip));
   903       /* obtain (fractional) value of x[j] in optimal basic solution
   904          to LP relaxation of the current subproblem */
   905       xassert(1 <= j && j <= n);
   906       beta = mip->col[j]->prim;
   907       /* since the value of x[j] is fractional, it is basic; compute
   908          corresponding row of the simplex table */
   909       len = lpx_eval_tab_row(mip, m+j, ind, val);
   910       /* kase < 0 means down-branch; kase > 0 means up-branch */
   911       for (kase = -1; kase <= +1; kase += 2)
   912       {  /* for down-branch we introduce new upper bound floor(beta)
   913             for x[j]; similarly, for up-branch we introduce new lower
   914             bound ceil(beta) for x[j]; in the current basis this new
   915             upper/lower bound is violated, so in the adjacent basis
   916             x[j] will leave the basis and go to its new upper/lower
   917             bound; we need to know which non-basic variable x[k] should
   918             enter the basis to keep dual feasibility */
   919 #if 0 /* 23/XI-2009 */
   920          k = lpx_dual_ratio_test(mip, len, ind, val, kase, 1e-7);
   921 #else
   922          k = lpx_dual_ratio_test(mip, len, ind, val, kase, 1e-9);
   923 #endif
   924          /* if no variable has been chosen, current basis being primal
   925             infeasible due to the new upper/lower bound of x[j] is dual
   926             unbounded, therefore, LP relaxation to corresponding branch
   927             has no primal feasible solution */
   928          if (k == 0)
   929          {  if (mip->dir == GLP_MIN)
   930             {  if (kase < 0)
   931                   *dn = +DBL_MAX;
   932                else
   933                   *up = +DBL_MAX;
   934             }
   935             else if (mip->dir == GLP_MAX)
   936             {  if (kase < 0)
   937                   *dn = -DBL_MAX;
   938                else
   939                   *up = -DBL_MAX;
   940             }
   941             else
   942                xassert(mip != mip);
   943             continue;
   944          }
   945          xassert(1 <= k && k <= m+n);
   946          /* row of the simplex table corresponding to specified basic
   947             variable x[j] is the following:
   948                x[j] = ... + alfa * x[k] + ... ;
   949             we need to know influence coefficient, alfa, at non-basic
   950             variable x[k] chosen with the dual ratio test */
   951          for (t = 1; t <= len; t++)
   952             if (ind[t] == k) break;
   953          xassert(1 <= t && t <= len);
   954          alfa = val[t];
   955          /* determine status and reduced cost of variable x[k] */
   956          if (k <= m)
   957          {  stat = mip->row[k]->stat;
   958             gamma = mip->row[k]->dual;
   959          }
   960          else
   961          {  stat = mip->col[k-m]->stat;
   962             gamma = mip->col[k-m]->dual;
   963          }
   964          /* x[k] cannot be basic or fixed non-basic */
   965          xassert(stat == GLP_NL || stat == GLP_NU || stat == GLP_NF);
   966          /* if the current basis is dual degenerative, some reduced
   967             costs, which are close to zero, may have wrong sign due to
   968             round-off errors, so correct the sign of gamma */
   969          if (mip->dir == GLP_MIN)
   970          {  if (stat == GLP_NL && gamma < 0.0 ||
   971                 stat == GLP_NU && gamma > 0.0 ||
   972                 stat == GLP_NF) gamma = 0.0;
   973          }
   974          else if (mip->dir == GLP_MAX)
   975          {  if (stat == GLP_NL && gamma > 0.0 ||
   976                 stat == GLP_NU && gamma < 0.0 ||
   977                 stat == GLP_NF) gamma = 0.0;
   978          }
   979          else
   980             xassert(mip != mip);
   981          /* determine the change of x[j] in the adjacent basis:
   982             delta x[j] = new x[j] - old x[j] */
   983          delta = (kase < 0 ? floor(beta) : ceil(beta)) - beta;
   984          /* compute the change of x[k] in the adjacent basis:
   985             delta x[k] = new x[k] - old x[k] = delta x[j] / alfa */
   986          delta /= alfa;
   987          /* compute the change of the objective in the adjacent basis:
   988             delta z = new z - old z = gamma * delta x[k] */
   989          dz = gamma * delta;
   990          if (mip->dir == GLP_MIN)
   991             xassert(dz >= 0.0);
   992          else if (mip->dir == GLP_MAX)
   993             xassert(dz <= 0.0);
   994          else
   995             xassert(mip != mip);
   996          /* compute the new objective value in the adjacent basis:
   997             new z = old z + delta z */
   998          if (kase < 0)
   999             *dn = mip->obj_val + dz;
  1000          else
  1001             *up = mip->obj_val + dz;
  1002       }
  1003       /*xprintf("obj = %g; dn = %g; up = %g\n",
  1004          mip->obj_val, *dn, *up);*/
  1005       return;
  1006 }
  1007 
  1008 /***********************************************************************
  1009 *  NAME
  1010 *
  1011 *  ios_round_bound - improve local bound by rounding
  1012 *
  1013 *  SYNOPSIS
  1014 *
  1015 *  #include "glpios.h"
  1016 *  double ios_round_bound(glp_tree *tree, double bound);
  1017 *
  1018 *  RETURNS
  1019 *
  1020 *  For the given local bound for any integer feasible solution to the
  1021 *  current subproblem the routine ios_round_bound returns an improved
  1022 *  local bound for the same integer feasible solution.
  1023 *
  1024 *  BACKGROUND
  1025 *
  1026 *  Let the current subproblem has the following objective function:
  1027 *
  1028 *     z =   sum  c[j] * x[j] + s >= b,                               (1)
  1029 *         j in J
  1030 *
  1031 *  where J = {j: c[j] is non-zero and integer, x[j] is integer}, s is
  1032 *  the sum of terms corresponding to fixed variables, b is an initial
  1033 *  local bound (minimization).
  1034 *
  1035 *  From (1) it follows that:
  1036 *
  1037 *     d *  sum  (c[j] / d) * x[j] + s >= b,                          (2)
  1038 *        j in J
  1039 *
  1040 *  or, equivalently,
  1041 *
  1042 *     sum  (c[j] / d) * x[j] >= (b - s) / d = h,                     (3)
  1043 *   j in J
  1044 *
  1045 *  where d = gcd(c[j]). Since the left-hand side of (3) is integer,
  1046 *  h = (b - s) / d can be rounded up to the nearest integer:
  1047 *
  1048 *     h' = ceil(h) = (b' - s) / d,                                   (4)
  1049 *
  1050 *  that gives an rounded, improved local bound:
  1051 *
  1052 *     b' = d * h' + s.                                               (5)
  1053 *
  1054 *  In case of maximization '>=' in (1) should be replaced by '<=' that
  1055 *  leads to the following formula:
  1056 *
  1057 *     h' = floor(h) = (b' - s) / d,                                  (6)
  1058 *
  1059 *  which should used in the same way as (4).
  1060 *
  1061 *  NOTE: If b is a valid local bound for a child of the current
  1062 *        subproblem, b' is also valid for that child subproblem. */
  1063 
  1064 double ios_round_bound(glp_tree *tree, double bound)
  1065 {     glp_prob *mip = tree->mip;
  1066       int n = mip->n;
  1067       int d, j, nn, *c = tree->iwrk;
  1068       double s, h;
  1069       /* determine c[j] and compute s */
  1070       nn = 0, s = mip->c0, d = 0;
  1071       for (j = 1; j <= n; j++)
  1072       {  GLPCOL *col = mip->col[j];
  1073          if (col->coef == 0.0) continue;
  1074          if (col->type == GLP_FX)
  1075          {  /* fixed variable */
  1076             s += col->coef * col->prim;
  1077          }
  1078          else
  1079          {  /* non-fixed variable */
  1080             if (col->kind != GLP_IV) goto skip;
  1081             if (col->coef != floor(col->coef)) goto skip;
  1082             if (fabs(col->coef) <= (double)INT_MAX)
  1083                c[++nn] = (int)fabs(col->coef);
  1084             else
  1085                d = 1;
  1086          }
  1087       }
  1088       /* compute d = gcd(c[1],...c[nn]) */
  1089       if (d == 0)
  1090       {  if (nn == 0) goto skip;
  1091          d = gcdn(nn, c);
  1092       }
  1093       xassert(d > 0);
  1094       /* compute new local bound */
  1095       if (mip->dir == GLP_MIN)
  1096       {  if (bound != +DBL_MAX)
  1097          {  h = (bound - s) / (double)d;
  1098             if (h >= floor(h) + 0.001)
  1099             {  /* round up */
  1100                h = ceil(h);
  1101                /*xprintf("d = %d; old = %g; ", d, bound);*/
  1102                bound = (double)d * h + s;
  1103                /*xprintf("new = %g\n", bound);*/
  1104             }
  1105          }
  1106       }
  1107       else if (mip->dir == GLP_MAX)
  1108       {  if (bound != -DBL_MAX)
  1109          {  h = (bound - s) / (double)d;
  1110             if (h <= ceil(h) - 0.001)
  1111             {  /* round down */
  1112                h = floor(h);
  1113                bound = (double)d * h + s;
  1114             }
  1115          }
  1116       }
  1117       else
  1118          xassert(mip != mip);
  1119 skip: return bound;
  1120 }
  1121 
  1122 /***********************************************************************
  1123 *  NAME
  1124 *
  1125 *  ios_is_hopeful - check if subproblem is hopeful
  1126 *
  1127 *  SYNOPSIS
  1128 *
  1129 *  #include "glpios.h"
  1130 *  int ios_is_hopeful(glp_tree *tree, double bound);
  1131 *
  1132 *  DESCRIPTION
  1133 *
  1134 *  Given the local bound of a subproblem the routine ios_is_hopeful
  1135 *  checks if the subproblem can have an integer optimal solution which
  1136 *  is better than the best one currently known.
  1137 *
  1138 *  RETURNS
  1139 *
  1140 *  If the subproblem can have a better integer optimal solution, the
  1141 *  routine returns non-zero; otherwise, if the corresponding branch can
  1142 *  be pruned, the routine returns zero. */
  1143 
  1144 int ios_is_hopeful(glp_tree *tree, double bound)
  1145 {     glp_prob *mip = tree->mip;
  1146       int ret = 1;
  1147       double eps;
  1148       if (mip->mip_stat == GLP_FEAS)
  1149       {  eps = tree->parm->tol_obj * (1.0 + fabs(mip->mip_obj));
  1150          switch (mip->dir)
  1151          {  case GLP_MIN:
  1152                if (bound >= mip->mip_obj - eps) ret = 0;
  1153                break;
  1154             case GLP_MAX:
  1155                if (bound <= mip->mip_obj + eps) ret = 0;
  1156                break;
  1157             default:
  1158                xassert(mip != mip);
  1159          }
  1160       }
  1161       else
  1162       {  switch (mip->dir)
  1163          {  case GLP_MIN:
  1164                if (bound == +DBL_MAX) ret = 0;
  1165                break;
  1166             case GLP_MAX:
  1167                if (bound == -DBL_MAX) ret = 0;
  1168                break;
  1169             default:
  1170                xassert(mip != mip);
  1171          }
  1172       }
  1173       return ret;
  1174 }
  1175 
  1176 /***********************************************************************
  1177 *  NAME
  1178 *
  1179 *  ios_best_node - find active node with best local bound
  1180 *
  1181 *  SYNOPSIS
  1182 *
  1183 *  #include "glpios.h"
  1184 *  int ios_best_node(glp_tree *tree);
  1185 *
  1186 *  DESCRIPTION
  1187 *
  1188 *  The routine ios_best_node finds an active node whose local bound is
  1189 *  best among other active nodes.
  1190 *
  1191 *  It is understood that the integer optimal solution of the original
  1192 *  mip problem cannot be better than the best bound, so the best bound
  1193 *  is an lower (minimization) or upper (maximization) global bound for
  1194 *  the original problem.
  1195 *
  1196 *  RETURNS
  1197 *
  1198 *  The routine ios_best_node returns the subproblem reference number
  1199 *  for the best node. However, if the tree is empty, it returns zero. */
  1200 
  1201 int ios_best_node(glp_tree *tree)
  1202 {     IOSNPD *node, *best = NULL;
  1203       switch (tree->mip->dir)
  1204       {  case GLP_MIN:
  1205             /* minimization */
  1206             for (node = tree->head; node != NULL; node = node->next)
  1207                if (best == NULL || best->bound > node->bound)
  1208                   best = node;
  1209             break;
  1210          case GLP_MAX:
  1211             /* maximization */
  1212             for (node = tree->head; node != NULL; node = node->next)
  1213                if (best == NULL || best->bound < node->bound)
  1214                   best = node;
  1215             break;
  1216          default:
  1217             xassert(tree != tree);
  1218       }
  1219       return best == NULL ? 0 : best->p;
  1220 }
  1221 
  1222 /***********************************************************************
  1223 *  NAME
  1224 *
  1225 *  ios_relative_gap - compute relative mip gap
  1226 *
  1227 *  SYNOPSIS
  1228 *
  1229 *  #include "glpios.h"
  1230 *  double ios_relative_gap(glp_tree *tree);
  1231 *
  1232 *  DESCRIPTION
  1233 *
  1234 *  The routine ios_relative_gap computes the relative mip gap using the
  1235 *  formula:
  1236 *
  1237 *     gap = |best_mip - best_bnd| / (|best_mip| + DBL_EPSILON),
  1238 *
  1239 *  where best_mip is the best integer feasible solution found so far,
  1240 *  best_bnd is the best (global) bound. If no integer feasible solution
  1241 *  has been found yet, rel_gap is set to DBL_MAX.
  1242 *
  1243 *  RETURNS
  1244 *
  1245 *  The routine ios_relative_gap returns the relative mip gap. */
  1246 
  1247 double ios_relative_gap(glp_tree *tree)
  1248 {     glp_prob *mip = tree->mip;
  1249       int p;
  1250       double best_mip, best_bnd, gap;
  1251       if (mip->mip_stat == GLP_FEAS)
  1252       {  best_mip = mip->mip_obj;
  1253          p = ios_best_node(tree);
  1254          if (p == 0)
  1255          {  /* the tree is empty */
  1256             gap = 0.0;
  1257          }
  1258          else
  1259          {  best_bnd = tree->slot[p].node->bound;
  1260             gap = fabs(best_mip - best_bnd) / (fabs(best_mip) +
  1261                DBL_EPSILON);
  1262          }
  1263       }
  1264       else
  1265       {  /* no integer feasible solution has been found yet */
  1266          gap = DBL_MAX;
  1267       }
  1268       return gap;
  1269 }
  1270 
  1271 /***********************************************************************
  1272 *  NAME
  1273 *
  1274 *  ios_solve_node - solve LP relaxation of current subproblem
  1275 *
  1276 *  SYNOPSIS
  1277 *
  1278 *  #include "glpios.h"
  1279 *  int ios_solve_node(glp_tree *tree);
  1280 *
  1281 *  DESCRIPTION
  1282 *
  1283 *  The routine ios_solve_node re-optimizes LP relaxation of the current
  1284 *  subproblem using the dual simplex method.
  1285 *
  1286 *  RETURNS
  1287 *
  1288 *  The routine returns the code which is reported by glp_simplex. */
  1289 
  1290 int ios_solve_node(glp_tree *tree)
  1291 {     glp_prob *mip = tree->mip;
  1292       glp_smcp parm;
  1293       int ret;
  1294       /* the current subproblem must exist */
  1295       xassert(tree->curr != NULL);
  1296       /* set some control parameters */
  1297       glp_init_smcp(&parm);
  1298       switch (tree->parm->msg_lev)
  1299       {  case GLP_MSG_OFF:
  1300             parm.msg_lev = GLP_MSG_OFF; break;
  1301          case GLP_MSG_ERR:
  1302             parm.msg_lev = GLP_MSG_ERR; break;
  1303          case GLP_MSG_ON:
  1304          case GLP_MSG_ALL:
  1305             parm.msg_lev = GLP_MSG_ON; break;
  1306          case GLP_MSG_DBG:
  1307             parm.msg_lev = GLP_MSG_ALL; break;
  1308          default:
  1309             xassert(tree != tree);
  1310       }
  1311       parm.meth = GLP_DUALP;
  1312       if (tree->parm->msg_lev < GLP_MSG_DBG)
  1313          parm.out_dly = tree->parm->out_dly;
  1314       else
  1315          parm.out_dly = 0;
  1316       /* if the incumbent objective value is already known, use it to
  1317          prematurely terminate the dual simplex search */
  1318       if (mip->mip_stat == GLP_FEAS)
  1319       {  switch (tree->mip->dir)
  1320          {  case GLP_MIN:
  1321                parm.obj_ul = mip->mip_obj;
  1322                break;
  1323             case GLP_MAX:
  1324                parm.obj_ll = mip->mip_obj;
  1325                break;
  1326             default:
  1327                xassert(mip != mip);
  1328          }
  1329       }
  1330       /* try to solve/re-optimize the LP relaxation */
  1331       ret = glp_simplex(mip, &parm);
  1332       tree->curr->solved++;
  1333 #if 0
  1334       xprintf("ret = %d; status = %d; pbs = %d; dbs = %d; some = %d\n",
  1335          ret, glp_get_status(mip), mip->pbs_stat, mip->dbs_stat,
  1336          mip->some);
  1337       lpx_print_sol(mip, "sol");
  1338 #endif
  1339       return ret;
  1340 }
  1341 
  1342 /**********************************************************************/
  1343 
  1344 IOSPOOL *ios_create_pool(glp_tree *tree)
  1345 {     /* create cut pool */
  1346       IOSPOOL *pool;
  1347 #if 0
  1348       pool = dmp_get_atom(tree->pool, sizeof(IOSPOOL));
  1349 #else
  1350       xassert(tree == tree);
  1351       pool = xmalloc(sizeof(IOSPOOL));
  1352 #endif
  1353       pool->size = 0;
  1354       pool->head = pool->tail = NULL;
  1355       pool->ord = 0, pool->curr = NULL;
  1356       return pool;
  1357 }
  1358 
  1359 int ios_add_row(glp_tree *tree, IOSPOOL *pool,
  1360       const char *name, int klass, int flags, int len, const int ind[],
  1361       const double val[], int type, double rhs)
  1362 {     /* add row (constraint) to the cut pool */
  1363       IOSCUT *cut;
  1364       IOSAIJ *aij;
  1365       int k;
  1366       xassert(pool != NULL);
  1367       cut = dmp_get_atom(tree->pool, sizeof(IOSCUT));
  1368       if (name == NULL || name[0] == '\0')
  1369          cut->name = NULL;
  1370       else
  1371       {  for (k = 0; name[k] != '\0'; k++)
  1372          {  if (k == 256)
  1373                xerror("glp_ios_add_row: cut name too long\n");
  1374             if (iscntrl((unsigned char)name[k]))
  1375                xerror("glp_ios_add_row: cut name contains invalid chara"
  1376                   "cter(s)\n");
  1377          }
  1378          cut->name = dmp_get_atom(tree->pool, strlen(name)+1);
  1379          strcpy(cut->name, name);
  1380       }
  1381       if (!(0 <= klass && klass <= 255))
  1382          xerror("glp_ios_add_row: klass = %d; invalid cut class\n",
  1383             klass);
  1384       cut->klass = (unsigned char)klass;
  1385       if (flags != 0)
  1386          xerror("glp_ios_add_row: flags = %d; invalid cut flags\n",
  1387             flags);
  1388       cut->ptr = NULL;
  1389       if (!(0 <= len && len <= tree->n))
  1390          xerror("glp_ios_add_row: len = %d; invalid cut length\n",
  1391             len);
  1392       for (k = 1; k <= len; k++)
  1393       {  aij = dmp_get_atom(tree->pool, sizeof(IOSAIJ));
  1394          if (!(1 <= ind[k] && ind[k] <= tree->n))
  1395             xerror("glp_ios_add_row: ind[%d] = %d; column index out of "
  1396                "range\n", k, ind[k]);
  1397          aij->j = ind[k];
  1398          aij->val = val[k];
  1399          aij->next = cut->ptr;
  1400          cut->ptr = aij;
  1401       }
  1402       if (!(type == GLP_LO || type == GLP_UP || type == GLP_FX))
  1403          xerror("glp_ios_add_row: type = %d; invalid cut type\n",
  1404             type);
  1405       cut->type = (unsigned char)type;
  1406       cut->rhs = rhs;
  1407       cut->prev = pool->tail;
  1408       cut->next = NULL;
  1409       if (cut->prev == NULL)
  1410          pool->head = cut;
  1411       else
  1412          cut->prev->next = cut;
  1413       pool->tail = cut;
  1414       pool->size++;
  1415       return pool->size;
  1416 }
  1417 
  1418 IOSCUT *ios_find_row(IOSPOOL *pool, int i)
  1419 {     /* find row (constraint) in the cut pool */
  1420       /* (smart linear search) */
  1421       xassert(pool != NULL);
  1422       xassert(1 <= i && i <= pool->size);
  1423       if (pool->ord == 0)
  1424       {  xassert(pool->curr == NULL);
  1425          pool->ord = 1;
  1426          pool->curr = pool->head;
  1427       }
  1428       xassert(pool->curr != NULL);
  1429       if (i < pool->ord)
  1430       {  if (i < pool->ord - i)
  1431          {  pool->ord = 1;
  1432             pool->curr = pool->head;
  1433             while (pool->ord != i)
  1434             {  pool->ord++;
  1435                xassert(pool->curr != NULL);
  1436                pool->curr = pool->curr->next;
  1437             }
  1438          }
  1439          else
  1440          {  while (pool->ord != i)
  1441             {  pool->ord--;
  1442                xassert(pool->curr != NULL);
  1443                pool->curr = pool->curr->prev;
  1444             }
  1445          }
  1446       }
  1447       else if (i > pool->ord)
  1448       {  if (i - pool->ord < pool->size - i)
  1449          {  while (pool->ord != i)
  1450             {  pool->ord++;
  1451                xassert(pool->curr != NULL);
  1452                pool->curr = pool->curr->next;
  1453             }
  1454          }
  1455          else
  1456          {  pool->ord = pool->size;
  1457             pool->curr = pool->tail;
  1458             while (pool->ord != i)
  1459             {  pool->ord--;
  1460                xassert(pool->curr != NULL);
  1461                pool->curr = pool->curr->prev;
  1462             }
  1463          }
  1464       }
  1465       xassert(pool->ord == i);
  1466       xassert(pool->curr != NULL);
  1467       return pool->curr;
  1468 }
  1469 
  1470 void ios_del_row(glp_tree *tree, IOSPOOL *pool, int i)
  1471 {     /* remove row (constraint) from the cut pool */
  1472       IOSCUT *cut;
  1473       IOSAIJ *aij;
  1474       xassert(pool != NULL);
  1475       if (!(1 <= i && i <= pool->size))
  1476          xerror("glp_ios_del_row: i = %d; cut number out of range\n",
  1477             i);
  1478       cut = ios_find_row(pool, i);
  1479       xassert(pool->curr == cut);
  1480       if (cut->next != NULL)
  1481          pool->curr = cut->next;
  1482       else if (cut->prev != NULL)
  1483          pool->ord--, pool->curr = cut->prev;
  1484       else
  1485          pool->ord = 0, pool->curr = NULL;
  1486       if (cut->name != NULL)
  1487          dmp_free_atom(tree->pool, cut->name, strlen(cut->name)+1);
  1488       if (cut->prev == NULL)
  1489       {  xassert(pool->head == cut);
  1490          pool->head = cut->next;
  1491       }
  1492       else
  1493       {  xassert(cut->prev->next == cut);
  1494          cut->prev->next = cut->next;
  1495       }
  1496       if (cut->next == NULL)
  1497       {  xassert(pool->tail == cut);
  1498          pool->tail = cut->prev;
  1499       }
  1500       else
  1501       {  xassert(cut->next->prev == cut);
  1502          cut->next->prev = cut->prev;
  1503       }
  1504       while (cut->ptr != NULL)
  1505       {  aij = cut->ptr;
  1506          cut->ptr = aij->next;
  1507          dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ));
  1508       }
  1509       dmp_free_atom(tree->pool, cut, sizeof(IOSCUT));
  1510       pool->size--;
  1511       return;
  1512 }
  1513 
  1514 void ios_clear_pool(glp_tree *tree, IOSPOOL *pool)
  1515 {     /* remove all rows (constraints) from the cut pool */
  1516       xassert(pool != NULL);
  1517       while (pool->head != NULL)
  1518       {  IOSCUT *cut = pool->head;
  1519          pool->head = cut->next;
  1520          if (cut->name != NULL)
  1521             dmp_free_atom(tree->pool, cut->name, strlen(cut->name)+1);
  1522          while (cut->ptr != NULL)
  1523          {  IOSAIJ *aij = cut->ptr;
  1524             cut->ptr = aij->next;
  1525             dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ));
  1526          }
  1527          dmp_free_atom(tree->pool, cut, sizeof(IOSCUT));
  1528       }
  1529       pool->size = 0;
  1530       pool->head = pool->tail = NULL;
  1531       pool->ord = 0, pool->curr = NULL;
  1532       return;
  1533 }
  1534 
  1535 void ios_delete_pool(glp_tree *tree, IOSPOOL *pool)
  1536 {     /* delete cut pool */
  1537       xassert(pool != NULL);
  1538       ios_clear_pool(tree, pool);
  1539       xfree(pool);
  1540       return;
  1541 }
  1542 
  1543 /**********************************************************************/
  1544 
  1545 #if 0
  1546 static int refer_to_node(glp_tree *tree, int j)
  1547 {     /* determine node number corresponding to binary variable x[j] or
  1548          its complement */
  1549       glp_prob *mip = tree->mip;
  1550       int n = mip->n;
  1551       int *ref;
  1552       if (j > 0)
  1553          ref = tree->n_ref;
  1554       else
  1555          ref = tree->c_ref, j = - j;
  1556       xassert(1 <= j && j <= n);
  1557       if (ref[j] == 0)
  1558       {  /* new node is needed */
  1559          SCG *g = tree->g;
  1560          int n_max = g->n_max;
  1561          ref[j] = scg_add_nodes(g, 1);
  1562          if (g->n_max > n_max)
  1563          {  int *save = tree->j_ref;
  1564             tree->j_ref = xcalloc(1+g->n_max, sizeof(int));
  1565             memcpy(&tree->j_ref[1], &save[1], g->n * sizeof(int));
  1566             xfree(save);
  1567          }
  1568          xassert(ref[j] == g->n);
  1569          tree->j_ref[ref[j]] = j;
  1570          xassert(tree->curr != NULL);
  1571          if (tree->curr->level > 0) tree->curr->own_nn++;
  1572       }
  1573       return ref[j];
  1574 }
  1575 #endif
  1576 
  1577 #if 0
  1578 void ios_add_edge(glp_tree *tree, int j1, int j2)
  1579 {     /* add new edge to the conflict graph */
  1580       glp_prob *mip = tree->mip;
  1581       int n = mip->n;
  1582       SCGRIB *e;
  1583       int first, i1, i2;
  1584       xassert(-n <= j1 && j1 <= +n && j1 != 0);
  1585       xassert(-n <= j2 && j2 <= +n && j2 != 0);
  1586       xassert(j1 != j2);
  1587       /* determine number of the first node, which was added for the
  1588          current subproblem */
  1589       xassert(tree->curr != NULL);
  1590       first = tree->g->n - tree->curr->own_nn + 1;
  1591       /* determine node numbers for both endpoints */
  1592       i1 = refer_to_node(tree, j1);
  1593       i2 = refer_to_node(tree, j2);
  1594       /* add edge (i1,i2) to the conflict graph */
  1595       e = scg_add_edge(tree->g, i1, i2);
  1596       /* if the current subproblem is not the root and both endpoints
  1597          were created on some previous levels, save the edge */
  1598       if (tree->curr->level > 0 && i1 < first && i2 < first)
  1599       {  IOSRIB *rib;
  1600          rib = dmp_get_atom(tree->pool, sizeof(IOSRIB));
  1601          rib->j1 = j1;
  1602          rib->j2 = j2;
  1603          rib->e = e;
  1604          rib->next = tree->curr->e_ptr;
  1605          tree->curr->e_ptr = rib;
  1606       }
  1607       return;
  1608 }
  1609 #endif
  1610 
  1611 /* eof */