1 /*===========================================================================*/
2 /* */
3 /* This file is part of the SYMPHONY MILP Solver Framework. */
4 /* */
5 /* SYMPHONY was jointly developed by Ted Ralphs (ted@lehigh.edu) and */
6 /* Laci Ladanyi (ladanyi@us.ibm.com). */
7 /* */
8 /* (c) Copyright 2000-2019 Ted Ralphs. All Rights Reserved. */
9 /* */
10 /* This software is licensed under the Eclipse Public License. Please see */
11 /* accompanying file for terms. */
12 /* */
13 /*===========================================================================*/
14
15 #include <stdlib.h>
16 #include <memory.h>
17 #include <string.h>
18 #include <math.h>
19
20 #include "sym_lp.h"
21 #include "sym_proccomm.h"
22 #include "sym_types.h"
23 #include "sym_macros.h"
24 #include "sym_constants.h"
25
26 /*===========================================================================*/
27
28 /*===========================================================================*\
29 * This file contains LP functions related to column operations.
30 \*===========================================================================*/
31
add_col_set(lp_prob * p,our_col_set * new_cols)32 void add_col_set(lp_prob *p, our_col_set *new_cols)
33 {
34 LPdata *lp_data = p->lp_data;
35 var_desc *evar, **extra, **vars = lp_data->vars;
36
37 char *status = lp_data->status;
38
39 int new_vars = new_cols->num_vars;
40 int i, j, oldn;
41 char *where_to_move;
42
43 int cnt = 0;
44 int *index;
45 char *lu;
46 double *bd;
47
48 int to_lb_num, *to_lb_ind, to_ub_num, *to_ub_ind;
49
50 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
51
52 colind_sort_extra(p);
53
54 if (new_cols->dual_feas == NOT_TDF){
55 to_lb_num = new_cols->rel_ub;
56 to_lb_ind = new_cols->rel_ub_ind;
57 to_ub_num = new_cols->rel_lb;
58 to_ub_ind = new_cols->rel_lb_ind;
59 }else{
60 to_ub_num = new_cols->rel_ub;
61 to_ub_ind = new_cols->rel_ub_ind;
62 to_lb_num = new_cols->rel_lb;
63 to_lb_ind = new_cols->rel_lb_ind;
64 }
65
66 if (new_vars){
67 size_lp_arrays(lp_data, TRUE, FALSE, 0, new_vars, new_cols->nzcnt);
68 }
69
70 lu = lp_data->tmp.c; /* n (max(n, new_vars), but already resized) */
71 index = lp_data->tmp.i1; /* n */
72 bd = lp_data->tmp.d; /* 2*n (MAX(n, 2*new_vars), but already resized) */
73
74 if (to_ub_num > 0){
75 memset(lu, 'U', to_ub_num);
76 /* Put the branching variables and base variable to the end
77 * of the list and the extra variables to the beginning */
78 for (i = to_ub_num - 1; i >= 0; i--){
79 j = to_ub_ind[i];
80 release_var(lp_data, j, MOVE_TO_UB); /* empty function for cplex */
81 status[j] = NOT_FIXED | (status[j] & NOT_REMOVABLE);
82 bd[cnt] = vars[j]->ub;
83 index[cnt++] = j;
84 }
85 }
86
87 if (to_lb_num > 0){
88 memset(lu + cnt, 'L', to_lb_num);
89 for (i = to_lb_num - 1; i >= 0; i--){
90 j = to_lb_ind[i];
91 release_var(lp_data, j, MOVE_TO_LB); /* empty function for cplex */
92 status[j] = NOT_FIXED | (status[j] & NOT_REMOVABLE);
93 bd[cnt] = vars[j]->lb;
94 index[cnt++] = j;
95 }
96 }
97
98 if (cnt > 0)
99 change_bounds(lp_data, cnt, index, lu, bd);
100
101 if (! new_vars)
102 return;
103
104 where_to_move = lp_data->tmp.c; /* new_vars */
105 /* In the current implementation everything not in the matrix was at its
106 * lower bound (0), therefore to restore dual feasibility they have to be
107 * moved to their upper bounds, while when we just add them to the problem
108 * they have to be moved to their lower bound */
109 memset(where_to_move,
110 new_cols->dual_feas == NOT_TDF ? MOVE_TO_UB : MOVE_TO_LB, new_vars);
111
112 add_cols(lp_data, new_vars, new_cols->nzcnt, new_cols->objx,
113 new_cols->matbeg, new_cols->matind, new_cols->matval,
114 new_cols->lb, new_cols->ub, where_to_move);
115 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
116 lp_data->col_set_changed = TRUE;
117 p->colset_changed = TRUE;
118 lp_data->ordering = COLIND_ORDERED;
119
120 /* update the extra parts of vars */
121 oldn = lp_data->n - new_vars;
122 extra = lp_data->vars + oldn;
123 for (i = new_vars - 1; i >= 0; i--){
124 evar = extra[i];
125 evar->userind = new_cols->userind[i];
126 evar->colind = oldn + i;
127 evar->lb = new_cols->lb[i];
128 evar->ub = new_cols->ub[i];
129 }
130
131 /* zero out x, i.e., set it to the LB */
132 memset(lp_data->x + oldn, 0, new_vars * DSIZE);
133 /* set status of the new vars to NOT_FIXED */
134 //memset(lp_data->status + oldn, NOT_FIXED, new_vars);
135 //TODO: char vs int
136 for (i = oldn; i<oldn+new_vars; i++) {
137 lp_data->status[i] = NOT_FIXED;
138 }
139 }
140
141 /*===========================================================================*/
142
143 /*===========================================================================*\
144 * Try to tighten bounds based on reduced cost and logical fixing
145 \*===========================================================================*/
146
tighten_bounds(lp_prob * p)147 void tighten_bounds(lp_prob *p)
148 {
149 LPdata *lp_data = p->lp_data;
150 double *dj = lp_data->dj;
151 //double *x = lp_data->x;
152 char *status = lp_data->status;
153 var_desc **vars = lp_data->vars;
154 int n = lp_data->n;
155 double lpetol = lp_data->lpetol;
156
157 double gap = 0.0, max_change;
158 int i, vars_recently_fixed_to_ub = 0;
159 int did_logical_fixing = FALSE, did_reduced_cost_fixing = FALSE;
160 int lb_vars = 0, perm_lb_vars = 0, ub_vars = 0, perm_ub_vars = 0;
161 int del_vars = 0, *delstat = NULL;
162
163 //char not_fixed__lb__switch, not_fixed__ub__switch;
164 int *ind = 0;
165 char *lu = 0;
166 double *bd = 0, *ub = 0, *lb = 0;
167 int cnt = 0;
168
169 colind_sort_extra(p);
170
171 check_ub(p);
172 if (p->has_ub){
173 gap = p->ub - lp_data->objval - p->par.granularity;
174 }
175
176 /*========================================================================*\
177 * Here is the reduced cost fixing.
178 *
179 * If the gap is negative that means that we are above the limit, so don't
180 * do anything.
181 * Otherwise we do regular rc tightening if one of the following holds:
182 * -- if we have done rc fixing before then the gap must have decreased
183 * significantly
184 * -- if we haven't done rc tightening before, then the gap must be relatively
185 * small compared to the upper bound
186 \*=======================================================================*/
187 if (p->par.do_reduced_cost_fixing && p->has_ub && gap > 0){
188 if (p->last_gap == 0.0 ?
189 (gap < p->par.gap_as_ub_frac * p->ub) :
190 (gap < p->par.gap_as_last_gap_frac * p->last_gap)){
191 /* Tighten the upper/lower bounds on the variables,
192 prepare to delete them and do some statistics. */
193 delstat = lp_data->tmp.i1; /* 2*n */
194 ind = lp_data->tmp.i1 + n;
195 lu = lp_data->tmp.c; /* n */
196 bd = lp_data->tmp.d; /* n */
197
198 get_bounds(lp_data);
199 ub = lp_data->ub;
200 lb = lp_data->lb;
201
202 p->vars_deletable = 0;
203 memset((char *)delstat, 0, n * ISIZE);
204 lb_vars = perm_lb_vars = ub_vars = perm_ub_vars = 0;
205 for (cnt = 0, i = n-1; i >= 0; i--){
206 if (fabs(dj[i]) < lpetol || !vars[i]->is_int){
207 continue;
208 }
209 max_change = gap/dj[i];
210 if (max_change > 0 && max_change < ub[i] - lb[i]){
211 if (lp_data->nf_status & NF_CHECK_NOTHING){
212 status[i] ^= NOT_FIXED__PERM_LB__SWITCH;
213 perm_lb_vars++;
214 }else{
215 status[i] ^= NOT_FIXED__TEMP_LB__SWITCH;
216 lb_vars++;
217 }
218 ind[cnt] = i;
219 lu[cnt] = 'U';
220 bd[cnt++] = vars[i]->is_int ? floor(lb[i] + max_change) :
221 lb[i] + max_change;
222 vars[i]->new_ub = bd[cnt-1];
223 p->bound_changes_in_iter++;
224 if (! (status[i] & NOT_REMOVABLE) && lb[i] == 0 &&
225 lb[i] == ub[i]){
226 p->vars_deletable++;
227 delstat[i] = 1;
228 }
229 }else if (max_change < 0 && max_change > lb[i] - ub[i]){
230 if (lp_data->nf_status & NF_CHECK_NOTHING){
231 status[i] ^= NOT_FIXED__PERM_UB__SWITCH;
232 perm_ub_vars++;
233 }else{
234 status[i] ^= NOT_FIXED__TEMP_UB__SWITCH;
235 ub_vars++;
236 }
237 ind[cnt] = i;
238 lu[cnt] = 'L';
239 bd[cnt++] = vars[i]->is_int ? ceil(ub[i] + max_change) :
240 ub[i] + max_change;
241 vars[i]->new_lb = bd[cnt-1];
242 p->bound_changes_in_iter++;
243 if (! (status[i] & NOT_REMOVABLE) && lb[i] == 0 &&
244 lb[i] == ub[i]){
245 p->vars_deletable++;
246 delstat[i] = 1;
247 }
248 vars_recently_fixed_to_ub++;
249 }
250 did_reduced_cost_fixing = TRUE;
251 }
252 p->vars_recently_fixed_to_ub += vars_recently_fixed_to_ub;
253 }
254 }
255
256 #ifdef COMPILE_IN_LP
257 if (p->bc_level==0 && p->par.do_reduced_cost_fixing) {
258 /* we are root node. we will save the reduced costs after each round of
259 * cuts. whenever ub is updated, we can come back and update bounds in
260 * the root
261 */
262 if (p->tm->par.tighten_root_bounds){
263 save_root_reduced_costs(p);
264 }
265 }
266 #endif
267
268 if (cnt > 0){
269 change_bounds(lp_data, cnt, ind, lu, bd);
270 }
271
272 /*========================================================================*\
273 * Logical fixing is done only if the number of variables recently fixed
274 * to upper bound reaches a given constant AND is at least a certain
275 * fraction of the total number of variables.
276 \*=======================================================================*/
277
278 if ((p->par.do_logical_fixing) &&
279 (p->vars_recently_fixed_to_ub >
280 p->par.fixed_to_ub_before_logical_fixing) &&
281 (p->vars_recently_fixed_to_ub >
282 p->par.fixed_to_ub_frac_before_logical_fixing * n)){
283 logical_fixing_u(p);
284 did_logical_fixing = TRUE;
285 }
286
287 if (! did_reduced_cost_fixing && ! did_logical_fixing)
288 return;
289
290 if (did_reduced_cost_fixing)
291 p->last_gap = gap;
292 if (did_logical_fixing)
293 p->vars_recently_fixed_to_ub = 0;
294
295 if (p->par.verbosity > 3){
296 if (ub_vars)
297 printf("total of %i variables with temp adjusted UB ...\n",ub_vars);
298 if (perm_ub_vars)
299 printf("total of %i variables with perm adjusted UB ...\n",perm_ub_vars);
300 if (lb_vars)
301 printf("total of %i variables with temp adjusted LB ...\n",lb_vars);
302 if (perm_lb_vars)
303 printf("total of %i variables with perm adjusted LB ...\n",perm_lb_vars);
304 }
305 p->vars_at_lb = lb_vars;
306 p->vars_at_ub = ub_vars;
307
308 /* if enough variables have been fixed, then physically compress the matrix,
309 * eliminating the columns that are fixed to zero */
310 if (p->vars_deletable > p->par.mat_col_compress_num &&
311 p->vars_deletable > n * p->par.mat_col_compress_ratio){
312
313 PRINT(p->par.verbosity,3, ("Compressing constraint matrix (col) ...\n"));
314 del_vars = delete_cols(lp_data, p->vars_deletable, delstat);
315 if (del_vars > 0){
316 lp_data->lp_is_modified = LP_HAS_BEEN_MODIFIED;
317 lp_data->col_set_changed = TRUE;
318 }
319 if (p->vars_deletable > del_vars){
320 PRINT(p->par.verbosity, 3,
321 ("%i vars were not removed because they were basic ...\n",
322 p->vars_deletable - del_vars));
323 }
324 if (del_vars > 0){
325 p->vars_deletable -= del_vars;
326 PRINT(p->par.verbosity, 3,
327 ("%i vars successfully removed from the problem ...\n",
328 del_vars));
329 for (i = p->base.varnum; i < n; i++){
330 if (delstat[i] != -1){
331 *(vars[delstat[i]]) = *(vars[i]);
332 vars[delstat[i]]->colind = delstat[i];
333 }
334 }
335 }
336 }
337 }
338
339 /*===========================================================================*\
340 *===========================================================================*
341 *
342 * IMPORTANT:
343 * No matter whether this routine was called with a primal feasible tableau
344 * or not, if everything prices out, that means that if we were to add every
345 * variable right now, we would still have a dual feasible tableau, therefore
346 * the dual simplex could be continued, with the dual obj value only
347 * increasing.
348 *
349 * NOTE: The tableau we know of is always dual feasible.
350 *
351 * This routine starts by collecting ALL (known and not known)
352 * variables having reduced cost between 0 and 'gap' into
353 * 'new_cols'. As soon as a non-dual-feasible variable is encountered,
354 * the routine switches to collect the non-dual-feasibles into 'new_cols'.
355 *
356 * At the end, the result is:
357 * -- if dual_feas is TDF_HAS_ALL, then new_cols contains the set of
358 * non-fixable variables not in the matrix.
359 * -- if dual_feas is NOT_TDF, then new_cols contains all (well, at most
360 * max_ndf_vars) not dual feasible variables, which could be added.
361 * -- if dual_feas is TDF_NOT_ALL then the whole tableau is dual feas, just
362 * we have ran out of space for storing the non-fixable variables
363 *
364 \*===========================================================================*/
365
366 /*===========================================================================*\
367 * THIS FUNCTION CAN BE CALLED ONLY IF
368 * p->lp_data->nf_status[0] != NF_CHECK_NOTHING !!!
369 \*===========================================================================*/
370
price_all_vars(lp_prob * p)371 our_col_set *price_all_vars(lp_prob *p)
372 {
373 LPdata *lp_data = p->lp_data;
374 double lpetol = lp_data->lpetol;
375 int m = lp_data->m, n = lp_data->n;
376 char *status = lp_data->status;
377 double *dj = lp_data->dj;
378 double *dual = lp_data->dualsol;
379
380 int bvarnum = p->base.varnum;
381 int extranum = lp_data->n - bvarnum;
382 var_desc **vars = lp_data->vars;
383 var_desc **extra = vars + bvarnum;
384
385 int next_not_fixed, not_fixed_num = lp_data->not_fixed_num;
386 int *not_fixed = lp_data->not_fixed;
387 char new_nf_status = NF_CHECK_UNTIL_LAST;
388 int nf_status = lp_data->nf_status;
389 int tmp_not_fixed_num = 0, *tmp_not_fixed, *itmp;
390
391 int cutnum;
392 cut_data **cuts;
393 row_data *rows;
394
395 char basevar = TRUE; /* just to keep gcc quiet */
396 int ind, prevind, curind, nextind = -1; /* just to keep gcc quiet */
397 double gap, red_cost;
398
399 char must_add;
400 int dual_feas;
401 int termcode = p->lp_data->termcode;
402
403 our_col_set *new_cols = (our_col_set *) calloc(1, sizeof(our_col_set));
404 int max_ndf_vars, max_nfix_vars, max_new_vars, max_new_nzcnt;
405 int new_vars=0, new_nzcnt=0, rel_lb=0, rel_ub=0;
406
407 int collen, *colind;
408 double obj, lb, ub, *colval;
409
410 int i, j, k;
411
412 #ifdef STATISTICS
413 int nfix = 0;
414 #endif
415
416 /* Compute how many non-dual-feasible cols we are willing to add.
417 * Also compute how many non-fixable cols we are willing to add.
418 * Then, to start with, set max_new_vars to be the latter. */
419 max_ndf_vars = (int) (n * p->par.max_non_dual_feas_to_add_frac);
420 if (max_ndf_vars < p->par.max_non_dual_feas_to_add_min)
421 max_ndf_vars = p->par.max_non_dual_feas_to_add_min;
422 if (max_ndf_vars > p->par.max_non_dual_feas_to_add_max)
423 max_ndf_vars = p->par.max_non_dual_feas_to_add_max;
424
425 if (termcode != LP_D_UNBOUNDED){
426 max_nfix_vars = 0;
427 }else{
428 max_nfix_vars = (int) (n * p->par.max_not_fixable_to_add_frac);
429 if (max_nfix_vars < p->par.max_not_fixable_to_add_min)
430 max_nfix_vars = p->par.max_not_fixable_to_add_min;
431 if (max_nfix_vars > p->par.max_not_fixable_to_add_max)
432 max_nfix_vars = p->par.max_not_fixable_to_add_max;
433 }
434
435 tmp_not_fixed = (int *) malloc(p->par.not_fixed_storage_size * ISIZE);
436
437 max_new_vars = MAX(max_nfix_vars, max_ndf_vars);
438 max_new_nzcnt = m*max_new_vars;
439
440 new_cols->rel_lb_ind = p->vars_at_lb ?
441 (int *) malloc(p->vars_at_lb * ISIZE) : NULL;
442 new_cols->rel_ub_ind = p->vars_at_ub ?
443 (int *) malloc(p->vars_at_ub * ISIZE) : NULL;
444 new_cols->objx = (double *) malloc(max_new_vars * DSIZE);
445 new_cols->lb = (double *) malloc(max_new_vars * DSIZE);
446 new_cols->ub = (double *) malloc(max_new_vars * DSIZE);
447 new_cols->matbeg = (int *) malloc((max_new_vars+1) * ISIZE);
448 new_cols->matbeg[0] = 0;
449 new_cols->matind = (int *) malloc(max_new_nzcnt * ISIZE);
450 new_cols->matval = (double *) malloc(max_new_nzcnt * DSIZE);
451 new_cols->userind = (int *) malloc(max_new_vars * ISIZE);
452
453 userind_sort_extra(p);
454
455 /* Collect the non-base lpcuts */
456 cutnum = m - p->base.cutnum;
457 cuts = (cut_data **) lp_data->tmp.p1; /* m (actually, cutnum < m) */
458 rows = lp_data->rows + p->base.cutnum;
459 for (i = cutnum - 1; i >= 0; i--)
460 cuts[i] = rows[i].cut;
461
462 colind = lp_data->tmp.i1; /* m */
463 colval = lp_data->tmp.d; /* 2*m */
464
465 must_add = FALSE;
466 dual_feas = TDF_HAS_ALL;
467 check_ub(p);
468 gap = p->has_ub ? p->ub - p->par.granularity - lp_data->objval :
469 SYM_INFINITY;
470
471 /*========================================================================*\
472 * Now we loop through every single variable, and get those to be
473 * added into the new_cols structure.
474 *
475 * In the loop 'i' runs on the base vars, 'j' on the extra vars and 'k'
476 * on the not_fixed ones.
477 * In every iteration we compute nextind, the next variable we know of,
478 * based on i and j.
479 * -- If we have run out of non-fixables and
480 * - nf_status == NF_CHECK_UNTIL_LAST then, we know of all non-fixable
481 * non-base variable and we are past them, then only base vars can be
482 * left to be checked. And nextind is processed.
483 * - nf_status == NF_CHECK_AFTER_LAST then there are more non-fixables,
484 * but we don't know what they are, then the user has to tell about the
485 * next variable.
486 * -- If nextind is smaller than the next non-fixable we know of then
487 * nextind is processed.
488 * -- If it is greater then the next non-fixable has to be processed.
489 \*=======================================================================*/
490
491 curind = prevind = -1;
492
493 for (i = 0, j = 0, k = 0; TRUE; prevind = curind){
494 switch ((i < bvarnum ? 1 : 0) + (j < extranum ? 2 : 0)){
495 case 0: /* none left */
496 nextind = -1; basevar = FALSE; break;
497 case 1: /* only base vars left */
498 nextind = vars[i]->userind; basevar = TRUE; break;
499 case 2: /* only extra vars left */
500 nextind = extra[j]->userind; basevar = FALSE; break;
501 case 3: /* both of them left */
502 if (vars[i]->userind < extra[j]->userind){
503 nextind = vars[i]->userind; basevar = TRUE;
504 }else{
505 nextind = extra[j]->userind; basevar = FALSE;
506 }
507 break;
508 }
509
510 /*=====================================================================*\
511 * If we have a chance to prove TDF or
512 * If we proved NOT_TDF but still have enough space to add new cols,
513 * then the user the user generates the next col (or says that the next
514 * col is what we have next), otherwise the next col is what we have
515 * next.
516 \*=====================================================================*/
517
518 if ((dual_feas != NOT_TDF) ||
519 (dual_feas == NOT_TDF && new_vars < max_ndf_vars)){
520 if (k < not_fixed_num){
521 /* If anything is left on the not_fixed list then compare it to
522 * the next in matrix (nextind) and get the smaller one */
523 next_not_fixed = not_fixed[k];
524 if (nextind == -1 || nextind > next_not_fixed){
525 /* FIXME: Catch the error message for this function */
526 curind = generate_column_u(p, cutnum, cuts,
527 prevind, next_not_fixed,
528 GENERATE_NEXTIND,
529 colval, colind, &collen, &obj,
530 &lb, &ub);
531 k++;
532 }else{
533 if (nextind == next_not_fixed)
534 k++;
535 curind = nextind;
536 }
537 }else{
538 /*If nothing is left in not_fixed then things depend on nf_status*/
539 if (nf_status == NF_CHECK_UNTIL_LAST){
540 curind = nextind;
541 }else{ /* NF_CHECK_AFTER_LAST */
542 /* FIXME: Catch the error message for this function */
543 curind = generate_column_u(p, cutnum, cuts,
544 prevind, nextind,
545 GENERATE_REAL_NEXTIND,
546 colval, colind, &collen, &obj,
547 &lb, &ub);
548 }
549 }
550 }else{
551 curind = nextind;
552 }
553
554 /* Now curind is the one we work with. If it's negative then there are
555 * no more cols. */
556 if (curind < 0)
557 break;
558
559 if (curind == nextind){
560 /* The next col is what we have next, i.e., it is in the matrix.
561 * We've got to check it only if it is fixed to LB or UB. */
562 ind = basevar ? i : j + bvarnum;
563 red_cost = dj[ind];
564 if (status[ind] & TEMP_FIXED_TO_LB || status[ind] & TEMP_FIXED_TO_UB){
565 if (status[ind] & TEMP_FIXED_TO_LB){
566 if (red_cost < -lpetol){
567 if (dual_feas != NOT_TDF){
568 dual_feas = NOT_TDF;
569 rel_lb = rel_ub = new_vars = new_nzcnt = 0;
570 }
571 new_cols->rel_lb_ind[rel_lb++] = ind;
572 }else if (dual_feas != NOT_TDF){
573 if (red_cost < gap){
574 new_cols->rel_lb_ind[rel_lb++] = ind;
575 }else{
576 new_cols->rel_lb_ind[rel_lb++] = - ind - 1;
577 }
578 }
579 }else{ /* TEMP_FIXED_TO_UB */
580 if (red_cost > lpetol){
581 if (dual_feas != NOT_TDF){
582 dual_feas = NOT_TDF;
583 rel_lb = rel_ub = new_vars = new_nzcnt = 0;
584 }
585 new_cols->rel_ub_ind[rel_ub++] = ind;
586 }else if (dual_feas != NOT_TDF){
587 if (red_cost > -gap){
588 new_cols->rel_ub_ind[rel_ub++] = ind;
589 }else{
590 new_cols->rel_ub_ind[rel_ub++] = - ind - 1;
591 }
592 }
593 }
594 }
595 if (basevar)
596 i++;
597 else
598 j++;
599
600 }else{ /* the next col is not in the matrix */
601 red_cost = obj - dot_product(colval, colind, collen, dual);
602 if (red_cost < -lpetol){
603 if (dual_feas != NOT_TDF){
604 dual_feas = NOT_TDF;
605 rel_lb = rel_ub = new_vars = new_nzcnt = 0;
606 }
607 must_add = TRUE;
608 }else if (dual_feas != NOT_TDF && red_cost < gap){
609 if (new_vars == max_nfix_vars){
610 /* Run out of space!! */
611 dual_feas = TDF_NOT_ALL;
612 #ifdef STATISTICS
613 if (!nfix)
614 nfix = new_vars;
615 nfix++;
616 #endif
617 }else{
618 must_add = TRUE;
619 }
620 }
621 if (must_add){
622 new_cols->objx[new_vars] = obj;
623 new_cols->lb[new_vars] = lb;
624 new_cols->ub[new_vars] = ub;
625 new_cols->matbeg[new_vars + 1] = new_cols->matbeg[new_vars]+collen;
626 memcpy((char *)(new_cols->matind+new_cols->matbeg[new_vars]),
627 (char *)colind, collen * ISIZE);
628 memcpy((char *)(new_cols->matval+new_cols->matbeg[new_vars]),
629 (char *)colval, collen * DSIZE);
630 new_nzcnt += collen;
631 new_cols->userind[new_vars++] = curind;
632 must_add = FALSE;
633 }
634 basevar = FALSE;
635 }
636 /* Add this variable to the not_fixed list if it cannot be permanently
637 fixed*/
638 if (red_cost > -gap && red_cost < gap && !basevar){
639 if (tmp_not_fixed_num < p->par.not_fixed_storage_size){
640 tmp_not_fixed[tmp_not_fixed_num++] = curind;
641 }else{
642 new_nf_status = NF_CHECK_AFTER_LAST;
643 }
644 }
645 }
646
647 new_cols->num_vars = new_vars;
648 new_cols->nzcnt = new_nzcnt;
649 new_cols->rel_lb = rel_lb;
650 new_cols->rel_ub = rel_ub;
651
652 switch (new_cols->dual_feas = dual_feas){
653 case NOT_TDF:
654 PRINT(p->par.verbosity, 5, ("price_all_vars() : NOT_TDF. [ %i ]\n",
655 new_vars + rel_ub + rel_lb));
656 p->vars_at_lb -= rel_lb;
657 p->vars_at_ub -= rel_ub;
658 for (i = rel_lb - 1; i >= 0; i--){
659 if (! (status[new_cols->rel_lb_ind[i]] & NOT_REMOVABLE))
660 p->vars_deletable--;
661 }
662 break;
663
664 case TDF_NOT_ALL:
665 PRINT(p->par.verbosity, 5, ("price_all_vars() : TDF_NOT_ALL.\n"));
666 #ifdef STATISTICS
667 PRINT(p->par.verbosity, 5, (" ( nonfix / maxnonfix : %i / %i )\n",
668 nfix, max_nfix_vars));
669 #endif
670
671 /*=====================================================================*\
672 * If total dual feasibility has been proved but there are too
673 * many vars to add, then (although we don't want to release
674 * those already fixed but not permanently fixable) we can
675 * permanently fix those which are marked such way.
676 \*=====================================================================*/
677
678 for (k = 0, i = 0; i < rel_lb; i++){
679 if ((j = new_cols->rel_lb_ind[i]) < 0){
680 j = -j-1;
681 status[j] ^= TEMP_PERM_LB__SWITCH;
682 }else{
683 new_cols->rel_lb_ind[k++] = j;
684 }
685 if (! (status[j] & NOT_REMOVABLE))
686 p->vars_deletable--;
687 }
688 new_cols->rel_lb = k;
689 for (k = 0, i = 0; i < rel_ub; i++){
690 if ((j = new_cols->rel_ub_ind[i]) < 0){
691 status[-j-1] ^= TEMP_PERM_UB__SWITCH;
692 }else{
693 new_cols->rel_ub_ind[k++] = j;
694 }
695 }
696 new_cols->rel_ub = k;
697
698 /*=====================================================================*\
699 * If we come here, there should only be variables to add in the case
700 * the LP is infeasible. Since we can't add all the variables in, we may
701 * as well just wait for the call to restore_feasibility and only add in
702 * the one (if there is one) that destroys the proof of infeasibility.
703 \*=====================================================================*/
704
705 new_vars = rel_ub = rel_lb = 0;
706
707 /* Update which variables have to be priced out. */
708 lp_data->nf_status = new_nf_status;
709 lp_data->not_fixed_num = tmp_not_fixed_num;
710 itmp = lp_data->not_fixed;
711 lp_data->not_fixed = tmp_not_fixed;
712 tmp_not_fixed = itmp;
713
714 break;
715
716 case TDF_HAS_ALL:
717 /* Now if total dual feasibility is proved and there aren't too
718 * many extra variables to be added... */
719 lp_data->nf_status = NF_CHECK_NOTHING;
720 lp_data->not_fixed_num = 0;
721 for (k = 0, i = 0; i < rel_lb; i++){
722 if ((j = new_cols->rel_lb_ind[i]) < 0){
723 j = -j-1;
724 status[j] ^= TEMP_PERM_LB__SWITCH;
725 }else{
726 new_cols->rel_lb_ind[k++] = j;
727 }
728 if (! (status[j] & NOT_REMOVABLE))
729 p->vars_deletable--;
730 }
731 new_cols->rel_lb = rel_lb = k;
732 for (k = 0, i = 0; i < rel_ub; i++){
733 if ((j = new_cols->rel_ub_ind[i]) < 0){
734 status[-j-1] ^= TEMP_PERM_UB__SWITCH;
735 }else{
736 new_cols->rel_ub_ind[k++] = j;
737 }
738 }
739 new_cols->rel_ub = rel_ub = k;
740
741 p->vars_at_lb = 0; /* they are either permanently fixed or released */
742 p->vars_at_ub = 0; /* they are either permanently fixed or released */
743
744 PRINT(p->par.verbosity, 5, ("price_all_vars() : TDF_HAS_ALL. [ %i ]\n",
745 new_vars + rel_ub + rel_lb));
746 break;
747 }
748
749 FREE(tmp_not_fixed);
750
751 if (rel_lb || rel_ub)
752 PRINT(p->par.verbosity, 1,
753 ("Released %i 0-variables and %i 1-variables.\n", rel_lb, rel_ub));
754
755 if (new_vars || rel_lb || rel_ub)
756 add_col_set(p, new_cols);
757
758 return(new_cols);
759 }
760
761 /*===========================================================================*/
762
763 /*===========================================================================*\
764 * This function is called with an infeasible problem in lp_data, and after
765 * we found that the lp is total dual feasible but there were too many
766 * variables to add (i.e., TDF_NOT_ALL).
767 * Therefore in new_cols there will be a bunch of vars wich are not permanently
768 * fixable so we should start to check those first.
769 *
770 * The final goal here is to find a variable which destroys a proof of primal
771 * infeasibility
772 \*===========================================================================*/
773
restore_lp_feasibility(lp_prob * p,our_col_set * new_cols)774 int restore_lp_feasibility(lp_prob *p, our_col_set *new_cols)
775 {
776 LPdata *lp_data = p->lp_data;
777 double lpetol = lp_data->lpetol;
778 char *status = lp_data->status;
779 double *dual = lp_data->dualsol;
780
781 int bvarnum = p->base.varnum;
782 int extranum = lp_data->n - bvarnum;
783 var_desc **vars = lp_data->vars;
784 var_desc **extra = vars + bvarnum;
785
786 int next_not_fixed, *not_fixed = lp_data->not_fixed;
787 int nf_status = lp_data->nf_status;
788 int not_fixed_num = lp_data->not_fixed_num;
789
790 int cutnum;
791 cut_data **cuts;
792
793 int infind, violation;
794
795 int collen, *colind;
796 double obj, lb = 0, ub = 0, *colval, *binvrow;
797
798 double gap, red_cost, prod;
799
800 char basevar = TRUE; /* just to keep gcc quiet */
801 int prevind, curind, nextind = -1; /* just to keep gcc quiet */
802
803 int i, j, k;
804
805 /* Get a proof of infeasibility and get the corresponding row of the basis
806 * inverse */
807
808 violation = get_proof_of_infeas(lp_data, &infind);
809
810 /*========================================================================*\
811 * Collect the non-base lpcuts. We would have to do the same as in
812 * price_all_vars(), but this function is called right after that, and thus
813 * lp_data->tmp.p1 still has the pointers to the cuts :-).
814 * And, price_all_vars did NOT resize the tmp arrays as it has failed if
815 * we came to this function.
816 * Also itmpm and dtmpm were reallocated there. (for big enough)
817 \*========================================================================*/
818
819 cutnum = lp_data->m - p->base.cutnum;
820 cuts = (cut_data **) lp_data->tmp.p1;
821 colind = lp_data->tmp.i1;
822 colval = lp_data->tmp.d;
823 binvrow = lp_data->tmp.d + lp_data->m;
824
825 get_binvrow(lp_data, infind, binvrow);
826
827 check_ub(p);
828 gap = p->has_ub ? p->ub - p->par.granularity - lp_data->objval :
829 SYM_INFINITY;
830
831 /* First check those released from their lower bound in price_all_vars(),
832 * and see if they destroy the proof of infeas */
833 for (i=new_cols->rel_lb-1; i>=0; i--){
834 j = new_cols->rel_lb_ind[i];
835 get_column(lp_data, j, colval, colind, &collen, &obj);
836 prod = dot_product(colval, colind, collen, binvrow);
837 if ((violation == LOWER_THAN_LB && prod < -lpetol) ||
838 (violation == HIGHER_THAN_UB && prod > lpetol)){
839 /* OK, just release this variable */
840 PRINT(p->par.verbosity, 2,
841 ("RELEASED_FROM_LB while restoring feasibility.\n"));
842 new_cols->num_vars = new_cols->rel_ub = new_cols->rel_lb = 0;
843 change_ub(lp_data, j, lp_data->vars[j]->ub);
844 status[j] ^= NOT_FIXED__TEMP_LB__SWITCH;
845 release_var(lp_data, j, MOVE_TO_LB);
846 return(TRUE);
847 }
848 }
849 new_cols->rel_lb = 0; /*We don't need these anymore*/
850
851 /* Now check those released from their upper bound */
852 for (i=new_cols->rel_ub-1; i>=0; i--){
853 j = new_cols->rel_ub_ind[i];
854 get_column(lp_data, j, colval, colind, &collen, &obj);
855 prod = dot_product(colval, colind, collen, binvrow);
856 if ((violation == LOWER_THAN_LB && prod > lpetol) ||
857 (violation == HIGHER_THAN_UB && prod < -lpetol)){
858 /* OK, just release this variable */
859 PRINT(p->par.verbosity, 2,
860 ("RELEASED_FROM_UB while restoring feasibility.\n"));
861 new_cols->num_vars = new_cols->rel_ub = new_cols->rel_lb = 0;
862 change_lb(lp_data, j, lp_data->vars[j]->lb);
863 status[j] ^= NOT_FIXED__TEMP_UB__SWITCH;
864 release_var(lp_data, j, MOVE_TO_UB);
865 return(TRUE);
866 }
867 }
868 new_cols->rel_ub = 0; /*We don't need these anymore*/
869
870 /* Now check the ones described in the new_vars part of new_cols.
871 * These are either already added, or we got that far */
872 for (i=0; i<new_cols->num_vars; i++){
873 prod = dot_product(new_cols->matval + new_cols->matbeg[i],
874 new_cols->matind + new_cols->matbeg[i],
875 new_cols->matbeg[i+1] - new_cols->matbeg[i],
876 binvrow);
877 if ((violation == LOWER_THAN_LB && prod < -lpetol) ||
878 (violation == HIGHER_THAN_UB && prod > lpetol)){
879 PRINT(p->par.verbosity, 2,
880 ("1 variable added while restoring feasibility.\n"));
881 new_cols->rel_ub = new_cols->rel_lb = 0;
882 new_cols->num_vars = 1;
883 if (i > 0){
884 new_cols->userind[0] = new_cols->userind[i];
885 new_cols->objx[0] = new_cols->objx[i];
886 new_cols->lb[0] = lb;
887 new_cols->ub[0] = ub;
888 memmove(new_cols->matind, new_cols->matind + new_cols->matbeg[i],
889 new_cols->nzcnt * ISIZE);
890 memmove(new_cols->matval, new_cols->matval + new_cols->matbeg[i],
891 new_cols->nzcnt * DSIZE);
892 new_cols->matbeg[1] = new_cols->nzcnt;
893 }
894 new_cols->nzcnt = new_cols->matbeg[i+1] - new_cols->matbeg[i];
895 add_col_set(p, new_cols);
896 return(TRUE);
897 }
898 }
899
900 /* OK, we are out of the if, so we still didn't get rid of the proof. */
901
902 userind_sort_extra(p);
903
904 /* Just to avoid copying */
905 colind = new_cols->matind;
906 colval = new_cols->matval;
907
908 /*========================================================================*\
909 * Go through all the columns not in the matrix starting from prevind,
910 * i.e., we start where price_all_vars gave up collecting not fixables.
911 * Those in the matrix are already tested; they were listed in
912 * rel_{lb,ub}_ind.
913 * To do this first get the right i,j,k for a loop awfully similar to the
914 * one in price_all_vars.
915 \*========================================================================*/
916
917 prevind = new_cols->userind[new_cols->num_vars-1];
918 i = bvarnum > 0 ? bfind(prevind, p->base.userind, bvarnum) + 1 : 0;
919 for (j = extranum - 1; j >= 0 && extra[j]->userind > prevind; j--); j++;
920 k = not_fixed_num > 0 ? bfind(prevind, not_fixed, not_fixed_num) + 1 : 0;
921
922 for (; ; prevind = curind){
923 if (k == not_fixed_num && nf_status != NF_CHECK_AFTER_LAST)
924 /* nothing can help now... */
925 break;
926 switch ((i < bvarnum ? 1 : 0) + (j < extranum ? 2 : 0)){
927 case 0: /* none left */
928 nextind = -1; break;
929 case 1: /* only base vars left */
930 nextind = vars[i]->userind; basevar = TRUE; break;
931 case 2: /* only extra vars left */
932 nextind = extra[j]->userind; basevar = FALSE; break;
933 case 3: /* both of them left */
934 if (vars[i]->userind < extra[j]->userind){
935 nextind = vars[i]->userind; basevar = TRUE;
936 }else{
937 nextind = extra[j]->userind; basevar = FALSE;
938 }
939 break;
940 }
941 if (k < not_fixed_num){
942 next_not_fixed = not_fixed[k];
943 if (nextind == -1 || nextind > next_not_fixed){
944 /* FIXME: Catch the error message for this function */
945 curind = generate_column_u(p, cutnum, cuts,
946 prevind, next_not_fixed,
947 GENERATE_NEXTIND, colval, colind,
948 &collen, &obj, &lb, &ub);
949 k++;
950 }else{
951 if (nextind == next_not_fixed)
952 k++;
953 curind = nextind;
954 }
955 }else{ /* no we know that NF_CHECK_AFTER_LAST */
956 /* FIXME: Catch the error message for this function */
957 curind = generate_column_u(p, cutnum, cuts,
958 prevind, nextind, GENERATE_REAL_NEXTIND,
959 colval, colind, &collen, &obj, &lb, &ub);
960 }
961
962 /* Now curind is the one we work with. If it's negative then there are
963 * no more cols. */
964 if (curind < 0)
965 break;
966
967 if (curind == nextind){
968 /* no point in testing curind: it is either in the matrix or has
969 * already been tested */
970 if (basevar)
971 i++;
972 else
973 j++;
974 }else{
975 prod = dot_product(colval, colind, collen, binvrow);
976 if ((violation == LOWER_THAN_LB && prod < -lpetol) ||
977 (violation == HIGHER_THAN_UB && prod > lpetol)){
978 red_cost = obj - dot_product(colval, colind, collen, dual);
979 if (red_cost < gap){ /* It is at 0 level anyway and dual feas*/
980 PRINT(p->par.verbosity, 2,
981 ("1 variable added while restoring feasibility.\n"));
982 new_cols->num_vars = 1;
983 new_cols->userind[0] = curind;
984 new_cols->objx[0] = obj;
985 new_cols->matbeg[1] = collen;
986 new_cols->nzcnt = collen;
987 add_col_set(p, new_cols);
988 return(TRUE);
989 }
990 }
991 }
992 }
993
994 /* We came out of the loop ==> primal feasibility cannot be restored */
995 return(FALSE);
996 }
997
998 /*===========================================================================*/
999
userind_sort_extra(lp_prob * p)1000 void userind_sort_extra(lp_prob *p)
1001 {
1002 LPdata *lp_data = p->lp_data;
1003 if (lp_data->n > p->base.varnum + 1){
1004 if (lp_data->ordering == COLIND_ORDERED){
1005 qsort((char *)(lp_data->vars + p->base.varnum),
1006 lp_data->n - p->base.varnum,
1007 sizeof(var_desc *), var_uind_comp);
1008 lp_data->ordering = USERIND_ORDERED;
1009 }
1010 }else{
1011 lp_data->ordering = COLIND_AND_USERIND_ORDERED;
1012 }
1013 }
1014
1015 /*===========================================================================*/
1016
colind_sort_extra(lp_prob * p)1017 void colind_sort_extra(lp_prob *p)
1018 {
1019 LPdata *lp_data = p->lp_data;
1020 if (lp_data->n > p->base.varnum + 1){
1021 if (lp_data->ordering == USERIND_ORDERED){
1022 qsort((char *)(lp_data->vars + p->base.varnum),
1023 lp_data->n - p->base.varnum,
1024 sizeof(var_desc *), var_cind_comp);
1025 lp_data->ordering = COLIND_ORDERED;
1026 }
1027 }else{
1028 lp_data->ordering = COLIND_AND_USERIND_ORDERED;
1029 }
1030 }
1031
1032 /*===========================================================================*/
1033
var_uind_comp(const void * v0,const void * v1)1034 int var_uind_comp(const void *v0, const void *v1)
1035 {
1036 return((*(var_desc **)v0)->userind - (*(var_desc **)v1)->userind);
1037 }
1038
1039 /*===========================================================================*/
1040
var_cind_comp(const void * v0,const void * v1)1041 int var_cind_comp(const void *v0, const void *v1)
1042 {
1043 return((*(var_desc **)v0)->colind - (*(var_desc **)v1)->colind);
1044 }
1045
1046
1047 /*===========================================================================*/
1048 #ifdef COMPILE_IN_LP
save_root_reduced_costs(lp_prob * p)1049 int save_root_reduced_costs(lp_prob *p)
1050 {
1051 int *indices;
1052 double *values, *lb, *ub;
1053 tm_prob *tm = p->tm;
1054 rc_desc *rc = NULL;
1055 int pos;
1056 int *tind = p->lp_data->tmp.i1;
1057 int cnt = 0, i, j;
1058 int n = p->lp_data->n;
1059 var_desc **vars = p->lp_data->vars;
1060 double lpetol = p->lp_data->lpetol;
1061 double *lp_lb, *lp_ub;
1062 double *dj = p->lp_data->dj;
1063
1064 get_bounds(p->lp_data);
1065 lp_lb = p->lp_data->lb;
1066 lp_ub = p->lp_data->ub;
1067 for (i = 0; i < n; i++){
1068 if (vars[i]->is_int && lp_ub[i]-lp_lb[i]>lpetol &&
1069 (dj[i] > lpetol || dj[i] < -lpetol)){
1070 tind[cnt] = i;
1071 cnt++;
1072 }
1073 }
1074 PRINT(p->par.verbosity, 5, ("there are %d non zero reduced costs for "
1075 "integer vars\n", cnt));
1076
1077 if (cnt==0) {
1078 return 0;
1079 }
1080 indices = (int *)malloc(cnt*ISIZE);
1081 values = (double *)malloc(cnt*DSIZE);
1082 lb = (double *)malloc(cnt*DSIZE);
1083 ub = (double *)malloc(cnt*DSIZE);
1084
1085 for (i = 0; i < cnt; i++){
1086 j = tind[i];
1087 indices[i] = vars[j]->userind;
1088 values[i] = dj[j];
1089 lb[i] = lp_lb[j];
1090 ub[i] = lp_ub[j];
1091 }
1092 /*
1093 for (int i=0;i<cnt;i++) {
1094 printf("var %d, %20.10f\n",indices[i],values[i]);
1095 }
1096 printf("\n\n");
1097 */
1098
1099 if (!tm->reduced_costs) {
1100 tm->reduced_costs = (rc_desc *) malloc(sizeof(rc_desc));
1101 rc = tm->reduced_costs;
1102 rc->size = 10;
1103 rc->num_rcs = 0;
1104 rc->indices = (int **)calloc(rc->size,sizeof(int *));
1105 rc->values = (double **)calloc(rc->size,sizeof(double *));
1106 rc->lb = (double **)calloc(rc->size,sizeof(double *));
1107 rc->ub = (double **)calloc(rc->size,sizeof(double *));
1108 rc->obj = (double *)malloc(rc->size*DSIZE);
1109 rc->cnt = (int *)calloc(rc->size,ISIZE);
1110 } else {
1111 rc = tm->reduced_costs;
1112 }
1113
1114 pos = rc->num_rcs%rc->size;
1115 if (rc->size==rc->num_rcs) {
1116 /* replace the oldest one with the new one */
1117 FREE(rc->indices[pos]);
1118 FREE(rc->values[pos]);
1119 FREE(rc->lb[pos]);
1120 FREE(rc->ub[pos]);
1121 }
1122 rc->indices[pos] = indices;
1123 rc->values[pos] = values;
1124 rc->lb[pos] = lb;
1125 rc->ub[pos] = ub;
1126 rc->cnt[pos] = cnt;
1127 rc->obj[pos] = p->lp_data->objval;
1128 if (rc->num_rcs < rc->size) {
1129 rc->num_rcs++;
1130 }
1131 return 0;
1132 }
1133
1134 /*===========================================================================*/
1135 #if 0
1136 int tighten_root_bounds(lp_prob *p)
1137 {
1138 /*
1139 * using the reduced costs that are saved from the root node, try to
1140 * improve variable bounds.
1141 * should be called whenever ub is updated.
1142 * change only bounds for root. not for the current node. the bounds for
1143 * current node are updated in tighten_bounds()
1144 */
1145 int i, j, k, l;
1146 rc_desc *rc = p->tm->reduced_costs;
1147 double gap, max_change;
1148 double *dj, *lb, *ub;
1149 int *saved_ind;
1150 int *ind = p->lp_data->tmp.i1;
1151 double *bd = p->lp_data->tmp.d;
1152 char *lu = p->lp_data->tmp.c;
1153 int cnt, total_changes = 0;
1154 double lpetol = p->lp_data->lpetol;
1155 bounds_change_desc *bnd_change;
1156 int *new_ind;
1157 int num_new_bounds;
1158 int verbosity = p->par.verbosity;
1159 int *oldindex;
1160 double *oldvalue;
1161 char *oldlu;
1162
1163 if (!rc) {
1164 return 0;
1165 }
1166
1167 if (!p->has_ub) {
1168 PRINT(verbosity, -1, ("tighten_root_bounds: cant tighten bounds if ub "
1169 "does not exist!\n"));
1170 return 0;
1171 }
1172
1173 new_ind = (int *)malloc(p->mip->n*ISIZE);
1174 for (i=0; i<rc->num_rcs;i++) {
1175 gap = p->ub - rc->obj[i] - p->par.granularity;
1176 if (gap<=lpetol) {
1177 continue;
1178 }
1179 saved_ind = rc->indices[i];
1180 dj = rc->values[i];
1181 lb = rc->lb[i];
1182 ub = rc->ub[i];
1183 cnt = 0;
1184 for (j=0; j<rc->cnt[i]; j++) {
1185 max_change = gap/dj[j];
1186 if (max_change > 0 && max_change < ub[j]-lb[j]){
1187 ind[cnt] = saved_ind[j];
1188 lu[cnt] = 'U';
1189 bd[cnt++] = floor(lb[j] + max_change);
1190 }else if (max_change < 0 && max_change > lb[j] - ub[j]){
1191 ind[cnt] = saved_ind[j];
1192 lu[cnt] = 'L';
1193 bd[cnt++] = ceil(ub[j] + max_change);
1194 }
1195 }
1196 PRINT(verbosity, 5, ("tighten_root_bounds: at node %d, tightening %d "
1197 "bounds in root\n", p->bc_index, cnt));
1198 if (cnt == 0) {
1199 continue;
1200 }
1201 /* add these changes to root node */
1202 if (p->tm->rootnode->desc.bnd_change) {
1203 bnd_change = p->tm->rootnode->desc.bnd_change;
1204 } else {
1205 p->tm->rootnode->desc.bnd_change = bnd_change =
1206 (bounds_change_desc *)calloc(1,sizeof(bounds_change_desc));
1207 }
1208 if (bnd_change->num_changes>0) {
1209 /*
1210 * update existing changes and store the new ones in a separate array
1211 */
1212 num_new_bounds=0;
1213 oldvalue = bnd_change->value;
1214 oldindex = bnd_change->index;
1215 oldlu = bnd_change->lbub;
1216 for (k=0; k<cnt; k++) {
1217 for (j=0; j<bnd_change->num_changes; j++) {
1218 if (oldindex[j]==ind[k] && oldlu[j]==lu[k]){
1219 if (lu[k]=='L' && oldvalue[j]<bd[k]) {
1220 oldvalue[j]=bd[k];
1221 total_changes++;
1222 } else if (lu[k]=='U' && oldvalue[j]>bd[k]) {
1223 oldvalue[j]=bd[k];
1224 total_changes++;
1225 }
1226 break;
1227 }
1228 }
1229 if (j>=bnd_change->num_changes) {
1230 new_ind[num_new_bounds] = k;
1231 num_new_bounds++;
1232 }
1233 }
1234 /* those changes that dint already have an entry and stored now */
1235 if (num_new_bounds) {
1236 int new_cnt = num_new_bounds+bnd_change->num_changes;
1237 bnd_change->index = (int *)realloc(bnd_change->index,
1238 ISIZE*new_cnt);
1239 bnd_change->lbub = (char *)realloc(bnd_change->lbub,
1240 CSIZE*new_cnt);
1241 bnd_change->value = (double *)realloc(bnd_change->value,
1242 DSIZE*new_cnt);
1243 oldvalue = bnd_change->value;
1244 oldindex = bnd_change->index;
1245 oldlu = bnd_change->lbub;
1246 l = bnd_change->num_changes;
1247 for (j=0; j<num_new_bounds; j++) {
1248 total_changes++;
1249 k = new_ind[j];
1250 oldindex[l] = ind[k];
1251 oldlu[l] = lu[k];
1252 oldvalue[l] = bd[k];
1253 bnd_change->num_changes++;
1254 l++;
1255 }
1256 }
1257 } else {
1258 bnd_change->index = (int *)malloc(cnt*ISIZE);
1259 bnd_change->lbub = (char *)malloc(cnt*CSIZE);
1260 bnd_change->value = (double *)malloc(cnt*DSIZE);
1261 bnd_change->index = (int *) memcpy(bnd_change->index, ind, ISIZE*cnt);
1262 bnd_change->lbub = (char *) memcpy(bnd_change->lbub, lu, CSIZE*cnt);
1263 bnd_change->value = (double *) memcpy(bnd_change->value, bd,
1264 DSIZE*cnt);
1265 bnd_change->num_changes = cnt;
1266 }
1267 }
1268 if (verbosity>5 && p->tm->rootnode->desc.bnd_change!=NULL) {
1269 printf("tighten_root_bounds: root now has %d changes\n",
1270 p->tm->rootnode->desc.bnd_change->num_changes);
1271 }
1272 FREE(new_ind);
1273 return 0;
1274 }
1275 #endif
1276 #endif
1277 /*===========================================================================*/
1278 /*===========================================================================*/
1279