/* glpnpp03.c */

 /***********************************************************************

 *  This code is part of GLPK (GNU Linear Programming Kit).

 *

 *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,

 *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,

 *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.

 *  E-mail: <mao@gnu.org>.

 *

 *  GLPK is free software: you can redistribute it and/or modify it

 *  under the terms of the GNU General Public License as published by

 *  the Free Software Foundation, either version 3 of the License, or

 *  (at your option) any later version.

 *

 *  GLPK is distributed in the hope that it will be useful, but WITHOUT

 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

 *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

 *  License for more details.

 *

 *  You should have received a copy of the GNU General Public License

 *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.

 ***********************************************************************/

 #include "glpnpp.h"

 /***********************************************************************

 *  NAME

 *

 *  npp_empty_row - process empty row

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  int npp_empty_row(NPP *npp, NPPROW *p);

 *

 *  DESCRIPTION

 *

 *  The routine npp_empty_row processes row p, which is empty, i.e.

 *  coefficients at all columns in this row are zero:

 *

 *     L[p] <= sum 0 x[j] <= U[p],                                    (1)

 *

 *  where L[p] <= U[p].

 *

 *  RETURNS

 *

 *  0 - success;

 *

 *  1 - problem has no primal feasible solution.

 *

 *  PROBLEM TRANSFORMATION

 *

 *  If the following conditions hold:

 *

 *     L[p] <= +eps,  U[p] >= -eps,                                   (2)

 *

 *  where eps is an absolute tolerance for row value, the row p is

 *  redundant. In this case it can be replaced by equivalent redundant

 *  row, which is free (unbounded), and then removed from the problem.

 *  Otherwise, the row p is infeasible and, thus, the problem has no

 *  primal feasible solution.

 *

 *  RECOVERING BASIC SOLUTION

 *

 *  See the routine npp_free_row.

 *

 *  RECOVERING INTERIOR-POINT SOLUTION

 *

 *  See the routine npp_free_row.

 *

 *  RECOVERING MIP SOLUTION

 *

 *  None needed. */

 int npp_empty_row(NPP *npp, NPPROW *p)

 {     /* process empty row */

       double eps = 1e-3;

       /* the row must be empty */

       xassert(p->ptr == NULL);

       /* check primal feasibility */

       if (p->lb > +eps || p->ub < -eps)

          return 1;

       /* replace the row by equivalent free (unbounded) row */

       p->lb = -DBL_MAX, p->ub = +DBL_MAX;

       /* and process it */

       npp_free_row(npp, p);

       return 0;

 }

 /***********************************************************************

 *  NAME

 *

 *  npp_empty_col - process empty column

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  int npp_empty_col(NPP *npp, NPPCOL *q);

 *

 *  DESCRIPTION

 *

 *  The routine npp_empty_col processes column q:

 *

 *     l[q] <= x[q] <= u[q],                                          (1)

 *

 *  where l[q] <= u[q], which is empty, i.e. has zero coefficients in

 *  all constraint rows.

 *

 *  RETURNS

 *

 *  0 - success;

 *

 *  1 - problem has no dual feasible solution.

 *

 *  PROBLEM TRANSFORMATION

 *

 *  The row of the dual system corresponding to the empty column is the

 *  following:

 *

 *     sum 0 pi[i] + lambda[q] = c[q],                                (2)

 *      i

 *

 *  from which it follows that:

 *

 *     lambda[q] = c[q].                                              (3)

 *

 *  If the following condition holds:

 *

 *     c[q] < - eps,                                                  (4)

 *

 *  where eps is an absolute tolerance for column multiplier, the lower

 *  column bound l[q] must be active to provide dual feasibility (note

 *  that being preprocessed the problem is always minimization). In this

 *  case the column can be fixed on its lower bound and removed from the

 *  problem (if the column is integral, its bounds are also assumed to

 *  be integral). And if the column has no lower bound (l[q] = -oo), the

 *  problem has no dual feasible solution.

 *

 *  If the following condition holds:

 *

 *     c[q] > + eps,                                                  (5)

 *

 *  the upper column bound u[q] must be active to provide dual

 *  feasibility. In this case the column can be fixed on its upper bound

 *  and removed from the problem. And if the column has no upper bound

 *  (u[q] = +oo), the problem has no dual feasible solution.

 *

 *  Finally, if the following condition holds:

 *

 *     - eps <= c[q] <= +eps,                                         (6)

 *

 *  dual feasibility does not depend on a particular value of column q.

 *  In this case the column can be fixed either on its lower bound (if

 *  l[q] > -oo) or on its upper bound (if u[q] < +oo) or at zero (if the

 *  column is unbounded) and then removed from the problem.

 *

 *  RECOVERING BASIC SOLUTION

 *

 *  See the routine npp_fixed_col. Having been recovered the column

 *  is assigned status GLP_NS. However, if actually it is not fixed

 *  (l[q] < u[q]), its status should be changed to GLP_NL, GLP_NU, or

 *  GLP_NF depending on which bound it was fixed on transformation stage.

 *

 *  RECOVERING INTERIOR-POINT SOLUTION

 *

 *  See the routine npp_fixed_col.

 *

 *  RECOVERING MIP SOLUTION

 *

 *  See the routine npp_fixed_col. */

 struct empty_col

 {     /* empty column */

       int q;

       /* column reference number */

       char stat;

       /* status in basic solution */

 };

 static int rcv_empty_col(NPP *npp, void *info);

 int npp_empty_col(NPP *npp, NPPCOL *q)

 {     /* process empty column */

       struct empty_col *info;

       double eps = 1e-3;

       /* the column must be empty */

       xassert(q->ptr == NULL);

       /* check dual feasibility */

       if (q->coef > +eps && q->lb == -DBL_MAX)

          return 1;

       if (q->coef < -eps && q->ub == +DBL_MAX)

          return 1;

       /* create transformation stack entry */

       info = npp_push_tse(npp,

          rcv_empty_col, sizeof(struct empty_col));

       info->q = q->j;

       /* fix the column */

       if (q->lb == -DBL_MAX && q->ub == +DBL_MAX)

       {  /* free column */

          info->stat = GLP_NF;

          q->lb = q->ub = 0.0;

       }

       else if (q->ub == +DBL_MAX)

 lo:   {  /* column with lower bound */

          info->stat = GLP_NL;

          q->ub = q->lb;

       }

       else if (q->lb == -DBL_MAX)

 up:   {  /* column with upper bound */

          info->stat = GLP_NU;

          q->lb = q->ub;

       }

       else if (q->lb != q->ub)

       {  /* double-bounded column */

          if (q->coef >= +DBL_EPSILON) goto lo;

          if (q->coef <= -DBL_EPSILON) goto up;

          if (fabs(q->lb) <= fabs(q->ub)) goto lo; else goto up;

       }

       else

       {  /* fixed column */

          info->stat = GLP_NS;

       }

       /* process fixed column */

       npp_fixed_col(npp, q);

       return 0;

 }

 static int rcv_empty_col(NPP *npp, void *_info)

 {     /* recover empty column */

       struct empty_col *info = _info;

       if (npp->sol == GLP_SOL)

          npp->c_stat[info->q] = info->stat;

       return 0;

 }

 /***********************************************************************

 *  NAME

 *

 *  npp_implied_value - process implied column value

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  int npp_implied_value(NPP *npp, NPPCOL *q, double s);

 *

 *  DESCRIPTION

 *

 *  For column q:

 *

 *     l[q] <= x[q] <= u[q],                                          (1)

 *

 *  where l[q] < u[q], the routine npp_implied_value processes its

 *  implied value s[q]. If this implied value satisfies to the current

 *  column bounds and integrality condition, the routine fixes column q

 *  at the given point. Note that the column is kept in the problem in

 *  any case.

 *

 *  RETURNS

 *

 *  0 - column has been fixed;

 *

 *  1 - implied value violates to current column bounds;

 *

 *  2 - implied value violates integrality condition.

 *

 *  ALGORITHM

 *

 *  Implied column value s[q] satisfies to the current column bounds if

 *  the following condition holds:

 *

 *     l[q] - eps <= s[q] <= u[q] + eps,                              (2)

 *

 *  where eps is an absolute tolerance for column value. If the column

 *  is integral, the following condition also must hold:

 *

 *     |s[q] - floor(s[q]+0.5)| <= eps,                               (3)

 *

 *  where floor(s[q]+0.5) is the nearest integer to s[q].

 *

 *  If both condition (2) and (3) are satisfied, the column can be fixed

 *  at the value s[q], or, if it is integral, at floor(s[q]+0.5).

 *  Otherwise, if s[q] violates (2) or (3), the problem has no feasible

 *  solution.

 *

 *  Note: If s[q] is close to l[q] or u[q], it seems to be reasonable to

 *  fix the column at its lower or upper bound, resp. rather than at the

 *  implied value. */

 int npp_implied_value(NPP *npp, NPPCOL *q, double s)

 {     /* process implied column value */

       double eps, nint;

       xassert(npp == npp);

       /* column must not be fixed */

       xassert(q->lb < q->ub);

       /* check integrality */

       if (q->is_int)

       {  nint = floor(s + 0.5);

          if (fabs(s - nint) <= 1e-5)

             s = nint;

          else

             return 2;

       }

       /* check current column lower bound */

       if (q->lb != -DBL_MAX)

       {  eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb));

          if (s < q->lb - eps) return 1;

          /* if s[q] is close to l[q], fix column at its lower bound

             rather than at the implied value */

          if (s < q->lb + 1e-3 * eps)

          {  q->ub = q->lb;

             return 0;

          }

       }

       /* check current column upper bound */

       if (q->ub != +DBL_MAX)

       {  eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub));

          if (s > q->ub + eps) return 1;

          /* if s[q] is close to u[q], fix column at its upper bound

             rather than at the implied value */

          if (s > q->ub - 1e-3 * eps)

          {  q->lb = q->ub;

             return 0;

          }

       }

       /* fix column at the implied value */

       q->lb = q->ub = s;

       return 0;

 }

 /***********************************************************************

 *  NAME

 *

 *  npp_eq_singlet - process row singleton (equality constraint)

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  int npp_eq_singlet(NPP *npp, NPPROW *p);

 *

 *  DESCRIPTION

 *

 *  The routine npp_eq_singlet processes row p, which is equiality

 *  constraint having the only non-zero coefficient:

 *

 *     a[p,q] x[q] = b.                                               (1)

 *

 *  RETURNS

 *

 *  0 - success;

 *

 *  1 - problem has no primal feasible solution;

 *

 *  2 - problem has no integer feasible solution.

 *

 *  PROBLEM TRANSFORMATION

 *

 *  The equality constraint defines implied value of column q:

 *

 *     x[q] = s[q] = b / a[p,q].                                      (2)

 *

 *  If the implied value s[q] satisfies to the column bounds (see the

 *  routine npp_implied_value), the column can be fixed at s[q] and

 *  removed from the problem. In this case row p becomes redundant, so

 *  it can be replaced by equivalent free row and also removed from the

 *  problem.

 *

 *  Note that the routine removes from the problem only row p. Column q

 *  becomes fixed, however, it is kept in the problem.

 *

 *  RECOVERING BASIC SOLUTION

 *

 *  In solution to the original problem row p is assigned status GLP_NS

 *  (active equality constraint), and column q is assigned status GLP_BS

 *  (basic column).

 *

 *  Multiplier for row p can be computed as follows. In the dual system

 *  of the original problem column q corresponds to the following row:

 *

 *     sum a[i,q] pi[i] + lambda[q] = c[q]  ==>

 *      i

 *

 *     sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q].

 *     i!=p

 *

 *  Therefore:

 *

 *               1

 *     pi[p] = ------ (c[q] - lambda[q] - sum a[i,q] pi[i]),          (3)

 *             a[p,q]                     i!=q

 *

 *  where lambda[q] = 0 (since column[q] is basic), and pi[i] for all

 *  i != p are known in solution to the transformed problem.

 *

 *  Value of column q in solution to the original problem is assigned

 *  its implied value s[q].

 *

 *  RECOVERING INTERIOR-POINT SOLUTION

 *

 *  Multiplier for row p is computed with formula (3). Value of column

 *  q is assigned its implied value s[q].

 *

 *  RECOVERING MIP SOLUTION

 *

 *  Value of column q is assigned its implied value s[q]. */

 struct eq_singlet

 {     /* row singleton (equality constraint) */

       int p;

       /* row reference number */

       int q;

       /* column reference number */

       double apq;

       /* constraint coefficient a[p,q] */

       double c;

       /* objective coefficient at x[q] */

       NPPLFE *ptr;

       /* list of non-zero coefficients a[i,q], i != p */

 };

 static int rcv_eq_singlet(NPP *npp, void *info);

 int npp_eq_singlet(NPP *npp, NPPROW *p)

 {     /* process row singleton (equality constraint) */

       struct eq_singlet *info;

       NPPCOL *q;

       NPPAIJ *aij;

       NPPLFE *lfe;

       int ret;

       double s;

       /* the row must be singleton equality constraint */

       xassert(p->lb == p->ub);

       xassert(p->ptr != NULL && p->ptr->r_next == NULL);

       /* compute and process implied column value */

       aij = p->ptr;

       q = aij->col;

       s = p->lb / aij->val;

       ret = npp_implied_value(npp, q, s);

       xassert(0 <= ret && ret <= 2);

       if (ret != 0) return ret;

       /* create transformation stack entry */

       info = npp_push_tse(npp,

          rcv_eq_singlet, sizeof(struct eq_singlet));

       info->p = p->i;

       info->q = q->j;

       info->apq = aij->val;

       info->c = q->coef;

       info->ptr = NULL;

       /* save column coefficients a[i,q], i != p (not needed for MIP

          solution) */

       if (npp->sol != GLP_MIP)

       {  for (aij = q->ptr; aij != NULL; aij = aij->c_next)

          {  if (aij->row == p) continue; /* skip a[p,q] */

             lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));

             lfe->ref = aij->row->i;

             lfe->val = aij->val;

             lfe->next = info->ptr;

             info->ptr = lfe;

          }

       }

       /* remove the row from the problem */

       npp_del_row(npp, p);

       return 0;

 }

 static int rcv_eq_singlet(NPP *npp, void *_info)

 {     /* recover row singleton (equality constraint) */

       struct eq_singlet *info = _info;

       NPPLFE *lfe;

       double temp;

       if (npp->sol == GLP_SOL)

       {  /* column q must be already recovered as GLP_NS */

          if (npp->c_stat[info->q] != GLP_NS)

          {  npp_error();

             return 1;

          }

          npp->r_stat[info->p] = GLP_NS;

          npp->c_stat[info->q] = GLP_BS;

       }

       if (npp->sol != GLP_MIP)

       {  /* compute multiplier for row p with formula (3) */

          temp = info->c;

          for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)

             temp -= lfe->val * npp->r_pi[lfe->ref];

          npp->r_pi[info->p] = temp / info->apq;

       }

       return 0;

 }

 /***********************************************************************

 *  NAME

 *

 *  npp_implied_lower - process implied column lower bound

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  int npp_implied_lower(NPP *npp, NPPCOL *q, double l);

 *

 *  DESCRIPTION

 *

 *  For column q:

 *

 *     l[q] <= x[q] <= u[q],                                          (1)

 *

 *  where l[q] < u[q], the routine npp_implied_lower processes its

 *  implied lower bound l'[q]. As the result the current column lower

 *  bound may increase. Note that the column is kept in the problem in

 *  any case.

 *

 *  RETURNS

 *

 *  0 - current column lower bound has not changed;

 *

 *  1 - current column lower bound has changed, but not significantly;

 *

 *  2 - current column lower bound has significantly changed;

 *

 *  3 - column has been fixed on its upper bound;

 *

 *  4 - implied lower bound violates current column upper bound.

 *

 *  ALGORITHM

 *

 *  If column q is integral, before processing its implied lower bound

 *  should be rounded up:

 *

 *              ( floor(l'[q]+0.5), if |l'[q] - floor(l'[q]+0.5)| <= eps

 *     l'[q] := <                                                     (2)

 *              ( ceil(l'[q]),      otherwise

 *

 *  where floor(l'[q]+0.5) is the nearest integer to l'[q], ceil(l'[q])

 *  is smallest integer not less than l'[q], and eps is an absolute

 *  tolerance for column value.

 *

 *  Processing implied column lower bound l'[q] includes the following

 *  cases:

 *

 *  1) if l'[q] < l[q] + eps, implied lower bound is redundant;

 *

 *  2) if l[q] + eps <= l[q] <= u[q] + eps, current column lower bound

 *     l[q] can be strengthened by replacing it with l'[q]. If in this

 *     case new column lower bound becomes close to current column upper

 *     bound u[q], the column can be fixed on its upper bound;

 *

 *  3) if l'[q] > u[q] + eps, implied lower bound violates current

 *     column upper bound u[q], in which case the problem has no primal

 *     feasible solution. */

 int npp_implied_lower(NPP *npp, NPPCOL *q, double l)

 {     /* process implied column lower bound */

       int ret;

       double eps, nint;

       xassert(npp == npp);

       /* column must not be fixed */

       xassert(q->lb < q->ub);

       /* implied lower bound must be finite */

       xassert(l != -DBL_MAX);

       /* if column is integral, round up l'[q] */

       if (q->is_int)

       {  nint = floor(l + 0.5);

          if (fabs(l - nint) <= 1e-5)

             l = nint;

          else

             l = ceil(l);

       }

       /* check current column lower bound */

       if (q->lb != -DBL_MAX)

       {  eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->lb));

          if (l < q->lb + eps)

          {  ret = 0; /* redundant */

             goto done;

          }

       }

       /* check current column upper bound */

       if (q->ub != +DBL_MAX)

       {  eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub));

          if (l > q->ub + eps)

          {  ret = 4; /* infeasible */

             goto done;

          }

          /* if l'[q] is close to u[q], fix column at its upper bound */

          if (l > q->ub - 1e-3 * eps)

          {  q->lb = q->ub;

             ret = 3; /* fixed */

             goto done;

          }

       }

       /* check if column lower bound changes significantly */

       if (q->lb == -DBL_MAX)

          ret = 2; /* significantly */

       else if (q->is_int && l > q->lb + 0.5)

          ret = 2; /* significantly */

       else if (l > q->lb + 0.30 * (1.0 + fabs(q->lb)))

          ret = 2; /* significantly */

       else

          ret = 1; /* not significantly */

       /* set new column lower bound */

       q->lb = l;

 done: return ret;

 }

 /***********************************************************************

 *  NAME

 *

 *  npp_implied_upper - process implied column upper bound

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  int npp_implied_upper(NPP *npp, NPPCOL *q, double u);

 *

 *  DESCRIPTION

 *

 *  For column q:

 *

 *     l[q] <= x[q] <= u[q],                                          (1)

 *

 *  where l[q] < u[q], the routine npp_implied_upper processes its

 *  implied upper bound u'[q]. As the result the current column upper

 *  bound may decrease. Note that the column is kept in the problem in

 *  any case.

 *

 *  RETURNS

 *

 *  0 - current column upper bound has not changed;

 *

 *  1 - current column upper bound has changed, but not significantly;

 *

 *  2 - current column upper bound has significantly changed;

 *

 *  3 - column has been fixed on its lower bound;

 *

 *  4 - implied upper bound violates current column lower bound.

 *

 *  ALGORITHM

 *

 *  If column q is integral, before processing its implied upper bound

 *  should be rounded down:

 *

 *              ( floor(u'[q]+0.5), if |u'[q] - floor(l'[q]+0.5)| <= eps

 *     u'[q] := <                                                     (2)

 *              ( floor(l'[q]),     otherwise

 *

 *  where floor(u'[q]+0.5) is the nearest integer to u'[q],

 *  floor(u'[q]) is largest integer not greater than u'[q], and eps is

 *  an absolute tolerance for column value.

 *

 *  Processing implied column upper bound u'[q] includes the following

 *  cases:

 *

 *  1) if u'[q] > u[q] - eps, implied upper bound is redundant;

 *

 *  2) if l[q] - eps <= u[q] <= u[q] - eps, current column upper bound

 *     u[q] can be strengthened by replacing it with u'[q]. If in this

 *     case new column upper bound becomes close to current column lower

 *     bound, the column can be fixed on its lower bound;

 *

 *  3) if u'[q] < l[q] - eps, implied upper bound violates current

 *     column lower bound l[q], in which case the problem has no primal

 *     feasible solution. */

 int npp_implied_upper(NPP *npp, NPPCOL *q, double u)

 {     int ret;

       double eps, nint;

       xassert(npp == npp);

       /* column must not be fixed */

       xassert(q->lb < q->ub);

       /* implied upper bound must be finite */

       xassert(u != +DBL_MAX);

       /* if column is integral, round down u'[q] */

       if (q->is_int)

       {  nint = floor(u + 0.5);

          if (fabs(u - nint) <= 1e-5)

             u = nint;

          else

             u = floor(u);

       }

       /* check current column upper bound */

       if (q->ub != +DBL_MAX)

       {  eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->ub));

          if (u > q->ub - eps)

          {  ret = 0; /* redundant */

             goto done;

          }

       }

       /* check current column lower bound */

       if (q->lb != -DBL_MAX)

       {  eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb));

          if (u < q->lb - eps)

          {  ret = 4; /* infeasible */

             goto done;

          }

          /* if u'[q] is close to l[q], fix column at its lower bound */

          if (u < q->lb + 1e-3 * eps)

          {  q->ub = q->lb;

             ret = 3; /* fixed */

             goto done;

          }

       }

       /* check if column upper bound changes significantly */

       if (q->ub == +DBL_MAX)

          ret = 2; /* significantly */

       else if (q->is_int && u < q->ub - 0.5)

          ret = 2; /* significantly */

       else if (u < q->ub - 0.30 * (1.0 + fabs(q->ub)))

          ret = 2; /* significantly */

       else

          ret = 1; /* not significantly */

       /* set new column upper bound */

       q->ub = u;

 done: return ret;

 }

 /***********************************************************************

 *  NAME

 *

 *  npp_ineq_singlet - process row singleton (inequality constraint)

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  int npp_ineq_singlet(NPP *npp, NPPROW *p);

 *

 *  DESCRIPTION

 *

 *  The routine npp_ineq_singlet processes row p, which is inequality

 *  constraint having the only non-zero coefficient:

 *

 *     L[p] <= a[p,q] * x[q] <= U[p],                                 (1)

 *

 *  where L[p] < U[p], L[p] > -oo and/or U[p] < +oo.

 *

 *  RETURNS

 *

 *  0 - current column bounds have not changed;

 *

 *  1 - current column bounds have changed, but not significantly;

 *

 *  2 - current column bounds have significantly changed;

 *

 *  3 - column has been fixed on its lower or upper bound;

 *

 *  4 - problem has no primal feasible solution.

 *

 *  PROBLEM TRANSFORMATION

 *

 *  Inequality constraint (1) defines implied bounds of column q:

 *

 *             (  L[p] / a[p,q],  if a[p,q] > 0

 *     l'[q] = <                                                      (2)

 *             (  U[p] / a[p,q],  if a[p,q] < 0

 *

 *             (  U[p] / a[p,q],  if a[p,q] > 0

 *     u'[q] = <                                                      (3)

 *             (  L[p] / a[p,q],  if a[p,q] < 0

 *

 *  If these implied bounds do not violate current bounds of column q:

 *

 *     l[q] <= x[q] <= u[q],                                          (4)

 *

 *  they can be used to strengthen the current column bounds:

 *

 *     l[q] := max(l[q], l'[q]),                                      (5)

 *

 *     u[q] := min(u[q], u'[q]).                                      (6)

 *

 *  (See the routines npp_implied_lower and npp_implied_upper.)

 *

 *  Once bounds of row p (1) have been carried over column q, the row

 *  becomes redundant, so it can be replaced by equivalent free row and

 *  removed from the problem.

 *

 *  Note that the routine removes from the problem only row p. Column q,

 *  even it has been fixed, is kept in the problem.

 *

 *  RECOVERING BASIC SOLUTION

 *

 *  Note that the row in the dual system corresponding to column q is

 *  the following:

 *

 *     sum a[i,q] pi[i] + lambda[q] = c[q]  ==>

 *      i

 *                                                                    (7)

 *     sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q],

 *     i!=p

 *

 *  where pi[i] for all i != p are known in solution to the transformed

 *  problem. Row p does not exist in the transformed problem, so it has

 *  zero multiplier there. This allows computing multiplier for column q

 *  in solution to the transformed problem:

 *

 *     lambda~[q] = c[q] - sum a[i,q] pi[i].                          (8)

 *                         i!=p

 *

 *  Let in solution to the transformed problem column q be non-basic

 *  with lower bound active (GLP_NL, lambda~[q] >= 0), and this lower

 *  bound be implied one l'[q]. From the original problem's standpoint

 *  this then means that actually the original column lower bound l[q]

 *  is inactive, and active is that row bound L[p] or U[p] that defines

 *  the implied bound l'[q] (2). In this case in solution to the

 *  original problem column q is assigned status GLP_BS while row p is

 *  assigned status GLP_NL (if a[p,q] > 0) or GLP_NU (if a[p,q] < 0).

 *  Since now column q is basic, its multiplier lambda[q] is zero. This

 *  allows using (7) and (8) to find multiplier for row p in solution to

 *  the original problem:

 *

 *               1

 *     pi[p] = ------ (c[q] - sum a[i,q] pi[i]) = lambda~[q] / a[p,q] (9)

 *             a[p,q]         i!=p

 *

 *  Now let in solution to the transformed problem column q be non-basic

 *  with upper bound active (GLP_NU, lambda~[q] <= 0), and this upper

 *  bound be implied one u'[q]. As in the previous case this then means

 *  that from the original problem's standpoint actually the original

 *  column upper bound u[q] is inactive, and active is that row bound

 *  L[p] or U[p] that defines the implied bound u'[q] (3). In this case

 *  in solution to the original problem column q is assigned status

 *  GLP_BS, row p is assigned status GLP_NU (if a[p,q] > 0) or GLP_NL

 *  (if a[p,q] < 0), and its multiplier is computed with formula (9).

 *

 *  Strengthening bounds of column q according to (5) and (6) may make

 *  it fixed. Thus, if in solution to the transformed problem column q is

 *  non-basic and fixed (GLP_NS), we can suppose that if lambda~[q] > 0,

 *  column q has active lower bound (GLP_NL), and if lambda~[q] < 0,

 *  column q has active upper bound (GLP_NU), reducing this case to two

 *  previous ones. If, however, lambda~[q] is close to zero or

 *  corresponding bound of row p does not exist (this may happen if

 *  lambda~[q] has wrong sign due to round-off errors, in which case it

 *  is expected to be close to zero, since solution is assumed to be dual

 *  feasible), column q can be assigned status GLP_BS (basic), and row p

 *  can be made active on its existing bound. In the latter case row

 *  multiplier pi[p] computed with formula (9) will be also close to

 *  zero, and dual feasibility will be kept.

 *

 *  In all other cases, namely, if in solution to the transformed

 *  problem column q is basic (GLP_BS), or non-basic with original lower

 *  bound l[q] active (GLP_NL), or non-basic with original upper bound

 *  u[q] active (GLP_NU), constraint (1) is inactive. So in solution to

 *  the original problem status of column q remains unchanged, row p is

 *  assigned status GLP_BS, and its multiplier pi[p] is assigned zero

 *  value.

 *

 *  RECOVERING INTERIOR-POINT SOLUTION

 *

 *  First, value of multiplier for column q in solution to the original

 *  problem is computed with formula (8). If lambda~[q] > 0 and column q

 *  has implied lower bound, or if lambda~[q] < 0 and column q has

 *  implied upper bound, this means that from the original problem's

 *  standpoint actually row p has corresponding active bound, in which

 *  case its multiplier pi[p] is computed with formula (9). In other

 *  cases, when the sign of lambda~[q] corresponds to original bound of

 *  column q, or when lambda~[q] =~ 0, value of row multiplier pi[p] is

 *  assigned zero value.

 *

 *  RECOVERING MIP SOLUTION

 *

 *  None needed. */

 struct ineq_singlet

 {     /* row singleton (inequality constraint) */

       int p;

       /* row reference number */

       int q;

       /* column reference number */

       double apq;

       /* constraint coefficient a[p,q] */

       double c;

       /* objective coefficient at x[q] */

       double lb;

       /* row lower bound */

       double ub;

       /* row upper bound */

       char lb_changed;

       /* this flag is set if column lower bound was changed */

       char ub_changed;

       /* this flag is set if column upper bound was changed */

       NPPLFE *ptr;

       /* list of non-zero coefficients a[i,q], i != p */

 };

 static int rcv_ineq_singlet(NPP *npp, void *info);

 int npp_ineq_singlet(NPP *npp, NPPROW *p)

 {     /* process row singleton (inequality constraint) */

       struct ineq_singlet *info;

       NPPCOL *q;

       NPPAIJ *apq, *aij;

       NPPLFE *lfe;

       int lb_changed, ub_changed;

       double ll, uu;

       /* the row must be singleton inequality constraint */

       xassert(p->lb != -DBL_MAX || p->ub != +DBL_MAX);

       xassert(p->lb < p->ub);

       xassert(p->ptr != NULL && p->ptr->r_next == NULL);

       /* compute implied column bounds */

       apq = p->ptr;

       q = apq->col;

       xassert(q->lb < q->ub);

       if (apq->val > 0.0)

       {  ll = (p->lb == -DBL_MAX ? -DBL_MAX : p->lb / apq->val);

          uu = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub / apq->val);

       }

       else

       {  ll = (p->ub == +DBL_MAX ? -DBL_MAX : p->ub / apq->val);

          uu = (p->lb == -DBL_MAX ? +DBL_MAX : p->lb / apq->val);

       }

       /* process implied column lower bound */

       if (ll == -DBL_MAX)

          lb_changed = 0;

       else

       {  lb_changed = npp_implied_lower(npp, q, ll);

          xassert(0 <= lb_changed && lb_changed <= 4);

          if (lb_changed == 4) return 4; /* infeasible */

       }

       /* process implied column upper bound */

       if (uu == +DBL_MAX)

          ub_changed = 0;

       else if (lb_changed == 3)

       {  /* column was fixed on its upper bound due to l'[q] = u[q] */

          /* note that L[p] < U[p], so l'[q] = u[q] < u'[q] */

          ub_changed = 0;

       }

       else

       {  ub_changed = npp_implied_upper(npp, q, uu);

          xassert(0 <= ub_changed && ub_changed <= 4);

          if (ub_changed == 4) return 4; /* infeasible */

       }

       /* if neither lower nor upper column bound was changed, the row

          is originally redundant and can be replaced by free row */

       if (!lb_changed && !ub_changed)

       {  p->lb = -DBL_MAX, p->ub = +DBL_MAX;

          npp_free_row(npp, p);

          return 0;

       }

       /* create transformation stack entry */

       info = npp_push_tse(npp,

          rcv_ineq_singlet, sizeof(struct ineq_singlet));

       info->p = p->i;

       info->q = q->j;

       info->apq = apq->val;

       info->c = q->coef;

       info->lb = p->lb;

       info->ub = p->ub;

       info->lb_changed = (char)lb_changed;

       info->ub_changed = (char)ub_changed;

       info->ptr = NULL;

       /* save column coefficients a[i,q], i != p (not needed for MIP

          solution) */

       if (npp->sol != GLP_MIP)

       {  for (aij = q->ptr; aij != NULL; aij = aij->c_next)

          {  if (aij == apq) continue; /* skip a[p,q] */

             lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));

             lfe->ref = aij->row->i;

             lfe->val = aij->val;

             lfe->next = info->ptr;

             info->ptr = lfe;

          }

       }

       /* remove the row from the problem */

       npp_del_row(npp, p);

       return lb_changed >= ub_changed ? lb_changed : ub_changed;

 }

 static int rcv_ineq_singlet(NPP *npp, void *_info)

 {     /* recover row singleton (inequality constraint) */

       struct ineq_singlet *info = _info;

       NPPLFE *lfe;

       double lambda;

       if (npp->sol == GLP_MIP) goto done;

       /* compute lambda~[q] in solution to the transformed problem

          with formula (8) */

       lambda = info->c;

       for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)

          lambda -= lfe->val * npp->r_pi[lfe->ref];

       if (npp->sol == GLP_SOL)

       {  /* recover basic solution */

          if (npp->c_stat[info->q] == GLP_BS)

          {  /* column q is basic, so row p is inactive */

             npp->r_stat[info->p] = GLP_BS;

             npp->r_pi[info->p] = 0.0;

          }

          else if (npp->c_stat[info->q] == GLP_NL)

 nl:      {  /* column q is non-basic with lower bound active */

             if (info->lb_changed)

             {  /* it is implied bound, so actually row p is active

                   while column q is basic */

                npp->r_stat[info->p] =

                   (char)(info->apq > 0.0 ? GLP_NL : GLP_NU);

                npp->c_stat[info->q] = GLP_BS;

                npp->r_pi[info->p] = lambda / info->apq;

             }

             else

             {  /* it is original bound, so row p is inactive */

                npp->r_stat[info->p] = GLP_BS;

                npp->r_pi[info->p] = 0.0;

             }

          }

          else if (npp->c_stat[info->q] == GLP_NU)

 nu:      {  /* column q is non-basic with upper bound active */

             if (info->ub_changed)

             {  /* it is implied bound, so actually row p is active

                   while column q is basic */

                npp->r_stat[info->p] =

                   (char)(info->apq > 0.0 ? GLP_NU : GLP_NL);

                npp->c_stat[info->q] = GLP_BS;

                npp->r_pi[info->p] = lambda / info->apq;

             }

             else

             {  /* it is original bound, so row p is inactive */

                npp->r_stat[info->p] = GLP_BS;

                npp->r_pi[info->p] = 0.0;

             }

          }

          else if (npp->c_stat[info->q] == GLP_NS)

          {  /* column q is non-basic and fixed; note, however, that in

                in the original problem it is non-fixed */

             if (lambda > +1e-7)

             {  if (info->apq > 0.0 && info->lb != -DBL_MAX ||

                    info->apq < 0.0 && info->ub != +DBL_MAX ||

                   !info->lb_changed)

                {  /* either corresponding bound of row p exists or

                      column q remains non-basic with its original lower

                      bound active */

                   npp->c_stat[info->q] = GLP_NL;

                   goto nl;

                }

             }

             if (lambda < -1e-7)

             {  if (info->apq > 0.0 && info->ub != +DBL_MAX ||

                    info->apq < 0.0 && info->lb != -DBL_MAX ||

                   !info->ub_changed)

                {  /* either corresponding bound of row p exists or

                      column q remains non-basic with its original upper

                      bound active */

                   npp->c_stat[info->q] = GLP_NU;

                   goto nu;

                }

             }

             /* either lambda~[q] is close to zero, or corresponding

                bound of row p does not exist, because lambda~[q] has

                wrong sign due to round-off errors; in the latter case

                lambda~[q] is also assumed to be close to zero; so, we

                can make row p active on its existing bound and column q

                basic; pi[p] will have wrong sign, but it also will be

                close to zero (rarus casus of dual degeneracy) */

             if (info->lb != -DBL_MAX && info->ub == +DBL_MAX)

             {  /* row lower bound exists, but upper bound doesn't */

                npp->r_stat[info->p] = GLP_NL;

             }

             else if (info->lb == -DBL_MAX && info->ub != +DBL_MAX)

             {  /* row upper bound exists, but lower bound doesn't */

                npp->r_stat[info->p] = GLP_NU;

             }

             else if (info->lb != -DBL_MAX && info->ub != +DBL_MAX)

             {  /* both row lower and upper bounds exist */

                /* to choose proper active row bound we should not use

                   lambda~[q], because its value being close to zero is

                   unreliable; so we choose that bound which provides

                   primal feasibility for original constraint (1) */

                if (info->apq * npp->c_value[info->q] <=

                    0.5 * (info->lb + info->ub))

                   npp->r_stat[info->p] = GLP_NL;

                else

                   npp->r_stat[info->p] = GLP_NU;

             }

             else

             {  npp_error();

                return 1;

             }

             npp->c_stat[info->q] = GLP_BS;

             npp->r_pi[info->p] = lambda / info->apq;

          }

          else

          {  npp_error();

             return 1;

          }

       }

       if (npp->sol == GLP_IPT)

       {  /* recover interior-point solution */

          if (lambda > +DBL_EPSILON && info->lb_changed ||

              lambda < -DBL_EPSILON && info->ub_changed)

          {  /* actually row p has corresponding active bound */

             npp->r_pi[info->p] = lambda / info->apq;

          }

          else

          {  /* either bounds of column q are both inactive or its

                original bound is active */

             npp->r_pi[info->p] = 0.0;

          }

       }

 done: return 0;

 }

 /***********************************************************************

 *  NAME

 *

 *  npp_implied_slack - process column singleton (implied slack variable)

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  void npp_implied_slack(NPP *npp, NPPCOL *q);

 *

 *  DESCRIPTION

 *

 *  The routine npp_implied_slack processes column q:

 *

 *     l[q] <= x[q] <= u[q],                                          (1)

 *

 *  where l[q] < u[q], having the only non-zero coefficient in row p,

 *  which is equality constraint:

 *

 *     sum a[p,j] x[j] + a[p,q] x[q] = b.                             (2)

 *     j!=q

 *

 *  PROBLEM TRANSFORMATION

 *

 *  (If x[q] is integral, this transformation must not be used.)

 *

 *  The term a[p,q] x[q] in constraint (2) can be considered as a slack

 *  variable that allows to carry bounds of column q over row p and then

 *  remove column q from the problem.

 *

 *  Constraint (2) can be written as follows:

 *

 *     sum a[p,j] x[j] = b - a[p,q] x[q].                             (3)

 *     j!=q

 *

 *  According to (1) constraint (3) is equivalent to the following

 *  inequality constraint:

 *

 *     L[p] <= sum a[p,j] x[j] <= U[p],                               (4)

 *             j!=q

 *

 *  where