1 /* glpnpp03.c */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 *
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7 * 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9 * E-mail: <mao@gnu.org>.
10 *
11 * GLPK is free software: you can redistribute it and/or modify it
12 * under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19 * License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23 ***********************************************************************/
24
25 #include "glpnpp.h"
26
27 /***********************************************************************
28 * NAME
29 *
30 * npp_empty_row - process empty row
31 *
32 * SYNOPSIS
33 *
34 * #include "glpnpp.h"
35 * int npp_empty_row(NPP *npp, NPPROW *p);
36 *
37 * DESCRIPTION
38 *
39 * The routine npp_empty_row processes row p, which is empty, i.e.
40 * coefficients at all columns in this row are zero:
41 *
42 * L[p] <= sum 0 x[j] <= U[p], (1)
43 *
44 * where L[p] <= U[p].
45 *
46 * RETURNS
47 *
48 * 0 - success;
49 *
50 * 1 - problem has no primal feasible solution.
51 *
52 * PROBLEM TRANSFORMATION
53 *
54 * If the following conditions hold:
55 *
56 * L[p] <= +eps, U[p] >= -eps, (2)
57 *
58 * where eps is an absolute tolerance for row value, the row p is
59 * redundant. In this case it can be replaced by equivalent redundant
60 * row, which is free (unbounded), and then removed from the problem.
61 * Otherwise, the row p is infeasible and, thus, the problem has no
62 * primal feasible solution.
63 *
64 * RECOVERING BASIC SOLUTION
65 *
66 * See the routine npp_free_row.
67 *
68 * RECOVERING INTERIOR-POINT SOLUTION
69 *
70 * See the routine npp_free_row.
71 *
72 * RECOVERING MIP SOLUTION
73 *
74 * None needed. */
75
npp_empty_row(NPP * npp,NPPROW * p)76 int npp_empty_row(NPP *npp, NPPROW *p)
77 { /* process empty row */
78 double eps = 1e-3;
79 /* the row must be empty */
80 xassert(p->ptr == NULL);
81 /* check primal feasibility */
82 if (p->lb > +eps || p->ub < -eps)
83 return 1;
84 /* replace the row by equivalent free (unbounded) row */
85 p->lb = -DBL_MAX, p->ub = +DBL_MAX;
86 /* and process it */
87 npp_free_row(npp, p);
88 return 0;
89 }
90
91 /***********************************************************************
92 * NAME
93 *
94 * npp_empty_col - process empty column
95 *
96 * SYNOPSIS
97 *
98 * #include "glpnpp.h"
99 * int npp_empty_col(NPP *npp, NPPCOL *q);
100 *
101 * DESCRIPTION
102 *
103 * The routine npp_empty_col processes column q:
104 *
105 * l[q] <= x[q] <= u[q], (1)
106 *
107 * where l[q] <= u[q], which is empty, i.e. has zero coefficients in
108 * all constraint rows.
109 *
110 * RETURNS
111 *
112 * 0 - success;
113 *
114 * 1 - problem has no dual feasible solution.
115 *
116 * PROBLEM TRANSFORMATION
117 *
118 * The row of the dual system corresponding to the empty column is the
119 * following:
120 *
121 * sum 0 pi[i] + lambda[q] = c[q], (2)
122 * i
123 *
124 * from which it follows that:
125 *
126 * lambda[q] = c[q]. (3)
127 *
128 * If the following condition holds:
129 *
130 * c[q] < - eps, (4)
131 *
132 * where eps is an absolute tolerance for column multiplier, the lower
133 * column bound l[q] must be active to provide dual feasibility (note
134 * that being preprocessed the problem is always minimization). In this
135 * case the column can be fixed on its lower bound and removed from the
136 * problem (if the column is integral, its bounds are also assumed to
137 * be integral). And if the column has no lower bound (l[q] = -oo), the
138 * problem has no dual feasible solution.
139 *
140 * If the following condition holds:
141 *
142 * c[q] > + eps, (5)
143 *
144 * the upper column bound u[q] must be active to provide dual
145 * feasibility. In this case the column can be fixed on its upper bound
146 * and removed from the problem. And if the column has no upper bound
147 * (u[q] = +oo), the problem has no dual feasible solution.
148 *
149 * Finally, if the following condition holds:
150 *
151 * - eps <= c[q] <= +eps, (6)
152 *
153 * dual feasibility does not depend on a particular value of column q.
154 * In this case the column can be fixed either on its lower bound (if
155 * l[q] > -oo) or on its upper bound (if u[q] < +oo) or at zero (if the
156 * column is unbounded) and then removed from the problem.
157 *
158 * RECOVERING BASIC SOLUTION
159 *
160 * See the routine npp_fixed_col. Having been recovered the column
161 * is assigned status GLP_NS. However, if actually it is not fixed
162 * (l[q] < u[q]), its status should be changed to GLP_NL, GLP_NU, or
163 * GLP_NF depending on which bound it was fixed on transformation stage.
164 *
165 * RECOVERING INTERIOR-POINT SOLUTION
166 *
167 * See the routine npp_fixed_col.
168 *
169 * RECOVERING MIP SOLUTION
170 *
171 * See the routine npp_fixed_col. */
172
173 struct empty_col
174 { /* empty column */
175 int q;
176 /* column reference number */
177 char stat;
178 /* status in basic solution */
179 };
180
181 static int rcv_empty_col(NPP *npp, void *info);
182
npp_empty_col(NPP * npp,NPPCOL * q)183 int npp_empty_col(NPP *npp, NPPCOL *q)
184 { /* process empty column */
185 struct empty_col *info;
186 double eps = 1e-3;
187 /* the column must be empty */
188 xassert(q->ptr == NULL);
189 /* check dual feasibility */
190 if (q->coef > +eps && q->lb == -DBL_MAX)
191 return 1;
192 if (q->coef < -eps && q->ub == +DBL_MAX)
193 return 1;
194 /* create transformation stack entry */
195 info = npp_push_tse(npp,
196 rcv_empty_col, sizeof(struct empty_col));
197 info->q = q->j;
198 /* fix the column */
199 if (q->lb == -DBL_MAX && q->ub == +DBL_MAX)
200 { /* free column */
201 info->stat = GLP_NF;
202 q->lb = q->ub = 0.0;
203 }
204 else if (q->ub == +DBL_MAX)
205 lo: { /* column with lower bound */
206 info->stat = GLP_NL;
207 q->ub = q->lb;
208 }
209 else if (q->lb == -DBL_MAX)
210 up: { /* column with upper bound */
211 info->stat = GLP_NU;
212 q->lb = q->ub;
213 }
214 else if (q->lb != q->ub)
215 { /* double-bounded column */
216 if (q->coef >= +DBL_EPSILON) goto lo;
217 if (q->coef <= -DBL_EPSILON) goto up;
218 if (fabs(q->lb) <= fabs(q->ub)) goto lo; else goto up;
219 }
220 else
221 { /* fixed column */
222 info->stat = GLP_NS;
223 }
224 /* process fixed column */
225 npp_fixed_col(npp, q);
226 return 0;
227 }
228
rcv_empty_col(NPP * npp,void * _info)229 static int rcv_empty_col(NPP *npp, void *_info)
230 { /* recover empty column */
231 struct empty_col *info = _info;
232 if (npp->sol == GLP_SOL)
233 npp->c_stat[info->q] = info->stat;
234 return 0;
235 }
236
237 /***********************************************************************
238 * NAME
239 *
240 * npp_implied_value - process implied column value
241 *
242 * SYNOPSIS
243 *
244 * #include "glpnpp.h"
245 * int npp_implied_value(NPP *npp, NPPCOL *q, double s);
246 *
247 * DESCRIPTION
248 *
249 * For column q:
250 *
251 * l[q] <= x[q] <= u[q], (1)
252 *
253 * where l[q] < u[q], the routine npp_implied_value processes its
254 * implied value s[q]. If this implied value satisfies to the current
255 * column bounds and integrality condition, the routine fixes column q
256 * at the given point. Note that the column is kept in the problem in
257 * any case.
258 *
259 * RETURNS
260 *
261 * 0 - column has been fixed;
262 *
263 * 1 - implied value violates to current column bounds;
264 *
265 * 2 - implied value violates integrality condition.
266 *
267 * ALGORITHM
268 *
269 * Implied column value s[q] satisfies to the current column bounds if
270 * the following condition holds:
271 *
272 * l[q] - eps <= s[q] <= u[q] + eps, (2)
273 *
274 * where eps is an absolute tolerance for column value. If the column
275 * is integral, the following condition also must hold:
276 *
277 * |s[q] - floor(s[q]+0.5)| <= eps, (3)
278 *
279 * where floor(s[q]+0.5) is the nearest integer to s[q].
280 *
281 * If both condition (2) and (3) are satisfied, the column can be fixed
282 * at the value s[q], or, if it is integral, at floor(s[q]+0.5).
283 * Otherwise, if s[q] violates (2) or (3), the problem has no feasible
284 * solution.
285 *
286 * Note: If s[q] is close to l[q] or u[q], it seems to be reasonable to
287 * fix the column at its lower or upper bound, resp. rather than at the
288 * implied value. */
289
npp_implied_value(NPP * npp,NPPCOL * q,double s)290 int npp_implied_value(NPP *npp, NPPCOL *q, double s)
291 { /* process implied column value */
292 double eps, nint;
293 xassert(npp == npp);
294 /* column must not be fixed */
295 xassert(q->lb < q->ub);
296 /* check integrality */
297 if (q->is_int)
298 { nint = floor(s + 0.5);
299 if (fabs(s - nint) <= 1e-5)
300 s = nint;
301 else
302 return 2;
303 }
304 /* check current column lower bound */
305 if (q->lb != -DBL_MAX)
306 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb));
307 if (s < q->lb - eps) return 1;
308 /* if s[q] is close to l[q], fix column at its lower bound
309 rather than at the implied value */
310 if (s < q->lb + 1e-3 * eps)
311 { q->ub = q->lb;
312 return 0;
313 }
314 }
315 /* check current column upper bound */
316 if (q->ub != +DBL_MAX)
317 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub));
318 if (s > q->ub + eps) return 1;
319 /* if s[q] is close to u[q], fix column at its upper bound
320 rather than at the implied value */
321 if (s > q->ub - 1e-3 * eps)
322 { q->lb = q->ub;
323 return 0;
324 }
325 }
326 /* fix column at the implied value */
327 q->lb = q->ub = s;
328 return 0;
329 }
330
331 /***********************************************************************
332 * NAME
333 *
334 * npp_eq_singlet - process row singleton (equality constraint)
335 *
336 * SYNOPSIS
337 *
338 * #include "glpnpp.h"
339 * int npp_eq_singlet(NPP *npp, NPPROW *p);
340 *
341 * DESCRIPTION
342 *
343 * The routine npp_eq_singlet processes row p, which is equiality
344 * constraint having the only non-zero coefficient:
345 *
346 * a[p,q] x[q] = b. (1)
347 *
348 * RETURNS
349 *
350 * 0 - success;
351 *
352 * 1 - problem has no primal feasible solution;
353 *
354 * 2 - problem has no integer feasible solution.
355 *
356 * PROBLEM TRANSFORMATION
357 *
358 * The equality constraint defines implied value of column q:
359 *
360 * x[q] = s[q] = b / a[p,q]. (2)
361 *
362 * If the implied value s[q] satisfies to the column bounds (see the
363 * routine npp_implied_value), the column can be fixed at s[q] and
364 * removed from the problem. In this case row p becomes redundant, so
365 * it can be replaced by equivalent free row and also removed from the
366 * problem.
367 *
368 * Note that the routine removes from the problem only row p. Column q
369 * becomes fixed, however, it is kept in the problem.
370 *
371 * RECOVERING BASIC SOLUTION
372 *
373 * In solution to the original problem row p is assigned status GLP_NS
374 * (active equality constraint), and column q is assigned status GLP_BS
375 * (basic column).
376 *
377 * Multiplier for row p can be computed as follows. In the dual system
378 * of the original problem column q corresponds to the following row:
379 *
380 * sum a[i,q] pi[i] + lambda[q] = c[q] ==>
381 * i
382 *
383 * sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q].
384 * i!=p
385 *
386 * Therefore:
387 *
388 * 1
389 * pi[p] = ------ (c[q] - lambda[q] - sum a[i,q] pi[i]), (3)
390 * a[p,q] i!=q
391 *
392 * where lambda[q] = 0 (since column[q] is basic), and pi[i] for all
393 * i != p are known in solution to the transformed problem.
394 *
395 * Value of column q in solution to the original problem is assigned
396 * its implied value s[q].
397 *
398 * RECOVERING INTERIOR-POINT SOLUTION
399 *
400 * Multiplier for row p is computed with formula (3). Value of column
401 * q is assigned its implied value s[q].
402 *
403 * RECOVERING MIP SOLUTION
404 *
405 * Value of column q is assigned its implied value s[q]. */
406
407 struct eq_singlet
408 { /* row singleton (equality constraint) */
409 int p;
410 /* row reference number */
411 int q;
412 /* column reference number */
413 double apq;
414 /* constraint coefficient a[p,q] */
415 double c;
416 /* objective coefficient at x[q] */
417 NPPLFE *ptr;
418 /* list of non-zero coefficients a[i,q], i != p */
419 };
420
421 static int rcv_eq_singlet(NPP *npp, void *info);
422
npp_eq_singlet(NPP * npp,NPPROW * p)423 int npp_eq_singlet(NPP *npp, NPPROW *p)
424 { /* process row singleton (equality constraint) */
425 struct eq_singlet *info;
426 NPPCOL *q;
427 NPPAIJ *aij;
428 NPPLFE *lfe;
429 int ret;
430 double s;
431 /* the row must be singleton equality constraint */
432 xassert(p->lb == p->ub);
433 xassert(p->ptr != NULL && p->ptr->r_next == NULL);
434 /* compute and process implied column value */
435 aij = p->ptr;
436 q = aij->col;
437 s = p->lb / aij->val;
438 ret = npp_implied_value(npp, q, s);
439 xassert(0 <= ret && ret <= 2);
440 if (ret != 0) return ret;
441 /* create transformation stack entry */
442 info = npp_push_tse(npp,
443 rcv_eq_singlet, sizeof(struct eq_singlet));
444 info->p = p->i;
445 info->q = q->j;
446 info->apq = aij->val;
447 info->c = q->coef;
448 info->ptr = NULL;
449 /* save column coefficients a[i,q], i != p (not needed for MIP
450 solution) */
451 if (npp->sol != GLP_MIP)
452 { for (aij = q->ptr; aij != NULL; aij = aij->c_next)
453 { if (aij->row == p) continue; /* skip a[p,q] */
454 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
455 lfe->ref = aij->row->i;
456 lfe->val = aij->val;
457 lfe->next = info->ptr;
458 info->ptr = lfe;
459 }
460 }
461 /* remove the row from the problem */
462 npp_del_row(npp, p);
463 return 0;
464 }
465
rcv_eq_singlet(NPP * npp,void * _info)466 static int rcv_eq_singlet(NPP *npp, void *_info)
467 { /* recover row singleton (equality constraint) */
468 struct eq_singlet *info = _info;
469 NPPLFE *lfe;
470 double temp;
471 if (npp->sol == GLP_SOL)
472 { /* column q must be already recovered as GLP_NS */
473 if (npp->c_stat[info->q] != GLP_NS)
474 { npp_error();
475 return 1;
476 }
477 npp->r_stat[info->p] = GLP_NS;
478 npp->c_stat[info->q] = GLP_BS;
479 }
480 if (npp->sol != GLP_MIP)
481 { /* compute multiplier for row p with formula (3) */
482 temp = info->c;
483 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
484 temp -= lfe->val * npp->r_pi[lfe->ref];
485 npp->r_pi[info->p] = temp / info->apq;
486 }
487 return 0;
488 }
489
490 /***********************************************************************
491 * NAME
492 *
493 * npp_implied_lower - process implied column lower bound
494 *
495 * SYNOPSIS
496 *
497 * #include "glpnpp.h"
498 * int npp_implied_lower(NPP *npp, NPPCOL *q, double l);
499 *
500 * DESCRIPTION
501 *
502 * For column q:
503 *
504 * l[q] <= x[q] <= u[q], (1)
505 *
506 * where l[q] < u[q], the routine npp_implied_lower processes its
507 * implied lower bound l'[q]. As the result the current column lower
508 * bound may increase. Note that the column is kept in the problem in
509 * any case.
510 *
511 * RETURNS
512 *
513 * 0 - current column lower bound has not changed;
514 *
515 * 1 - current column lower bound has changed, but not significantly;
516 *
517 * 2 - current column lower bound has significantly changed;
518 *
519 * 3 - column has been fixed on its upper bound;
520 *
521 * 4 - implied lower bound violates current column upper bound.
522 *
523 * ALGORITHM
524 *
525 * If column q is integral, before processing its implied lower bound
526 * should be rounded up:
527 *
528 * ( floor(l'[q]+0.5), if |l'[q] - floor(l'[q]+0.5)| <= eps
529 * l'[q] := < (2)
530 * ( ceil(l'[q]), otherwise
531 *
532 * where floor(l'[q]+0.5) is the nearest integer to l'[q], ceil(l'[q])
533 * is smallest integer not less than l'[q], and eps is an absolute
534 * tolerance for column value.
535 *
536 * Processing implied column lower bound l'[q] includes the following
537 * cases:
538 *
539 * 1) if l'[q] < l[q] + eps, implied lower bound is redundant;
540 *
541 * 2) if l[q] + eps <= l[q] <= u[q] + eps, current column lower bound
542 * l[q] can be strengthened by replacing it with l'[q]. If in this
543 * case new column lower bound becomes close to current column upper
544 * bound u[q], the column can be fixed on its upper bound;
545 *
546 * 3) if l'[q] > u[q] + eps, implied lower bound violates current
547 * column upper bound u[q], in which case the problem has no primal
548 * feasible solution. */
549
npp_implied_lower(NPP * npp,NPPCOL * q,double l)550 int npp_implied_lower(NPP *npp, NPPCOL *q, double l)
551 { /* process implied column lower bound */
552 int ret;
553 double eps, nint;
554 xassert(npp == npp);
555 /* column must not be fixed */
556 xassert(q->lb < q->ub);
557 /* implied lower bound must be finite */
558 xassert(l != -DBL_MAX);
559 /* if column is integral, round up l'[q] */
560 if (q->is_int)
561 { nint = floor(l + 0.5);
562 if (fabs(l - nint) <= 1e-5)
563 l = nint;
564 else
565 l = ceil(l);
566 }
567 /* check current column lower bound */
568 if (q->lb != -DBL_MAX)
569 { eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->lb));
570 if (l < q->lb + eps)
571 { ret = 0; /* redundant */
572 goto done;
573 }
574 }
575 /* check current column upper bound */
576 if (q->ub != +DBL_MAX)
577 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub));
578 if (l > q->ub + eps)
579 { ret = 4; /* infeasible */
580 goto done;
581 }
582 /* if l'[q] is close to u[q], fix column at its upper bound */
583 if (l > q->ub - 1e-3 * eps)
584 { q->lb = q->ub;
585 ret = 3; /* fixed */
586 goto done;
587 }
588 }
589 /* check if column lower bound changes significantly */
590 if (q->lb == -DBL_MAX)
591 ret = 2; /* significantly */
592 else if (q->is_int && l > q->lb + 0.5)
593 ret = 2; /* significantly */
594 else if (l > q->lb + 0.30 * (1.0 + fabs(q->lb)))
595 ret = 2; /* significantly */
596 else
597 ret = 1; /* not significantly */
598 /* set new column lower bound */
599 q->lb = l;
600 done: return ret;
601 }
602
603 /***********************************************************************
604 * NAME
605 *
606 * npp_implied_upper - process implied column upper bound
607 *
608 * SYNOPSIS
609 *
610 * #include "glpnpp.h"
611 * int npp_implied_upper(NPP *npp, NPPCOL *q, double u);
612 *
613 * DESCRIPTION
614 *
615 * For column q:
616 *
617 * l[q] <= x[q] <= u[q], (1)
618 *
619 * where l[q] < u[q], the routine npp_implied_upper processes its
620 * implied upper bound u'[q]. As the result the current column upper
621 * bound may decrease. Note that the column is kept in the problem in
622 * any case.
623 *
624 * RETURNS
625 *
626 * 0 - current column upper bound has not changed;
627 *
628 * 1 - current column upper bound has changed, but not significantly;
629 *
630 * 2 - current column upper bound has significantly changed;
631 *
632 * 3 - column has been fixed on its lower bound;
633 *
634 * 4 - implied upper bound violates current column lower bound.
635 *
636 * ALGORITHM
637 *
638 * If column q is integral, before processing its implied upper bound
639 * should be rounded down:
640 *
641 * ( floor(u'[q]+0.5), if |u'[q] - floor(l'[q]+0.5)| <= eps
642 * u'[q] := < (2)
643 * ( floor(l'[q]), otherwise
644 *
645 * where floor(u'[q]+0.5) is the nearest integer to u'[q],
646 * floor(u'[q]) is largest integer not greater than u'[q], and eps is
647 * an absolute tolerance for column value.
648 *
649 * Processing implied column upper bound u'[q] includes the following
650 * cases:
651 *
652 * 1) if u'[q] > u[q] - eps, implied upper bound is redundant;
653 *
654 * 2) if l[q] - eps <= u[q] <= u[q] - eps, current column upper bound
655 * u[q] can be strengthened by replacing it with u'[q]. If in this
656 * case new column upper bound becomes close to current column lower
657 * bound, the column can be fixed on its lower bound;
658 *
659 * 3) if u'[q] < l[q] - eps, implied upper bound violates current
660 * column lower bound l[q], in which case the problem has no primal
661 * feasible solution. */
662
npp_implied_upper(NPP * npp,NPPCOL * q,double u)663 int npp_implied_upper(NPP *npp, NPPCOL *q, double u)
664 { int ret;
665 double eps, nint;
666 xassert(npp == npp);
667 /* column must not be fixed */
668 xassert(q->lb < q->ub);
669 /* implied upper bound must be finite */
670 xassert(u != +DBL_MAX);
671 /* if column is integral, round down u'[q] */
672 if (q->is_int)
673 { nint = floor(u + 0.5);
674 if (fabs(u - nint) <= 1e-5)
675 u = nint;
676 else
677 u = floor(u);
678 }
679 /* check current column upper bound */
680 if (q->ub != +DBL_MAX)
681 { eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->ub));
682 if (u > q->ub - eps)
683 { ret = 0; /* redundant */
684 goto done;
685 }
686 }
687 /* check current column lower bound */
688 if (q->lb != -DBL_MAX)
689 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb));
690 if (u < q->lb - eps)
691 { ret = 4; /* infeasible */
692 goto done;
693 }
694 /* if u'[q] is close to l[q], fix column at its lower bound */
695 if (u < q->lb + 1e-3 * eps)
696 { q->ub = q->lb;
697 ret = 3; /* fixed */
698 goto done;
699 }
700 }
701 /* check if column upper bound changes significantly */
702 if (q->ub == +DBL_MAX)
703 ret = 2; /* significantly */
704 else if (q->is_int && u < q->ub - 0.5)
705 ret = 2; /* significantly */
706 else if (u < q->ub - 0.30 * (1.0 + fabs(q->ub)))
707 ret = 2; /* significantly */
708 else
709 ret = 1; /* not significantly */
710 /* set new column upper bound */
711 q->ub = u;
712 done: return ret;
713 }
714
715 /***********************************************************************
716 * NAME
717 *
718 * npp_ineq_singlet - process row singleton (inequality constraint)
719 *
720 * SYNOPSIS
721 *
722 * #include "glpnpp.h"
723 * int npp_ineq_singlet(NPP *npp, NPPROW *p);
724 *
725 * DESCRIPTION
726 *
727 * The routine npp_ineq_singlet processes row p, which is inequality
728 * constraint having the only non-zero coefficient:
729 *
730 * L[p] <= a[p,q] * x[q] <= U[p], (1)
731 *
732 * where L[p] < U[p], L[p] > -oo and/or U[p] < +oo.
733 *
734 * RETURNS
735 *
736 * 0 - current column bounds have not changed;
737 *
738 * 1 - current column bounds have changed, but not significantly;
739 *
740 * 2 - current column bounds have significantly changed;
741 *
742 * 3 - column has been fixed on its lower or upper bound;
743 *
744 * 4 - problem has no primal feasible solution.
745 *
746 * PROBLEM TRANSFORMATION
747 *
748 * Inequality constraint (1) defines implied bounds of column q:
749 *
750 * ( L[p] / a[p,q], if a[p,q] > 0
751 * l'[q] = < (2)
752 * ( U[p] / a[p,q], if a[p,q] < 0
753 *
754 * ( U[p] / a[p,q], if a[p,q] > 0
755 * u'[q] = < (3)
756 * ( L[p] / a[p,q], if a[p,q] < 0
757 *
758 * If these implied bounds do not violate current bounds of column q:
759 *
760 * l[q] <= x[q] <= u[q], (4)
761 *
762 * they can be used to strengthen the current column bounds:
763 *
764 * l[q] := max(l[q], l'[q]), (5)
765 *
766 * u[q] := min(u[q], u'[q]). (6)
767 *
768 * (See the routines npp_implied_lower and npp_implied_upper.)
769 *
770 * Once bounds of row p (1) have been carried over column q, the row
771 * becomes redundant, so it can be replaced by equivalent free row and
772 * removed from the problem.
773 *
774 * Note that the routine removes from the problem only row p. Column q,
775 * even it has been fixed, is kept in the problem.
776 *
777 * RECOVERING BASIC SOLUTION
778 *
779 * Note that the row in the dual system corresponding to column q is
780 * the following:
781 *
782 * sum a[i,q] pi[i] + lambda[q] = c[q] ==>
783 * i
784 * (7)
785 * sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q],
786 * i!=p
787 *
788 * where pi[i] for all i != p are known in solution to the transformed
789 * problem. Row p does not exist in the transformed problem, so it has
790 * zero multiplier there. This allows computing multiplier for column q
791 * in solution to the transformed problem:
792 *
793 * lambda~[q] = c[q] - sum a[i,q] pi[i]. (8)
794 * i!=p
795 *
796 * Let in solution to the transformed problem column q be non-basic
797 * with lower bound active (GLP_NL, lambda~[q] >= 0), and this lower
798 * bound be implied one l'[q]. From the original problem's standpoint
799 * this then means that actually the original column lower bound l[q]
800 * is inactive, and active is that row bound L[p] or U[p] that defines
801 * the implied bound l'[q] (2). In this case in solution to the
802 * original problem column q is assigned status GLP_BS while row p is
803 * assigned status GLP_NL (if a[p,q] > 0) or GLP_NU (if a[p,q] < 0).
804 * Since now column q is basic, its multiplier lambda[q] is zero. This
805 * allows using (7) and (8) to find multiplier for row p in solution to
806 * the original problem:
807 *
808 * 1
809 * pi[p] = ------ (c[q] - sum a[i,q] pi[i]) = lambda~[q] / a[p,q] (9)
810 * a[p,q] i!=p
811 *
812 * Now let in solution to the transformed problem column q be non-basic
813 * with upper bound active (GLP_NU, lambda~[q] <= 0), and this upper
814 * bound be implied one u'[q]. As in the previous case this then means
815 * that from the original problem's standpoint actually the original
816 * column upper bound u[q] is inactive, and active is that row bound
817 * L[p] or U[p] that defines the implied bound u'[q] (3). In this case
818 * in solution to the original problem column q is assigned status
819 * GLP_BS, row p is assigned status GLP_NU (if a[p,q] > 0) or GLP_NL
820 * (if a[p,q] < 0), and its multiplier is computed with formula (9).
821 *
822 * Strengthening bounds of column q according to (5) and (6) may make
823 * it fixed. Thus, if in solution to the transformed problem column q is
824 * non-basic and fixed (GLP_NS), we can suppose that if lambda~[q] > 0,
825 * column q has active lower bound (GLP_NL), and if lambda~[q] < 0,
826 * column q has active upper bound (GLP_NU), reducing this case to two
827 * previous ones. If, however, lambda~[q] is close to zero or
828 * corresponding bound of row p does not exist (this may happen if
829 * lambda~[q] has wrong sign due to round-off errors, in which case it
830 * is expected to be close to zero, since solution is assumed to be dual
831 * feasible), column q can be assigned status GLP_BS (basic), and row p
832 * can be made active on its existing bound. In the latter case row
833 * multiplier pi[p] computed with formula (9) will be also close to
834 * zero, and dual feasibility will be kept.
835 *
836 * In all other cases, namely, if in solution to the transformed
837 * problem column q is basic (GLP_BS), or non-basic with original lower
838 * bound l[q] active (GLP_NL), or non-basic with original upper bound
839 * u[q] active (GLP_NU), constraint (1) is inactive. So in solution to
840 * the original problem status of column q remains unchanged, row p is
841 * assigned status GLP_BS, and its multiplier pi[p] is assigned zero
842 * value.
843 *
844 * RECOVERING INTERIOR-POINT SOLUTION
845 *
846 * First, value of multiplier for column q in solution to the original
847 * problem is computed with formula (8). If lambda~[q] > 0 and column q
848 * has implied lower bound, or if lambda~[q] < 0 and column q has
849 * implied upper bound, this means that from the original problem's
850 * standpoint actually row p has corresponding active bound, in which
851 * case its multiplier pi[p] is computed with formula (9). In other
852 * cases, when the sign of lambda~[q] corresponds to original bound of
853 * column q, or when lambda~[q] =~ 0, value of row multiplier pi[p] is
854 * assigned zero value.
855 *
856 * RECOVERING MIP SOLUTION
857 *
858 * None needed. */
859
860 struct ineq_singlet
861 { /* row singleton (inequality constraint) */
862 int p;
863 /* row reference number */
864 int q;
865 /* column reference number */
866 double apq;
867 /* constraint coefficient a[p,q] */
868 double c;
869 /* objective coefficient at x[q] */
870 double lb;
871 /* row lower bound */
872 double ub;
873 /* row upper bound */
874 char lb_changed;
875 /* this flag is set if column lower bound was changed */
876 char ub_changed;
877 /* this flag is set if column upper bound was changed */
878 NPPLFE *ptr;
879 /* list of non-zero coefficients a[i,q], i != p */
880 };
881
882 static int rcv_ineq_singlet(NPP *npp, void *info);
883
npp_ineq_singlet(NPP * npp,NPPROW * p)884 int npp_ineq_singlet(NPP *npp, NPPROW *p)
885 { /* process row singleton (inequality constraint) */
886 struct ineq_singlet *info;
887 NPPCOL *q;
888 NPPAIJ *apq, *aij;
889 NPPLFE *lfe;
890 int lb_changed, ub_changed;
891 double ll, uu;
892 /* the row must be singleton inequality constraint */
893 xassert(p->lb != -DBL_MAX || p->ub != +DBL_MAX);
894 xassert(p->lb < p->ub);
895 xassert(p->ptr != NULL && p->ptr->r_next == NULL);
896 /* compute implied column bounds */
897 apq = p->ptr;
898 q = apq->col;
899 xassert(q->lb < q->ub);
900 if (apq->val > 0.0)
901 { ll = (p->lb == -DBL_MAX ? -DBL_MAX : p->lb / apq->val);
902 uu = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub / apq->val);
903 }
904 else
905 { ll = (p->ub == +DBL_MAX ? -DBL_MAX : p->ub / apq->val);
906 uu = (p->lb == -DBL_MAX ? +DBL_MAX : p->lb / apq->val);
907 }
908 /* process implied column lower bound */
909 if (ll == -DBL_MAX)
910 lb_changed = 0;
911 else
912 { lb_changed = npp_implied_lower(npp, q, ll);
913 xassert(0 <= lb_changed && lb_changed <= 4);
914 if (lb_changed == 4) return 4; /* infeasible */
915 }
916 /* process implied column upper bound */
917 if (uu == +DBL_MAX)
918 ub_changed = 0;
919 else if (lb_changed == 3)
920 { /* column was fixed on its upper bound due to l'[q] = u[q] */
921 /* note that L[p] < U[p], so l'[q] = u[q] < u'[q] */
922 ub_changed = 0;
923 }
924 else
925 { ub_changed = npp_implied_upper(npp, q, uu);
926 xassert(0 <= ub_changed && ub_changed <= 4);
927 if (ub_changed == 4) return 4; /* infeasible */
928 }
929 /* if neither lower nor upper column bound was changed, the row
930 is originally redundant and can be replaced by free row */
931 if (!lb_changed && !ub_changed)
932 { p->lb = -DBL_MAX, p->ub = +DBL_MAX;
933 npp_free_row(npp, p);
934 return 0;
935 }
936 /* create transformation stack entry */
937 info = npp_push_tse(npp,
938 rcv_ineq_singlet, sizeof(struct ineq_singlet));
939 info->p = p->i;
940 info->q = q->j;
941 info->apq = apq->val;
942 info->c = q->coef;
943 info->lb = p->lb;
944 info->ub = p->ub;
945 info->lb_changed = (char)lb_changed;
946 info->ub_changed = (char)ub_changed;
947 info->ptr = NULL;
948 /* save column coefficients a[i,q], i != p (not needed for MIP
949 solution) */
950 if (npp->sol != GLP_MIP)
951 { for (aij = q->ptr; aij != NULL; aij = aij->c_next)
952 { if (aij == apq) continue; /* skip a[p,q] */
953 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
954 lfe->ref = aij->row->i;
955 lfe->val = aij->val;
956 lfe->next = info->ptr;
957 info->ptr = lfe;
958 }
959 }
960 /* remove the row from the problem */
961 npp_del_row(npp, p);
962 return lb_changed >= ub_changed ? lb_changed : ub_changed;
963 }
964
rcv_ineq_singlet(NPP * npp,void * _info)965 static int rcv_ineq_singlet(NPP *npp, void *_info)
966 { /* recover row singleton (inequality constraint) */
967 struct ineq_singlet *info = _info;
968 NPPLFE *lfe;
969 double lambda;
970 if (npp->sol == GLP_MIP) goto done;
971 /* compute lambda~[q] in solution to the transformed problem
972 with formula (8) */
973 lambda = info->c;
974 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
975 lambda -= lfe->val * npp->r_pi[lfe->ref];
976 if (npp->sol == GLP_SOL)
977 { /* recover basic solution */
978 if (npp->c_stat[info->q] == GLP_BS)
979 { /* column q is basic, so row p is inactive */
980 npp->r_stat[info->p] = GLP_BS;
981 npp->r_pi[info->p] = 0.0;
982 }
983 else if (npp->c_stat[info->q] == GLP_NL)
984 nl: { /* column q is non-basic with lower bound active */
985 if (info->lb_changed)
986 { /* it is implied bound, so actually row p is active
987 while column q is basic */
988 npp->r_stat[info->p] =
989 (char)(info->apq > 0.0 ? GLP_NL : GLP_NU);
990 npp->c_stat[info->q] = GLP_BS;
991 npp->r_pi[info->p] = lambda / info->apq;
992 }
993 else
994 { /* it is original bound, so row p is inactive */
995 npp->r_stat[info->p] = GLP_BS;
996 npp->r_pi[info->p] = 0.0;
997 }
998 }
999 else if (npp->c_stat[info->q] == GLP_NU)
1000 nu: { /* column q is non-basic with upper bound active */
1001 if (info->ub_changed)
1002 { /* it is implied bound, so actually row p is active
1003 while column q is basic */
1004 npp->r_stat[info->p] =
1005 (char)(info->apq > 0.0 ? GLP_NU : GLP_NL);
1006 npp->c_stat[info->q] = GLP_BS;
1007 npp->r_pi[info->p] = lambda / info->apq;
1008 }
1009 else
1010 { /* it is original bound, so row p is inactive */
1011 npp->r_stat[info->p] = GLP_BS;
1012 npp->r_pi[info->p] = 0.0;
1013 }
1014 }
1015 else if (npp->c_stat[info->q] == GLP_NS)
1016 { /* column q is non-basic and fixed; note, however, that in
1017 in the original problem it is non-fixed */
1018 if (lambda > +1e-7)
1019 { if (info->apq > 0.0 && info->lb != -DBL_MAX ||
1020 info->apq < 0.0 && info->ub != +DBL_MAX ||
1021 !info->lb_changed)
1022 { /* either corresponding bound of row p exists or
1023 column q remains non-basic with its original lower
1024 bound active */
1025 npp->c_stat[info->q] = GLP_NL;
1026 goto nl;
1027 }
1028 }
1029 if (lambda < -1e-7)
1030 { if (info->apq > 0.0 && info->ub != +DBL_MAX ||
1031 info->apq < 0.0 && info->lb != -DBL_MAX ||
1032 !info->ub_changed)
1033 { /* either corresponding bound of row p exists or
1034 column q remains non-basic with its original upper
1035 bound active */
1036 npp->c_stat[info->q] = GLP_NU;
1037 goto nu;
1038 }
1039 }
1040 /* either lambda~[q] is close to zero, or corresponding
1041 bound of row p does not exist, because lambda~[q] has
1042 wrong sign due to round-off errors; in the latter case
1043 lambda~[q] is also assumed to be close to zero; so, we
1044 can make row p active on its existing bound and column q
1045 basic; pi[p] will have wrong sign, but it also will be
1046 close to zero (rarus casus of dual degeneracy) */
1047 if (info->lb != -DBL_MAX && info->ub == +DBL_MAX)
1048 { /* row lower bound exists, but upper bound doesn't */
1049 npp->r_stat[info->p] = GLP_NL;
1050 }
1051 else if (info->lb == -DBL_MAX && info->ub != +DBL_MAX)
1052 { /* row upper bound exists, but lower bound doesn't */
1053 npp->r_stat[info->p] = GLP_NU;
1054 }
1055 else if (info->lb != -DBL_MAX && info->ub != +DBL_MAX)
1056 { /* both row lower and upper bounds exist */
1057 /* to choose proper active row bound we should not use
1058 lambda~[q], because its value being close to zero is
1059 unreliable; so we choose that bound which provides
1060 primal feasibility for original constraint (1) */
1061 if (info->apq * npp->c_value[info->q] <=
1062 0.5 * (info->lb + info->ub))
1063 npp->r_stat[info->p] = GLP_NL;
1064 else
1065 npp->r_stat[info->p] = GLP_NU;
1066 }
1067 else
1068 { npp_error();
1069 return 1;
1070 }
1071 npp->c_stat[info->q] = GLP_BS;
1072 npp->r_pi[info->p] = lambda / info->apq;
1073 }
1074 else
1075 { npp_error();
1076 return 1;
1077 }
1078 }
1079 if (npp->sol == GLP_IPT)
1080 { /* recover interior-point solution */
1081 if (lambda > +DBL_EPSILON && info->lb_changed ||
1082 lambda < -DBL_EPSILON && info->ub_changed)
1083 { /* actually row p has corresponding active bound */
1084 npp->r_pi[info->p] = lambda / info->apq;
1085 }
1086 else
1087 { /* either bounds of column q are both inactive or its
1088 original bound is active */
1089 npp->r_pi[info->p] = 0.0;
1090 }
1091 }
1092 done: return 0;
1093 }
1094
1095 /***********************************************************************
1096 * NAME
1097 *
1098 * npp_implied_slack - process column singleton (implied slack variable)
1099 *
1100 * SYNOPSIS
1101 *
1102 * #include "glpnpp.h"
1103 * void npp_implied_slack(NPP *npp, NPPCOL *q);
1104 *
1105 * DESCRIPTION
1106 *
1107 * The routine npp_implied_slack processes column q:
1108 *
1109 * l[q] <= x[q] <= u[q], (1)
1110 *
1111 * where l[q] < u[q], having the only non-zero coefficient in row p,
1112 * which is equality constraint:
1113 *
1114 * sum a[p,j] x[j] + a[p,q] x[q] = b. (2)
1115 * j!=q
1116 *
1117 * PROBLEM TRANSFORMATION
1118 *
1119 * (If x[q] is integral, this transformation must not be used.)
1120 *
1121 * The term a[p,q] x[q] in constraint (2) can be considered as a slack
1122 * variable that allows to carry bounds of column q over row p and then
1123 * remove column q from the problem.
1124 *
1125 * Constraint (2) can be written as follows:
1126 *
1127 * sum a[p,j] x[j] = b - a[p,q] x[q]. (3)
1128 * j!=q
1129 *
1130 * According to (1) constraint (3) is equivalent to the following
1131 * inequality constraint:
1132 *
1133 * L[p] <= sum a[p,j] x[j] <= U[p], (4)
1134 * j!=q
1135 *
1136 * where
1137 *
1138 * ( b - a[p,q] u[q], if a[p,q] > 0
1139 * L[p] = < (5)
1140 * ( b - a[p,q] l[q], if a[p,q] < 0
1141 *
1142 * ( b - a[p,q] l[q], if a[p,q] > 0
1143 * U[p] = < (6)
1144 * ( b - a[p,q] u[q], if a[p,q] < 0
1145 *
1146 * From (2) it follows that:
1147 *
1148 * 1
1149 * x[q] = ------ (b - sum a[p,j] x[j]). (7)
1150 * a[p,q] j!=q
1151 *
1152 * In order to eliminate x[q] from the objective row we substitute it
1153 * from (6) to that row:
1154 *
1155 * z = sum c[j] x[j] + c[q] x[q] + c[0] =
1156 * j!=q
1157 * 1
1158 * = sum c[j] x[j] + c[q] [------ (b - sum a[p,j] x[j])] + c0 =
1159 * j!=q a[p,q] j!=q
1160 *
1161 * = sum c~[j] x[j] + c~[0],
1162 * j!=q
1163 * a[p,j] b
1164 * c~[j] = c[j] - c[q] ------, c~0 = c0 - c[q] ------ (8)
1165 * a[p,q] a[p,q]
1166 *
1167 * are values of objective coefficients and constant term, resp., in
1168 * the transformed problem.
1169 *
1170 * Note that column q is column singleton, so in the dual system of the
1171 * original problem it corresponds to the following row singleton:
1172 *
1173 * a[p,q] pi[p] + lambda[q] = c[q]. (9)
1174 *
1175 * In the transformed problem row (9) would be the following:
1176 *
1177 * a[p,q] pi~[p] + lambda[q] = c~[q] = 0. (10)
1178 *
1179 * Subtracting (10) from (9) we have:
1180 *
1181 * a[p,q] (pi[p] - pi~[p]) = c[q]
1182 *
1183 * that gives the following formula to compute multiplier for row p in
1184 * solution to the original problem using its value in solution to the
1185 * transformed problem:
1186 *
1187 * pi[p] = pi~[p] + c[q] / a[p,q]. (11)
1188 *
1189 * RECOVERING BASIC SOLUTION
1190 *
1191 * Status of column q in solution to the original problem is defined
1192 * by status of row p in solution to the transformed problem and the
1193 * sign of coefficient a[p,q] in the original inequality constraint (2)
1194 * as follows:
1195 *
1196 * +-----------------------+---------+--------------------+
1197 * | Status of row p | Sign of | Status of column q |
1198 * | (transformed problem) | a[p,q] | (original problem) |
1199 * +-----------------------+---------+--------------------+
1200 * | GLP_BS | + / - | GLP_BS |
1201 * | GLP_NL | + | GLP_NU |
1202 * | GLP_NL | - | GLP_NL |
1203 * | GLP_NU | + | GLP_NL |
1204 * | GLP_NU | - | GLP_NU |
1205 * | GLP_NF | + / - | GLP_NF |
1206 * +-----------------------+---------+--------------------+
1207 *
1208 * Value of column q is computed with formula (7). Since originally row
1209 * p is equality constraint, its status is assigned GLP_NS, and value of
1210 * its multiplier pi[p] is computed with formula (11).
1211 *
1212 * RECOVERING INTERIOR-POINT SOLUTION
1213 *
1214 * Value of column q is computed with formula (7). Row multiplier value
1215 * pi[p] is computed with formula (11).
1216 *
1217 * RECOVERING MIP SOLUTION
1218 *
1219 * Value of column q is computed with formula (7). */
1220
1221 struct implied_slack
1222 { /* column singleton (implied slack variable) */
1223 int p;
1224 /* row reference number */
1225 int q;
1226 /* column reference number */
1227 double apq;
1228 /* constraint coefficient a[p,q] */
1229 double b;
1230 /* right-hand side of original equality constraint */
1231 double c;
1232 /* original objective coefficient at x[q] */
1233 NPPLFE *ptr;
1234 /* list of non-zero coefficients a[p,j], j != q */
1235 };
1236
1237 static int rcv_implied_slack(NPP *npp, void *info);
1238
npp_implied_slack(NPP * npp,NPPCOL * q)1239 void npp_implied_slack(NPP *npp, NPPCOL *q)
1240 { /* process column singleton (implied slack variable) */
1241 struct implied_slack *info;
1242 NPPROW *p;
1243 NPPAIJ *aij;
1244 NPPLFE *lfe;
1245 /* the column must be non-integral non-fixed singleton */
1246 xassert(!q->is_int);
1247 xassert(q->lb < q->ub);
1248 xassert(q->ptr != NULL && q->ptr->c_next == NULL);
1249 /* corresponding row must be equality constraint */
1250 aij = q->ptr;
1251 p = aij->row;
1252 xassert(p->lb == p->ub);
1253 /* create transformation stack entry */
1254 info = npp_push_tse(npp,
1255 rcv_implied_slack, sizeof(struct implied_slack));
1256 info->p = p->i;
1257 info->q = q->j;
1258 info->apq = aij->val;
1259 info->b = p->lb;
1260 info->c = q->coef;
1261 info->ptr = NULL;
1262 /* save row coefficients a[p,j], j != q, and substitute x[q]
1263 into the objective row */
1264 for (aij = p->ptr; aij != NULL; aij = aij->r_next)
1265 { if (aij->col == q) continue; /* skip a[p,q] */
1266 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
1267 lfe->ref = aij->col->j;
1268 lfe->val = aij->val;
1269 lfe->next = info->ptr;
1270 info->ptr = lfe;
1271 aij->col->coef -= info->c * (aij->val / info->apq);
1272 }
1273 npp->c0 += info->c * (info->b / info->apq);
1274 /* compute new row bounds */
1275 if (info->apq > 0.0)
1276 { p->lb = (q->ub == +DBL_MAX ?
1277 -DBL_MAX : info->b - info->apq * q->ub);
1278 p->ub = (q->lb == -DBL_MAX ?
1279 +DBL_MAX : info->b - info->apq * q->lb);
1280 }
1281 else
1282 { p->lb = (q->lb == -DBL_MAX ?
1283 -DBL_MAX : info->b - info->apq * q->lb);
1284 p->ub = (q->ub == +DBL_MAX ?
1285 +DBL_MAX : info->b - info->apq * q->ub);
1286 }
1287 /* remove the column from the problem */
1288 npp_del_col(npp, q);
1289 return;
1290 }
1291
rcv_implied_slack(NPP * npp,void * _info)1292 static int rcv_implied_slack(NPP *npp, void *_info)
1293 { /* recover column singleton (implied slack variable) */
1294 struct implied_slack *info = _info;
1295 NPPLFE *lfe;
1296 double temp;
1297 if (npp->sol == GLP_SOL)
1298 { /* assign statuses to row p and column q */
1299 if (npp->r_stat[info->p] == GLP_BS ||
1300 npp->r_stat[info->p] == GLP_NF)
1301 npp->c_stat[info->q] = npp->r_stat[info->p];
1302 else if (npp->r_stat[info->p] == GLP_NL)
1303 npp->c_stat[info->q] =
1304 (char)(info->apq > 0.0 ? GLP_NU : GLP_NL);
1305 else if (npp->r_stat[info->p] == GLP_NU)
1306 npp->c_stat[info->q] =
1307 (char)(info->apq > 0.0 ? GLP_NL : GLP_NU);
1308 else
1309 { npp_error();
1310 return 1;
1311 }
1312 npp->r_stat[info->p] = GLP_NS;
1313 }
1314 if (npp->sol != GLP_MIP)
1315 { /* compute multiplier for row p */
1316 npp->r_pi[info->p] += info->c / info->apq;
1317 }
1318 /* compute value of column q */
1319 temp = info->b;
1320 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
1321 temp -= lfe->val * npp->c_value[lfe->ref];
1322 npp->c_value[info->q] = temp / info->apq;
1323 return 0;
1324 }
1325
1326 /***********************************************************************
1327 * NAME
1328 *
1329 * npp_implied_free - process column singleton (implied free variable)
1330 *
1331 * SYNOPSIS
1332 *
1333 * #include "glpnpp.h"
1334 * int npp_implied_free(NPP *npp, NPPCOL *q);
1335 *
1336 * DESCRIPTION
1337 *
1338 * The routine npp_implied_free processes column q:
1339 *
1340 * l[q] <= x[q] <= u[q], (1)
1341 *
1342 * having non-zero coefficient in the only row p, which is inequality
1343 * constraint:
1344 *
1345 * L[p] <= sum a[p,j] x[j] + a[p,q] x[q] <= U[p], (2)
1346 * j!=q
1347 *
1348 * where l[q] < u[q], L[p] < U[p], L[p] > -oo and/or U[p] < +oo.
1349 *
1350 * RETURNS
1351 *
1352 * 0 - success;
1353 *
1354 * 1 - column lower and/or upper bound(s) can be active;
1355 *
1356 * 2 - problem has no dual feasible solution.
1357 *
1358 * PROBLEM TRANSFORMATION
1359 *
1360 * Constraint (2) can be written as follows:
1361 *
1362 * L[p] - sum a[p,j] x[j] <= a[p,q] x[q] <= U[p] - sum a[p,j] x[j],
1363 * j!=q j!=q
1364 *
1365 * from which it follows that:
1366 *
1367 * alfa <= a[p,q] x[q] <= beta, (3)
1368 *
1369 * where
1370 *
1371 * alfa = inf(L[p] - sum a[p,j] x[j]) =
1372 * j!=q
1373 *
1374 * = L[p] - sup sum a[p,j] x[j] = (4)
1375 * j!=q
1376 *
1377 * = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j],
1378 * j in Jp j in Jn
1379 *
1380 * beta = sup(L[p] - sum a[p,j] x[j]) =
1381 * j!=q
1382 *
1383 * = L[p] - inf sum a[p,j] x[j] = (5)
1384 * j!=q
1385 *
1386 * = L[p] - sum a[p,j] l[j] - sum a[p,j] u[j],
1387 * j in Jp j in Jn
1388 *
1389 * Jp = {j != q: a[p,j] > 0}, Jn = {j != q: a[p,j] < 0}. (6)
1390 *
1391 * Inequality (3) defines implied bounds of variable x[q]:
1392 *
1393 * l'[q] <= x[q] <= u'[q], (7)
1394 *
1395 * where
1396 *
1397 * ( alfa / a[p,q], if a[p,q] > 0
1398 * l'[q] = < (8a)
1399 * ( beta / a[p,q], if a[p,q] < 0
1400 *
1401 * ( beta / a[p,q], if a[p,q] > 0
1402 * u'[q] = < (8b)
1403 * ( alfa / a[p,q], if a[p,q] < 0
1404 *
1405 * Thus, if l'[q] > l[q] - eps and u'[q] < u[q] + eps, where eps is
1406 * an absolute tolerance for column value, column bounds (1) cannot be
1407 * active, in which case column q can be replaced by equivalent free
1408 * (unbounded) column.
1409 *
1410 * Note that column q is column singleton, so in the dual system of the
1411 * original problem it corresponds to the following row singleton:
1412 *
1413 * a[p,q] pi[p] + lambda[q] = c[q], (9)
1414 *
1415 * from which it follows that:
1416 *
1417 * pi[p] = (c[q] - lambda[q]) / a[p,q]. (10)
1418 *
1419 * Let x[q] be implied free (unbounded) variable. Then column q can be
1420 * only basic, so its multiplier lambda[q] is equal to zero, and from
1421 * (10) we have:
1422 *
1423 * pi[p] = c[q] / a[p,q]. (11)
1424 *
1425 * There are possible three cases:
1426 *
1427 * 1) pi[p] < -eps, where eps is an absolute tolerance for row
1428 * multiplier. In this case, to provide dual feasibility of the
1429 * original problem, row p must be active on its lower bound, and
1430 * if its lower bound does not exist (L[p] = -oo), the problem has
1431 * no dual feasible solution;
1432 *
1433 * 2) pi[p] > +eps. In this case row p must be active on its upper
1434 * bound, and if its upper bound does not exist (U[p] = +oo), the
1435 * problem has no dual feasible solution;
1436 *
1437 * 3) -eps <= pi[p] <= +eps. In this case any (either lower or upper)
1438 * bound of row p can be active, because this does not affect dual
1439 * feasibility.
1440 *
1441 * Thus, in all three cases original inequality constraint (2) can be
1442 * replaced by equality constraint, where the right-hand side is either
1443 * lower or upper bound of row p, and bounds of column q can be removed
1444 * that makes it free (unbounded). (May note that this transformation
1445 * can be followed by transformation "Column singleton (implied slack
1446 * variable)" performed by the routine npp_implied_slack.)
1447 *
1448 * RECOVERING BASIC SOLUTION
1449 *
1450 * Status of row p in solution to the original problem is determined
1451 * by its status in solution to the transformed problem and its bound,
1452 * which was choosen to be active:
1453 *
1454 * +-----------------------+--------+--------------------+
1455 * | Status of row p | Active | Status of row p |
1456 * | (transformed problem) | bound | (original problem) |
1457 * +-----------------------+--------+--------------------+
1458 * | GLP_BS | L[p] | GLP_BS |
1459 * | GLP_BS | U[p] | GLP_BS |
1460 * | GLP_NS | L[p] | GLP_NL |
1461 * | GLP_NS | U[p] | GLP_NU |
1462 * +-----------------------+--------+--------------------+
1463 *
1464 * Value of row multiplier pi[p] (as well as value of column q) in
1465 * solution to the original problem is the same as in solution to the
1466 * transformed problem.
1467 *
1468 * RECOVERING INTERIOR-POINT SOLUTION
1469 *
1470 * Value of row multiplier pi[p] in solution to the original problem is
1471 * the same as in solution to the transformed problem.
1472 *
1473 * RECOVERING MIP SOLUTION
1474 *
1475 * None needed. */
1476
1477 struct implied_free
1478 { /* column singleton (implied free variable) */
1479 int p;
1480 /* row reference number */
1481 char stat;
1482 /* row status:
1483 GLP_NL - active constraint on lower bound
1484 GLP_NU - active constraint on upper bound */
1485 };
1486
1487 static int rcv_implied_free(NPP *npp, void *info);
1488
npp_implied_free(NPP * npp,NPPCOL * q)1489 int npp_implied_free(NPP *npp, NPPCOL *q)
1490 { /* process column singleton (implied free variable) */
1491 struct implied_free *info;
1492 NPPROW *p;
1493 NPPAIJ *apq, *aij;
1494 double alfa, beta, l, u, pi, eps;
1495 /* the column must be non-fixed singleton */
1496 xassert(q->lb < q->ub);
1497 xassert(q->ptr != NULL && q->ptr->c_next == NULL);
1498 /* corresponding row must be inequality constraint */
1499 apq = q->ptr;
1500 p = apq->row;
1501 xassert(p->lb != -DBL_MAX || p->ub != +DBL_MAX);
1502 xassert(p->lb < p->ub);
1503 /* compute alfa */
1504 alfa = p->lb;
1505 if (alfa != -DBL_MAX)
1506 { for (aij = p->ptr; aij != NULL; aij = aij->r_next)
1507 { if (aij == apq) continue; /* skip a[p,q] */
1508 if (aij->val > 0.0)
1509 { if (aij->col->ub == +DBL_MAX)
1510 { alfa = -DBL_MAX;
1511 break;
1512 }
1513 alfa -= aij->val * aij->col->ub;
1514 }
1515 else /* < 0.0 */
1516 { if (aij->col->lb == -DBL_MAX)
1517 { alfa = -DBL_MAX;
1518 break;
1519 }
1520 alfa -= aij->val * aij->col->lb;
1521 }
1522 }
1523 }
1524 /* compute beta */
1525 beta = p->ub;
1526 if (beta != +DBL_MAX)
1527 { for (aij = p->ptr; aij != NULL; aij = aij->r_next)
1528 { if (aij == apq) continue; /* skip a[p,q] */
1529 if (aij->val > 0.0)
1530 { if (aij->col->lb == -DBL_MAX)
1531 { beta = +DBL_MAX;
1532 break;
1533 }
1534 beta -= aij->val * aij->col->lb;
1535 }
1536 else /* < 0.0 */
1537 { if (aij->col->ub == +DBL_MAX)
1538 { beta = +DBL_MAX;
1539 break;
1540 }
1541 beta -= aij->val * aij->col->ub;
1542 }
1543 }
1544 }
1545 /* compute implied column lower bound l'[q] */
1546 if (apq->val > 0.0)
1547 l = (alfa == -DBL_MAX ? -DBL_MAX : alfa / apq->val);
1548 else /* < 0.0 */
1549 l = (beta == +DBL_MAX ? -DBL_MAX : beta / apq->val);
1550 /* compute implied column upper bound u'[q] */
1551 if (apq->val > 0.0)
1552 u = (beta == +DBL_MAX ? +DBL_MAX : beta / apq->val);
1553 else
1554 u = (alfa == -DBL_MAX ? +DBL_MAX : alfa / apq->val);
1555 /* check if column lower bound l[q] can be active */
1556 if (q->lb != -DBL_MAX)
1557 { eps = 1e-9 + 1e-12 * fabs(q->lb);
1558 if (l < q->lb - eps) return 1; /* yes, it can */
1559 }
1560 /* check if column upper bound u[q] can be active */
1561 if (q->ub != +DBL_MAX)
1562 { eps = 1e-9 + 1e-12 * fabs(q->ub);
1563 if (u > q->ub + eps) return 1; /* yes, it can */
1564 }
1565 /* okay; make column q free (unbounded) */
1566 q->lb = -DBL_MAX, q->ub = +DBL_MAX;
1567 /* create transformation stack entry */
1568 info = npp_push_tse(npp,
1569 rcv_implied_free, sizeof(struct implied_free));
1570 info->p = p->i;
1571 info->stat = -1;
1572 /* compute row multiplier pi[p] */
1573 pi = q->coef / apq->val;
1574 /* check dual feasibility for row p */
1575 if (pi > +DBL_EPSILON)
1576 { /* lower bound L[p] must be active */
1577 if (p->lb != -DBL_MAX)
1578 nl: { info->stat = GLP_NL;
1579 p->ub = p->lb;
1580 }
1581 else
1582 { if (pi > +1e-5) return 2; /* dual infeasibility */
1583 /* take a chance on U[p] */
1584 xassert(p->ub != +DBL_MAX);
1585 goto nu;
1586 }
1587 }
1588 else if (pi < -DBL_EPSILON)
1589 { /* upper bound U[p] must be active */
1590 if (p->ub != +DBL_MAX)
1591 nu: { info->stat = GLP_NU;
1592 p->lb = p->ub;
1593 }
1594 else
1595 { if (pi < -1e-5) return 2; /* dual infeasibility */
1596 /* take a chance on L[p] */
1597 xassert(p->lb != -DBL_MAX);
1598 goto nl;
1599 }
1600 }
1601 else
1602 { /* any bound (either L[p] or U[p]) can be made active */
1603 if (p->ub == +DBL_MAX)
1604 { xassert(p->lb != -DBL_MAX);
1605 goto nl;
1606 }
1607 if (p->lb == -DBL_MAX)
1608 { xassert(p->ub != +DBL_MAX);
1609 goto nu;
1610 }
1611 if (fabs(p->lb) <= fabs(p->ub)) goto nl; else goto nu;
1612 }
1613 return 0;
1614 }
1615
rcv_implied_free(NPP * npp,void * _info)1616 static int rcv_implied_free(NPP *npp, void *_info)
1617 { /* recover column singleton (implied free variable) */
1618 struct implied_free *info = _info;
1619 if (npp->sol == GLP_SOL)
1620 { if (npp->r_stat[info->p] == GLP_BS)
1621 npp->r_stat[info->p] = GLP_BS;
1622 else if (npp->r_stat[info->p] == GLP_NS)
1623 { xassert(info->stat == GLP_NL || info->stat == GLP_NU);
1624 npp->r_stat[info->p] = info->stat;
1625 }
1626 else
1627 { npp_error();
1628 return 1;
1629 }
1630 }
1631 return 0;
1632 }
1633
1634 /***********************************************************************
1635 * NAME
1636 *
1637 * npp_eq_doublet - process row doubleton (equality constraint)
1638 *
1639 * SYNOPSIS
1640 *
1641 * #include "glpnpp.h"
1642 * NPPCOL *npp_eq_doublet(NPP *npp, NPPROW *p);
1643 *
1644 * DESCRIPTION
1645 *
1646 * The routine npp_eq_doublet processes row p, which is equality
1647 * constraint having exactly two non-zero coefficients:
1648 *
1649 * a[p,q] x[q] + a[p,r] x[r] = b. (1)
1650 *
1651 * As the result of processing one of columns q or r is eliminated from
1652 * all other rows and, thus, becomes column singleton of type "implied
1653 * slack variable". Row p is not changed and along with column q and r
1654 * remains in the problem.
1655 *
1656 * RETURNS
1657 *
1658 * The routine npp_eq_doublet returns pointer to the descriptor of that
1659 * column q or r which has been eliminated. If, due to some reason, the
1660 * elimination was not performed, the routine returns NULL.
1661 *
1662 * PROBLEM TRANSFORMATION
1663 *
1664 * First, we decide which column q or r will be eliminated. Let it be
1665 * column q. Consider i-th constraint row, where column q has non-zero
1666 * coefficient a[i,q] != 0:
1667 *
1668 * L[i] <= sum a[i,j] x[j] <= U[i]. (2)
1669 * j
1670 *
1671 * In order to eliminate column q from row (2) we subtract from it row
1672 * (1) multiplied by gamma[i] = a[i,q] / a[p,q], i.e. we replace in the
1673 * transformed problem row (2) by its linear combination with row (1).
1674 * This transformation changes only coefficients in columns q and r,
1675 * and bounds of row i as follows:
1676 *
1677 * a~[i,q] = a[i,q] - gamma[i] a[p,q] = 0, (3)
1678 *
1679 * a~[i,r] = a[i,r] - gamma[i] a[p,r], (4)
1680 *
1681 * L~[i] = L[i] - gamma[i] b, (5)
1682 *
1683 * U~[i] = U[i] - gamma[i] b. (6)
1684 *
1685 * RECOVERING BASIC SOLUTION
1686 *
1687 * The transformation of the primal system of the original problem:
1688 *
1689 * L <= A x <= U (7)
1690 *
1691 * is equivalent to multiplying from the left a transformation matrix F
1692 * by components of this primal system, which in the transformed problem
1693 * becomes the following:
1694 *
1695 * F L <= F A x <= F U ==> L~ <= A~x <= U~. (8)
1696 *
1697 * The matrix F has the following structure:
1698 *
1699 * ( 1 -gamma[1] )
1700 * ( )
1701 * ( 1 -gamma[2] )
1702 * ( )
1703 * ( ... ... )
1704 * ( )
1705 * F = ( 1 -gamma[p-1] ) (9)
1706 * ( )
1707 * ( 1 )
1708 * ( )
1709 * ( -gamma[p+1] 1 )
1710 * ( )
1711 * ( ... ... )
1712 *
1713 * where its column containing elements -gamma[i] corresponds to row p
1714 * of the primal system.
1715 *
1716 * From (8) it follows that the dual system of the original problem:
1717 *
1718 * A'pi + lambda = c, (10)
1719 *
1720 * in the transformed problem becomes the following:
1721 *
1722 * A'F'inv(F')pi + lambda = c ==> (A~)'pi~ + lambda = c, (11)
1723 *
1724 * where:
1725 *
1726 * pi~ = inv(F')pi (12)
1727 *
1728 * is the vector of row multipliers in the transformed problem. Thus:
1729 *
1730 * pi = F'pi~. (13)
1731 *
1732 * Therefore, as it follows from (13), value of multiplier for row p in
1733 * solution to the original problem can be computed as follows:
1734 *
1735 * pi[p] = pi~[p] - sum gamma[i] pi~[i], (14)
1736 * i
1737 *
1738 * where pi~[i] = pi[i] is multiplier for row i (i != p).
1739 *
1740 * Note that the statuses of all rows and columns are not changed.
1741 *
1742 * RECOVERING INTERIOR-POINT SOLUTION
1743 *
1744 * Multiplier for row p in solution to the original problem is computed
1745 * with formula (14).
1746 *
1747 * RECOVERING MIP SOLUTION
1748 *
1749 * None needed. */
1750
1751 struct eq_doublet
1752 { /* row doubleton (equality constraint) */
1753 int p;
1754 /* row reference number */
1755 double apq;
1756 /* constraint coefficient a[p,q] */
1757 NPPLFE *ptr;
1758 /* list of non-zero coefficients a[i,q], i != p */
1759 };
1760
1761 static int rcv_eq_doublet(NPP *npp, void *info);
1762
npp_eq_doublet(NPP * npp,NPPROW * p)1763 NPPCOL *npp_eq_doublet(NPP *npp, NPPROW *p)
1764 { /* process row doubleton (equality constraint) */
1765 struct eq_doublet *info;
1766 NPPROW *i;
1767 NPPCOL *q, *r;
1768 NPPAIJ *apq, *apr, *aiq, *air, *next;
1769 NPPLFE *lfe;
1770 double gamma;
1771 /* the row must be doubleton equality constraint */
1772 xassert(p->lb == p->ub);
1773 xassert(p->ptr != NULL && p->ptr->r_next != NULL &&
1774 p->ptr->r_next->r_next == NULL);
1775 /* choose column to be eliminated */
1776 { NPPAIJ *a1, *a2;
1777 a1 = p->ptr, a2 = a1->r_next;
1778 if (fabs(a2->val) < 0.001 * fabs(a1->val))
1779 { /* only first column can be eliminated, because second one
1780 has too small constraint coefficient */
1781 apq = a1, apr = a2;
1782 }
1783 else if (fabs(a1->val) < 0.001 * fabs(a2->val))
1784 { /* only second column can be eliminated, because first one
1785 has too small constraint coefficient */
1786 apq = a2, apr = a1;
1787 }
1788 else
1789 { /* both columns are appropriate; choose that one which is
1790 shorter to minimize fill-in */
1791 if (npp_col_nnz(npp, a1->col) <= npp_col_nnz(npp, a2->col))
1792 { /* first column is shorter */
1793 apq = a1, apr = a2;
1794 }
1795 else
1796 { /* second column is shorter */
1797 apq = a2, apr = a1;
1798 }
1799 }
1800 }
1801 /* now columns q and r have been chosen */
1802 q = apq->col, r = apr->col;
1803 /* create transformation stack entry */
1804 info = npp_push_tse(npp,
1805 rcv_eq_doublet, sizeof(struct eq_doublet));
1806 info->p = p->i;
1807 info->apq = apq->val;
1808 info->ptr = NULL;
1809 /* transform each row i (i != p), where a[i,q] != 0, to eliminate
1810 column q */
1811 for (aiq = q->ptr; aiq != NULL; aiq = next)
1812 { next = aiq->c_next;
1813 if (aiq == apq) continue; /* skip row p */
1814 i = aiq->row; /* row i to be transformed */
1815 /* save constraint coefficient a[i,q] */
1816 if (npp->sol != GLP_MIP)
1817 { lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
1818 lfe->ref = i->i;
1819 lfe->val = aiq->val;
1820 lfe->next = info->ptr;
1821 info->ptr = lfe;
1822 }
1823 /* find coefficient a[i,r] in row i */
1824 for (air = i->ptr; air != NULL; air = air->r_next)
1825 if (air->col == r) break;
1826 /* if a[i,r] does not exist, create a[i,r] = 0 */
1827 if (air == NULL)
1828 air = npp_add_aij(npp, i, r, 0.0);
1829 /* compute gamma[i] = a[i,q] / a[p,q] */
1830 gamma = aiq->val / apq->val;
1831 /* (row i) := (row i) - gamma[i] * (row p); see (3)-(6) */
1832 /* new a[i,q] is exact zero due to elimnation; remove it from
1833 row i */
1834 npp_del_aij(npp, aiq);
1835 /* compute new a[i,r] */
1836 air->val -= gamma * apr->val;
1837 /* if new a[i,r] is close to zero due to numeric cancelation,
1838 remove it from row i */
1839 if (fabs(air->val) <= 1e-10)
1840 npp_del_aij(npp, air);
1841 /* compute new lower and upper bounds of row i */
1842 if (i->lb == i->ub)
1843 i->lb = i->ub = (i->lb - gamma * p->lb);
1844 else
1845 { if (i->lb != -DBL_MAX)
1846 i->lb -= gamma * p->lb;
1847 if (i->ub != +DBL_MAX)
1848 i->ub -= gamma * p->lb;
1849 }
1850 }
1851 return q;
1852 }
1853
rcv_eq_doublet(NPP * npp,void * _info)1854 static int rcv_eq_doublet(NPP *npp, void *_info)
1855 { /* recover row doubleton (equality constraint) */
1856 struct eq_doublet *info = _info;
1857 NPPLFE *lfe;
1858 double gamma, temp;
1859 /* we assume that processing row p is followed by processing
1860 column q as singleton of type "implied slack variable", in
1861 which case row p must always be active equality constraint */
1862 if (npp->sol == GLP_SOL)
1863 { if (npp->r_stat[info->p] != GLP_NS)
1864 { npp_error();
1865 return 1;
1866 }
1867 }
1868 if (npp->sol != GLP_MIP)
1869 { /* compute value of multiplier for row p; see (14) */
1870 temp = npp->r_pi[info->p];
1871 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
1872 { gamma = lfe->val / info->apq; /* a[i,q] / a[p,q] */
1873 temp -= gamma * npp->r_pi[lfe->ref];
1874 }
1875 npp->r_pi[info->p] = temp;
1876 }
1877 return 0;
1878 }
1879
1880 /***********************************************************************
1881 * NAME
1882 *
1883 * npp_forcing_row - process forcing row
1884 *
1885 * SYNOPSIS
1886 *
1887 * #include "glpnpp.h"
1888 * int npp_forcing_row(NPP *npp, NPPROW *p, int at);
1889 *
1890 * DESCRIPTION
1891 *
1892 * The routine npp_forcing row processes row p of general format:
1893 *
1894 * L[p] <= sum a[p,j] x[j] <= U[p], (1)
1895 * j
1896 *
1897 * l[j] <= x[j] <= u[j], (2)
1898 *
1899 * where L[p] <= U[p] and l[j] < u[j] for all a[p,j] != 0. It is also
1900 * assumed that:
1901 *
1902 * 1) if at = 0 then |L[p] - U'[p]| <= eps, where U'[p] is implied
1903 * row upper bound (see below), eps is an absolute tolerance for row
1904 * value;
1905 *
1906 * 2) if at = 1 then |U[p] - L'[p]| <= eps, where L'[p] is implied
1907 * row lower bound (see below).
1908 *
1909 * RETURNS
1910 *
1911 * 0 - success;
1912 *
1913 * 1 - cannot fix columns due to too small constraint coefficients.
1914 *
1915 * PROBLEM TRANSFORMATION
1916 *
1917 * Implied lower and upper bounds of row (1) are determined by bounds
1918 * of corresponding columns (variables) as follows:
1919 *
1920 * L'[p] = inf sum a[p,j] x[j] =
1921 * j
1922 * (3)
1923 * = sum a[p,j] l[j] + sum a[p,j] u[j],
1924 * j in Jp j in Jn
1925 *
1926 * U'[p] = sup sum a[p,j] x[j] =
1927 * (4)
1928 * = sum a[p,j] u[j] + sum a[p,j] l[j],
1929 * j in Jp j in Jn
1930 *
1931 * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)
1932 *
1933 * If L[p] =~ U'[p] (at = 0), solution can be primal feasible only when
1934 * all variables take their boundary values as defined by (4):
1935 *
1936 * ( u[j], if j in Jp
1937 * x[j] = < (6)
1938 * ( l[j], if j in Jn
1939 *
1940 * Similarly, if U[p] =~ L'[p] (at = 1), solution can be primal feasible
1941 * only when all variables take their boundary values as defined by (3):
1942 *
1943 * ( l[j], if j in Jp
1944 * x[j] = < (7)
1945 * ( u[j], if j in Jn
1946 *
1947 * Condition (6) or (7) allows fixing all columns (variables x[j])
1948 * in row (1) on their bounds and then removing them from the problem
1949 * (see the routine npp_fixed_col). Due to this row p becomes redundant,
1950 * so it can be replaced by equivalent free (unbounded) row and also
1951 * removed from the problem (see the routine npp_free_row).
1952 *
1953 * 1. To apply this transformation row (1) should not have coefficients
1954 * whose magnitude is too small, i.e. all a[p,j] should satisfy to
1955 * the following condition:
1956 *
1957 * |a[p,j]| >= eps * max(1, |a[p,k]|), (8)
1958 * k
1959 * where eps is a relative tolerance for constraint coefficients.
1960 * Otherwise, fixing columns may be numerically unreliable and may
1961 * lead to wrong solution.
1962 *
1963 * 2. The routine fixes columns and remove bounds of row p, however,
1964 * it does not remove the row and columns from the problem.
1965 *
1966 * RECOVERING BASIC SOLUTION
1967 *
1968 * In the transformed problem row p being inactive constraint is
1969 * assigned status GLP_BS (as the result of transformation of free
1970 * row), and all columns in this row are assigned status GLP_NS (as the
1971 * result of transformation of fixed columns).
1972 *
1973 * Note that in the dual system of the transformed (as well as original)
1974 * problem every column j in row p corresponds to the following row:
1975 *
1976 * sum a[i,j] pi[i] + a[p,j] pi[p] + lambda[j] = c[j], (9)
1977 * i!=p
1978 *
1979 * from which it follows that:
1980 *
1981 * lambda[j] = c[j] - sum a[i,j] pi[i] - a[p,j] pi[p]. (10)
1982 * i!=p
1983 *
1984 * In the transformed problem values of all multipliers pi[i] are known
1985 * (including pi[i], whose value is zero, since row p is inactive).
1986 * Thus, using formula (10) it is possible to compute values of
1987 * multipliers lambda[j] for all columns in row p.
1988 *
1989 * Note also that in the original problem all columns in row p are
1990 * bounded, not fixed. So status GLP_NS assigned to every such column
1991 * must be changed to GLP_NL or GLP_NU depending on which bound the
1992 * corresponding column has been fixed. This status change may lead to
1993 * dual feasibility violation for solution of the original problem,
1994 * because now column multipliers must satisfy to the following
1995 * condition:
1996 *
1997 * ( >= 0, if status of column j is GLP_NL,
1998 * lambda[j] < (11)
1999 * ( <= 0, if status of column j is GLP_NU.
2000 *
2001 * If this condition holds, solution to the original problem is the
2002 * same as to the transformed problem. Otherwise, we have to perform
2003 * one degenerate pivoting step of the primal simplex method to obtain
2004 * dual feasible (hence, optimal) solution to the original problem as
2005 * follows. If, on problem transformation, row p was made active on its
2006 * lower bound (case at = 0), we change its status to GLP_NL (or GLP_NS)
2007 * and start increasing its multiplier pi[p]. Otherwise, if row p was
2008 * made active on its upper bound (case at = 1), we change its status
2009 * to GLP_NU (or GLP_NS) and start decreasing pi[p]. From (10) it
2010 * follows that:
2011 *
2012 * delta lambda[j] = - a[p,j] * delta pi[p] = - a[p,j] pi[p]. (12)
2013 *
2014 * Simple analysis of formulae (3)-(5) shows that changing pi[p] in the
2015 * specified direction causes increasing lambda[j] for every column j
2016 * assigned status GLP_NL (delta lambda[j] > 0) and decreasing lambda[j]
2017 * for every column j assigned status GLP_NU (delta lambda[j] < 0). It
2018 * is understood that once the last lambda[q], which violates condition
2019 * (11), has reached zero, multipliers lambda[j] for all columns get
2020 * valid signs. Such column q can be determined as follows. Let d[j] be
2021 * initial value of lambda[j] (i.e. reduced cost of column j) in the
2022 * transformed problem computed with formula (10) when pi[p] = 0. Then
2023 * lambda[j] = d[j] + delta lambda[j], and from (12) it follows that
2024 * lambda[j] becomes zero if:
2025 *
2026 * delta lambda[j] = - a[p,j] pi[p] = - d[j] ==>
2027 * (13)
2028 * pi[p] = d[j] / a[p,j].
2029 *
2030 * Therefore, the last column q, for which lambda[q] becomes zero, can
2031 * be determined from the following condition:
2032 *
2033 * |d[q] / a[p,q]| = max |pi[p]| = max |d[j] / a[p,j]|, (14)
2034 * j in D j in D
2035 *
2036 * where D is a set of columns j whose, reduced costs d[j] have invalid
2037 * signs, i.e. violate condition (11). (Thus, if D is empty, solution
2038 * to the original problem is the same as solution to the transformed
2039 * problem, and no correction is needed as was noticed above.) In
2040 * solution to the original problem column q is assigned status GLP_BS,
2041 * since it replaces column of auxiliary variable of row p (becoming
2042 * active) in the basis, and multiplier for row p is assigned its new
2043 * value, which is pi[p] = d[q] / a[p,q]. Note that due to primal
2044 * degeneracy values of all columns having non-zero coefficients in row
2045 * p remain unchanged.
2046 *
2047 * RECOVERING INTERIOR-POINT SOLUTION
2048 *
2049 * Value of multiplier pi[p] in solution to the original problem is
2050 * corrected in the same way as for basic solution. Values of all
2051 * columns having non-zero coefficients in row p remain unchanged.
2052 *
2053 * RECOVERING MIP SOLUTION
2054 *
2055 * None needed. */
2056
2057 struct forcing_col
2058 { /* column fixed on its bound by forcing row */
2059 int j;
2060 /* column reference number */
2061 char stat;
2062 /* original column status:
2063 GLP_NL - fixed on lower bound
2064 GLP_NU - fixed on upper bound */
2065 double a;
2066 /* constraint coefficient a[p,j] */
2067 double c;
2068 /* objective coefficient c[j] */
2069 NPPLFE *ptr;
2070 /* list of non-zero coefficients a[i,j], i != p */
2071 struct forcing_col *next;
2072 /* pointer to another column fixed by forcing row */
2073 };
2074
2075 struct forcing_row
2076 { /* forcing row */
2077 int p;
2078 /* row reference number */
2079 char stat;
2080 /* status assigned to the row if it becomes active:
2081 GLP_NS - active equality constraint
2082 GLP_NL - inequality constraint with lower bound active
2083 GLP_NU - inequality constraint with upper bound active */
2084 struct forcing_col *ptr;
2085 /* list of all columns having non-zero constraint coefficient
2086 a[p,j] in the forcing row */
2087 };
2088
2089 static int rcv_forcing_row(NPP *npp, void *info);
2090
npp_forcing_row(NPP * npp,NPPROW * p,int at)2091 int npp_forcing_row(NPP *npp, NPPROW *p, int at)
2092 { /* process forcing row */
2093 struct forcing_row *info;
2094 struct forcing_col *col = NULL;
2095 NPPCOL *j;
2096 NPPAIJ *apj, *aij;
2097 NPPLFE *lfe;
2098 double big;
2099 xassert(at == 0 || at == 1);
2100 /* determine maximal magnitude of the row coefficients */
2101 big = 1.0;
2102 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
2103 if (big < fabs(apj->val)) big = fabs(apj->val);
2104 /* if there are too small coefficients in the row, transformation
2105 should not be applied */
2106 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
2107 if (fabs(apj->val) < 1e-7 * big) return 1;
2108 /* create transformation stack entry */
2109 info = npp_push_tse(npp,
2110 rcv_forcing_row, sizeof(struct forcing_row));
2111 info->p = p->i;
2112 if (p->lb == p->ub)
2113 { /* equality constraint */
2114 info->stat = GLP_NS;
2115 }
2116 else if (at == 0)
2117 { /* inequality constraint; case L[p] = U'[p] */
2118 info->stat = GLP_NL;
2119 xassert(p->lb != -DBL_MAX);
2120 }
2121 else /* at == 1 */
2122 { /* inequality constraint; case U[p] = L'[p] */
2123 info->stat = GLP_NU;
2124 xassert(p->ub != +DBL_MAX);
2125 }
2126 info->ptr = NULL;
2127 /* scan the forcing row, fix columns at corresponding bounds, and
2128 save column information (the latter is not needed for MIP) */
2129 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
2130 { /* column j has non-zero coefficient in the forcing row */
2131 j = apj->col;
2132 /* it must be non-fixed */
2133 xassert(j->lb < j->ub);
2134 /* allocate stack entry to save column information */
2135 if (npp->sol != GLP_MIP)
2136 { col = dmp_get_atom(npp->stack, sizeof(struct forcing_col));
2137 col->j = j->j;
2138 col->stat = -1; /* will be set below */
2139 col->a = apj->val;
2140 col->c = j->coef;
2141 col->ptr = NULL;
2142 col->next = info->ptr;
2143 info->ptr = col;
2144 }
2145 /* fix column j */
2146 if (at == 0 && apj->val < 0.0 || at != 0 && apj->val > 0.0)
2147 { /* at its lower bound */
2148 if (npp->sol != GLP_MIP)
2149 col->stat = GLP_NL;
2150 xassert(j->lb != -DBL_MAX);
2151 j->ub = j->lb;
2152 }
2153 else
2154 { /* at its upper bound */
2155 if (npp->sol != GLP_MIP)
2156 col->stat = GLP_NU;
2157 xassert(j->ub != +DBL_MAX);
2158 j->lb = j->ub;
2159 }
2160 /* save column coefficients a[i,j], i != p */
2161 if (npp->sol != GLP_MIP)
2162 { for (aij = j->ptr; aij != NULL; aij = aij->c_next)
2163 { if (aij == apj) continue; /* skip a[p,j] */
2164 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
2165 lfe->ref = aij->row->i;
2166 lfe->val = aij->val;
2167 lfe->next = col->ptr;
2168 col->ptr = lfe;
2169 }
2170 }
2171 }
2172 /* make the row free (unbounded) */
2173 p->lb = -DBL_MAX, p->ub = +DBL_MAX;
2174 return 0;
2175 }
2176
rcv_forcing_row(NPP * npp,void * _info)2177 static int rcv_forcing_row(NPP *npp, void *_info)
2178 { /* recover forcing row */
2179 struct forcing_row *info = _info;
2180 struct forcing_col *col, *piv;
2181 NPPLFE *lfe;
2182 double d, big, temp;
2183 if (npp->sol == GLP_MIP) goto done;
2184 /* initially solution to the original problem is the same as
2185 to the transformed problem, where row p is inactive constraint
2186 with pi[p] = 0, and all columns are non-basic */
2187 if (npp->sol == GLP_SOL)
2188 { if (npp->r_stat[info->p] != GLP_BS)
2189 { npp_error();
2190 return 1;
2191 }
2192 for (col = info->ptr; col != NULL; col = col->next)
2193 { if (npp->c_stat[col->j] != GLP_NS)
2194 { npp_error();
2195 return 1;
2196 }
2197 npp->c_stat[col->j] = col->stat; /* original status */
2198 }
2199 }
2200 /* compute reduced costs d[j] for all columns with formula (10)
2201 and store them in col.c instead objective coefficients */
2202 for (col = info->ptr; col != NULL; col = col->next)
2203 { d = col->c;
2204 for (lfe = col->ptr; lfe != NULL; lfe = lfe->next)
2205 d -= lfe->val * npp->r_pi[lfe->ref];
2206 col->c = d;
2207 }
2208 /* consider columns j, whose multipliers lambda[j] has wrong
2209 sign in solution to the transformed problem (where lambda[j] =
2210 d[j]), and choose column q, whose multipler lambda[q] reaches
2211 zero last on changing row multiplier pi[p]; see (14) */
2212 piv = NULL, big = 0.0;
2213 for (col = info->ptr; col != NULL; col = col->next)
2214 { d = col->c; /* d[j] */
2215 temp = fabs(d / col->a);
2216 if (col->stat == GLP_NL)
2217 { /* column j has active lower bound */
2218 if (d < 0.0 && big < temp)
2219 piv = col, big = temp;
2220 }
2221 else if (col->stat == GLP_NU)
2222 { /* column j has active upper bound */
2223 if (d > 0.0 && big < temp)
2224 piv = col, big = temp;
2225 }
2226 else
2227 { npp_error();
2228 return 1;
2229 }
2230 }
2231 /* if column q does not exist, no correction is needed */
2232 if (piv != NULL)
2233 { /* correct solution; row p becomes active constraint while
2234 column q becomes basic */
2235 if (npp->sol == GLP_SOL)
2236 { npp->r_stat[info->p] = info->stat;
2237 npp->c_stat[piv->j] = GLP_BS;
2238 }
2239 /* assign new value to row multiplier pi[p] = d[p] / a[p,q] */
2240 npp->r_pi[info->p] = piv->c / piv->a;
2241 }
2242 done: return 0;
2243 }
2244
2245 /***********************************************************************
2246 * NAME
2247 *
2248 * npp_analyze_row - perform general row analysis
2249 *
2250 * SYNOPSIS
2251 *
2252 * #include "glpnpp.h"
2253 * int npp_analyze_row(NPP *npp, NPPROW *p);
2254 *
2255 * DESCRIPTION
2256 *
2257 * The routine npp_analyze_row performs analysis of row p of general
2258 * format:
2259 *
2260 * L[p] <= sum a[p,j] x[j] <= U[p], (1)
2261 * j
2262 *
2263 * l[j] <= x[j] <= u[j], (2)
2264 *
2265 * where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0.
2266 *
2267 * RETURNS
2268 *
2269 * 0x?0 - row lower bound does not exist or is redundant;
2270 *
2271 * 0x?1 - row lower bound can be active;
2272 *
2273 * 0x?2 - row lower bound is a forcing bound;
2274 *
2275 * 0x0? - row upper bound does not exist or is redundant;
2276 *
2277 * 0x1? - row upper bound can be active;
2278 *
2279 * 0x2? - row upper bound is a forcing bound;
2280 *
2281 * 0x33 - row bounds are inconsistent with column bounds.
2282 *
2283 * ALGORITHM
2284 *
2285 * Analysis of row (1) is based on analysis of its implied lower and
2286 * upper bounds, which are determined by bounds of corresponding columns
2287 * (variables) as follows:
2288 *
2289 * L'[p] = inf sum a[p,j] x[j] =
2290 * j
2291 * (3)
2292 * = sum a[p,j] l[j] + sum a[p,j] u[j],
2293 * j in Jp j in Jn
2294 *
2295 * U'[p] = sup sum a[p,j] x[j] =
2296 * (4)
2297 * = sum a[p,j] u[j] + sum a[p,j] l[j],
2298 * j in Jp j in Jn
2299 *
2300 * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)
2301 *
2302 * (Note that bounds of all columns in row p are assumed to be correct,
2303 * so L'[p] <= U'[p].)
2304 *
2305 * Analysis of row lower bound L[p] includes the following cases:
2306 *
2307 * 1) if L[p] > U'[p] + eps, where eps is an absolute tolerance for row
2308 * value, row lower bound L[p] and implied row upper bound U'[p] are
2309 * inconsistent, ergo, the problem has no primal feasible solution;
2310 *
2311 * 2) if U'[p] - eps <= L[p] <= U'[p] + eps, i.e. if L[p] =~ U'[p],
2312 * the row is a forcing row on its lower bound (see description of
2313 * the routine npp_forcing_row);
2314 *
2315 * 3) if L[p] > L'[p] + eps, row lower bound L[p] can be active (this
2316 * conclusion does not account other rows in the problem);
2317 *
2318 * 4) if L[p] <= L'[p] + eps, row lower bound L[p] cannot be active, so
2319 * it is redundant and can be removed (replaced by -oo).
2320 *
2321 * Analysis of row upper bound U[p] is performed in a similar way and
2322 * includes the following cases:
2323 *
2324 * 1) if U[p] < L'[p] - eps, row upper bound U[p] and implied row lower
2325 * bound L'[p] are inconsistent, ergo the problem has no primal
2326 * feasible solution;
2327 *
2328 * 2) if L'[p] - eps <= U[p] <= L'[p] + eps, i.e. if U[p] =~ L'[p],
2329 * the row is a forcing row on its upper bound (see description of
2330 * the routine npp_forcing_row);
2331 *
2332 * 3) if U[p] < U'[p] - eps, row upper bound U[p] can be active (this
2333 * conclusion does not account other rows in the problem);
2334 *
2335 * 4) if U[p] >= U'[p] - eps, row upper bound U[p] cannot be active, so
2336 * it is redundant and can be removed (replaced by +oo). */
2337
npp_analyze_row(NPP * npp,NPPROW * p)2338 int npp_analyze_row(NPP *npp, NPPROW *p)
2339 { /* perform general row analysis */
2340 NPPAIJ *aij;
2341 int ret = 0x00;
2342 double l, u, eps;
2343 xassert(npp == npp);
2344 /* compute implied lower bound L'[p]; see (3) */
2345 l = 0.0;
2346 for (aij = p->ptr; aij != NULL; aij = aij->r_next)
2347 { if (aij->val > 0.0)
2348 { if (aij->col->lb == -DBL_MAX)
2349 { l = -DBL_MAX;
2350 break;
2351 }
2352 l += aij->val * aij->col->lb;
2353 }
2354 else /* aij->val < 0.0 */
2355 { if (aij->col->ub == +DBL_MAX)
2356 { l = -DBL_MAX;
2357 break;
2358 }
2359 l += aij->val * aij->col->ub;
2360 }
2361 }
2362 /* compute implied upper bound U'[p]; see (4) */
2363 u = 0.0;
2364 for (aij = p->ptr; aij != NULL; aij = aij->r_next)
2365 { if (aij->val > 0.0)
2366 { if (aij->col->ub == +DBL_MAX)
2367 { u = +DBL_MAX;
2368 break;
2369 }
2370 u += aij->val * aij->col->ub;
2371 }
2372 else /* aij->val < 0.0 */
2373 { if (aij->col->lb == -DBL_MAX)
2374 { u = +DBL_MAX;
2375 break;
2376 }
2377 u += aij->val * aij->col->lb;
2378 }
2379 }
2380 /* column bounds are assumed correct, so L'[p] <= U'[p] */
2381 /* check if row lower bound is consistent */
2382 if (p->lb != -DBL_MAX)
2383 { eps = 1e-3 + 1e-6 * fabs(p->lb);
2384 if (p->lb - eps > u)
2385 { ret = 0x33;
2386 goto done;
2387 }
2388 }
2389 /* check if row upper bound is consistent */
2390 if (p->ub != +DBL_MAX)
2391 { eps = 1e-3 + 1e-6 * fabs(p->ub);
2392 if (p->ub + eps < l)
2393 { ret = 0x33;
2394 goto done;
2395 }
2396 }
2397 /* check if row lower bound can be active/forcing */
2398 if (p->lb != -DBL_MAX)
2399 { eps = 1e-9 + 1e-12 * fabs(p->lb);
2400 if (p->lb - eps > l)
2401 { if (p->lb + eps <= u)
2402 ret |= 0x01;
2403 else
2404 ret |= 0x02;
2405 }
2406 }
2407 /* check if row upper bound can be active/forcing */
2408 if (p->ub != +DBL_MAX)
2409 { eps = 1e-9 + 1e-12 * fabs(p->ub);
2410 if (p->ub + eps < u)
2411 { /* check if the upper bound is forcing */
2412 if (p->ub - eps >= l)
2413 ret |= 0x10;
2414 else
2415 ret |= 0x20;
2416 }
2417 }
2418 done: return ret;
2419 }
2420
2421 /***********************************************************************
2422 * NAME
2423 *
2424 * npp_inactive_bound - remove row lower/upper inactive bound
2425 *
2426 * SYNOPSIS
2427 *
2428 * #include "glpnpp.h"
2429 * void npp_inactive_bound(NPP *npp, NPPROW *p, int which);
2430 *
2431 * DESCRIPTION
2432 *
2433 * The routine npp_inactive_bound removes lower (if which = 0) or upper
2434 * (if which = 1) bound of row p:
2435 *
2436 * L[p] <= sum a[p,j] x[j] <= U[p],
2437 *
2438 * which (bound) is assumed to be redundant.
2439 *
2440 * PROBLEM TRANSFORMATION
2441 *
2442 * If which = 0, current lower bound L[p] of row p is assigned -oo.
2443 * If which = 1, current upper bound U[p] of row p is assigned +oo.
2444 *
2445 * RECOVERING BASIC SOLUTION
2446 *
2447 * If in solution to the transformed problem row p is inactive
2448 * constraint (GLP_BS), its status is not changed in solution to the
2449 * original problem. Otherwise, status of row p in solution to the
2450 * original problem is defined by its type before transformation and
2451 * its status in solution to the transformed problem as follows:
2452 *
2453 * +---------------------+-------+---------------+---------------+
2454 * | Row | Flag | Row status in | Row status in |
2455 * | type | which | transfmd soln | original soln |
2456 * +---------------------+-------+---------------+---------------+
2457 * | sum >= L[p] | 0 | GLP_NF | GLP_NL |
2458 * | sum <= U[p] | 1 | GLP_NF | GLP_NU |
2459 * | L[p] <= sum <= U[p] | 0 | GLP_NU | GLP_NU |
2460 * | L[p] <= sum <= U[p] | 1 | GLP_NL | GLP_NL |
2461 * | sum = L[p] = U[p] | 0 | GLP_NU | GLP_NS |
2462 * | sum = L[p] = U[p] | 1 | GLP_NL | GLP_NS |
2463 * +---------------------+-------+---------------+---------------+
2464 *
2465 * RECOVERING INTERIOR-POINT SOLUTION
2466 *
2467 * None needed.
2468 *
2469 * RECOVERING MIP SOLUTION
2470 *
2471 * None needed. */
2472
2473 struct inactive_bound
2474 { /* row inactive bound */
2475 int p;
2476 /* row reference number */
2477 char stat;
2478 /* row status (if active constraint) */
2479 };
2480
2481 static int rcv_inactive_bound(NPP *npp, void *info);
2482
npp_inactive_bound(NPP * npp,NPPROW * p,int which)2483 void npp_inactive_bound(NPP *npp, NPPROW *p, int which)
2484 { /* remove row lower/upper inactive bound */
2485 struct inactive_bound *info;
2486 if (npp->sol == GLP_SOL)
2487 { /* create transformation stack entry */
2488 info = npp_push_tse(npp,
2489 rcv_inactive_bound, sizeof(struct inactive_bound));
2490 info->p = p->i;
2491 if (p->ub == +DBL_MAX)
2492 info->stat = GLP_NL;
2493 else if (p->lb == -DBL_MAX)
2494 info->stat = GLP_NU;
2495 else if (p->lb != p->ub)
2496 info->stat = (char)(which == 0 ? GLP_NU : GLP_NL);
2497 else
2498 info->stat = GLP_NS;
2499 }
2500 /* remove row inactive bound */
2501 if (which == 0)
2502 { xassert(p->lb != -DBL_MAX);
2503 p->lb = -DBL_MAX;
2504 }
2505 else if (which == 1)
2506 { xassert(p->ub != +DBL_MAX);
2507 p->ub = +DBL_MAX;
2508 }
2509 else
2510 xassert(which != which);
2511 return;
2512 }
2513
rcv_inactive_bound(NPP * npp,void * _info)2514 static int rcv_inactive_bound(NPP *npp, void *_info)
2515 { /* recover row status */
2516 struct inactive_bound *info = _info;
2517 if (npp->sol != GLP_SOL)
2518 { npp_error();
2519 return 1;
2520 }
2521 if (npp->r_stat[info->p] == GLP_BS)
2522 npp->r_stat[info->p] = GLP_BS;
2523 else
2524 npp->r_stat[info->p] = info->stat;
2525 return 0;
2526 }
2527
2528 /***********************************************************************
2529 * NAME
2530 *
2531 * npp_implied_bounds - determine implied column bounds
2532 *
2533 * SYNOPSIS
2534 *
2535 * #include "glpnpp.h"
2536 * void npp_implied_bounds(NPP *npp, NPPROW *p);
2537 *
2538 * DESCRIPTION
2539 *
2540 * The routine npp_implied_bounds inspects general row (constraint) p:
2541 *
2542 * L[p] <= sum a[p,j] x[j] <= U[p], (1)
2543 *
2544 * l[j] <= x[j] <= u[j], (2)
2545 *
2546 * where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0, to compute
2547 * implied bounds of columns (variables x[j]) in this row.
2548 *
2549 * The routine stores implied column bounds l'[j] and u'[j] in column
2550 * descriptors (NPPCOL); it does not change current column bounds l[j]
2551 * and u[j]. (Implied column bounds can be then used to strengthen the
2552 * current column bounds; see the routines npp_implied_lower and
2553 * npp_implied_upper).
2554 *
2555 * ALGORITHM
2556 *
2557 * Current column bounds (2) define implied lower and upper bounds of
2558 * row (1) as follows:
2559 *
2560 * L'[p] = inf sum a[p,j] x[j] =
2561 * j
2562 * (3)
2563 * = sum a[p,j] l[j] + sum a[p,j] u[j],
2564 * j in Jp j in Jn
2565 *
2566 * U'[p] = sup sum a[p,j] x[j] =
2567 * (4)
2568 * = sum a[p,j] u[j] + sum a[p,j] l[j],
2569 * j in Jp j in Jn
2570 *
2571 * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)
2572 *
2573 * (Note that bounds of all columns in row p are assumed to be correct,
2574 * so L'[p] <= U'[p].)
2575 *
2576 * If L[p] > L'[p] and/or U[p] < U'[p], the lower and/or upper bound of
2577 * row (1) can be active, in which case such row defines implied bounds
2578 * of its variables.
2579 *
2580 * Let x[k] be some variable having in row (1) coefficient a[p,k] != 0.
2581 * Consider a case when row lower bound can be active (L[p] > L'[p]):
2582 *
2583 * sum a[p,j] x[j] >= L[p] ==>
2584 * j
2585 *
2586 * sum a[p,j] x[j] + a[p,k] x[k] >= L[p] ==>
2587 * j!=k
2588 * (6)
2589 * a[p,k] x[k] >= L[p] - sum a[p,j] x[j] ==>
2590 * j!=k
2591 *
2592 * a[p,k] x[k] >= L[p,k],
2593 *
2594 * where
2595 *
2596 * L[p,k] = inf(L[p] - sum a[p,j] x[j]) =
2597 * j!=k
2598 *
2599 * = L[p] - sup sum a[p,j] x[j] = (7)
2600 * j!=k
2601 *
2602 * = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j].
2603 * j in Jp\{k} j in Jn\{k}
2604 *
2605 * Thus:
2606 *
2607 * x[k] >= l'[k] = L[p,k] / a[p,k], if a[p,k] > 0, (8)
2608 *
2609 * x[k] <= u'[k] = L[p,k] / a[p,k], if a[p,k] < 0. (9)
2610 *
2611 * where l'[k] and u'[k] are implied lower and upper bounds of variable
2612 * x[k], resp.
2613 *
2614 * Now consider a similar case when row upper bound can be active
2615 * (U[p] < U'[p]):
2616 *
2617 * sum a[p,j] x[j] <= U[p] ==>
2618 * j
2619 *
2620 * sum a[p,j] x[j] + a[p,k] x[k] <= U[p] ==>
2621 * j!=k
2622 * (10)
2623 * a[p,k] x[k] <= U[p] - sum a[p,j] x[j] ==>
2624 * j!=k
2625 *
2626 * a[p,k] x[k] <= U[p,k],
2627 *
2628 * where:
2629 *
2630 * U[p,k] = sup(U[p] - sum a[p,j] x[j]) =
2631 * j!=k
2632 *
2633 * = U[p] - inf sum a[p,j] x[j] = (11)
2634 * j!=k
2635 *
2636 * = U[p] - sum a[p,j] l[j] - sum a[p,j] u[j].
2637 * j in Jp\{k} j in Jn\{k}
2638 *
2639 * Thus:
2640 *
2641 * x[k] <= u'[k] = U[p,k] / a[p,k], if a[p,k] > 0, (12)
2642 *
2643 * x[k] >= l'[k] = U[p,k] / a[p,k], if a[p,k] < 0. (13)
2644 *
2645 * Note that in formulae (8), (9), (12), and (13) coefficient a[p,k]
2646 * must not be too small in magnitude relatively to other non-zero
2647 * coefficients in row (1), i.e. the following condition must hold:
2648 *
2649 * |a[p,k]| >= eps * max(1, |a[p,j]|), (14)
2650 * j
2651 *
2652 * where eps is a relative tolerance for constraint coefficients.
2653 * Otherwise the implied column bounds can be numerical inreliable. For
2654 * example, using formula (8) for the following inequality constraint:
2655 *
2656 * 1e-12 x1 - x2 - x3 >= 0,
2657 *
2658 * where x1 >= -1, x2, x3, >= 0, may lead to numerically unreliable
2659 * conclusion that x1 >= 0.
2660 *
2661 * Using formulae (8), (9), (12), and (13) to compute implied bounds
2662 * for one variable requires |J| operations, where J = {j: a[p,j] != 0},
2663 * because this needs computing L[p,k] and U[p,k]. Thus, computing
2664 * implied bounds for all variables in row (1) would require |J|^2
2665 * operations, that is not a good technique. However, the total number
2666 * of operations can be reduced to |J| as follows.
2667 *
2668 * Let a[p,k] > 0. Then from (7) and (11) we have:
2669 *
2670 * L[p,k] = L[p] - (U'[p] - a[p,k] u[k]) =
2671 *
2672 * = L[p] - U'[p] + a[p,k] u[k],
2673 *
2674 * U[p,k] = U[p] - (L'[p] - a[p,k] l[k]) =
2675 *
2676 * = U[p] - L'[p] + a[p,k] l[k],
2677 *
2678 * where L'[p] and U'[p] are implied row lower and upper bounds defined
2679 * by formulae (3) and (4). Substituting these expressions into (8) and
2680 * (12) gives:
2681 *
2682 * l'[k] = L[p,k] / a[p,k] = u[k] + (L[p] - U'[p]) / a[p,k], (15)
2683 *
2684 * u'[k] = U[p,k] / a[p,k] = l[k] + (U[p] - L'[p]) / a[p,k]. (16)
2685 *
2686 * Similarly, if a[p,k] < 0, according to (7) and (11) we have:
2687 *
2688 * L[p,k] = L[p] - (U'[p] - a[p,k] l[k]) =
2689 *
2690 * = L[p] - U'[p] + a[p,k] l[k],
2691 *
2692 * U[p,k] = U[p] - (L'[p] - a[p,k] u[k]) =
2693 *
2694 * = U[p] - L'[p] + a[p,k] u[k],
2695 *
2696 * and substituting these expressions into (8) and (12) gives:
2697 *
2698 * l'[k] = U[p,k] / a[p,k] = u[k] + (U[p] - L'[p]) / a[p,k], (17)
2699 *
2700 * u'[k] = L[p,k] / a[p,k] = l[k] + (L[p] - U'[p]) / a[p,k]. (18)
2701 *
2702 * Note that formulae (15)-(18) can be used only if L'[p] and U'[p]
2703 * exist. However, if for some variable x[j] it happens that l[j] = -oo
2704 * and/or u[j] = +oo, values of L'[p] (if a[p,j] > 0) and/or U'[p] (if
2705 * a[p,j] < 0) are undefined. Consider, therefore, the most general
2706 * situation, when some column bounds (2) may not exist.
2707 *
2708 * Let:
2709 *
2710 * J' = {j : (a[p,j] > 0 and l[j] = -oo) or
2711 * (19)
2712 * (a[p,j] < 0 and u[j] = +oo)}.
2713 *
2714 * Then (assuming that row upper bound U[p] can be active) the following
2715 * three cases are possible:
2716 *
2717 * 1) |J'| = 0. In this case L'[p] exists, thus, for all variables x[j]
2718 * in row (1) we can use formulae (16) and (17);
2719 *
2720 * 2) J' = {k}. In this case L'[p] = -oo, however, U[p,k] (11) exists,
2721 * so for variable x[k] we can use formulae (12) and (13). Note that
2722 * for all other variables x[j] (j != k) l'[j] = -oo (if a[p,j] < 0)
2723 * or u'[j] = +oo (if a[p,j] > 0);
2724 *
2725 * 3) |J'| > 1. In this case for all variables x[j] in row [1] we have
2726 * l'[j] = -oo (if a[p,j] < 0) or u'[j] = +oo (if a[p,j] > 0).
2727 *
2728 * Similarly, let:
2729 *
2730 * J'' = {j : (a[p,j] > 0 and u[j] = +oo) or
2731 * (20)
2732 * (a[p,j] < 0 and l[j] = -oo)}.
2733 *
2734 * Then (assuming that row lower bound L[p] can be active) the following
2735 * three cases are possible:
2736 *
2737 * 1) |J''| = 0. In this case U'[p] exists, thus, for all variables x[j]
2738 * in row (1) we can use formulae (15) and (18);
2739 *
2740 * 2) J'' = {k}. In this case U'[p] = +oo, however, L[p,k] (7) exists,
2741 * so for variable x[k] we can use formulae (8) and (9). Note that
2742 * for all other variables x[j] (j != k) l'[j] = -oo (if a[p,j] > 0)
2743 * or u'[j] = +oo (if a[p,j] < 0);
2744 *
2745 * 3) |J''| > 1. In this case for all variables x[j] in row (1) we have
2746 * l'[j] = -oo (if a[p,j] > 0) or u'[j] = +oo (if a[p,j] < 0). */
2747
npp_implied_bounds(NPP * npp,NPPROW * p)2748 void npp_implied_bounds(NPP *npp, NPPROW *p)
2749 { NPPAIJ *apj, *apk;
2750 double big, eps, temp;
2751 xassert(npp == npp);
2752 /* initialize implied bounds for all variables and determine
2753 maximal magnitude of row coefficients a[p,j] */
2754 big = 1.0;
2755 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
2756 { apj->col->ll.ll = -DBL_MAX, apj->col->uu.uu = +DBL_MAX;
2757 if (big < fabs(apj->val)) big = fabs(apj->val);
2758 }
2759 eps = 1e-6 * big;
2760 /* process row lower bound (assuming that it can be active) */
2761 if (p->lb != -DBL_MAX)
2762 { apk = NULL;
2763 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
2764 { if (apj->val > 0.0 && apj->col->ub == +DBL_MAX ||
2765 apj->val < 0.0 && apj->col->lb == -DBL_MAX)
2766 { if (apk == NULL)
2767 apk = apj;
2768 else
2769 goto skip1;
2770 }
2771 }
2772 /* if a[p,k] = NULL then |J'| = 0 else J' = { k } */
2773 temp = p->lb;
2774 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
2775 { if (apj == apk)
2776 /* skip a[p,k] */;
2777 else if (apj->val > 0.0)
2778 temp -= apj->val * apj->col->ub;
2779 else /* apj->val < 0.0 */
2780 temp -= apj->val * apj->col->lb;
2781 }
2782 /* compute column implied bounds */
2783 if (apk == NULL)
2784 { /* temp = L[p] - U'[p] */
2785 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
2786 { if (apj->val >= +eps)
2787 { /* l'[j] := u[j] + (L[p] - U'[p]) / a[p,j] */
2788 apj->col->ll.ll = apj->col->ub + temp / apj->val;
2789 }
2790 else if (apj->val <= -eps)
2791 { /* u'[j] := l[j] + (L[p] - U'[p]) / a[p,j] */
2792 apj->col->uu.uu = apj->col->lb + temp / apj->val;
2793 }
2794 }
2795 }
2796 else
2797 { /* temp = L[p,k] */
2798 if (apk->val >= +eps)
2799 { /* l'[k] := L[p,k] / a[p,k] */
2800 apk->col->ll.ll = temp / apk->val;
2801 }
2802 else if (apk->val <= -eps)
2803 { /* u'[k] := L[p,k] / a[p,k] */
2804 apk->col->uu.uu = temp / apk->val;
2805 }
2806 }
2807 skip1: ;
2808 }
2809 /* process row upper bound (assuming that it can be active) */
2810 if (p->ub != +DBL_MAX)
2811 { apk = NULL;
2812 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
2813 { if (apj->val > 0.0 && apj->col->lb == -DBL_MAX ||
2814 apj->val < 0.0 && apj->col->ub == +DBL_MAX)
2815 { if (apk == NULL)
2816 apk = apj;
2817 else
2818 goto skip2;
2819 }
2820 }
2821 /* if a[p,k] = NULL then |J''| = 0 else J'' = { k } */
2822 temp = p->ub;
2823 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
2824 { if (apj == apk)
2825 /* skip a[p,k] */;
2826 else if (apj->val > 0.0)
2827 temp -= apj->val * apj->col->lb;
2828 else /* apj->val < 0.0 */
2829 temp -= apj->val * apj->col->ub;
2830 }
2831 /* compute column implied bounds */
2832 if (apk == NULL)
2833 { /* temp = U[p] - L'[p] */
2834 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
2835 { if (apj->val >= +eps)
2836 { /* u'[j] := l[j] + (U[p] - L'[p]) / a[p,j] */
2837 apj->col->uu.uu = apj->col->lb + temp / apj->val;
2838 }
2839 else if (apj->val <= -eps)
2840 { /* l'[j] := u[j] + (U[p] - L'[p]) / a[p,j] */
2841 apj->col->ll.ll = apj->col->ub + temp / apj->val;
2842 }
2843 }
2844 }
2845 else
2846 { /* temp = U[p,k] */
2847 if (apk->val >= +eps)
2848 { /* u'[k] := U[p,k] / a[p,k] */
2849 apk->col->uu.uu = temp / apk->val;
2850 }
2851 else if (apk->val <= -eps)
2852 { /* l'[k] := U[p,k] / a[p,k] */
2853 apk->col->ll.ll = temp / apk->val;
2854 }
2855 }
2856 skip2: ;
2857 }
2858 return;
2859 }
2860
2861 /* eof */
2862