lemon-project-template-glpk

comparison deps/glpk/src/glpios06.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:92ea9d295659
1 /* glpios06.c (MIR cut generator) */
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, 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 ***********************************************************************/
24
25 #include "glpios.h"
26
27 #define _MIR_DEBUG 0
28
29 #define MAXAGGR 5
30 /* maximal number of rows which can be aggregated */
31
32 struct MIR
33 { /* MIR cut generator working area */
34 /*--------------------------------------------------------------*/
35 /* global information valid for the root subproblem */
36 int m;
37 /* number of rows (in the root subproblem) */
38 int n;
39 /* number of columns */
40 char *skip; /* char skip[1+m]; */
41 /* skip[i], 1 <= i <= m, is a flag that means that row i should
42 not be used because (1) it is not suitable, or (2) because it
43 has been used in the aggregated constraint */
44 char *isint; /* char isint[1+m+n]; */
45 /* isint[k], 1 <= k <= m+n, is a flag that means that variable
46 x[k] is integer (otherwise, continuous) */
47 double *lb; /* double lb[1+m+n]; */
48 /* lb[k], 1 <= k <= m+n, is lower bound of x[k]; -DBL_MAX means
49 that x[k] has no lower bound */
50 int *vlb; /* int vlb[1+m+n]; */
51 /* vlb[k] = k', 1 <= k <= m+n, is the number of integer variable,
52 which defines variable lower bound x[k] >= lb[k] * x[k']; zero
53 means that x[k] has simple lower bound */
54 double *ub; /* double ub[1+m+n]; */
55 /* ub[k], 1 <= k <= m+n, is upper bound of x[k]; +DBL_MAX means
56 that x[k] has no upper bound */
57 int *vub; /* int vub[1+m+n]; */
58 /* vub[k] = k', 1 <= k <= m+n, is the number of integer variable,
59 which defines variable upper bound x[k] <= ub[k] * x[k']; zero
60 means that x[k] has simple upper bound */
61 /*--------------------------------------------------------------*/
62 /* current (fractional) point to be separated */
63 double *x; /* double x[1+m+n]; */
64 /* x[k] is current value of auxiliary (1 <= k <= m) or structural
65 (m+1 <= k <= m+n) variable */
66 /*--------------------------------------------------------------*/
67 /* aggregated constraint sum a[k] * x[k] = b, which is a linear
68 combination of original constraints transformed to equalities
69 by introducing auxiliary variables */
70 int agg_cnt;
71 /* number of rows (original constraints) used to build aggregated
72 constraint, 1 <= agg_cnt <= MAXAGGR */
73 int *agg_row; /* int agg_row[1+MAXAGGR]; */
74 /* agg_row[k], 1 <= k <= agg_cnt, is the row number used to build
75 aggregated constraint */
76 IOSVEC *agg_vec; /* IOSVEC agg_vec[1:m+n]; */
77 /* sparse vector of aggregated constraint coefficients, a[k] */
78 double agg_rhs;
79 /* right-hand side of the aggregated constraint, b */
80 /*--------------------------------------------------------------*/
81 /* bound substitution flags for modified constraint */
82 char *subst; /* char subst[1+m+n]; */
83 /* subst[k], 1 <= k <= m+n, is a bound substitution flag used for
84 variable x[k]:
85 '?' - x[k] is missing in modified constraint
86 'L' - x[k] = (lower bound) + x'[k]
87 'U' - x[k] = (upper bound) - x'[k] */
88 /*--------------------------------------------------------------*/
89 /* modified constraint sum a'[k] * x'[k] = b', where x'[k] >= 0,
90 derived from aggregated constraint by substituting bounds;
91 note that due to substitution of variable bounds there may be
92 additional terms in the modified constraint */
93 IOSVEC *mod_vec; /* IOSVEC mod_vec[1:m+n]; */
94 /* sparse vector of modified constraint coefficients, a'[k] */
95 double mod_rhs;
96 /* right-hand side of the modified constraint, b' */
97 /*--------------------------------------------------------------*/
98 /* cutting plane sum alpha[k] * x[k] <= beta */
99 IOSVEC *cut_vec; /* IOSVEC cut_vec[1:m+n]; */
100 /* sparse vector of cutting plane coefficients, alpha[k] */
101 double cut_rhs;
102 /* right-hand size of the cutting plane, beta */
103 };
104
105 /***********************************************************************
106 * NAME
107 *
108 * ios_mir_init - initialize MIR cut generator
109 *
110 * SYNOPSIS
111 *
112 * #include "glpios.h"
113 * void *ios_mir_init(glp_tree *tree);
114 *
115 * DESCRIPTION
116 *
117 * The routine ios_mir_init initializes the MIR cut generator assuming
118 * that the current subproblem is the root subproblem.
119 *
120 * RETURNS
121 *
122 * The routine ios_mir_init returns a pointer to the MIR cut generator
123 * working area. */
124
125 static void set_row_attrib(glp_tree *tree, struct MIR *mir)
126 { /* set global row attributes */
127 glp_prob *mip = tree->mip;
128 int m = mir->m;
129 int k;
130 for (k = 1; k <= m; k++)
131 { GLPROW *row = mip->row[k];
132 mir->skip[k] = 0;
133 mir->isint[k] = 0;
134 switch (row->type)
135 { case GLP_FR:
136 mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
137 case GLP_LO:
138 mir->lb[k] = row->lb, mir->ub[k] = +DBL_MAX; break;
139 case GLP_UP:
140 mir->lb[k] = -DBL_MAX, mir->ub[k] = row->ub; break;
141 case GLP_DB:
142 mir->lb[k] = row->lb, mir->ub[k] = row->ub; break;
143 case GLP_FX:
144 mir->lb[k] = mir->ub[k] = row->lb; break;
145 default:
146 xassert(row != row);
147 }
148 mir->vlb[k] = mir->vub[k] = 0;
149 }
150 return;
151 }
152
153 static void set_col_attrib(glp_tree *tree, struct MIR *mir)
154 { /* set global column attributes */
155 glp_prob *mip = tree->mip;
156 int m = mir->m;
157 int n = mir->n;
158 int k;
159 for (k = m+1; k <= m+n; k++)
160 { GLPCOL *col = mip->col[k-m];
161 switch (col->kind)
162 { case GLP_CV:
163 mir->isint[k] = 0; break;
164 case GLP_IV:
165 mir->isint[k] = 1; break;
166 default:
167 xassert(col != col);
168 }
169 switch (col->type)
170 { case GLP_FR:
171 mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
172 case GLP_LO:
173 mir->lb[k] = col->lb, mir->ub[k] = +DBL_MAX; break;
174 case GLP_UP:
175 mir->lb[k] = -DBL_MAX, mir->ub[k] = col->ub; break;
176 case GLP_DB:
177 mir->lb[k] = col->lb, mir->ub[k] = col->ub; break;
178 case GLP_FX:
179 mir->lb[k] = mir->ub[k] = col->lb; break;
180 default:
181 xassert(col != col);
182 }
183 mir->vlb[k] = mir->vub[k] = 0;
184 }
185 return;
186 }
187
188 static void set_var_bounds(glp_tree *tree, struct MIR *mir)
189 { /* set variable bounds */
190 glp_prob *mip = tree->mip;
191 int m = mir->m;
192 GLPAIJ *aij;
193 int i, k1, k2;
194 double a1, a2;
195 for (i = 1; i <= m; i++)
196 { /* we need the row to be '>= 0' or '<= 0' */
197 if (!(mir->lb[i] == 0.0 && mir->ub[i] == +DBL_MAX ||
198 mir->lb[i] == -DBL_MAX && mir->ub[i] == 0.0)) continue;
199 /* take first term */
200 aij = mip->row[i]->ptr;
201 if (aij == NULL) continue;
202 k1 = m + aij->col->j, a1 = aij->val;
203 /* take second term */
204 aij = aij->r_next;
205 if (aij == NULL) continue;
206 k2 = m + aij->col->j, a2 = aij->val;
207 /* there must be only two terms */
208 if (aij->r_next != NULL) continue;
209 /* interchange terms, if needed */
210 if (!mir->isint[k1] && mir->isint[k2])
211 ;
212 else if (mir->isint[k1] && !mir->isint[k2])
213 { k2 = k1, a2 = a1;
214 k1 = m + aij->col->j, a1 = aij->val;
215 }
216 else
217 { /* both terms are either continuous or integer */
218 continue;
219 }
220 /* x[k2] should be double-bounded */
221 if (mir->lb[k2] == -DBL_MAX || mir->ub[k2] == +DBL_MAX ||
222 mir->lb[k2] == mir->ub[k2]) continue;
223 /* change signs, if necessary */
224 if (mir->ub[i] == 0.0) a1 = - a1, a2 = - a2;
225 /* now the row has the form a1 * x1 + a2 * x2 >= 0, where x1
226 is continuous, x2 is integer */
227 if (a1 > 0.0)
228 { /* x1 >= - (a2 / a1) * x2 */
229 if (mir->vlb[k1] == 0)
230 { /* set variable lower bound for x1 */
231 mir->lb[k1] = - a2 / a1;
232 mir->vlb[k1] = k2;
233 /* the row should not be used */
234 mir->skip[i] = 1;
235 }
236 }
237 else /* a1 < 0.0 */
238 { /* x1 <= - (a2 / a1) * x2 */
239 if (mir->vub[k1] == 0)
240 { /* set variable upper bound for x1 */
241 mir->ub[k1] = - a2 / a1;
242 mir->vub[k1] = k2;
243 /* the row should not be used */
244 mir->skip[i] = 1;
245 }
246 }
247 }
248 return;
249 }
250
251 static void mark_useless_rows(glp_tree *tree, struct MIR *mir)
252 { /* mark rows which should not be used */
253 glp_prob *mip = tree->mip;
254 int m = mir->m;
255 GLPAIJ *aij;
256 int i, k, nv;
257 for (i = 1; i <= m; i++)
258 { /* free rows should not be used */
259 if (mir->lb[i] == -DBL_MAX && mir->ub[i] == +DBL_MAX)
260 { mir->skip[i] = 1;
261 continue;
262 }
263 nv = 0;
264 for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
265 { k = m + aij->col->j;
266 /* rows with free variables should not be used */
267 if (mir->lb[k] == -DBL_MAX && mir->ub[k] == +DBL_MAX)
268 { mir->skip[i] = 1;
269 break;
270 }
271 /* rows with integer variables having infinite (lower or
272 upper) bound should not be used */
273 if (mir->isint[k] && mir->lb[k] == -DBL_MAX ||
274 mir->isint[k] && mir->ub[k] == +DBL_MAX)
275 { mir->skip[i] = 1;
276 break;
277 }
278 /* count non-fixed variables */
279 if (!(mir->vlb[k] == 0 && mir->vub[k] == 0 &&
280 mir->lb[k] == mir->ub[k])) nv++;
281 }
282 /* rows with all variables fixed should not be used */
283 if (nv == 0)
284 { mir->skip[i] = 1;
285 continue;
286 }
287 }
288 return;
289 }
290
291 void *ios_mir_init(glp_tree *tree)
292 { /* initialize MIR cut generator */
293 glp_prob *mip = tree->mip;
294 int m = mip->m;
295 int n = mip->n;
296 struct MIR *mir;
297 #if _MIR_DEBUG
298 xprintf("ios_mir_init: warning: debug mode enabled\n");
299 #endif
300 /* allocate working area */
301 mir = xmalloc(sizeof(struct MIR));
302 mir->m = m;
303 mir->n = n;
304 mir->skip = xcalloc(1+m, sizeof(char));
305 mir->isint = xcalloc(1+m+n, sizeof(char));
306 mir->lb = xcalloc(1+m+n, sizeof(double));
307 mir->vlb = xcalloc(1+m+n, sizeof(int));
308 mir->ub = xcalloc(1+m+n, sizeof(double));
309 mir->vub = xcalloc(1+m+n, sizeof(int));
310 mir->x = xcalloc(1+m+n, sizeof(double));
311 mir->agg_row = xcalloc(1+MAXAGGR, sizeof(int));
312 mir->agg_vec = ios_create_vec(m+n);
313 mir->subst = xcalloc(1+m+n, sizeof(char));
314 mir->mod_vec = ios_create_vec(m+n);
315 mir->cut_vec = ios_create_vec(m+n);
316 /* set global row attributes */
317 set_row_attrib(tree, mir);
318 /* set global column attributes */
319 set_col_attrib(tree, mir);
320 /* set variable bounds */
321 set_var_bounds(tree, mir);
322 /* mark rows which should not be used */
323 mark_useless_rows(tree, mir);
324 return mir;
325 }
326
327 /***********************************************************************
328 * NAME
329 *
330 * ios_mir_gen - generate MIR cuts
331 *
332 * SYNOPSIS
333 *
334 * #include "glpios.h"
335 * void ios_mir_gen(glp_tree *tree, void *gen, IOSPOOL *pool);
336 *
337 * DESCRIPTION
338 *
339 * The routine ios_mir_gen generates MIR cuts for the current point and
340 * adds them to the cut pool. */
341
342 static void get_current_point(glp_tree *tree, struct MIR *mir)
343 { /* obtain current point */
344 glp_prob *mip = tree->mip;
345 int m = mir->m;
346 int n = mir->n;
347 int k;
348 for (k = 1; k <= m; k++)
349 mir->x[k] = mip->row[k]->prim;
350 for (k = m+1; k <= m+n; k++)
351 mir->x[k] = mip->col[k-m]->prim;
352 return;
353 }
354
355 #if _MIR_DEBUG
356 static void check_current_point(struct MIR *mir)
357 { /* check current point */
358 int m = mir->m;
359 int n = mir->n;
360 int k, kk;
361 double lb, ub, eps;
362 for (k = 1; k <= m+n; k++)
363 { /* determine lower bound */
364 lb = mir->lb[k];
365 kk = mir->vlb[k];
366 if (kk != 0)
367 { xassert(lb != -DBL_MAX);
368 xassert(!mir->isint[k]);
369 xassert(mir->isint[kk]);
370 lb *= mir->x[kk];
371 }
372 /* check lower bound */
373 if (lb != -DBL_MAX)
374 { eps = 1e-6 * (1.0 + fabs(lb));
375 xassert(mir->x[k] >= lb - eps);
376 }
377 /* determine upper bound */
378 ub = mir->ub[k];
379 kk = mir->vub[k];
380 if (kk != 0)
381 { xassert(ub != +DBL_MAX);
382 xassert(!mir->isint[k]);
383 xassert(mir->isint[kk]);
384 ub *= mir->x[kk];
385 }
386 /* check upper bound */
387 if (ub != +DBL_MAX)
388 { eps = 1e-6 * (1.0 + fabs(ub));
389 xassert(mir->x[k] <= ub + eps);
390 }
391 }
392 return;
393 }
394 #endif
395
396 static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i)
397 { /* use original i-th row as initial aggregated constraint */
398 glp_prob *mip = tree->mip;
399 int m = mir->m;
400 GLPAIJ *aij;
401 xassert(1 <= i && i <= m);
402 xassert(!mir->skip[i]);
403 /* mark i-th row in order not to use it in the same aggregated
404 constraint */
405 mir->skip[i] = 2;
406 mir->agg_cnt = 1;
407 mir->agg_row[1] = i;
408 /* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary
409 variable of row i, x[m+j] are structural variables */
410 ios_clear_vec(mir->agg_vec);
411 ios_set_vj(mir->agg_vec, i, 1.0);
412 for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
413 ios_set_vj(mir->agg_vec, m + aij->col->j, - aij->val);
414 mir->agg_rhs = 0.0;
415 #if _MIR_DEBUG
416 ios_check_vec(mir->agg_vec);
417 #endif
418 return;
419 }
420
421 #if _MIR_DEBUG
422 static void check_agg_row(struct MIR *mir)
423 { /* check aggregated constraint */
424 int m = mir->m;
425 int n = mir->n;
426 int j, k;
427 double r, big;
428 /* compute the residual r = sum a[k] * x[k] - b and determine
429 big = max(1, |a[k]|, |b|) */
430 r = 0.0, big = 1.0;
431 for (j = 1; j <= mir->agg_vec->nnz; j++)
432 { k = mir->agg_vec->ind[j];
433 xassert(1 <= k && k <= m+n);
434 r += mir->agg_vec->val[j] * mir->x[k];
435 if (big < fabs(mir->agg_vec->val[j]))
436 big = fabs(mir->agg_vec->val[j]);
437 }
438 r -= mir->agg_rhs;
439 if (big < fabs(mir->agg_rhs))
440 big = fabs(mir->agg_rhs);
441 /* the residual must be close to zero */
442 xassert(fabs(r) <= 1e-6 * big);
443 return;
444 }
445 #endif
446
447 static void subst_fixed_vars(struct MIR *mir)
448 { /* substitute fixed variables into aggregated constraint */
449 int m = mir->m;
450 int n = mir->n;
451 int j, k;
452 for (j = 1; j <= mir->agg_vec->nnz; j++)
453 { k = mir->agg_vec->ind[j];
454 xassert(1 <= k && k <= m+n);
455 if (mir->vlb[k] == 0 && mir->vub[k] == 0 &&
456 mir->lb[k] == mir->ub[k])
457 { /* x[k] is fixed */
458 mir->agg_rhs -= mir->agg_vec->val[j] * mir->lb[k];
459 mir->agg_vec->val[j] = 0.0;
460 }
461 }
462 /* remove terms corresponding to fixed variables */
463 ios_clean_vec(mir->agg_vec, DBL_EPSILON);
464 #if _MIR_DEBUG
465 ios_check_vec(mir->agg_vec);
466 #endif
467 return;
468 }
469
470 static void bound_subst_heur(struct MIR *mir)
471 { /* bound substitution heuristic */
472 int m = mir->m;
473 int n = mir->n;
474 int j, k, kk;
475 double d1, d2;
476 for (j = 1; j <= mir->agg_vec->nnz; j++)
477 { k = mir->agg_vec->ind[j];
478 xassert(1 <= k && k <= m+n);
479 if (mir->isint[k]) continue; /* skip integer variable */
480 /* compute distance from x[k] to its lower bound */
481 kk = mir->vlb[k];
482 if (kk == 0)
483 { if (mir->lb[k] == -DBL_MAX)
484 d1 = DBL_MAX;
485 else
486 d1 = mir->x[k] - mir->lb[k];
487 }
488 else
489 { xassert(1 <= kk && kk <= m+n);
490 xassert(mir->isint[kk]);
491 xassert(mir->lb[k] != -DBL_MAX);
492 d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
493 }
494 /* compute distance from x[k] to its upper bound */
495 kk = mir->vub[k];
496 if (kk == 0)
497 { if (mir->vub[k] == +DBL_MAX)
498 d2 = DBL_MAX;
499 else
500 d2 = mir->ub[k] - mir->x[k];
501 }
502 else
503 { xassert(1 <= kk && kk <= m+n);
504 xassert(mir->isint[kk]);
505 xassert(mir->ub[k] != +DBL_MAX);
506 d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
507 }
508 /* x[k] cannot be free */
509 xassert(d1 != DBL_MAX || d2 != DBL_MAX);
510 /* choose the bound which is closer to x[k] */
511 xassert(mir->subst[k] == '?');
512 if (d1 <= d2)
513 mir->subst[k] = 'L';
514 else
515 mir->subst[k] = 'U';
516 }
517 return;
518 }
519
520 static void build_mod_row(struct MIR *mir)
521 { /* substitute bounds and build modified constraint */
522 int m = mir->m;
523 int n = mir->n;
524 int j, jj, k, kk;
525 /* initially modified constraint is aggregated constraint */
526 ios_copy_vec(mir->mod_vec, mir->agg_vec);
527 mir->mod_rhs = mir->agg_rhs;
528 #if _MIR_DEBUG
529 ios_check_vec(mir->mod_vec);
530 #endif
531 /* substitute bounds for continuous variables; note that due to
532 substitution of variable bounds additional terms may appear in
533 modified constraint */
534 for (j = mir->mod_vec->nnz; j >= 1; j--)
535 { k = mir->mod_vec->ind[j];
536 xassert(1 <= k && k <= m+n);
537 if (mir->isint[k]) continue; /* skip integer variable */
538 if (mir->subst[k] == 'L')
539 { /* x[k] = (lower bound) + x'[k] */
540 xassert(mir->lb[k] != -DBL_MAX);
541 kk = mir->vlb[k];
542 if (kk == 0)
543 { /* x[k] = lb[k] + x'[k] */
544 mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
545 }
546 else
547 { /* x[k] = lb[k] * x[kk] + x'[k] */
548 xassert(mir->isint[kk]);
549 jj = mir->mod_vec->pos[kk];
550 if (jj == 0)
551 { ios_set_vj(mir->mod_vec, kk, 1.0);
552 jj = mir->mod_vec->pos[kk];
553 mir->mod_vec->val[jj] = 0.0;
554 }
555 mir->mod_vec->val[jj] +=
556 mir->mod_vec->val[j] * mir->lb[k];
557 }
558 }
559 else if (mir->subst[k] == 'U')
560 { /* x[k] = (upper bound) - x'[k] */
561 xassert(mir->ub[k] != +DBL_MAX);
562 kk = mir->vub[k];
563 if (kk == 0)
564 { /* x[k] = ub[k] - x'[k] */
565 mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
566 }
567 else
568 { /* x[k] = ub[k] * x[kk] - x'[k] */
569 xassert(mir->isint[kk]);
570 jj = mir->mod_vec->pos[kk];
571 if (jj == 0)
572 { ios_set_vj(mir->mod_vec, kk, 1.0);
573 jj = mir->mod_vec->pos[kk];
574 mir->mod_vec->val[jj] = 0.0;
575 }
576 mir->mod_vec->val[jj] +=
577 mir->mod_vec->val[j] * mir->ub[k];
578 }
579 mir->mod_vec->val[j] = - mir->mod_vec->val[j];
580 }
581 else
582 xassert(k != k);
583 }
584 #if _MIR_DEBUG
585 ios_check_vec(mir->mod_vec);
586 #endif
587 /* substitute bounds for integer variables */
588 for (j = 1; j <= mir->mod_vec->nnz; j++)
589 { k = mir->mod_vec->ind[j];
590 xassert(1 <= k && k <= m+n);
591 if (!mir->isint[k]) continue; /* skip continuous variable */
592 xassert(mir->subst[k] == '?');
593 xassert(mir->vlb[k] == 0 && mir->vub[k] == 0);
594 xassert(mir->lb[k] != -DBL_MAX && mir->ub[k] != +DBL_MAX);
595 if (fabs(mir->lb[k]) <= fabs(mir->ub[k]))
596 { /* x[k] = lb[k] + x'[k] */
597 mir->subst[k] = 'L';
598 mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
599 }
600 else
601 { /* x[k] = ub[k] - x'[k] */
602 mir->subst[k] = 'U';
603 mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
604 mir->mod_vec->val[j] = - mir->mod_vec->val[j];
605 }
606 }
607 #if _MIR_DEBUG
608 ios_check_vec(mir->mod_vec);
609 #endif
610 return;
611 }
612
613 #if _MIR_DEBUG
614 static void check_mod_row(struct MIR *mir)
615 { /* check modified constraint */
616 int m = mir->m;
617 int n = mir->n;
618 int j, k, kk;
619 double r, big, x;
620 /* compute the residual r = sum a'[k] * x'[k] - b' and determine
621 big = max(1, |a[k]|, |b|) */
622 r = 0.0, big = 1.0;
623 for (j = 1; j <= mir->mod_vec->nnz; j++)
624 { k = mir->mod_vec->ind[j];
625 xassert(1 <= k && k <= m+n);
626 if (mir->subst[k] == 'L')
627 { /* x'[k] = x[k] - (lower bound) */
628 xassert(mir->lb[k] != -DBL_MAX);
629 kk = mir->vlb[k];
630 if (kk == 0)
631 x = mir->x[k] - mir->lb[k];
632 else
633 x = mir->x[k] - mir->lb[k] * mir->x[kk];
634 }
635 else if (mir->subst[k] == 'U')
636 { /* x'[k] = (upper bound) - x[k] */
637 xassert(mir->ub[k] != +DBL_MAX);
638 kk = mir->vub[k];
639 if (kk == 0)
640 x = mir->ub[k] - mir->x[k];
641 else
642 x = mir->ub[k] * mir->x[kk] - mir->x[k];
643 }
644 else
645 xassert(k != k);
646 r += mir->mod_vec->val[j] * x;
647 if (big < fabs(mir->mod_vec->val[j]))
648 big = fabs(mir->mod_vec->val[j]);
649 }
650 r -= mir->mod_rhs;
651 if (big < fabs(mir->mod_rhs))
652 big = fabs(mir->mod_rhs);
653 /* the residual must be close to zero */
654 xassert(fabs(r) <= 1e-6 * big);
655 return;
656 }
657 #endif
658
659 /***********************************************************************
660 * mir_ineq - construct MIR inequality
661 *
662 * Given the single constraint mixed integer set
663 *
664 * |N|
665 * X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s},
666 * + + j in N
667 *
668 * this routine constructs the mixed integer rounding (MIR) inequality
669 *
670 * sum alpha[j] * x[j] <= beta + gamma * s,
671 * j in N
672 *
673 * which is valid for X.
674 *
675 * If the MIR inequality has been successfully constructed, the routine
676 * returns zero. Otherwise, if b is close to nearest integer, there may
677 * be numeric difficulties due to big coefficients; so in this case the
678 * routine returns non-zero. */
679
680 static int mir_ineq(const int n, const double a[], const double b,
681 double alpha[], double *beta, double *gamma)
682 { int j;
683 double f, t;
684 if (fabs(b - floor(b + .5)) < 0.01)
685 return 1;
686 f = b - floor(b);
687 for (j = 1; j <= n; j++)
688 { t = (a[j] - floor(a[j])) - f;
689 if (t <= 0.0)
690 alpha[j] = floor(a[j]);
691 else
692 alpha[j] = floor(a[j]) + t / (1.0 - f);
693 }
694 *beta = floor(b);
695 *gamma = 1.0 / (1.0 - f);
696 return 0;
697 }
698
699 /***********************************************************************
700 * cmir_ineq - construct c-MIR inequality
701 *
702 * Given the mixed knapsack set
703 *
704 * MK |N|
705 * X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s,
706 * + + j in N
707 *
708 * x[j] <= u[j]},
709 *
710 * a subset C of variables to be complemented, and a divisor delta > 0,
711 * this routine constructs the complemented MIR (c-MIR) inequality
712 *
713 * sum alpha[j] * x[j] <= beta + gamma * s,
714 * j in N
715 * MK
716 * which is valid for X .
717 *
718 * If the c-MIR inequality has been successfully constructed, the
719 * routine returns zero. Otherwise, if there is a risk of numerical
720 * difficulties due to big coefficients (see comments to the routine
721 * mir_ineq), the routine cmir_ineq returns non-zero. */
722
723 static int cmir_ineq(const int n, const double a[], const double b,
724 const double u[], const char cset[], const double delta,
725 double alpha[], double *beta, double *gamma)
726 { int j;
727 double *aa, bb;
728 aa = alpha, bb = b;
729 for (j = 1; j <= n; j++)
730 { aa[j] = a[j] / delta;
731 if (cset[j])
732 aa[j] = - aa[j], bb -= a[j] * u[j];
733 }
734 bb /= delta;
735 if (mir_ineq(n, aa, bb, alpha, beta, gamma)) return 1;
736 for (j = 1; j <= n; j++)
737 { if (cset[j])
738 alpha[j] = - alpha[j], *beta += alpha[j] * u[j];
739 }
740 *gamma /= delta;
741 return 0;
742 }
743
744 /***********************************************************************
745 * cmir_sep - c-MIR separation heuristic
746 *
747 * Given the mixed knapsack set
748 *
749 * MK |N|
750 * X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s,
751 * + + j in N
752 *
753 * x[j] <= u[j]}
754 *
755 * * *
756 * and a fractional point (x , s ), this routine tries to construct
757 * c-MIR inequality
758 *
759 * sum alpha[j] * x[j] <= beta + gamma * s,
760 * j in N
761 * MK
762 * which is valid for X and has (desirably maximal) violation at the
763 * fractional point given. This is attained by choosing an appropriate
764 * set C of variables to be complemented and a divisor delta > 0, which
765 * together define corresponding c-MIR inequality.
766 *
767 * If a violated c-MIR inequality has been successfully constructed,
768 * the routine returns its violation:
769 *
770 * * *
771 * sum alpha[j] * x [j] - beta - gamma * s ,
772 * j in N
773 *
774 * which is positive. In case of failure the routine returns zero. */
775
776 struct vset { int j; double v; };
777
778 static int cmir_cmp(const void *p1, const void *p2)
779 { const struct vset *v1 = p1, *v2 = p2;
780 if (v1->v < v2->v) return -1;
781 if (v1->v > v2->v) return +1;
782 return 0;
783 }
784
785 static double cmir_sep(const int n, const double a[], const double b,
786 const double u[], const double x[], const double s,
787 double alpha[], double *beta, double *gamma)
788 { int fail, j, k, nv, v;
789 double delta, eps, d_try[1+3], r, r_best;
790 char *cset;
791 struct vset *vset;
792 /* allocate working arrays */
793 cset = xcalloc(1+n, sizeof(char));
794 vset = xcalloc(1+n, sizeof(struct vset));
795 /* choose initial C */
796 for (j = 1; j <= n; j++)
797 cset[j] = (char)(x[j] >= 0.5 * u[j]);
798 /* choose initial delta */
799 r_best = delta = 0.0;
800 for (j = 1; j <= n; j++)
801 { xassert(a[j] != 0.0);
802 /* if x[j] is close to its bounds, skip it */
803 eps = 1e-9 * (1.0 + fabs(u[j]));
804 if (x[j] < eps || x[j] > u[j] - eps) continue;
805 /* try delta = |a[j]| to construct c-MIR inequality */
806 fail = cmir_ineq(n, a, b, u, cset, fabs(a[j]), alpha, beta,
807 gamma);
808 if (fail) continue;
809 /* compute violation */
810 r = - (*beta) - (*gamma) * s;
811 for (k = 1; k <= n; k++) r += alpha[k] * x[k];
812 if (r_best < r) r_best = r, delta = fabs(a[j]);
813 }
814 if (r_best < 0.001) r_best = 0.0;
815 if (r_best == 0.0) goto done;
816 xassert(delta > 0.0);
817 /* try to increase violation by dividing delta by 2, 4, and 8,
818 respectively */
819 d_try[1] = delta / 2.0;
820 d_try[2] = delta / 4.0;
821 d_try[3] = delta / 8.0;
822 for (j = 1; j <= 3; j++)
823 { /* construct c-MIR inequality */
824 fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,
825 gamma);
826 if (fail) continue;
827 /* compute violation */
828 r = - (*beta) - (*gamma) * s;
829 for (k = 1; k <= n; k++) r += alpha[k] * x[k];
830 if (r_best < r) r_best = r, delta = d_try[j];
831 }
832 /* build subset of variables lying strictly between their bounds
833 and order it by nondecreasing values of |x[j] - u[j]/2| */
834 nv = 0;
835 for (j = 1; j <= n; j++)
836 { /* if x[j] is close to its bounds, skip it */
837 eps = 1e-9 * (1.0 + fabs(u[j]));
838 if (x[j] < eps || x[j] > u[j] - eps) continue;
839 /* add x[j] to the subset */
840 nv++;
841 vset[nv].j = j;
842 vset[nv].v = fabs(x[j] - 0.5 * u[j]);
843 }
844 qsort(&vset[1], nv, sizeof(struct vset), cmir_cmp);
845 /* try to increase violation by successively complementing each
846 variable in the subset */
847 for (v = 1; v <= nv; v++)
848 { j = vset[v].j;
849 /* replace x[j] by its complement or vice versa */
850 cset[j] = (char)!cset[j];
851 /* construct c-MIR inequality */
852 fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
853 /* restore the variable */
854 cset[j] = (char)!cset[j];
855 /* do not replace the variable in case of failure */
856 if (fail) continue;
857 /* compute violation */
858 r = - (*beta) - (*gamma) * s;
859 for (k = 1; k <= n; k++) r += alpha[k] * x[k];
860 if (r_best < r) r_best = r, cset[j] = (char)!cset[j];
861 }
862 /* construct the best c-MIR inequality chosen */
863 fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
864 xassert(!fail);
865 done: /* free working arrays */
866 xfree(cset);
867 xfree(vset);
868 /* return to the calling routine */
869 return r_best;
870 }
871
872 static double generate(struct MIR *mir)
873 { /* try to generate violated c-MIR cut for modified constraint */
874 int m = mir->m;
875 int n = mir->n;
876 int j, k, kk, nint;
877 double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma;
878 ios_copy_vec(mir->cut_vec, mir->mod_vec);
879 mir->cut_rhs = mir->mod_rhs;
880 /* remove small terms, which can appear due to substitution of
881 variable bounds */
882 ios_clean_vec(mir->cut_vec, DBL_EPSILON);
883 #if _MIR_DEBUG
884 ios_check_vec(mir->cut_vec);
885 #endif
886 /* remove positive continuous terms to obtain MK relaxation */
887 for (j = 1; j <= mir->cut_vec->nnz; j++)
888 { k = mir->cut_vec->ind[j];
889 xassert(1 <= k && k <= m+n);
890 if (!mir->isint[k] && mir->cut_vec->val[j] > 0.0)
891 mir->cut_vec->val[j] = 0.0;
892 }
893 ios_clean_vec(mir->cut_vec, 0.0);
894 #if _MIR_DEBUG
895 ios_check_vec(mir->cut_vec);
896 #endif
897 /* move integer terms to the beginning of the sparse vector and
898 determine the number of integer variables */
899 nint = 0;
900 for (j = 1; j <= mir->cut_vec->nnz; j++)
901 { k = mir->cut_vec->ind[j];
902 xassert(1 <= k && k <= m+n);
903 if (mir->isint[k])
904 { double temp;
905 nint++;
906 /* interchange elements [nint] and [j] */
907 kk = mir->cut_vec->ind[nint];
908 mir->cut_vec->pos[k] = nint;
909 mir->cut_vec->pos[kk] = j;
910 mir->cut_vec->ind[nint] = k;
911 mir->cut_vec->ind[j] = kk;
912 temp = mir->cut_vec->val[nint];
913 mir->cut_vec->val[nint] = mir->cut_vec->val[j];
914 mir->cut_vec->val[j] = temp;
915 }
916 }
917 #if _MIR_DEBUG
918 ios_check_vec(mir->cut_vec);
919 #endif
920 /* if there is no integer variable, nothing to generate */
921 if (nint == 0) goto done;
922 /* allocate working arrays */
923 u = xcalloc(1+nint, sizeof(double));
924 x = xcalloc(1+nint, sizeof(double));
925 alpha = xcalloc(1+nint, sizeof(double));
926 /* determine u and x */
927 for (j = 1; j <= nint; j++)
928 { k = mir->cut_vec->ind[j];
929 xassert(m+1 <= k && k <= m+n);
930 xassert(mir->isint[k]);
931 u[j] = mir->ub[k] - mir->lb[k];
932 xassert(u[j] >= 1.0);
933 if (mir->subst[k] == 'L')
934 x[j] = mir->x[k] - mir->lb[k];
935 else if (mir->subst[k] == 'U')
936 x[j] = mir->ub[k] - mir->x[k];
937 else
938 xassert(k != k);
939 xassert(x[j] >= -0.001);
940 if (x[j] < 0.0) x[j] = 0.0;
941 }
942 /* compute s = - sum of continuous terms */
943 s = 0.0;
944 for (j = nint+1; j <= mir->cut_vec->nnz; j++)
945 { double x;
946 k = mir->cut_vec->ind[j];
947 xassert(1 <= k && k <= m+n);
948 /* must be continuous */
949 xassert(!mir->isint[k]);
950 if (mir->subst[k] == 'L')
951 { xassert(mir->lb[k] != -DBL_MAX);
952 kk = mir->vlb[k];
953 if (kk == 0)
954 x = mir->x[k] - mir->lb[k];
955 else
956 x = mir->x[k] - mir->lb[k] * mir->x[kk];
957 }
958 else if (mir->subst[k] == 'U')
959 { xassert(mir->ub[k] != +DBL_MAX);
960 kk = mir->vub[k];
961 if (kk == 0)
962 x = mir->ub[k] - mir->x[k];
963 else
964 x = mir->ub[k] * mir->x[kk] - mir->x[k];
965 }
966 else
967 xassert(k != k);
968 xassert(x >= -0.001);
969 if (x < 0.0) x = 0.0;
970 s -= mir->cut_vec->val[j] * x;
971 }
972 xassert(s >= 0.0);
973 /* apply heuristic to obtain most violated c-MIR inequality */
974 b = mir->cut_rhs;
975 r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,
976 &beta, &gamma);
977 if (r_best == 0.0) goto skip;
978 xassert(r_best > 0.0);
979 /* convert to raw cut */
980 /* sum alpha[j] * x[j] <= beta + gamma * s */
981 for (j = 1; j <= nint; j++)
982 mir->cut_vec->val[j] = alpha[j];
983 for (j = nint+1; j <= mir->cut_vec->nnz; j++)
984 { k = mir->cut_vec->ind[j];
985 if (k <= m+n) mir->cut_vec->val[j] *= gamma;
986 }
987 mir->cut_rhs = beta;
988 #if _MIR_DEBUG
989 ios_check_vec(mir->cut_vec);
990 #endif
991 skip: /* free working arrays */
992 xfree(u);
993 xfree(x);
994 xfree(alpha);
995 done: return r_best;
996 }
997
998 #if _MIR_DEBUG
999 static void check_raw_cut(struct MIR *mir, double r_best)
1000 { /* check raw cut before back bound substitution */
1001 int m = mir->m;
1002 int n = mir->n;
1003 int j, k, kk;
1004 double r, big, x;
1005 /* compute the residual r = sum a[k] * x[k] - b and determine
1006 big = max(1, |a[k]|, |b|) */
1007 r = 0.0, big = 1.0;
1008 for (j = 1; j <= mir->cut_vec->nnz; j++)
1009 { k = mir->cut_vec->ind[j];
1010 xassert(1 <= k && k <= m+n);
1011 if (mir->subst[k] == 'L')
1012 { xassert(mir->lb[k] != -DBL_MAX);
1013 kk = mir->vlb[k];
1014 if (kk == 0)
1015 x = mir->x[k] - mir->lb[k];
1016 else
1017 x = mir->x[k] - mir->lb[k] * mir->x[kk];
1018 }
1019 else if (mir->subst[k] == 'U')
1020 { xassert(mir->ub[k] != +DBL_MAX);
1021 kk = mir->vub[k];
1022 if (kk == 0)
1023 x = mir->ub[k] - mir->x[k];
1024 else
1025 x = mir->ub[k] * mir->x[kk] - mir->x[k];
1026 }
1027 else
1028 xassert(k != k);
1029 r += mir->cut_vec->val[j] * x;
1030 if (big < fabs(mir->cut_vec->val[j]))
1031 big = fabs(mir->cut_vec->val[j]);
1032 }
1033 r -= mir->cut_rhs;
1034 if (big < fabs(mir->cut_rhs))
1035 big = fabs(mir->cut_rhs);
1036 /* the residual must be close to r_best */
1037 xassert(fabs(r - r_best) <= 1e-6 * big);
1038 return;
1039 }
1040 #endif
1041
1042 static void back_subst(struct MIR *mir)
1043 { /* back substitution of original bounds */
1044 int m = mir->m;
1045 int n = mir->n;
1046 int j, jj, k, kk;
1047 /* at first, restore bounds of integer variables (because on
1048 restoring variable bounds of continuous variables we need
1049 original, not shifted, bounds of integer variables) */
1050 for (j = 1; j <= mir->cut_vec->nnz; j++)
1051 { k = mir->cut_vec->ind[j];
1052 xassert(1 <= k && k <= m+n);
1053 if (!mir->isint[k]) continue; /* skip continuous */
1054 if (mir->subst[k] == 'L')
1055 { /* x'[k] = x[k] - lb[k] */
1056 xassert(mir->lb[k] != -DBL_MAX);
1057 xassert(mir->vlb[k] == 0);
1058 mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
1059 }
1060 else if (mir->subst[k] == 'U')
1061 { /* x'[k] = ub[k] - x[k] */
1062 xassert(mir->ub[k] != +DBL_MAX);
1063 xassert(mir->vub[k] == 0);
1064 mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
1065 mir->cut_vec->val[j] = - mir->cut_vec->val[j];
1066 }
1067 else
1068 xassert(k != k);
1069 }
1070 /* now restore bounds of continuous variables */
1071 for (j = 1; j <= mir->cut_vec->nnz; j++)
1072 { k = mir->cut_vec->ind[j];
1073 xassert(1 <= k && k <= m+n);
1074 if (mir->isint[k]) continue; /* skip integer */
1075 if (mir->subst[k] == 'L')
1076 { /* x'[k] = x[k] - (lower bound) */
1077 xassert(mir->lb[k] != -DBL_MAX);
1078 kk = mir->vlb[k];
1079 if (kk == 0)
1080 { /* x'[k] = x[k] - lb[k] */
1081 mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
1082 }
1083 else
1084 { /* x'[k] = x[k] - lb[k] * x[kk] */
1085 jj = mir->cut_vec->pos[kk];
1086 #if 0
1087 xassert(jj != 0);
1088 #else
1089 if (jj == 0)
1090 { ios_set_vj(mir->cut_vec, kk, 1.0);
1091 jj = mir->cut_vec->pos[kk];
1092 xassert(jj != 0);
1093 mir->cut_vec->val[jj] = 0.0;
1094 }
1095 #endif
1096 mir->cut_vec->val[jj] -= mir->cut_vec->val[j] *
1097 mir->lb[k];
1098 }
1099 }
1100 else if (mir->subst[k] == 'U')
1101 { /* x'[k] = (upper bound) - x[k] */
1102 xassert(mir->ub[k] != +DBL_MAX);
1103 kk = mir->vub[k];
1104 if (kk == 0)
1105 { /* x'[k] = ub[k] - x[k] */
1106 mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
1107 }
1108 else
1109 { /* x'[k] = ub[k] * x[kk] - x[k] */
1110 jj = mir->cut_vec->pos[kk];
1111 if (jj == 0)
1112 { ios_set_vj(mir->cut_vec, kk, 1.0);
1113 jj = mir->cut_vec->pos[kk];
1114 xassert(jj != 0);
1115 mir->cut_vec->val[jj] = 0.0;
1116 }
1117 mir->cut_vec->val[jj] += mir->cut_vec->val[j] *
1118 mir->ub[k];
1119 }
1120 mir->cut_vec->val[j] = - mir->cut_vec->val[j];
1121 }
1122 else
1123 xassert(k != k);
1124 }
1125 #if _MIR_DEBUG
1126 ios_check_vec(mir->cut_vec);
1127 #endif
1128 return;
1129 }
1130
1131 #if _MIR_DEBUG
1132 static void check_cut_row(struct MIR *mir, double r_best)
1133 { /* check the cut after back bound substitution or elimination of
1134 auxiliary variables */
1135 int m = mir->m;
1136 int n = mir->n;
1137 int j, k;
1138 double r, big;
1139 /* compute the residual r = sum a[k] * x[k] - b and determine
1140 big = max(1, |a[k]|, |b|) */
1141 r = 0.0, big = 1.0;
1142 for (j = 1; j <= mir->cut_vec->nnz; j++)
1143 { k = mir->cut_vec->ind[j];
1144 xassert(1 <= k && k <= m+n);
1145 r += mir->cut_vec->val[j] * mir->x[k];
1146 if (big < fabs(mir->cut_vec->val[j]))
1147 big = fabs(mir->cut_vec->val[j]);
1148 }
1149 r -= mir->cut_rhs;
1150 if (big < fabs(mir->cut_rhs))
1151 big = fabs(mir->cut_rhs);
1152 /* the residual must be close to r_best */
1153 xassert(fabs(r - r_best) <= 1e-6 * big);
1154 return;
1155 }
1156 #endif
1157
1158 static void subst_aux_vars(glp_tree *tree, struct MIR *mir)
1159 { /* final substitution to eliminate auxiliary variables */
1160 glp_prob *mip = tree->mip;
1161 int m = mir->m;
1162 int n = mir->n;
1163 GLPAIJ *aij;
1164 int j, k, kk, jj;
1165 for (j = mir->cut_vec->nnz; j >= 1; j--)
1166 { k = mir->cut_vec->ind[j];
1167 xassert(1 <= k && k <= m+n);
1168 if (k > m) continue; /* skip structurals */
1169 for (aij = mip->row[k]->ptr; aij != NULL; aij = aij->r_next)
1170 { kk = m + aij->col->j; /* structural */
1171 jj = mir->cut_vec->pos[kk];
1172 if (jj == 0)
1173 { ios_set_vj(mir->cut_vec, kk, 1.0);
1174 jj = mir->cut_vec->pos[kk];
1175 mir->cut_vec->val[jj] = 0.0;
1176 }
1177 mir->cut_vec->val[jj] += mir->cut_vec->val[j] * aij->val;
1178 }
1179 mir->cut_vec->val[j] = 0.0;
1180 }
1181 ios_clean_vec(mir->cut_vec, 0.0);
1182 return;
1183 }
1184
1185 static void add_cut(glp_tree *tree, struct MIR *mir)
1186 { /* add constructed cut inequality to the cut pool */
1187 int m = mir->m;
1188 int n = mir->n;
1189 int j, k, len;
1190 int *ind = xcalloc(1+n, sizeof(int));
1191 double *val = xcalloc(1+n, sizeof(double));
1192 len = 0;
1193 for (j = mir->cut_vec->nnz; j >= 1; j--)
1194 { k = mir->cut_vec->ind[j];
1195 xassert(m+1 <= k && k <= m+n);
1196 len++, ind[len] = k - m, val[len] = mir->cut_vec->val[j];
1197 }
1198 #if 0
1199 ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP,
1200 mir->cut_rhs);
1201 #else
1202 glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
1203 mir->cut_rhs);
1204 #endif
1205 xfree(ind);
1206 xfree(val);
1207 return;
1208 }
1209
1210 static int aggregate_row(glp_tree *tree, struct MIR *mir)
1211 { /* try to aggregate another row */
1212 glp_prob *mip = tree->mip;
1213 int m = mir->m;
1214 int n = mir->n;
1215 GLPAIJ *aij;
1216 IOSVEC *v;
1217 int ii, j, jj, k, kk, kappa = 0, ret = 0;
1218 double d1, d2, d, d_max = 0.0;
1219 /* choose appropriate structural variable in the aggregated row
1220 to be substituted */
1221 for (j = 1; j <= mir->agg_vec->nnz; j++)
1222 { k = mir->agg_vec->ind[j];
1223 xassert(1 <= k && k <= m+n);
1224 if (k <= m) continue; /* skip auxiliary var */
1225 if (mir->isint[k]) continue; /* skip integer var */
1226 if (fabs(mir->agg_vec->val[j]) < 0.001) continue;
1227 /* compute distance from x[k] to its lower bound */
1228 kk = mir->vlb[k];
1229 if (kk == 0)
1230 { if (mir->lb[k] == -DBL_MAX)
1231 d1 = DBL_MAX;
1232 else
1233 d1 = mir->x[k] - mir->lb[k];
1234 }
1235 else
1236 { xassert(1 <= kk && kk <= m+n);
1237 xassert(mir->isint[kk]);
1238 xassert(mir->lb[k] != -DBL_MAX);
1239 d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
1240 }
1241 /* compute distance from x[k] to its upper bound */
1242 kk = mir->vub[k];
1243 if (kk == 0)
1244 { if (mir->vub[k] == +DBL_MAX)
1245 d2 = DBL_MAX;
1246 else
1247 d2 = mir->ub[k] - mir->x[k];
1248 }
1249 else
1250 { xassert(1 <= kk && kk <= m+n);
1251 xassert(mir->isint[kk]);
1252 xassert(mir->ub[k] != +DBL_MAX);
1253 d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
1254 }
1255 /* x[k] cannot be free */
1256 xassert(d1 != DBL_MAX || d2 != DBL_MAX);
1257 /* d = min(d1, d2) */
1258 d = (d1 <= d2 ? d1 : d2);
1259 xassert(d != DBL_MAX);
1260 /* should not be close to corresponding bound */
1261 if (d < 0.001) continue;
1262 if (d_max < d) d_max = d, kappa = k;
1263 }
1264 if (kappa == 0)
1265 { /* nothing chosen */
1266 ret = 1;
1267 goto done;
1268 }
1269 /* x[kappa] has been chosen */
1270 xassert(m+1 <= kappa && kappa <= m+n);
1271 xassert(!mir->isint[kappa]);
1272 /* find another row, which have not been used yet, to eliminate
1273 x[kappa] from the aggregated row */
1274 for (ii = 1; ii <= m; ii++)
1275 { if (mir->skip[ii]) continue;
1276 for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
1277 if (aij->col->j == kappa - m) break;
1278 if (aij != NULL && fabs(aij->val) >= 0.001) break;
1279 }
1280 if (ii > m)
1281 { /* nothing found */
1282 ret = 2;
1283 goto done;
1284 }
1285 /* row ii has been found; include it in the aggregated list */
1286 mir->agg_cnt++;
1287 xassert(mir->agg_cnt <= MAXAGGR);
1288 mir->agg_row[mir->agg_cnt] = ii;
1289 mir->skip[ii] = 2;
1290 /* v := new row */
1291 v = ios_create_vec(m+n);
1292 ios_set_vj(v, ii, 1.0);
1293 for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
1294 ios_set_vj(v, m + aij->col->j, - aij->val);
1295 #if _MIR_DEBUG
1296 ios_check_vec(v);
1297 #endif
1298 /* perform gaussian elimination to remove x[kappa] */
1299 j = mir->agg_vec->pos[kappa];
1300 xassert(j != 0);
1301 jj = v->pos[kappa];
1302 xassert(jj != 0);
1303 ios_linear_comb(mir->agg_vec,
1304 - mir->agg_vec->val[j] / v->val[jj], v);
1305 ios_delete_vec(v);
1306 ios_set_vj(mir->agg_vec, kappa, 0.0);
1307 #if _MIR_DEBUG
1308 ios_check_vec(mir->agg_vec);
1309 #endif
1310 done: return ret;
1311 }
1312
1313 void ios_mir_gen(glp_tree *tree, void *gen)
1314 { /* main routine to generate MIR cuts */
1315 glp_prob *mip = tree->mip;
1316 struct MIR *mir = gen;
1317 int m = mir->m;
1318 int n = mir->n;
1319 int i;
1320 double r_best;
1321 xassert(mip->m >= m);
1322 xassert(mip->n == n);
1323 /* obtain current point */
1324 get_current_point(tree, mir);
1325 #if _MIR_DEBUG
1326 /* check current point */
1327 check_current_point(mir);
1328 #endif
1329 /* reset bound substitution flags */
1330 memset(&mir->subst[1], '?', m+n);
1331 /* try to generate a set of violated MIR cuts */
1332 for (i = 1; i <= m; i++)
1333 { if (mir->skip[i]) continue;
1334 /* use original i-th row as initial aggregated constraint */
1335 initial_agg_row(tree, mir, i);
1336 loop: ;
1337 #if _MIR_DEBUG
1338 /* check aggregated row */
1339 check_agg_row(mir);
1340 #endif
1341 /* substitute fixed variables into aggregated constraint */
1342 subst_fixed_vars(mir);
1343 #if _MIR_DEBUG
1344 /* check aggregated row */
1345 check_agg_row(mir);
1346 #endif
1347 #if _MIR_DEBUG
1348 /* check bound substitution flags */
1349 { int k;
1350 for (k = 1; k <= m+n; k++)
1351 xassert(mir->subst[k] == '?');
1352 }
1353 #endif
1354 /* apply bound substitution heuristic */
1355 bound_subst_heur(mir);
1356 /* substitute bounds and build modified constraint */
1357 build_mod_row(mir);
1358 #if _MIR_DEBUG
1359 /* check modified row */
1360 check_mod_row(mir);
1361 #endif
1362 /* try to generate violated c-MIR cut for modified row */
1363 r_best = generate(mir);
1364 if (r_best > 0.0)
1365 { /* success */
1366 #if _MIR_DEBUG
1367 /* check raw cut before back bound substitution */
1368 check_raw_cut(mir, r_best);
1369 #endif
1370 /* back substitution of original bounds */
1371 back_subst(mir);
1372 #if _MIR_DEBUG
1373 /* check the cut after back bound substitution */
1374 check_cut_row(mir, r_best);
1375 #endif
1376 /* final substitution to eliminate auxiliary variables */
1377 subst_aux_vars(tree, mir);
1378 #if _MIR_DEBUG
1379 /* check the cut after elimination of auxiliaries */
1380 check_cut_row(mir, r_best);
1381 #endif
1382 /* add constructed cut inequality to the cut pool */
1383 add_cut(tree, mir);
1384 }
1385 /* reset bound substitution flags */
1386 { int j, k;
1387 for (j = 1; j <= mir->mod_vec->nnz; j++)
1388 { k = mir->mod_vec->ind[j];
1389 xassert(1 <= k && k <= m+n);
1390 xassert(mir->subst[k] != '?');
1391 mir->subst[k] = '?';
1392 }
1393 }
1394 if (r_best == 0.0)
1395 { /* failure */
1396 if (mir->agg_cnt < MAXAGGR)
1397 { /* try to aggregate another row */
1398 if (aggregate_row(tree, mir) == 0) goto loop;
1399 }
1400 }
1401 /* unmark rows used in the aggregated constraint */
1402 { int k, ii;
1403 for (k = 1; k <= mir->agg_cnt; k++)
1404 { ii = mir->agg_row[k];
1405 xassert(1 <= ii && ii <= m);
1406 xassert(mir->skip[ii] == 2);
1407 mir->skip[ii] = 0;
1408 }
1409 }
1410 }
1411 return;
1412 }
1413
1414 /***********************************************************************
1415 * NAME
1416 *
1417 * ios_mir_term - terminate MIR cut generator
1418 *
1419 * SYNOPSIS
1420 *
1421 * #include "glpios.h"
1422 * void ios_mir_term(void *gen);
1423 *
1424 * DESCRIPTION
1425 *
1426 * The routine ios_mir_term deletes the MIR cut generator working area
1427 * freeing all the memory allocated to it. */
1428
1429 void ios_mir_term(void *gen)
1430 { struct MIR *mir = gen;
1431 xfree(mir->skip);
1432 xfree(mir->isint);
1433 xfree(mir->lb);
1434 xfree(mir->vlb);
1435 xfree(mir->ub);
1436 xfree(mir->vub);
1437 xfree(mir->x);
1438 xfree(mir->agg_row);
1439 ios_delete_vec(mir->agg_vec);
1440 xfree(mir->subst);
1441 ios_delete_vec(mir->mod_vec);
1442 ios_delete_vec(mir->cut_vec);
1443 xfree(mir);
1444 return;
1445 }
1446
1447 /* eof */