lemon-project-template-glpk
view deps/glpk/src/glpios09.c @ 9:33de93886c88
Import GLPK 4.47
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 20:59:10 +0100 |
| parents | |
| children |
line source
1 /* glpios09.c (branching heuristics) */
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, 2011 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 ***********************************************************************/
25 #include "glpios.h"
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. */
47 static int branch_first(glp_tree *T, int *next);
48 static int branch_last(glp_tree *T, int *next);
49 static int branch_mostf(glp_tree *T, int *next);
50 static int branch_drtom(glp_tree *T, int *next);
52 int 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 }
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. */
89 static 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 }
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. */
116 static 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 }
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...".) */
146 static 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 }
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. */
193 static 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;
316 skip: /* 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 }
383 done: *_next = next;
384 return jj;
385 }
387 /**********************************************************************/
389 struct 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 };
407 void *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 }
423 static 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 }
485 void 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 }
524 skip: return;
525 }
527 void 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 }
540 static 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);
579 done: return psi;
580 }
582 static 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 }
597 int 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 }
655 done: *_next = sel;
656 return jjj;
657 }
659 /* eof */