COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpnpp03.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 96.5 KB
Line 
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
76int 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
173struct empty_col
174{     /* empty column */
175      int q;
176      /* column reference number */
177      char stat;
178      /* status in basic solution */
179};
180
181static int rcv_empty_col(NPP *npp, void *info);
182
183int 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)
205lo:   {  /* column with lower bound */
206         info->stat = GLP_NL;
207         q->ub = q->lb;
208      }
209      else if (q->lb == -DBL_MAX)
210up:   {  /* 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
229static 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
290int 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
407struct 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
421static int rcv_eq_singlet(NPP *npp, void *info);
422
423int 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
466static 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
550int 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;
600done: 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
663int 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;
712done: 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
860struct 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
882static int rcv_ineq_singlet(NPP *npp, void *info);
883
884int 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
965static 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)
984nl:      {  /* 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)
1000nu:      {  /* 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      }
1092done: 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
1221struct 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
1237static int rcv_implied_slack(NPP *npp, void *info);
1238
1239void 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
1292static 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
1477struct 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
1487static int rcv_implied_free(NPP *npp, void *info);
1488
1489int 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)
1578nl:      {  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)
1591nu:      {  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
1616static 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
1751struct 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
1761static int rcv_eq_doublet(NPP *npp, void *info);
1762
1763NPPCOL *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
1854static 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
2057struct 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
2075struct 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
2089static int rcv_forcing_row(NPP *npp, void *info);
2090
2091int 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
2177static 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      }
2242done: 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
2338int 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      }
2418done: 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
2473struct inactive_bound
2474{     /* row inactive bound */
2475      int p;
2476      /* row reference number */
2477      char stat;
2478      /* row status (if active constraint) */
2479};
2480
2481static int rcv_inactive_bound(NPP *npp, void *info);
2482
2483void 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
2514static 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
2748void 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         }
2807skip1:   ;
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         }
2856skip2:   ;
2857      }
2858      return;
2859}
2860
2861/* eof */
Note: See TracBrowser for help on using the repository browser.