COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpios09.c @ 2:4c8956a7bdf4

Last change on this file since 2:4c8956a7bdf4 was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 25.6 KB
Line 
1/* glpios09.c (branching heuristics) */
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_choose_var - select variable to branch on
31*
32*  SYNOPSIS
33*
34*  #include "glpios.h"
35*  int ios_choose_var(glp_tree *T, int *next);
36*
37*  The routine ios_choose_var chooses a variable from the candidate
38*  list to branch on. Additionally the routine provides a flag stored
39*  in the location next to suggests which of the child subproblems
40*  should be solved next.
41*
42*  RETURNS
43*
44*  The routine ios_choose_var returns the ordinal number of the column
45*  choosen. */
46
47static int branch_first(glp_tree *T, int *next);
48static int branch_last(glp_tree *T, int *next);
49static int branch_mostf(glp_tree *T, int *next);
50static int branch_drtom(glp_tree *T, int *next);
51
52int ios_choose_var(glp_tree *T, int *next)
53{     int j;
54      if (T->parm->br_tech == GLP_BR_FFV)
55      {  /* branch on first fractional variable */
56         j = branch_first(T, next);
57      }
58      else if (T->parm->br_tech == GLP_BR_LFV)
59      {  /* branch on last fractional variable */
60         j = branch_last(T, next);
61      }
62      else if (T->parm->br_tech == GLP_BR_MFV)
63      {  /* branch on most fractional variable */
64         j = branch_mostf(T, next);
65      }
66      else if (T->parm->br_tech == GLP_BR_DTH)
67      {  /* branch using the heuristic by Dreebeck and Tomlin */
68         j = branch_drtom(T, next);
69      }
70      else if (T->parm->br_tech == GLP_BR_PCH)
71      {  /* hybrid pseudocost heuristic */
72         j = ios_pcost_branch(T, next);
73      }
74      else
75         xassert(T != T);
76      return j;
77}
78
79/***********************************************************************
80*  branch_first - choose first branching variable
81*
82*  This routine looks up the list of structural variables and chooses
83*  the first one, which is of integer kind and has fractional value in
84*  optimal solution to the current LP relaxation.
85*
86*  This routine also selects the branch to be solved next where integer
87*  infeasibility of the chosen variable is less than in other one. */
88
89static int branch_first(glp_tree *T, int *_next)
90{     int j, next;
91      double beta;
92      /* choose the column to branch on */
93      for (j = 1; j <= T->n; j++)
94         if (T->non_int[j]) break;
95      xassert(1 <= j && j <= T->n);
96      /* select the branch to be solved next */
97      beta = glp_get_col_prim(T->mip, j);
98      if (beta - floor(beta) < ceil(beta) - beta)
99         next = GLP_DN_BRNCH;
100      else
101         next = GLP_UP_BRNCH;
102      *_next = next;
103      return j;
104}
105
106/***********************************************************************
107*  branch_last - choose last branching variable
108*
109*  This routine looks up the list of structural variables and chooses
110*  the last one, which is of integer kind and has fractional value in
111*  optimal solution to the current LP relaxation.
112*
113*  This routine also selects the branch to be solved next where integer
114*  infeasibility of the chosen variable is less than in other one. */
115
116static int branch_last(glp_tree *T, int *_next)
117{     int j, next;
118      double beta;
119      /* choose the column to branch on */
120      for (j = T->n; j >= 1; j--)
121         if (T->non_int[j]) break;
122      xassert(1 <= j && j <= T->n);
123      /* select the branch to be solved next */
124      beta = glp_get_col_prim(T->mip, j);
125      if (beta - floor(beta) < ceil(beta) - beta)
126         next = GLP_DN_BRNCH;
127      else
128         next = GLP_UP_BRNCH;
129      *_next = next;
130      return j;
131}
132
133/***********************************************************************
134*  branch_mostf - choose most fractional branching variable
135*
136*  This routine looks up the list of structural variables and chooses
137*  that one, which is of integer kind and has most fractional value in
138*  optimal solution to the current LP relaxation.
139*
140*  This routine also selects the branch to be solved next where integer
141*  infeasibility of the chosen variable is less than in other one.
142*
143*  (Alexander Martin notices that "...most infeasible is as good as
144*  random...".) */
145
146static int branch_mostf(glp_tree *T, int *_next)
147{     int j, jj, next;
148      double beta, most, temp;
149      /* choose the column to branch on */
150      jj = 0, most = DBL_MAX;
151      for (j = 1; j <= T->n; j++)
152      {  if (T->non_int[j])
153         {  beta = glp_get_col_prim(T->mip, j);
154            temp = floor(beta) + 0.5;
155            if (most > fabs(beta - temp))
156            {  jj = j, most = fabs(beta - temp);
157               if (beta < temp)
158                  next = GLP_DN_BRNCH;
159               else
160                  next = GLP_UP_BRNCH;
161            }
162         }
163      }
164      *_next = next;
165      return jj;
166}
167
168/***********************************************************************
169*  branch_drtom - choose branching var using Driebeck-Tomlin heuristic
170*
171*  This routine chooses a structural variable, which is required to be
172*  integral and has fractional value in optimal solution of the current
173*  LP relaxation, using a heuristic proposed by Driebeck and Tomlin.
174*
175*  The routine also selects the branch to be solved next, again due to
176*  Driebeck and Tomlin.
177*
178*  This routine is based on the heuristic proposed in:
179*
180*  Driebeck N.J. An algorithm for the solution of mixed-integer
181*  programming problems, Management Science, 12: 576-87 (1966);
182*
183*  and improved in:
184*
185*  Tomlin J.A. Branch and bound methods for integer and non-convex
186*  programming, in J.Abadie (ed.), Integer and Nonlinear Programming,
187*  North-Holland, Amsterdam, pp. 437-50 (1970).
188*
189*  Must note that this heuristic is time-expensive, because computing
190*  one-step degradation (see the routine below) requires one BTRAN for
191*  each fractional-valued structural variable. */
192
193static int branch_drtom(glp_tree *T, int *_next)
194{     glp_prob *mip = T->mip;
195      int m = mip->m;
196      int n = mip->n;
197      char *non_int = T->non_int;
198      int j, jj, k, t, next, kase, len, stat, *ind;
199      double x, dk, alfa, delta_j, delta_k, delta_z, dz_dn, dz_up,
200         dd_dn, dd_up, degrad, *val;
201      /* basic solution of LP relaxation must be optimal */
202      xassert(glp_get_status(mip) == GLP_OPT);
203      /* allocate working arrays */
204      ind = xcalloc(1+n, sizeof(int));
205      val = xcalloc(1+n, sizeof(double));
206      /* nothing has been chosen so far */
207      jj = 0, degrad = -1.0;
208      /* walk through the list of columns (structural variables) */
209      for (j = 1; j <= n; j++)
210      {  /* if j-th column is not marked as fractional, skip it */
211         if (!non_int[j]) continue;
212         /* obtain (fractional) value of j-th column in basic solution
213            of LP relaxation */
214         x = glp_get_col_prim(mip, j);
215         /* since the value of j-th column is fractional, the column is
216            basic; compute corresponding row of the simplex table */
217         len = glp_eval_tab_row(mip, m+j, ind, val);
218         /* the following fragment computes a change in the objective
219            function: delta Z = new Z - old Z, where old Z is the
220            objective value in the current optimal basis, and new Z is
221            the objective value in the adjacent basis, for two cases:
222            1) if new upper bound ub' = floor(x[j]) is introduced for
223               j-th column (down branch);
224            2) if new lower bound lb' = ceil(x[j]) is introduced for
225               j-th column (up branch);
226            since in both cases the solution remaining dual feasible
227            becomes primal infeasible, one implicit simplex iteration
228            is performed to determine the change delta Z;
229            it is obvious that new Z, which is never better than old Z,
230            is a lower (minimization) or upper (maximization) bound of
231            the objective function for down- and up-branches. */
232         for (kase = -1; kase <= +1; kase += 2)
233         {  /* if kase < 0, the new upper bound of x[j] is introduced;
234               in this case x[j] should decrease in order to leave the
235               basis and go to its new upper bound */
236            /* if kase > 0, the new lower bound of x[j] is introduced;
237               in this case x[j] should increase in order to leave the
238               basis and go to its new lower bound */
239            /* apply the dual ratio test in order to determine which
240               auxiliary or structural variable should enter the basis
241               to keep dual feasibility */
242            k = glp_dual_rtest(mip, len, ind, val, kase, 1e-9);
243            if (k != 0) k = ind[k];
244            /* if no non-basic variable has been chosen, LP relaxation
245               of corresponding branch being primal infeasible and dual
246               unbounded has no primal feasible solution; in this case
247               the change delta Z is formally set to infinity */
248            if (k == 0)
249            {  delta_z =
250                  (T->mip->dir == GLP_MIN ? +DBL_MAX : -DBL_MAX);
251               goto skip;
252            }
253            /* row of the simplex table that corresponds to non-basic
254               variable x[k] choosen by the dual ratio test is:
255                  x[j] = ... + alfa * x[k] + ...
256               where alfa is the influence coefficient (an element of
257               the simplex table row) */
258            /* determine the coefficient alfa */
259            for (t = 1; t <= len; t++) if (ind[t] == k) break;
260            xassert(1 <= t && t <= len);
261            alfa = val[t];
262            /* since in the adjacent basis the variable x[j] becomes
263               non-basic, knowing its value in the current basis we can
264               determine its change delta x[j] = new x[j] - old x[j] */
265            delta_j = (kase < 0 ? floor(x) : ceil(x)) - x;
266            /* and knowing the coefficient alfa we can determine the
267               corresponding change delta x[k] = new x[k] - old x[k],
268               where old x[k] is a value of x[k] in the current basis,
269               and new x[k] is a value of x[k] in the adjacent basis */
270            delta_k = delta_j / alfa;
271            /* Tomlin noticed that if the variable x[k] is of integer
272               kind, its change cannot be less (eventually) than one in
273               the magnitude */
274            if (k > m && glp_get_col_kind(mip, k-m) != GLP_CV)
275            {  /* x[k] is structural integer variable */
276               if (fabs(delta_k - floor(delta_k + 0.5)) > 1e-3)
277               {  if (delta_k > 0.0)
278                     delta_k = ceil(delta_k);  /* +3.14 -> +4 */
279                  else
280                     delta_k = floor(delta_k); /* -3.14 -> -4 */
281               }
282            }
283            /* now determine the status and reduced cost of x[k] in the
284               current basis */
285            if (k <= m)
286            {  stat = glp_get_row_stat(mip, k);
287               dk = glp_get_row_dual(mip, k);
288            }
289            else
290            {  stat = glp_get_col_stat(mip, k-m);
291               dk = glp_get_col_dual(mip, k-m);
292            }
293            /* if the current basis is dual degenerate, some reduced
294               costs which are close to zero may have wrong sign due to
295               round-off errors, so correct the sign of d[k] */
296            switch (T->mip->dir)
297            {  case GLP_MIN:
298                  if (stat == GLP_NL && dk < 0.0 ||
299                      stat == GLP_NU && dk > 0.0 ||
300                      stat == GLP_NF) dk = 0.0;
301                  break;
302               case GLP_MAX:
303                  if (stat == GLP_NL && dk > 0.0 ||
304                      stat == GLP_NU && dk < 0.0 ||
305                      stat == GLP_NF) dk = 0.0;
306                  break;
307               default:
308                  xassert(T != T);
309            }
310            /* now knowing the change of x[k] and its reduced cost d[k]
311               we can compute the corresponding change in the objective
312               function delta Z = new Z - old Z = d[k] * delta x[k];
313               note that due to Tomlin's modification new Z can be even
314               worse than in the adjacent basis */
315            delta_z = dk * delta_k;
316skip:       /* new Z is never better than old Z, therefore the change
317               delta Z is always non-negative (in case of minimization)
318               or non-positive (in case of maximization) */
319            switch (T->mip->dir)
320            {  case GLP_MIN: xassert(delta_z >= 0.0); break;
321               case GLP_MAX: xassert(delta_z <= 0.0); break;
322               default: xassert(T != T);
323            }
324            /* save the change in the objective fnction for down- and
325               up-branches, respectively */
326            if (kase < 0) dz_dn = delta_z; else dz_up = delta_z;
327         }
328         /* thus, in down-branch no integer feasible solution can be
329            better than Z + dz_dn, and in up-branch no integer feasible
330            solution can be better than Z + dz_up, where Z is value of
331            the objective function in the current basis */
332         /* following the heuristic by Driebeck and Tomlin we choose a
333            column (i.e. structural variable) which provides largest
334            degradation of the objective function in some of branches;
335            besides, we select the branch with smaller degradation to
336            be solved next and keep other branch with larger degradation
337            in the active list hoping to minimize the number of further
338            backtrackings */
339         if (degrad < fabs(dz_dn) || degrad < fabs(dz_up))
340         {  jj = j;
341            if (fabs(dz_dn) < fabs(dz_up))
342            {  /* select down branch to be solved next */
343               next = GLP_DN_BRNCH;
344               degrad = fabs(dz_up);
345            }
346            else
347            {  /* select up branch to be solved next */
348               next = GLP_UP_BRNCH;
349               degrad = fabs(dz_dn);
350            }
351            /* save the objective changes for printing */
352            dd_dn = dz_dn, dd_up = dz_up;
353            /* if down- or up-branch has no feasible solution, we does
354               not need to consider other candidates (in principle, the
355               corresponding branch could be pruned right now) */
356            if (degrad == DBL_MAX) break;
357         }
358      }
359      /* free working arrays */
360      xfree(ind);
361      xfree(val);
362      /* something must be chosen */
363      xassert(1 <= jj && jj <= n);
364#if 1 /* 02/XI-2009 */
365      if (degrad < 1e-6 * (1.0 + 0.001 * fabs(mip->obj_val)))
366      {  jj = branch_mostf(T, &next);
367         goto done;
368      }
369#endif
370      if (T->parm->msg_lev >= GLP_MSG_DBG)
371      {  xprintf("branch_drtom: column %d chosen to branch on\n", jj);
372         if (fabs(dd_dn) == DBL_MAX)
373            xprintf("branch_drtom: down-branch is infeasible\n");
374         else
375            xprintf("branch_drtom: down-branch bound is %.9e\n",
376               lpx_get_obj_val(mip) + dd_dn);
377         if (fabs(dd_up) == DBL_MAX)
378            xprintf("branch_drtom: up-branch   is infeasible\n");
379         else
380            xprintf("branch_drtom: up-branch   bound is %.9e\n",
381               lpx_get_obj_val(mip) + dd_up);
382      }
383done: *_next = next;
384      return jj;
385}
386
387/**********************************************************************/
388
389struct csa
390{     /* common storage area */
391      int *dn_cnt; /* int dn_cnt[1+n]; */
392      /* dn_cnt[j] is the number of subproblems, whose LP relaxations
393         have been solved and which are down-branches for variable x[j];
394         dn_cnt[j] = 0 means the down pseudocost is uninitialized */
395      double *dn_sum; /* double dn_sum[1+n]; */
396      /* dn_sum[j] is the sum of per unit degradations of the objective
397         over all dn_cnt[j] subproblems */
398      int *up_cnt; /* int up_cnt[1+n]; */
399      /* up_cnt[j] is the number of subproblems, whose LP relaxations
400         have been solved and which are up-branches for variable x[j];
401         up_cnt[j] = 0 means the up pseudocost is uninitialized */
402      double *up_sum; /* double up_sum[1+n]; */
403      /* up_sum[j] is the sum of per unit degradations of the objective
404         over all up_cnt[j] subproblems */
405};
406
407void *ios_pcost_init(glp_tree *tree)
408{     /* initialize working data used on pseudocost branching */
409      struct csa *csa;
410      int n = tree->n, j;
411      csa = xmalloc(sizeof(struct csa));
412      csa->dn_cnt = xcalloc(1+n, sizeof(int));
413      csa->dn_sum = xcalloc(1+n, sizeof(double));
414      csa->up_cnt = xcalloc(1+n, sizeof(int));
415      csa->up_sum = xcalloc(1+n, sizeof(double));
416      for (j = 1; j <= n; j++)
417      {  csa->dn_cnt[j] = csa->up_cnt[j] = 0;
418         csa->dn_sum[j] = csa->up_sum[j] = 0.0;
419      }
420      return csa;
421}
422
423static double eval_degrad(glp_prob *P, int j, double bnd)
424{     /* compute degradation of the objective on fixing x[j] at given
425         value with a limited number of dual simplex iterations */
426      /* this routine fixes column x[j] at specified value bnd,
427         solves resulting LP, and returns a lower bound to degradation
428         of the objective, degrad >= 0 */
429      glp_prob *lp;
430      glp_smcp parm;
431      int ret;
432      double degrad;
433      /* the current basis must be optimal */
434      xassert(glp_get_status(P) == GLP_OPT);
435      /* create a copy of P */
436      lp = glp_create_prob();
437      glp_copy_prob(lp, P, 0);
438      /* fix column x[j] at specified value */
439      glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd);
440      /* try to solve resulting LP */
441      glp_init_smcp(&parm);
442      parm.msg_lev = GLP_MSG_OFF;
443      parm.meth = GLP_DUAL;
444      parm.it_lim = 30;
445      parm.out_dly = 1000;
446      parm.meth = GLP_DUAL;
447      ret = glp_simplex(lp, &parm);
448      if (ret == 0 || ret == GLP_EITLIM)
449      {  if (glp_get_prim_stat(lp) == GLP_NOFEAS)
450         {  /* resulting LP has no primal feasible solution */
451            degrad = DBL_MAX;
452         }
453         else if (glp_get_dual_stat(lp) == GLP_FEAS)
454         {  /* resulting basis is optimal or at least dual feasible,
455               so we have the correct lower bound to degradation */
456            if (P->dir == GLP_MIN)
457               degrad = lp->obj_val - P->obj_val;
458            else if (P->dir == GLP_MAX)
459               degrad = P->obj_val - lp->obj_val;
460            else
461               xassert(P != P);
462            /* degradation cannot be negative by definition */
463            /* note that the lower bound to degradation may be close
464               to zero even if its exact value is zero due to round-off
465               errors on computing the objective value */
466            if (degrad < 1e-6 * (1.0 + 0.001 * fabs(P->obj_val)))
467               degrad = 0.0;
468         }
469         else
470         {  /* the final basis reported by the simplex solver is dual
471               infeasible, so we cannot determine a non-trivial lower
472               bound to degradation */
473            degrad = 0.0;
474         }
475      }
476      else
477      {  /* the simplex solver failed */
478         degrad = 0.0;
479      }
480      /* delete the copy of P */
481      glp_delete_prob(lp);
482      return degrad;
483}
484
485void ios_pcost_update(glp_tree *tree)
486{     /* update history information for pseudocost branching */
487      /* this routine is called every time when LP relaxation of the
488         current subproblem has been solved to optimality with all lazy
489         and cutting plane constraints included */
490      int j;
491      double dx, dz, psi;
492      struct csa *csa = tree->pcost;
493      xassert(csa != NULL);
494      xassert(tree->curr != NULL);
495      /* if the current subproblem is the root, skip updating */
496      if (tree->curr->up == NULL) goto skip;
497      /* determine branching variable x[j], which was used in the
498         parent subproblem to create the current subproblem */
499      j = tree->curr->up->br_var;
500      xassert(1 <= j && j <= tree->n);
501      /* determine the change dx[j] = new x[j] - old x[j],
502         where new x[j] is a value of x[j] in optimal solution to LP
503         relaxation of the current subproblem, old x[j] is a value of
504         x[j] in optimal solution to LP relaxation of the parent
505         subproblem */
506      dx = tree->mip->col[j]->prim - tree->curr->up->br_val;
507      xassert(dx != 0.0);
508      /* determine corresponding change dz = new dz - old dz in the
509         objective function value */
510      dz = tree->mip->obj_val - tree->curr->up->lp_obj;
511      /* determine per unit degradation of the objective function */
512      psi = fabs(dz / dx);
513      /* update history information */
514      if (dx < 0.0)
515      {  /* the current subproblem is down-branch */
516         csa->dn_cnt[j]++;
517         csa->dn_sum[j] += psi;
518      }
519      else /* dx > 0.0 */
520      {  /* the current subproblem is up-branch */
521         csa->up_cnt[j]++;
522         csa->up_sum[j] += psi;
523      }
524skip: return;
525}
526
527void ios_pcost_free(glp_tree *tree)
528{     /* free working area used on pseudocost branching */
529      struct csa *csa = tree->pcost;
530      xassert(csa != NULL);
531      xfree(csa->dn_cnt);
532      xfree(csa->dn_sum);
533      xfree(csa->up_cnt);
534      xfree(csa->up_sum);
535      xfree(csa);
536      tree->pcost = NULL;
537      return;
538}
539
540static double eval_psi(glp_tree *T, int j, int brnch)
541{     /* compute estimation of pseudocost of variable x[j] for down-
542         or up-branch */
543      struct csa *csa = T->pcost;
544      double beta, degrad, psi;
545      xassert(csa != NULL);
546      xassert(1 <= j && j <= T->n);
547      if (brnch == GLP_DN_BRNCH)
548      {  /* down-branch */
549         if (csa->dn_cnt[j] == 0)
550         {  /* initialize down pseudocost */
551            beta = T->mip->col[j]->prim;
552            degrad = eval_degrad(T->mip, j, floor(beta));
553            if (degrad == DBL_MAX)
554            {  psi = DBL_MAX;
555               goto done;
556            }
557            csa->dn_cnt[j] = 1;
558            csa->dn_sum[j] = degrad / (beta - floor(beta));
559         }
560         psi = csa->dn_sum[j] / (double)csa->dn_cnt[j];
561      }
562      else if (brnch == GLP_UP_BRNCH)
563      {  /* up-branch */
564         if (csa->up_cnt[j] == 0)
565         {  /* initialize up pseudocost */
566            beta = T->mip->col[j]->prim;
567            degrad = eval_degrad(T->mip, j, ceil(beta));
568            if (degrad == DBL_MAX)
569            {  psi = DBL_MAX;
570               goto done;
571            }
572            csa->up_cnt[j] = 1;
573            csa->up_sum[j] = degrad / (ceil(beta) - beta);
574         }
575         psi = csa->up_sum[j] / (double)csa->up_cnt[j];
576      }
577      else
578         xassert(brnch != brnch);
579done: return psi;
580}
581
582static void progress(glp_tree *T)
583{     /* display progress of pseudocost initialization */
584      struct csa *csa = T->pcost;
585      int j, nv = 0, ni = 0;
586      for (j = 1; j <= T->n; j++)
587      {  if (glp_ios_can_branch(T, j))
588         {  nv++;
589            if (csa->dn_cnt[j] > 0 && csa->up_cnt[j] > 0) ni++;
590         }
591      }
592      xprintf("Pseudocosts initialized for %d of %d variables\n",
593         ni, nv);
594      return;
595}
596
597int ios_pcost_branch(glp_tree *T, int *_next)
598{     /* choose branching variable with pseudocost branching */
599      glp_long t = xtime();
600      int j, jjj, sel;
601      double beta, psi, d1, d2, d, dmax;
602      /* initialize the working arrays */
603      if (T->pcost == NULL)
604         T->pcost = ios_pcost_init(T);
605      /* nothing has been chosen so far */
606      jjj = 0, dmax = -1.0;
607      /* go through the list of branching candidates */
608      for (j = 1; j <= T->n; j++)
609      {  if (!glp_ios_can_branch(T, j)) continue;
610         /* determine primal value of x[j] in optimal solution to LP
611            relaxation of the current subproblem */
612         beta = T->mip->col[j]->prim;
613         /* estimate pseudocost of x[j] for down-branch */
614         psi = eval_psi(T, j, GLP_DN_BRNCH);
615         if (psi == DBL_MAX)
616         {  /* down-branch has no primal feasible solution */
617            jjj = j, sel = GLP_DN_BRNCH;
618            goto done;
619         }
620         /* estimate degradation of the objective for down-branch */
621         d1 = psi * (beta - floor(beta));
622         /* estimate pseudocost of x[j] for up-branch */
623         psi = eval_psi(T, j, GLP_UP_BRNCH);
624         if (psi == DBL_MAX)
625         {  /* up-branch has no primal feasible solution */
626            jjj = j, sel = GLP_UP_BRNCH;
627            goto done;
628         }
629         /* estimate degradation of the objective for up-branch */
630         d2 = psi * (ceil(beta) - beta);
631         /* determine d = max(d1, d2) */
632         d = (d1 > d2 ? d1 : d2);
633         /* choose x[j] which provides maximal estimated degradation of
634            the objective either in down- or up-branch */
635         if (dmax < d)
636         {  dmax = d;
637            jjj = j;
638            /* continue the search from a subproblem, where degradation
639               is less than in other one */
640            sel = (d1 <= d2 ? GLP_DN_BRNCH : GLP_UP_BRNCH);
641         }
642         /* display progress of pseudocost initialization */
643         if (T->parm->msg_lev >= GLP_ON)
644         {  if (xdifftime(xtime(), t) >= 10.0)
645            {  progress(T);
646               t = xtime();
647            }
648         }
649      }
650      if (dmax == 0.0)
651      {  /* no degradation is indicated; choose a variable having most
652            fractional value */
653         jjj = branch_mostf(T, &sel);
654      }
655done: *_next = sel;
656      return jjj;
657}
658
659/* eof */
Note: See TracBrowser for help on using the repository browser.