lemon-project-template-glpk
view deps/glpk/src/glpnpp04.c @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 21:43:29 +0100 |
| parents | |
| children |
line source
1 /* glpnpp04.c */
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 "glpnpp.h"
27 /***********************************************************************
28 * NAME
29 *
30 * npp_binarize_prob - binarize MIP problem
31 *
32 * SYNOPSIS
33 *
34 * #include "glpnpp.h"
35 * int npp_binarize_prob(NPP *npp);
36 *
37 * DESCRIPTION
38 *
39 * The routine npp_binarize_prob replaces in the original MIP problem
40 * every integer variable:
41 *
42 * l[q] <= x[q] <= u[q], (1)
43 *
44 * where l[q] < u[q], by an equivalent sum of binary variables.
45 *
46 * RETURNS
47 *
48 * The routine returns the number of integer variables for which the
49 * transformation failed, because u[q] - l[q] > d_max.
50 *
51 * PROBLEM TRANSFORMATION
52 *
53 * If variable x[q] has non-zero lower bound, it is first processed
54 * with the routine npp_lbnd_col. Thus, we can assume that:
55 *
56 * 0 <= x[q] <= u[q]. (2)
57 *
58 * If u[q] = 1, variable x[q] is already binary, so further processing
59 * is not needed. Let, therefore, that 2 <= u[q] <= d_max, and n be a
60 * smallest integer such that u[q] <= 2^n - 1 (n >= 2, since u[q] >= 2).
61 * Then variable x[q] can be replaced by the following sum:
62 *
63 * n-1
64 * x[q] = sum 2^k x[k], (3)
65 * k=0
66 *
67 * where x[k] are binary columns (variables). If u[q] < 2^n - 1, the
68 * following additional inequality constraint must be also included in
69 * the transformed problem:
70 *
71 * n-1
72 * sum 2^k x[k] <= u[q]. (4)
73 * k=0
74 *
75 * Note: Assuming that in the transformed problem x[q] becomes binary
76 * variable x[0], this transformation causes new n-1 binary variables
77 * to appear.
78 *
79 * Substituting x[q] from (3) to the objective row gives:
80 *
81 * z = sum c[j] x[j] + c[0] =
82 * j
83 *
84 * = sum c[j] x[j] + c[q] x[q] + c[0] =
85 * j!=q
86 * n-1
87 * = sum c[j] x[j] + c[q] sum 2^k x[k] + c[0] =
88 * j!=q k=0
89 * n-1
90 * = sum c[j] x[j] + sum c[k] x[k] + c[0],
91 * j!=q k=0
92 *
93 * where:
94 *
95 * c[k] = 2^k c[q], k = 0, ..., n-1. (5)
96 *
97 * And substituting x[q] from (3) to i-th constraint row i gives:
98 *
99 * L[i] <= sum a[i,j] x[j] <= U[i] ==>
100 * j
101 *
102 * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==>
103 * j!=q
104 * n-1
105 * L[i] <= sum a[i,j] x[j] + a[i,q] sum 2^k x[k] <= U[i] ==>
106 * j!=q k=0
107 * n-1
108 * L[i] <= sum a[i,j] x[j] + sum a[i,k] x[k] <= U[i],
109 * j!=q k=0
110 *
111 * where:
112 *
113 * a[i,k] = 2^k a[i,q], k = 0, ..., n-1. (6)
114 *
115 * RECOVERING SOLUTION
116 *
117 * Value of variable x[q] is computed with formula (3). */
119 struct binarize
120 { int q;
121 /* column reference number for x[q] = x[0] */
122 int j;
123 /* column reference number for x[1]; x[2] has reference number
124 j+1, x[3] - j+2, etc. */
125 int n;
126 /* total number of binary variables, n >= 2 */
127 };
129 static int rcv_binarize_prob(NPP *npp, void *info);
131 int npp_binarize_prob(NPP *npp)
132 { /* binarize MIP problem */
133 struct binarize *info;
134 NPPROW *row;
135 NPPCOL *col, *bin;
136 NPPAIJ *aij;
137 int u, n, k, temp, nfails, nvars, nbins, nrows;
138 /* new variables will be added to the end of the column list, so
139 we go from the end to beginning of the column list */
140 nfails = nvars = nbins = nrows = 0;
141 for (col = npp->c_tail; col != NULL; col = col->prev)
142 { /* skip continuous variable */
143 if (!col->is_int) continue;
144 /* skip fixed variable */
145 if (col->lb == col->ub) continue;
146 /* skip binary variable */
147 if (col->lb == 0.0 && col->ub == 1.0) continue;
148 /* check if the transformation is applicable */
149 if (col->lb < -1e6 || col->ub > +1e6 ||
150 col->ub - col->lb > 4095.0)
151 { /* unfortunately, not */
152 nfails++;
153 continue;
154 }
155 /* process integer non-binary variable x[q] */
156 nvars++;
157 /* make x[q] non-negative, if its lower bound is non-zero */
158 if (col->lb != 0.0)
159 npp_lbnd_col(npp, col);
160 /* now 0 <= x[q] <= u[q] */
161 xassert(col->lb == 0.0);
162 u = (int)col->ub;
163 xassert(col->ub == (double)u);
164 /* if x[q] is binary, further processing is not needed */
165 if (u == 1) continue;
166 /* determine smallest n such that u <= 2^n - 1 (thus, n is the
167 number of binary variables needed) */
168 n = 2, temp = 4;
169 while (u >= temp)
170 n++, temp += temp;
171 nbins += n;
172 /* create transformation stack entry */
173 info = npp_push_tse(npp,
174 rcv_binarize_prob, sizeof(struct binarize));
175 info->q = col->j;
176 info->j = 0; /* will be set below */
177 info->n = n;
178 /* if u < 2^n - 1, we need one additional row for (4) */
179 if (u < temp - 1)
180 { row = npp_add_row(npp), nrows++;
181 row->lb = -DBL_MAX, row->ub = u;
182 }
183 else
184 row = NULL;
185 /* in the transformed problem variable x[q] becomes binary
186 variable x[0], so its objective and constraint coefficients
187 are not changed */
188 col->ub = 1.0;
189 /* include x[0] into constraint (4) */
190 if (row != NULL)
191 npp_add_aij(npp, row, col, 1.0);
192 /* add other binary variables x[1], ..., x[n-1] */
193 for (k = 1, temp = 2; k < n; k++, temp += temp)
194 { /* add new binary variable x[k] */
195 bin = npp_add_col(npp);
196 bin->is_int = 1;
197 bin->lb = 0.0, bin->ub = 1.0;
198 bin->coef = (double)temp * col->coef;
199 /* store column reference number for x[1] */
200 if (info->j == 0)
201 info->j = bin->j;
202 else
203 xassert(info->j + (k-1) == bin->j);
204 /* duplicate constraint coefficients for x[k]; this also
205 automatically includes x[k] into constraint (4) */
206 for (aij = col->ptr; aij != NULL; aij = aij->c_next)
207 npp_add_aij(npp, aij->row, bin, (double)temp * aij->val);
208 }
209 }
210 if (nvars > 0)
211 xprintf("%d integer variable(s) were replaced by %d binary one"
212 "s\n", nvars, nbins);
213 if (nrows > 0)
214 xprintf("%d row(s) were added due to binarization\n", nrows);
215 if (nfails > 0)
216 xprintf("Binarization failed for %d integer variable(s)\n",
217 nfails);
218 return nfails;
219 }
221 static int rcv_binarize_prob(NPP *npp, void *_info)
222 { /* recovery binarized variable */
223 struct binarize *info = _info;
224 int k, temp;
225 double sum;
226 /* compute value of x[q]; see formula (3) */
227 sum = npp->c_value[info->q];
228 for (k = 1, temp = 2; k < info->n; k++, temp += temp)
229 sum += (double)temp * npp->c_value[info->j + (k-1)];
230 npp->c_value[info->q] = sum;
231 return 0;
232 }
234 /**********************************************************************/
236 struct elem
237 { /* linear form element a[j] x[j] */
238 double aj;
239 /* non-zero coefficient value */
240 NPPCOL *xj;
241 /* pointer to variable (column) */
242 struct elem *next;
243 /* pointer to another term */
244 };
246 static struct elem *copy_form(NPP *npp, NPPROW *row, double s)
247 { /* copy linear form */
248 NPPAIJ *aij;
249 struct elem *ptr, *e;
250 ptr = NULL;
251 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
252 { e = dmp_get_atom(npp->pool, sizeof(struct elem));
253 e->aj = s * aij->val;
254 e->xj = aij->col;
255 e->next = ptr;
256 ptr = e;
257 }
258 return ptr;
259 }
261 static void drop_form(NPP *npp, struct elem *ptr)
262 { /* drop linear form */
263 struct elem *e;
264 while (ptr != NULL)
265 { e = ptr;
266 ptr = e->next;
267 dmp_free_atom(npp->pool, e, sizeof(struct elem));
268 }
269 return;
270 }
272 /***********************************************************************
273 * NAME
274 *
275 * npp_is_packing - test if constraint is packing inequality
276 *
277 * SYNOPSIS
278 *
279 * #include "glpnpp.h"
280 * int npp_is_packing(NPP *npp, NPPROW *row);
281 *
282 * RETURNS
283 *
284 * If the specified row (constraint) is packing inequality (see below),
285 * the routine npp_is_packing returns non-zero. Otherwise, it returns
286 * zero.
287 *
288 * PACKING INEQUALITIES
289 *
290 * In canonical format the packing inequality is the following:
291 *
292 * sum x[j] <= 1, (1)
293 * j in J
294 *
295 * where all variables x[j] are binary. This inequality expresses the
296 * condition that in any integer feasible solution at most one variable
297 * from set J can take non-zero (unity) value while other variables
298 * must be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because
299 * if J is empty or |J| = 1, the inequality (1) is redundant.
300 *
301 * In general case the packing inequality may include original variables
302 * x[j] as well as their complements x~[j]:
303 *
304 * sum x[j] + sum x~[j] <= 1, (2)
305 * j in Jp j in Jn
306 *
307 * where Jp and Jn are not intersected. Therefore, using substitution
308 * x~[j] = 1 - x[j] gives the packing inequality in generalized format:
309 *
310 * sum x[j] - sum x[j] <= 1 - |Jn|. (3)
311 * j in Jp j in Jn */
313 int npp_is_packing(NPP *npp, NPPROW *row)
314 { /* test if constraint is packing inequality */
315 NPPCOL *col;
316 NPPAIJ *aij;
317 int b;
318 xassert(npp == npp);
319 if (!(row->lb == -DBL_MAX && row->ub != +DBL_MAX))
320 return 0;
321 b = 1;
322 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
323 { col = aij->col;
324 if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
325 return 0;
326 if (aij->val == +1.0)
327 ;
328 else if (aij->val == -1.0)
329 b--;
330 else
331 return 0;
332 }
333 if (row->ub != (double)b) return 0;
334 return 1;
335 }
337 /***********************************************************************
338 * NAME
339 *
340 * npp_hidden_packing - identify hidden packing inequality
341 *
342 * SYNOPSIS
343 *
344 * #include "glpnpp.h"
345 * int npp_hidden_packing(NPP *npp, NPPROW *row);
346 *
347 * DESCRIPTION
348 *
349 * The routine npp_hidden_packing processes specified inequality
350 * constraint, which includes only binary variables, and the number of
351 * the variables is not less than two. If the original inequality is
352 * equivalent to a packing inequality, the routine replaces it by this
353 * equivalent inequality. If the original constraint is double-sided
354 * inequality, it is replaced by a pair of single-sided inequalities,
355 * if necessary.
356 *
357 * RETURNS
358 *
359 * If the original inequality constraint was replaced by equivalent
360 * packing inequality, the routine npp_hidden_packing returns non-zero.
361 * Otherwise, it returns zero.
362 *
363 * PROBLEM TRANSFORMATION
364 *
365 * Consider an inequality constraint:
366 *
367 * sum a[j] x[j] <= b, (1)
368 * j in J
369 *
370 * where all variables x[j] are binary, and |J| >= 2. (In case of '>='
371 * inequality it can be transformed to '<=' format by multiplying both
372 * its sides by -1.)
373 *
374 * Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution
375 * x[j] = 1 - x~[j] for all j in Jn, we have:
376 *
377 * sum a[j] x[j] <= b ==>
378 * j in J
379 *
380 * sum a[j] x[j] + sum a[j] x[j] <= b ==>
381 * j in Jp j in Jn
382 *
383 * sum a[j] x[j] + sum a[j] (1 - x~[j]) <= b ==>
384 * j in Jp j in Jn
385 *
386 * sum a[j] x[j] - sum a[j] x~[j] <= b - sum a[j].
387 * j in Jp j in Jn j in Jn
388 *
389 * Thus, meaning the transformation above, we can assume that in
390 * inequality (1) all coefficients a[j] are positive. Moreover, we can
391 * assume that a[j] <= b. In fact, let a[j] > b; then the following
392 * three cases are possible:
393 *
394 * 1) b < 0. In this case inequality (1) is infeasible, so the problem
395 * has no feasible solution (see the routine npp_analyze_row);
396 *
397 * 2) b = 0. In this case inequality (1) is a forcing inequality on its
398 * upper bound (see the routine npp_forcing row), from which it
399 * follows that all variables x[j] should be fixed at zero;
400 *
401 * 3) b > 0. In this case inequality (1) defines an implied zero upper
402 * bound for variable x[j] (see the routine npp_implied_bounds), from
403 * which it follows that x[j] should be fixed at zero.
404 *
405 * It is assumed that all three cases listed above have been recognized
406 * by the routine npp_process_prob, which performs basic MIP processing
407 * prior to a call the routine npp_hidden_packing. So, if one of these
408 * cases occurs, we should just skip processing such constraint.
409 *
410 * Thus, let 0 < a[j] <= b. Then it is obvious that constraint (1) is
411 * equivalent to packing inquality only if:
412 *
413 * a[j] + a[k] > b + eps (2)
414 *
415 * for all j, k in J, j != k, where eps is an absolute tolerance for
416 * row (linear form) value. Checking the condition (2) for all j and k,
417 * j != k, requires time O(|J|^2). However, this time can be reduced to
418 * O(|J|), if use minimal a[j] and a[k], in which case it is sufficient
419 * to check the condition (2) only once.
420 *
421 * Once the original inequality (1) is replaced by equivalent packing
422 * inequality, we need to perform back substitution x~[j] = 1 - x[j] for
423 * all j in Jn (see above).
424 *
425 * RECOVERING SOLUTION
426 *
427 * None needed. */
429 static int hidden_packing(NPP *npp, struct elem *ptr, double *_b)
430 { /* process inequality constraint: sum a[j] x[j] <= b;
431 0 - specified row is NOT hidden packing inequality;
432 1 - specified row is packing inequality;
433 2 - specified row is hidden packing inequality. */
434 struct elem *e, *ej, *ek;
435 int neg;
436 double b = *_b, eps;
437 xassert(npp == npp);
438 /* a[j] must be non-zero, x[j] must be binary, for all j in J */
439 for (e = ptr; e != NULL; e = e->next)
440 { xassert(e->aj != 0.0);
441 xassert(e->xj->is_int);
442 xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0);
443 }
444 /* check if the specified inequality constraint already has the
445 form of packing inequality */
446 neg = 0; /* neg is |Jn| */
447 for (e = ptr; e != NULL; e = e->next)
448 { if (e->aj == +1.0)
449 ;
450 else if (e->aj == -1.0)
451 neg++;
452 else
453 break;
454 }
455 if (e == NULL)
456 { /* all coefficients a[j] are +1 or -1; check rhs b */
457 if (b == (double)(1 - neg))
458 { /* it is packing inequality; no processing is needed */
459 return 1;
460 }
461 }
462 /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j]
463 positive; the result is a~[j] = |a[j]| and new rhs b */
464 for (e = ptr; e != NULL; e = e->next)
465 if (e->aj < 0) b -= e->aj;
466 /* now a[j] > 0 for all j in J (actually |a[j]| are used) */
467 /* if a[j] > b, skip processing--this case must not appear */
468 for (e = ptr; e != NULL; e = e->next)
469 if (fabs(e->aj) > b) return 0;
470 /* now 0 < a[j] <= b for all j in J */
471 /* find two minimal coefficients a[j] and a[k], j != k */
472 ej = NULL;
473 for (e = ptr; e != NULL; e = e->next)
474 if (ej == NULL || fabs(ej->aj) > fabs(e->aj)) ej = e;
475 xassert(ej != NULL);
476 ek = NULL;
477 for (e = ptr; e != NULL; e = e->next)
478 if (e != ej)
479 if (ek == NULL || fabs(ek->aj) > fabs(e->aj)) ek = e;
480 xassert(ek != NULL);
481 /* the specified constraint is equivalent to packing inequality
482 iff a[j] + a[k] > b + eps */
483 eps = 1e-3 + 1e-6 * fabs(b);
484 if (fabs(ej->aj) + fabs(ek->aj) <= b + eps) return 0;
485 /* perform back substitution x~[j] = 1 - x[j] and construct the
486 final equivalent packing inequality in generalized format */
487 b = 1.0;
488 for (e = ptr; e != NULL; e = e->next)
489 { if (e->aj > 0.0)
490 e->aj = +1.0;
491 else /* e->aj < 0.0 */
492 e->aj = -1.0, b -= 1.0;
493 }
494 *_b = b;
495 return 2;
496 }
498 int npp_hidden_packing(NPP *npp, NPPROW *row)
499 { /* identify hidden packing inequality */
500 NPPROW *copy;
501 NPPAIJ *aij;
502 struct elem *ptr, *e;
503 int kase, ret, count = 0;
504 double b;
505 /* the row must be inequality constraint */
506 xassert(row->lb < row->ub);
507 for (kase = 0; kase <= 1; kase++)
508 { if (kase == 0)
509 { /* process row upper bound */
510 if (row->ub == +DBL_MAX) continue;
511 ptr = copy_form(npp, row, +1.0);
512 b = + row->ub;
513 }
514 else
515 { /* process row lower bound */
516 if (row->lb == -DBL_MAX) continue;
517 ptr = copy_form(npp, row, -1.0);
518 b = - row->lb;
519 }
520 /* now the inequality has the form "sum a[j] x[j] <= b" */
521 ret = hidden_packing(npp, ptr, &b);
522 xassert(0 <= ret && ret <= 2);
523 if (kase == 1 && ret == 1 || ret == 2)
524 { /* the original inequality has been identified as hidden
525 packing inequality */
526 count++;
527 #ifdef GLP_DEBUG
528 xprintf("Original constraint:\n");
529 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
530 xprintf(" %+g x%d", aij->val, aij->col->j);
531 if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb);
532 if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub);
533 xprintf("\n");
534 xprintf("Equivalent packing inequality:\n");
535 for (e = ptr; e != NULL; e = e->next)
536 xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j);
537 xprintf(", <= %g\n", b);
538 #endif
539 if (row->lb == -DBL_MAX || row->ub == +DBL_MAX)
540 { /* the original row is single-sided inequality; no copy
541 is needed */
542 copy = NULL;
543 }
544 else
545 { /* the original row is double-sided inequality; we need
546 to create its copy for other bound before replacing it
547 with the equivalent inequality */
548 copy = npp_add_row(npp);
549 if (kase == 0)
550 { /* the copy is for lower bound */
551 copy->lb = row->lb, copy->ub = +DBL_MAX;
552 }
553 else
554 { /* the copy is for upper bound */
555 copy->lb = -DBL_MAX, copy->ub = row->ub;
556 }
557 /* copy original row coefficients */
558 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
559 npp_add_aij(npp, copy, aij->col, aij->val);
560 }
561 /* replace the original inequality by equivalent one */
562 npp_erase_row(npp, row);
563 row->lb = -DBL_MAX, row->ub = b;
564 for (e = ptr; e != NULL; e = e->next)
565 npp_add_aij(npp, row, e->xj, e->aj);
566 /* continue processing lower bound for the copy */
567 if (copy != NULL) row = copy;
568 }
569 drop_form(npp, ptr);
570 }
571 return count;
572 }
574 /***********************************************************************
575 * NAME
576 *
577 * npp_implied_packing - identify implied packing inequality
578 *
579 * SYNOPSIS
580 *
581 * #include "glpnpp.h"
582 * int npp_implied_packing(NPP *npp, NPPROW *row, int which,
583 * NPPCOL *var[], char set[]);
584 *
585 * DESCRIPTION
586 *
587 * The routine npp_implied_packing processes specified row (constraint)
588 * of general format:
589 *
590 * L <= sum a[j] x[j] <= U. (1)
591 * j
592 *
593 * If which = 0, only lower bound L, which must exist, is considered,
594 * while upper bound U is ignored. Similarly, if which = 1, only upper
595 * bound U, which must exist, is considered, while lower bound L is
596 * ignored. Thus, if the specified row is a double-sided inequality or
597 * equality constraint, this routine should be called twice for both
598 * lower and upper bounds.
599 *
600 * The routine npp_implied_packing attempts to find a non-trivial (i.e.
601 * having not less than two binary variables) packing inequality:
602 *
603 * sum x[j] - sum x[j] <= 1 - |Jn|, (2)
604 * j in Jp j in Jn
605 *
606 * which is relaxation of the constraint (1) in the sense that any
607 * solution satisfying to that constraint also satisfies to the packing
608 * inequality (2). If such relaxation exists, the routine stores
609 * pointers to descriptors of corresponding binary variables and their
610 * flags, resp., to locations var[1], var[2], ..., var[len] and set[1],
611 * set[2], ..., set[len], where set[j] = 0 means that j in Jp and
612 * set[j] = 1 means that j in Jn.
613 *
614 * RETURNS
615 *
616 * The routine npp_implied_packing returns len, which is the total
617 * number of binary variables in the packing inequality found, len >= 2.
618 * However, if the relaxation does not exist, the routine returns zero.
619 *
620 * ALGORITHM
621 *
622 * If which = 0, the constraint coefficients (1) are multiplied by -1
623 * and b is assigned -L; if which = 1, the constraint coefficients (1)
624 * are not changed and b is assigned +U. In both cases the specified
625 * constraint gets the following format:
626 *
627 * sum a[j] x[j] <= b. (3)
628 * j
629 *
630 * (Note that (3) is a relaxation of (1), because one of bounds L or U
631 * is ignored.)
632 *
633 * Let J be set of binary variables, Kp be set of non-binary (integer
634 * or continuous) variables with a[j] > 0, and Kn be set of non-binary
635 * variables with a[j] < 0. Then the inequality (3) can be written as
636 * follows:
637 *
638 * sum a[j] x[j] <= b - sum a[j] x[j] - sum a[j] x[j]. (4)
639 * j in J j in Kp j in Kn
640 *
641 * To get rid of non-binary variables we can replace the inequality (4)
642 * by the following relaxed inequality:
643 *
644 * sum a[j] x[j] <= b~, (5)
645 * j in J
646 *
647 * where:
648 *
649 * b~ = sup(b - sum a[j] x[j] - sum a[j] x[j]) =
650 * j in Kp j in Kn
651 *
652 * = b - inf sum a[j] x[j] - inf sum a[j] x[j] = (6)
653 * j in Kp j in Kn
654 *
655 * = b - sum a[j] l[j] - sum a[j] u[j].
656 * j in Kp j in Kn
657 *
658 * Note that if lower bound l[j] (if j in Kp) or upper bound u[j]
659 * (if j in Kn) of some non-binary variable x[j] does not exist, then
660 * formally b = +oo, in which case further analysis is not performed.
661 *
662 * Let Bp = {j in J: a[j] > 0}, Bn = {j in J: a[j] < 0}. To make all
663 * the inequality coefficients in (5) positive, we replace all x[j] in
664 * Bn by their complementaries, substituting x[j] = 1 - x~[j] for all
665 * j in Bn, that gives:
666 *
667 * sum a[j] x[j] - sum a[j] x~[j] <= b~ - sum a[j]. (7)
668 * j in Bp j in Bn j in Bn
669 *
670 * This inequality is a relaxation of the original constraint (1), and
671 * it is a binary knapsack inequality. Writing it in the standard format
672 * we have:
673 *
674 * sum alfa[j] z[j] <= beta, (8)
675 * j in J
676 *
677 * where:
678 * ( + a[j], if j in Bp,
679 * alfa[j] = < (9)
680 * ( - a[j], if j in Bn,
681 *
682 * ( x[j], if j in Bp,
683 * z[j] = < (10)
684 * ( 1 - x[j], if j in Bn,
685 *
686 * beta = b~ - sum a[j]. (11)
687 * j in Bn
688 *
689 * In the inequality (8) all coefficients are positive, therefore, the
690 * packing relaxation to be found for this inequality is the following:
691 *
692 * sum z[j] <= 1. (12)
693 * j in P
694 *
695 * It is obvious that set P within J, which we would like to find, must
696 * satisfy to the following condition:
697 *
698 * alfa[j] + alfa[k] > beta + eps for all j, k in P, j != k, (13)
699 *
700 * where eps is an absolute tolerance for value of the linear form.
701 * Thus, it is natural to take P = {j: alpha[j] > (beta + eps) / 2}.
702 * Moreover, if in the equality (8) there exist coefficients alfa[k],
703 * for which alfa[k] <= (beta + eps) / 2, but which, nevertheless,
704 * satisfies to the condition (13) for all j in P, *one* corresponding
705 * variable z[k] (having, for example, maximal coefficient alfa[k]) can
706 * be included in set P, that allows increasing the number of binary
707 * variables in (12) by one.
708 *
709 * Once the set P has been built, for the inequality (12) we need to
710 * perform back substitution according to (10) in order to express it
711 * through the original binary variables. As the result of such back
712 * substitution the relaxed packing inequality get its final format (2),
713 * where Jp = J intersect Bp, and Jn = J intersect Bn. */
715 int npp_implied_packing(NPP *npp, NPPROW *row, int which,
716 NPPCOL *var[], char set[])
717 { struct elem *ptr, *e, *i, *k;
718 int len = 0;
719 double b, eps;
720 /* build inequality (3) */
721 if (which == 0)
722 { ptr = copy_form(npp, row, -1.0);
723 xassert(row->lb != -DBL_MAX);
724 b = - row->lb;
725 }
726 else if (which == 1)
727 { ptr = copy_form(npp, row, +1.0);
728 xassert(row->ub != +DBL_MAX);
729 b = + row->ub;
730 }
731 /* remove non-binary variables to build relaxed inequality (5);
732 compute its right-hand side b~ with formula (6) */
733 for (e = ptr; e != NULL; e = e->next)
734 { if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0))
735 { /* x[j] is non-binary variable */
736 if (e->aj > 0.0)
737 { if (e->xj->lb == -DBL_MAX) goto done;
738 b -= e->aj * e->xj->lb;
739 }
740 else /* e->aj < 0.0 */
741 { if (e->xj->ub == +DBL_MAX) goto done;
742 b -= e->aj * e->xj->ub;
743 }
744 /* a[j] = 0 means that variable x[j] is removed */
745 e->aj = 0.0;
746 }
747 }
748 /* substitute x[j] = 1 - x~[j] to build knapsack inequality (8);
749 compute its right-hand side beta with formula (11) */
750 for (e = ptr; e != NULL; e = e->next)
751 if (e->aj < 0.0) b -= e->aj;
752 /* if beta is close to zero, the knapsack inequality is either
753 infeasible or forcing inequality; this must never happen, so
754 we skip further analysis */
755 if (b < 1e-3) goto done;
756 /* build set P as well as sets Jp and Jn, and determine x[k] as
757 explained above in comments to the routine */
758 eps = 1e-3 + 1e-6 * b;
759 i = k = NULL;
760 for (e = ptr; e != NULL; e = e->next)
761 { /* note that alfa[j] = |a[j]| */
762 if (fabs(e->aj) > 0.5 * (b + eps))
763 { /* alfa[j] > (b + eps) / 2; include x[j] in set P, i.e. in
764 set Jp or Jn */
765 var[++len] = e->xj;
766 set[len] = (char)(e->aj > 0.0 ? 0 : 1);
767 /* alfa[i] = min alfa[j] over all j included in set P */
768 if (i == NULL || fabs(i->aj) > fabs(e->aj)) i = e;
769 }
770 else if (fabs(e->aj) >= 1e-3)
771 { /* alfa[k] = max alfa[j] over all j not included in set P;
772 we skip coefficient a[j] if it is close to zero to avoid
773 numerically unreliable results */
774 if (k == NULL || fabs(k->aj) < fabs(e->aj)) k = e;
775 }
776 }
777 /* if alfa[k] satisfies to condition (13) for all j in P, include
778 x[k] in P */
779 if (i != NULL && k != NULL && fabs(i->aj) + fabs(k->aj) > b + eps)
780 { var[++len] = k->xj;
781 set[len] = (char)(k->aj > 0.0 ? 0 : 1);
782 }
783 /* trivial packing inequality being redundant must never appear,
784 so we just ignore it */
785 if (len < 2) len = 0;
786 done: drop_form(npp, ptr);
787 return len;
788 }
790 /***********************************************************************
791 * NAME
792 *
793 * npp_is_covering - test if constraint is covering inequality
794 *
795 * SYNOPSIS
796 *
797 * #include "glpnpp.h"
798 * int npp_is_covering(NPP *npp, NPPROW *row);
799 *
800 * RETURNS
801 *
802 * If the specified row (constraint) is covering inequality (see below),
803 * the routine npp_is_covering returns non-zero. Otherwise, it returns
804 * zero.
805 *
806 * COVERING INEQUALITIES
807 *
808 * In canonical format the covering inequality is the following:
809 *
810 * sum x[j] >= 1, (1)
811 * j in J
812 *
813 * where all variables x[j] are binary. This inequality expresses the
814 * condition that in any integer feasible solution variables in set J
815 * cannot be all equal to zero at the same time, i.e. at least one
816 * variable must take non-zero (unity) value. W.l.o.g. it is assumed
817 * that |J| >= 2, because if J is empty, the inequality (1) is
818 * infeasible, and if |J| = 1, the inequality (1) is a forcing row.
819 *
820 * In general case the covering inequality may include original
821 * variables x[j] as well as their complements x~[j]:
822 *
823 * sum x[j] + sum x~[j] >= 1, (2)
824 * j in Jp j in Jn
825 *
826 * where Jp and Jn are not intersected. Therefore, using substitution
827 * x~[j] = 1 - x[j] gives the packing inequality in generalized format:
828 *
829 * sum x[j] - sum x[j] >= 1 - |Jn|. (3)
830 * j in Jp j in Jn
831 *
832 * (May note that the inequality (3) cuts off infeasible solutions,
833 * where x[j] = 0 for all j in Jp and x[j] = 1 for all j in Jn.)
834 *
835 * NOTE: If |J| = 2, the inequality (3) is equivalent to packing
836 * inequality (see the routine npp_is_packing). */
838 int npp_is_covering(NPP *npp, NPPROW *row)
839 { /* test if constraint is covering inequality */
840 NPPCOL *col;
841 NPPAIJ *aij;
842 int b;
843 xassert(npp == npp);
844 if (!(row->lb != -DBL_MAX && row->ub == +DBL_MAX))
845 return 0;
846 b = 1;
847 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
848 { col = aij->col;
849 if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
850 return 0;
851 if (aij->val == +1.0)
852 ;
853 else if (aij->val == -1.0)
854 b--;
855 else
856 return 0;
857 }
858 if (row->lb != (double)b) return 0;
859 return 1;
860 }
862 /***********************************************************************
863 * NAME
864 *
865 * npp_hidden_covering - identify hidden covering inequality
866 *
867 * SYNOPSIS
868 *
869 * #include "glpnpp.h"
870 * int npp_hidden_covering(NPP *npp, NPPROW *row);
871 *
872 * DESCRIPTION
873 *
874 * The routine npp_hidden_covering processes specified inequality
875 * constraint, which includes only binary variables, and the number of
876 * the variables is not less than three. If the original inequality is
877 * equivalent to a covering inequality (see below), the routine
878 * replaces it by the equivalent inequality. If the original constraint
879 * is double-sided inequality, it is replaced by a pair of single-sided
880 * inequalities, if necessary.
881 *
882 * RETURNS
883 *
884 * If the original inequality constraint was replaced by equivalent
885 * covering inequality, the routine npp_hidden_covering returns
886 * non-zero. Otherwise, it returns zero.
887 *
888 * PROBLEM TRANSFORMATION
889 *
890 * Consider an inequality constraint:
891 *
892 * sum a[j] x[j] >= b, (1)
893 * j in J
894 *
895 * where all variables x[j] are binary, and |J| >= 3. (In case of '<='
896 * inequality it can be transformed to '>=' format by multiplying both
897 * its sides by -1.)
898 *
899 * Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution
900 * x[j] = 1 - x~[j] for all j in Jn, we have:
901 *
902 * sum a[j] x[j] >= b ==>
903 * j in J
904 *
905 * sum a[j] x[j] + sum a[j] x[j] >= b ==>
906 * j in Jp j in Jn
907 *
908 * sum a[j] x[j] + sum a[j] (1 - x~[j]) >= b ==>
909 * j in Jp j in Jn
910 *
911 * sum m a[j] x[j] - sum a[j] x~[j] >= b - sum a[j].
912 * j in Jp j in Jn j in Jn
913 *
914 * Thus, meaning the transformation above, we can assume that in
915 * inequality (1) all coefficients a[j] are positive. Moreover, we can
916 * assume that b > 0, because otherwise the inequality (1) would be
917 * redundant (see the routine npp_analyze_row). It is then obvious that
918 * constraint (1) is equivalent to covering inequality only if:
919 *
920 * a[j] >= b, (2)
921 *
922 * for all j in J.
923 *
924 * Once the original inequality (1) is replaced by equivalent covering
925 * inequality, we need to perform back substitution x~[j] = 1 - x[j] for
926 * all j in Jn (see above).
927 *
928 * RECOVERING SOLUTION
929 *
930 * None needed. */
932 static int hidden_covering(NPP *npp, struct elem *ptr, double *_b)
933 { /* process inequality constraint: sum a[j] x[j] >= b;
934 0 - specified row is NOT hidden covering inequality;
935 1 - specified row is covering inequality;
936 2 - specified row is hidden covering inequality. */
937 struct elem *e;
938 int neg;
939 double b = *_b, eps;
940 xassert(npp == npp);
941 /* a[j] must be non-zero, x[j] must be binary, for all j in J */
942 for (e = ptr; e != NULL; e = e->next)
943 { xassert(e->aj != 0.0);
944 xassert(e->xj->is_int);
945 xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0);
946 }
947 /* check if the specified inequality constraint already has the
948 form of covering inequality */
949 neg = 0; /* neg is |Jn| */
950 for (e = ptr; e != NULL; e = e->next)
951 { if (e->aj == +1.0)
952 ;
953 else if (e->aj == -1.0)
954 neg++;
955 else
956 break;
957 }
958 if (e == NULL)
959 { /* all coefficients a[j] are +1 or -1; check rhs b */
960 if (b == (double)(1 - neg))
961 { /* it is covering inequality; no processing is needed */
962 return 1;
963 }
964 }
965 /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j]
966 positive; the result is a~[j] = |a[j]| and new rhs b */
967 for (e = ptr; e != NULL; e = e->next)
968 if (e->aj < 0) b -= e->aj;
969 /* now a[j] > 0 for all j in J (actually |a[j]| are used) */
970 /* if b <= 0, skip processing--this case must not appear */
971 if (b < 1e-3) return 0;
972 /* now a[j] > 0 for all j in J, and b > 0 */
973 /* the specified constraint is equivalent to covering inequality
974 iff a[j] >= b for all j in J */
975 eps = 1e-9 + 1e-12 * fabs(b);
976 for (e = ptr; e != NULL; e = e->next)
977 if (fabs(e->aj) < b - eps) return 0;
978 /* perform back substitution x~[j] = 1 - x[j] and construct the
979 final equivalent covering inequality in generalized format */
980 b = 1.0;
981 for (e = ptr; e != NULL; e = e->next)
982 { if (e->aj > 0.0)
983 e->aj = +1.0;
984 else /* e->aj < 0.0 */
985 e->aj = -1.0, b -= 1.0;
986 }
987 *_b = b;
988 return 2;
989 }
991 int npp_hidden_covering(NPP *npp, NPPROW *row)
992 { /* identify hidden covering inequality */
993 NPPROW *copy;
994 NPPAIJ *aij;
995 struct elem *ptr, *e;
996 int kase, ret, count = 0;
997 double b;
998 /* the row must be inequality constraint */
999 xassert(row->lb < row->ub);
1000 for (kase = 0; kase <= 1; kase++)
1001 { if (kase == 0)
1002 { /* process row lower bound */
1003 if (row->lb == -DBL_MAX) continue;
1004 ptr = copy_form(npp, row, +1.0);
1005 b = + row->lb;
1006 }
1007 else
1008 { /* process row upper bound */
1009 if (row->ub == +DBL_MAX) continue;
1010 ptr = copy_form(npp, row, -1.0);
1011 b = - row->ub;
1012 }
1013 /* now the inequality has the form "sum a[j] x[j] >= b" */
1014 ret = hidden_covering(npp, ptr, &b);
1015 xassert(0 <= ret && ret <= 2);
1016 if (kase == 1 && ret == 1 || ret == 2)
1017 { /* the original inequality has been identified as hidden
1018 covering inequality */
1019 count++;
1020 #ifdef GLP_DEBUG
1021 xprintf("Original constraint:\n");
1022 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1023 xprintf(" %+g x%d", aij->val, aij->col->j);
1024 if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb);
1025 if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub);
1026 xprintf("\n");
1027 xprintf("Equivalent covering inequality:\n");
1028 for (e = ptr; e != NULL; e = e->next)
1029 xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j);
1030 xprintf(", >= %g\n", b);
1031 #endif
1032 if (row->lb == -DBL_MAX || row->ub == +DBL_MAX)
1033 { /* the original row is single-sided inequality; no copy
1034 is needed */
1035 copy = NULL;
1036 }
1037 else
1038 { /* the original row is double-sided inequality; we need
1039 to create its copy for other bound before replacing it
1040 with the equivalent inequality */
1041 copy = npp_add_row(npp);
1042 if (kase == 0)
1043 { /* the copy is for upper bound */
1044 copy->lb = -DBL_MAX, copy->ub = row->ub;
1045 }
1046 else
1047 { /* the copy is for lower bound */
1048 copy->lb = row->lb, copy->ub = +DBL_MAX;
1049 }
1050 /* copy original row coefficients */
1051 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1052 npp_add_aij(npp, copy, aij->col, aij->val);
1053 }
1054 /* replace the original inequality by equivalent one */
1055 npp_erase_row(npp, row);
1056 row->lb = b, row->ub = +DBL_MAX;
1057 for (e = ptr; e != NULL; e = e->next)
1058 npp_add_aij(npp, row, e->xj, e->aj);
1059 /* continue processing upper bound for the copy */
1060 if (copy != NULL) row = copy;
1061 }
1062 drop_form(npp, ptr);
1063 }
1064 return count;
1065 }
1067 /***********************************************************************
1068 * NAME
1069 *
1070 * npp_is_partitioning - test if constraint is partitioning equality
1071 *
1072 * SYNOPSIS
1073 *
1074 * #include "glpnpp.h"
1075 * int npp_is_partitioning(NPP *npp, NPPROW *row);
1076 *
1077 * RETURNS
1078 *
1079 * If the specified row (constraint) is partitioning equality (see
1080 * below), the routine npp_is_partitioning returns non-zero. Otherwise,
1081 * it returns zero.
1082 *
1083 * PARTITIONING EQUALITIES
1084 *
1085 * In canonical format the partitioning equality is the following:
1086 *
1087 * sum x[j] = 1, (1)
1088 * j in J
1089 *
1090 * where all variables x[j] are binary. This equality expresses the
1091 * condition that in any integer feasible solution exactly one variable
1092 * in set J must take non-zero (unity) value while other variables must
1093 * be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because if
1094 * J is empty, the inequality (1) is infeasible, and if |J| = 1, the
1095 * inequality (1) is a fixing row.
1096 *
1097 * In general case the partitioning equality may include original
1098 * variables x[j] as well as their complements x~[j]:
1099 *
1100 * sum x[j] + sum x~[j] = 1, (2)
1101 * j in Jp j in Jn
1102 *
1103 * where Jp and Jn are not intersected. Therefore, using substitution
1104 * x~[j] = 1 - x[j] leads to the partitioning equality in generalized
1105 * format:
1106 *
1107 * sum x[j] - sum x[j] = 1 - |Jn|. (3)
1108 * j in Jp j in Jn */
1110 int npp_is_partitioning(NPP *npp, NPPROW *row)
1111 { /* test if constraint is partitioning equality */
1112 NPPCOL *col;
1113 NPPAIJ *aij;
1114 int b;
1115 xassert(npp == npp);
1116 if (row->lb != row->ub) return 0;
1117 b = 1;
1118 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1119 { col = aij->col;
1120 if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0))
1121 return 0;
1122 if (aij->val == +1.0)
1123 ;
1124 else if (aij->val == -1.0)
1125 b--;
1126 else
1127 return 0;
1128 }
1129 if (row->lb != (double)b) return 0;
1130 return 1;
1131 }
1133 /***********************************************************************
1134 * NAME
1135 *
1136 * npp_reduce_ineq_coef - reduce inequality constraint coefficients
1137 *
1138 * SYNOPSIS
1139 *
1140 * #include "glpnpp.h"
1141 * int npp_reduce_ineq_coef(NPP *npp, NPPROW *row);
1142 *
1143 * DESCRIPTION
1144 *
1145 * The routine npp_reduce_ineq_coef processes specified inequality
1146 * constraint attempting to replace it by an equivalent constraint,
1147 * where magnitude of coefficients at binary variables is smaller than
1148 * in the original constraint. If the inequality is double-sided, it is
1149 * replaced by a pair of single-sided inequalities, if necessary.
1150 *
1151 * RETURNS
1152 *
1153 * The routine npp_reduce_ineq_coef returns the number of coefficients
1154 * reduced.
1155 *
1156 * BACKGROUND
1157 *
1158 * Consider an inequality constraint:
1159 *
1160 * sum a[j] x[j] >= b. (1)
1161 * j in J
1162 *
1163 * (In case of '<=' inequality it can be transformed to '>=' format by
1164 * multiplying both its sides by -1.) Let x[k] be a binary variable;
1165 * other variables can be integer as well as continuous. We can write
1166 * constraint (1) as follows:
1167 *
1168 * a[k] x[k] + t[k] >= b, (2)
1169 *
1170 * where:
1171 *
1172 * t[k] = sum a[j] x[j]. (3)
1173 * j in J\{k}
1174 *
1175 * Since x[k] is binary, constraint (2) is equivalent to disjunction of
1176 * the following two constraints:
1177 *
1178 * x[k] = 0, t[k] >= b (4)
1179 *
1180 * OR
1181 *
1182 * x[k] = 1, t[k] >= b - a[k]. (5)
1183 *
1184 * Let also that for the partial sum t[k] be known some its implied
1185 * lower bound inf t[k].
1186 *
1187 * Case a[k] > 0. Let inf t[k] < b, since otherwise both constraints
1188 * (4) and (5) and therefore constraint (2) are redundant.
1189 * If inf t[k] > b - a[k], only constraint (5) is redundant, in which
1190 * case it can be replaced with the following redundant and therefore
1191 * equivalent constraint:
1192 *
1193 * t[k] >= b - a'[k] = inf t[k], (6)
1194 *
1195 * where:
1196 *
1197 * a'[k] = b - inf t[k]. (7)
1198 *
1199 * Thus, the original constraint (2) is equivalent to the following
1200 * constraint with coefficient at variable x[k] changed:
1201 *
1202 * a'[k] x[k] + t[k] >= b. (8)
1203 *
1204 * From inf t[k] < b it follows that a'[k] > 0, i.e. the coefficient
1205 * at x[k] keeps its sign. And from inf t[k] > b - a[k] it follows that
1206 * a'[k] < a[k], i.e. the coefficient reduces in magnitude.
1207 *
1208 * Case a[k] < 0. Let inf t[k] < b - a[k], since otherwise both
1209 * constraints (4) and (5) and therefore constraint (2) are redundant.
1210 * If inf t[k] > b, only constraint (4) is redundant, in which case it
1211 * can be replaced with the following redundant and therefore equivalent
1212 * constraint:
1213 *
1214 * t[k] >= b' = inf t[k]. (9)
1215 *
1216 * Rewriting constraint (5) as follows:
1217 *
1218 * t[k] >= b - a[k] = b' - a'[k], (10)
1219 *
1220 * where:
1221 *
1222 * a'[k] = a[k] + b' - b = a[k] + inf t[k] - b, (11)
1223 *
1224 * we can see that disjunction of constraint (9) and (10) is equivalent
1225 * to disjunction of constraint (4) and (5), from which it follows that
1226 * the original constraint (2) is equivalent to the following constraint
1227 * with both coefficient at variable x[k] and right-hand side changed:
1228 *
1229 * a'[k] x[k] + t[k] >= b'. (12)
1230 *
1231 * From inf t[k] < b - a[k] it follows that a'[k] < 0, i.e. the
1232 * coefficient at x[k] keeps its sign. And from inf t[k] > b it follows
1233 * that a'[k] > a[k], i.e. the coefficient reduces in magnitude.
1234 *
1235 * PROBLEM TRANSFORMATION
1236 *
1237 * In the routine npp_reduce_ineq_coef the following implied lower
1238 * bound of the partial sum (3) is used:
1239 *
1240 * inf t[k] = sum a[j] l[j] + sum a[j] u[j], (13)
1241 * j in Jp\{k} k in Jn\{k}
1242 *
1243 * where Jp = {j : a[j] > 0}, Jn = {j : a[j] < 0}, l[j] and u[j] are
1244 * lower and upper bounds, resp., of variable x[j].
1245 *
1246 * In order to compute inf t[k] more efficiently, the following formula,
1247 * which is equivalent to (13), is actually used:
1248 *
1249 * ( h - a[k] l[k] = h, if a[k] > 0,
1250 * inf t[k] = < (14)
1251 * ( h - a[k] u[k] = h - a[k], if a[k] < 0,
1252 *
1253 * where:
1254 *
1255 * h = sum a[j] l[j] + sum a[j] u[j] (15)
1256 * j in Jp j in Jn
1257 *
1258 * is the implied lower bound of row (1).
1259 *
1260 * Reduction of positive coefficient (a[k] > 0) does not change value
1261 * of h, since l[k] = 0. In case of reduction of negative coefficient
1262 * (a[k] < 0) from (11) it follows that:
1263 *
1264 * delta a[k] = a'[k] - a[k] = inf t[k] - b (> 0), (16)
1265 *
1266 * so new value of h (accounting that u[k] = 1) can be computed as
1267 * follows:
1268 *
1269 * h := h + delta a[k] = h + (inf t[k] - b). (17)
1270 *
1271 * RECOVERING SOLUTION
1272 *
1273 * None needed. */
1275 static int reduce_ineq_coef(NPP *npp, struct elem *ptr, double *_b)
1276 { /* process inequality constraint: sum a[j] x[j] >= b */
1277 /* returns: the number of coefficients reduced */
1278 struct elem *e;
1279 int count = 0;
1280 double h, inf_t, new_a, b = *_b;
1281 xassert(npp == npp);
1282 /* compute h; see (15) */
1283 h = 0.0;
1284 for (e = ptr; e != NULL; e = e->next)
1285 { if (e->aj > 0.0)
1286 { if (e->xj->lb == -DBL_MAX) goto done;
1287 h += e->aj * e->xj->lb;
1288 }
1289 else /* e->aj < 0.0 */
1290 { if (e->xj->ub == +DBL_MAX) goto done;
1291 h += e->aj * e->xj->ub;
1292 }
1293 }
1294 /* perform reduction of coefficients at binary variables */
1295 for (e = ptr; e != NULL; e = e->next)
1296 { /* skip non-binary variable */
1297 if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0))
1298 continue;
1299 if (e->aj > 0.0)
1300 { /* compute inf t[k]; see (14) */
1301 inf_t = h;
1302 if (b - e->aj < inf_t && inf_t < b)
1303 { /* compute reduced coefficient a'[k]; see (7) */
1304 new_a = b - inf_t;
1305 if (new_a >= +1e-3 &&
1306 e->aj - new_a >= 0.01 * (1.0 + e->aj))
1307 { /* accept a'[k] */
1308 #ifdef GLP_DEBUG
1309 xprintf("+");
1310 #endif
1311 e->aj = new_a;
1312 count++;
1313 }
1314 }
1315 }
1316 else /* e->aj < 0.0 */
1317 { /* compute inf t[k]; see (14) */
1318 inf_t = h - e->aj;
1319 if (b < inf_t && inf_t < b - e->aj)
1320 { /* compute reduced coefficient a'[k]; see (11) */
1321 new_a = e->aj + (inf_t - b);
1322 if (new_a <= -1e-3 &&
1323 new_a - e->aj >= 0.01 * (1.0 - e->aj))
1324 { /* accept a'[k] */
1325 #ifdef GLP_DEBUG
1326 xprintf("-");
1327 #endif
1328 e->aj = new_a;
1329 /* update h; see (17) */
1330 h += (inf_t - b);
1331 /* compute b'; see (9) */
1332 b = inf_t;
1333 count++;
1334 }
1335 }
1336 }
1337 }
1338 *_b = b;
1339 done: return count;
1340 }
1342 int npp_reduce_ineq_coef(NPP *npp, NPPROW *row)
1343 { /* reduce inequality constraint coefficients */
1344 NPPROW *copy;
1345 NPPAIJ *aij;
1346 struct elem *ptr, *e;
1347 int kase, count[2];
1348 double b;
1349 /* the row must be inequality constraint */
1350 xassert(row->lb < row->ub);
1351 count[0] = count[1] = 0;
1352 for (kase = 0; kase <= 1; kase++)
1353 { if (kase == 0)
1354 { /* process row lower bound */
1355 if (row->lb == -DBL_MAX) continue;
1356 #ifdef GLP_DEBUG
1357 xprintf("L");
1358 #endif
1359 ptr = copy_form(npp, row, +1.0);
1360 b = + row->lb;
1361 }
1362 else
1363 { /* process row upper bound */
1364 if (row->ub == +DBL_MAX) continue;
1365 #ifdef GLP_DEBUG
1366 xprintf("U");
1367 #endif
1368 ptr = copy_form(npp, row, -1.0);
1369 b = - row->ub;
1370 }
1371 /* now the inequality has the form "sum a[j] x[j] >= b" */
1372 count[kase] = reduce_ineq_coef(npp, ptr, &b);
1373 if (count[kase] > 0)
1374 { /* the original inequality has been replaced by equivalent
1375 one with coefficients reduced */
1376 if (row->lb == -DBL_MAX || row->ub == +DBL_MAX)
1377 { /* the original row is single-sided inequality; no copy
1378 is needed */
1379 copy = NULL;
1380 }
1381 else
1382 { /* the original row is double-sided inequality; we need
1383 to create its copy for other bound before replacing it
1384 with the equivalent inequality */
1385 #ifdef GLP_DEBUG
1386 xprintf("*");
1387 #endif
1388 copy = npp_add_row(npp);
1389 if (kase == 0)
1390 { /* the copy is for upper bound */
1391 copy->lb = -DBL_MAX, copy->ub = row->ub;
1392 }
1393 else
1394 { /* the copy is for lower bound */
1395 copy->lb = row->lb, copy->ub = +DBL_MAX;
1396 }
1397 /* copy original row coefficients */
1398 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1399 npp_add_aij(npp, copy, aij->col, aij->val);
1400 }
1401 /* replace the original inequality by equivalent one */
1402 npp_erase_row(npp, row);
1403 row->lb = b, row->ub = +DBL_MAX;
1404 for (e = ptr; e != NULL; e = e->next)
1405 npp_add_aij(npp, row, e->xj, e->aj);
1406 /* continue processing upper bound for the copy */
1407 if (copy != NULL) row = copy;
1408 }
1409 drop_form(npp, ptr);
1410 }
1411 return count[0] + count[1];
1412 }
1414 /* eof */