/* 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

 *

 *            ( b - a[p,q] u[q],  if a[p,q] > 0

 *     L[p] = <                                                       (5)

 *            ( b - a[p,q] l[q],  if a[p,q] < 0

 *

 *            ( b - a[p,q] l[q],  if a[p,q] > 0

 *     U[p] = <                                                       (6)

 *            ( b - a[p,q] u[q],  if a[p,q] < 0

 *

 *  From (2) it follows that:

 *

 *              1

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

 *            a[p,q]      j!=q

 *

 *  In order to eliminate x[q] from the objective row we substitute it

 *  from (6) to that row:

 *

 *     z = sum c[j] x[j] + c[q] x[q] + c[0] =

 *         j!=q

 *                                 1

 *       = sum c[j] x[j] + c[q] [------ (b - sum a[p,j] x[j])] + c0 =

 *         j!=q                  a[p,q]      j!=q

 *

 *       = sum c~[j] x[j] + c~[0],

 *         j!=q

 *                         a[p,j]                     b

 *     c~[j] = c[j] - c[q] ------,  c~0 = c0 - c[q] ------            (8)

 *                         a[p,q]                   a[p,q]

 *

 *  are values of objective coefficients and constant term, resp., in

 *  the transformed problem.

 *

 *  Note that column q is column singleton, so in the dual system of the

 *  original problem it corresponds to the following row singleton:

 *

 *     a[p,q] pi[p] + lambda[q] = c[q].                               (9)

 *

 *  In the transformed problem row (9) would be the following:

 *

 *     a[p,q] pi~[p] + lambda[q] = c~[q] = 0.                        (10)

 *

 *  Subtracting (10) from (9) we have:

 *

 *     a[p,q] (pi[p] - pi~[p]) = c[q]

 *

 *  that gives the following formula to compute multiplier for row p in

 *  solution to the original problem using its value in solution to the

 *  transformed problem:

 *

 *     pi[p] = pi~[p] + c[q] / a[p,q].                               (11)

 *

 *  RECOVERING BASIC SOLUTION

 *

 *  Status of column q in solution to the original problem is defined

 *  by status of row p in solution to the transformed problem and the

 *  sign of coefficient a[p,q] in the original inequality constraint (2)

 *  as follows:

 *

 *     +-----------------------+---------+--------------------+

 *     |    Status of row p    | Sign of | Status of column q |

 *     | (transformed problem) | a[p,q]  | (original problem) |

 *     +-----------------------+---------+--------------------+

 *     |        GLP_BS         |  + / -  |       GLP_BS       |

 *     |        GLP_NL         |    +    |       GLP_NU       |

 *     |        GLP_NL         |    -    |       GLP_NL       |

 *     |        GLP_NU         |    +    |       GLP_NL       |

 *     |        GLP_NU         |    -    |       GLP_NU       |

 *     |        GLP_NF         |  + / -  |       GLP_NF       |

 *     +-----------------------+---------+--------------------+

 *

 *  Value of column q is computed with formula (7). Since originally row

 *  p is equality constraint, its status is assigned GLP_NS, and value of

 *  its multiplier pi[p] is computed with formula (11).

 *

 *  RECOVERING INTERIOR-POINT SOLUTION

 *

 *  Value of column q is computed with formula (7). Row multiplier value

 *  pi[p] is computed with formula (11).

 *

 *  RECOVERING MIP SOLUTION

 *

 *  Value of column q is computed with formula (7). */

 struct implied_slack

 {     /* column singleton (implied slack variable) */

       int p;

       /* row reference number */

       int q;

       /* column reference number */

       double apq;

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

       double b;

       /* right-hand side of original equality constraint */

       double c;

       /* original objective coefficient at x[q] */

       NPPLFE *ptr;

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

 };

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

 void npp_implied_slack(NPP *npp, NPPCOL *q)

 {     /* process column singleton (implied slack variable) */

       struct implied_slack *info;

       NPPROW *p;

       NPPAIJ *aij;

       NPPLFE *lfe;

       /* the column must be non-integral non-fixed singleton */

       xassert(!q->is_int);

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

       xassert(q->ptr != NULL && q->ptr->c_next == NULL);

       /* corresponding row must be equality constraint */

       aij = q->ptr;

       p = aij->row;

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

       /* create transformation stack entry */

       info = npp_push_tse(npp,

          rcv_implied_slack, sizeof(struct implied_slack));

       info->p = p->i;

       info->q = q->j;

       info->apq = aij->val;

       info->b = p->lb;

       info->c = q->coef;

       info->ptr = NULL;

       /* save row coefficients a[p,j], j != q, and substitute x[q]

          into the objective row */

       for (aij = p->ptr; aij != NULL; aij = aij->r_next)

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

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

          lfe->ref = aij->col->j;

          lfe->val = aij->val;

          lfe->next = info->ptr;

          info->ptr = lfe;

          aij->col->coef -= info->c * (aij->val / info->apq);

       }

       npp->c0 += info->c * (info->b / info->apq);

       /* compute new row bounds */

       if (info->apq > 0.0)

       {  p->lb = (q->ub == +DBL_MAX ?

             -DBL_MAX : info->b - info->apq * q->ub);

          p->ub = (q->lb == -DBL_MAX ?

             +DBL_MAX : info->b - info->apq * q->lb);

       }

       else

       {  p->lb = (q->lb == -DBL_MAX ?

             -DBL_MAX : info->b - info->apq * q->lb);

          p->ub = (q->ub == +DBL_MAX ?

             +DBL_MAX : info->b - info->apq * q->ub);

       }

       /* remove the column from the problem */

       npp_del_col(npp, q);

       return;

 }

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

 {     /* recover column singleton (implied slack variable) */

       struct implied_slack *info = _info;

       NPPLFE *lfe;

       double temp;

       if (npp->sol == GLP_SOL)

       {  /* assign statuses to row p and column q */

          if (npp->r_stat[info->p] == GLP_BS ||

              npp->r_stat[info->p] == GLP_NF)

             npp->c_stat[info->q] = npp->r_stat[info->p];

          else if (npp->r_stat[info->p] == GLP_NL)

             npp->c_stat[info->q] =

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

          else if (npp->r_stat[info->p] == GLP_NU)

             npp->c_stat[info->q] =

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

          else

          {  npp_error();

             return 1;

          }

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

       }

       if (npp->sol != GLP_MIP)

       {  /* compute multiplier for row p */

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

       }

       /* compute value of column q */

       temp = info->b;

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

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

       npp->c_value[info->q] = temp / info->apq;

       return 0;

 }

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

 *  NAME

 *

 *  npp_implied_free - process column singleton (implied free variable)

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

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

 *

 *  DESCRIPTION

 *

 *  The routine npp_implied_free processes column q:

 *

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

 *

 *  having non-zero coefficient in the only row p, which is inequality

 *  constraint:

 *

 *     L[p] <= sum a[p,j] x[j] + a[p,q] x[q] <= U[p],                 (2)

 *             j!=q

 *

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

 *

 *  RETURNS

 *

 *  0 - success;

 *

 *  1 - column lower and/or upper bound(s) can be active;

 *

 *  2 - problem has no dual feasible solution.

 *

 *  PROBLEM TRANSFORMATION

 *

 *  Constraint (2) can be written as follows:

 *

 *     L[p] - sum a[p,j] x[j] <= a[p,q] x[q] <= U[p] - sum a[p,j] x[j],

 *            j!=q                                     j!=q

 *

 *  from which it follows that:

 *

 *     alfa <= a[p,q] x[q] <= beta,                                   (3)

 *

 *  where

 *

 *     alfa = inf(L[p] - sum a[p,j] x[j]) =

 *                       j!=q

 *

 *          = L[p] - sup sum a[p,j] x[j] =                            (4)

 *                       j!=q

 *

 *          = L[p] -  sum  a[p,j] u[j] -  sum  a[p,j] l[j],

 *                  j in Jp             j in Jn

 *

 *     beta = sup(L[p] - sum a[p,j] x[j]) =

 *                       j!=q

 *

 *          = L[p] - inf sum a[p,j] x[j] =                            (5)

 *                       j!=q

 *

 *          = L[p] -  sum  a[p,j] l[j] -  sum  a[p,j] u[j],

 *                  j in Jp             j in Jn

 *

 *     Jp = {j != q: a[p,j] > 0},  Jn = {j != q: a[p,j] < 0}.         (6)

 *

 *  Inequality (3) defines implied bounds of variable x[q]:

 *

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

 *

 *  where

 *

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

 *     l'[q] = <                                                     (8a)

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

 *

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

 *     u'[q] = <                                                     (8b)

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

 *

 *  Thus, if l'[q] > l[q] - eps and u'[q] < u[q] + eps, where eps is

 *  an absolute tolerance for column value, column bounds (1) cannot be

 *  active, in which case column q can be replaced by equivalent free

 *  (unbounded) column.

 *

 *  Note that column q is column singleton, so in the dual system of the

 *  original problem it corresponds to the following row singleton:

 *

 *     a[p,q] pi[p] + lambda[q] = c[q],                               (9)

 *

 *  from which it follows that:

 *

 *     pi[p] = (c[q] - lambda[q]) / a[p,q].                          (10)

 *

 *  Let x[q] be implied free (unbounded) variable. Then column q can be

 *  only basic, so its multiplier lambda[q] is equal to zero, and from

 *  (10) we have:

 *

 *     pi[p] = c[q] / a[p,q].                                        (11)

 *

 *  There are possible three cases:

 *

 *  1) pi[p] < -eps, where eps is an absolute tolerance for row

 *     multiplier. In this case, to provide dual feasibility of the

 *     original problem, row p must be active on its lower bound, and

 *     if its lower bound does not exist (L[p] = -oo), the problem has

 *     no dual feasible solution;

 *

 *  2) pi[p] > +eps. In this case row p must be active on its upper

 *     bound, and if its upper bound does not exist (U[p] = +oo), the

 *     problem has no dual feasible solution;

 *

 *  3) -eps <= pi[p] <= +eps. In this case any (either lower or upper)

 *     bound of row p can be active, because this does not affect dual

 *     feasibility.

 *

 *  Thus, in all three cases original inequality constraint (2) can be

 *  replaced by equality constraint, where the right-hand side is either

 *  lower or upper bound of row p, and bounds of column q can be removed

 *  that makes it free (unbounded). (May note that this transformation

 *  can be followed by transformation "Column singleton (implied slack

 *  variable)" performed by the routine npp_implied_slack.)

 *

 *  RECOVERING BASIC SOLUTION

 *

 *  Status of row p in solution to the original problem is determined

 *  by its status in solution to the transformed problem and its bound,

 *  which was choosen to be active:

 *

 *     +-----------------------+--------+--------------------+

 *     |    Status of row p    | Active | Status of row p    |

 *     | (transformed problem) | bound  | (original problem) |

 *     +-----------------------+--------+--------------------+

 *     |        GLP_BS         |  L[p]  |       GLP_BS       |

 *     |        GLP_BS         |  U[p]  |       GLP_BS       |

 *     |        GLP_NS         |  L[p]  |       GLP_NL       |

 *     |        GLP_NS         |  U[p]  |       GLP_NU       |

 *     +-----------------------+--------+--------------------+

 *

 *  Value of row multiplier pi[p] (as well as value of column q) in

 *  solution to the original problem is the same as in solution to the

 *  transformed problem.

 *

 *  RECOVERING INTERIOR-POINT SOLUTION

 *

 *  Value of row multiplier pi[p] in solution to the original problem is

 *  the same as in solution to the transformed problem.

 *

 *  RECOVERING MIP SOLUTION

 *

 *  None needed. */

 struct implied_free

 {     /* column singleton (implied free variable) */

       int p;

       /* row reference number */

       char stat;

       /* row status:

          GLP_NL - active constraint on lower bound

          GLP_NU - active constraint on upper bound */

 };

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

 int npp_implied_free(NPP *npp, NPPCOL *q)

 {     /* process column singleton (implied free variable) */

       struct implied_free *info;

       NPPROW *p;

       NPPAIJ *apq, *aij;

       double alfa, beta, l, u, pi, eps;

       /* the column must be non-fixed singleton */

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

       xassert(q->ptr != NULL && q->ptr->c_next == NULL);

       /* corresponding row must be inequality constraint */

       apq = q->ptr;

       p = apq->row;

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

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

       /* compute alfa */

       alfa = p->lb;

       if (alfa != -DBL_MAX)

       {  for (aij = p->ptr; aij != NULL; aij = aij->r_next)

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

             if (aij->val > 0.0)

             {  if (aij->col->ub == +DBL_MAX)

                {  alfa = -DBL_MAX;

                   break;

                }

                alfa -= aij->val * aij->col->ub;

             }

             else /* < 0.0 */

             {  if (aij->col->lb == -DBL_MAX)

                {  alfa = -DBL_MAX;

                   break;

                }

                alfa -= aij->val * aij->col->lb;

             }

          }

       }

       /* compute beta */

       beta = p->ub;

       if (beta != +DBL_MAX)

       {  for (aij = p->ptr; aij != NULL; aij = aij->r_next)

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

             if (aij->val > 0.0)

             {  if (aij->col->lb == -DBL_MAX)

                {  beta = +DBL_MAX;

                   break;

                }

                beta -= aij->val * aij->col->lb;

             }

             else /* < 0.0 */

             {  if (aij->col->ub == +DBL_MAX)

                {  beta = +DBL_MAX;

                   break;

                }

                beta -= aij->val * aij->col->ub;

             }

          }

       }

       /* compute implied column lower bound l'[q] */

       if (apq->val > 0.0)

          l = (alfa == -DBL_MAX ? -DBL_MAX : alfa / apq->val);

       else /* < 0.0 */

          l = (beta == +DBL_MAX ? -DBL_MAX : beta / apq->val);

       /* compute implied column upper bound u'[q] */

       if (apq->val > 0.0)

          u = (beta == +DBL_MAX ? +DBL_MAX : beta / apq->val);

       else

          u = (alfa == -DBL_MAX ? +DBL_MAX : alfa / apq->val);

       /* check if column lower bound l[q] can be active */

       if (q->lb != -DBL_MAX)

       {  eps = 1e-9 + 1e-12 * fabs(q->lb);

          if (l < q->lb - eps) return 1; /* yes, it can */

       }

       /* check if column upper bound u[q] can be active */

       if (q->ub != +DBL_MAX)

       {  eps = 1e-9 + 1e-12 * fabs(q->ub);

          if (u > q->ub + eps) return 1; /* yes, it can */

       }

       /* okay; make column q free (unbounded) */

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

       /* create transformation stack entry */

       info = npp_push_tse(npp,

          rcv_implied_free, sizeof(struct implied_free));

       info->p = p->i;

       info->stat = -1;

       /* compute row multiplier pi[p] */

       pi = q->coef / apq->val;

       /* check dual feasibility for row p */

       if (pi > +DBL_EPSILON)

       {  /* lower bound L[p] must be active */

          if (p->lb != -DBL_MAX)

 nl:      {  info->stat = GLP_NL;

             p->ub = p->lb;

          }

          else

          {  if (pi > +1e-5) return 2; /* dual infeasibility */

             /* take a chance on U[p] */

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

             goto nu;

          }

       }

       else if (pi < -DBL_EPSILON)

       {  /* upper bound U[p] must be active */

          if (p->ub != +DBL_MAX)

 nu:      {  info->stat = GLP_NU;

             p->lb = p->ub;

          }

          else

          {  if (pi < -1e-5) return 2; /* dual infeasibility */

             /* take a chance on L[p] */

             xassert(p->lb != -DBL_MAX);

             goto nl;

          }

       }

       else

       {  /* any bound (either L[p] or U[p]) can be made active  */

          if (p->ub == +DBL_MAX)

          {  xassert(p->lb != -DBL_MAX);

             goto nl;

          }

          if (p->lb == -DBL_MAX)

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

             goto nu;

          }

          if (fabs(p->lb) <= fabs(p->ub)) goto nl; else goto nu;

       }

       return 0;

 }

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

 {     /* recover column singleton (implied free variable) */

       struct implied_free *info = _info;

       if (npp->sol == GLP_SOL)

       {  if (npp->r_stat[info->p] == GLP_BS)

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

          else if (npp->r_stat[info->p] == GLP_NS)

          {  xassert(info->stat == GLP_NL || info->stat == GLP_NU);

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

          }

          else

          {  npp_error();

             return 1;

          }

       }

       return 0;

 }

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

 *  NAME

 *

 *  npp_eq_doublet - process row doubleton (equality constraint)

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  NPPCOL *npp_eq_doublet(NPP *npp, NPPROW *p);

 *

 *  DESCRIPTION

 *

 *  The routine npp_eq_doublet processes row p, which is equality

 *  constraint having exactly two non-zero coefficients:

 *

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

 *

 *  As the result of processing one of columns q or r is eliminated from

 *  all other rows and, thus, becomes column singleton of type "implied

 *  slack variable". Row p is not changed and along with column q and r

 *  remains in the problem.

 *

 *  RETURNS

 *

 *  The routine npp_eq_doublet returns pointer to the descriptor of that

 *  column q or r which has been eliminated. If, due to some reason, the

 *  elimination was not performed, the routine returns NULL.

 *

 *  PROBLEM TRANSFORMATION

 *

 *  First, we decide which column q or r will be eliminated. Let it be

 *  column q. Consider i-th constraint row, where column q has non-zero

 *  coefficient a[i,q] != 0:

 *

 *     L[i] <= sum a[i,j] x[j] <= U[i].                               (2)

 *              j

 *

 *  In order to eliminate column q from row (2) we subtract from it row

 *  (1) multiplied by gamma[i] = a[i,q] / a[p,q], i.e. we replace in the

 *  transformed problem row (2) by its linear combination with row (1).

 *  This transformation changes only coefficients in columns q and r,

 *  and bounds of row i as follows:

 *

 *     a~[i,q] = a[i,q] - gamma[i] a[p,q] = 0,                        (3)

 *

 *     a~[i,r] = a[i,r] - gamma[i] a[p,r],                            (4)

 *

 *       L~[i] = L[i] - gamma[i] b,                                   (5)

 *

 *       U~[i] = U[i] - gamma[i] b.                                   (6)

 *

 *  RECOVERING BASIC SOLUTION

 *

 *  The transformation of the primal system of the original problem:

 *

 *     L <= A x <= U                                                  (7)

 *

 *  is equivalent to multiplying from the left a transformation matrix F

 *  by components of this primal system, which in the transformed problem

 *  becomes the following:

 *

 *     F L <= F A x <= F U  ==>  L~ <= A~x <= U~.                     (8)

 *

 *  The matrix F has the following structure:

 *

 *         ( 1           -gamma[1]            )

 *         (                                  )

 *         (    1        -gamma[2]            )

 *         (                                  )

 *         (      ...       ...               )

 *         (                                  )

 *     F = (          1  -gamma[p-1]          )                       (9)

 *         (                                  )

 *         (                 1                )

 *         (                                  )

 *         (             -gamma[p+1]  1       )

 *         (                                  )

 *         (                ...          ...  )

 *

 *  where its column containing elements -gamma[i] corresponds to row p

 *  of the primal system.

 *

 *  From (8) it follows that the dual system of the original problem:

 *

 *     A'pi + lambda = c,                                            (10)

 *

 *  in the transformed problem becomes the following:

 *

 *     A'F'inv(F')pi + lambda = c  ==>  (A~)'pi~ + lambda = c,       (11)

 *

 *  where:

 *

 *     pi~ = inv(F')pi                                               (12)

 *

 *  is the vector of row multipliers in the transformed problem. Thus:

 *

 *     pi = F'pi~.                                                   (13)

 *

 *  Therefore, as it follows from (13), value of multiplier for row p in

 *  solution to the original problem can be computed as follows:

 *

 *     pi[p] = pi~[p] - sum gamma[i] pi~[i],                         (14)

 *                       i

 *

 *  where pi~[i] = pi[i] is multiplier for row i (i != p).

 *

 *  Note that the statuses of all rows and columns are not changed.

 *

 *  RECOVERING INTERIOR-POINT SOLUTION

 *

 *  Multiplier for row p in solution to the original problem is computed

 *  with formula (14).

 *

 *  RECOVERING MIP SOLUTION

 *

 *  None needed. */

 struct eq_doublet

 {     /* row doubleton (equality constraint) */

       int p;

       /* row reference number */

       double apq;

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

       NPPLFE *ptr;

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

 };

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

 NPPCOL *npp_eq_doublet(NPP *npp, NPPROW *p)

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

       struct eq_doublet *info;

       NPPROW *i;

       NPPCOL *q, *r;

       NPPAIJ *apq, *apr, *aiq, *air, *next;

       NPPLFE *lfe;

       double gamma;

       /* the row must be doubleton equality constraint */

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

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

               p->ptr->r_next->r_next == NULL);

       /* choose column to be eliminated */

       {  NPPAIJ *a1, *a2;

          a1 = p->ptr, a2 = a1->r_next;

          if (fabs(a2->val) < 0.001 * fabs(a1->val))

          {  /* only first column can be eliminated, because second one

                has too small constraint coefficient */

             apq = a1, apr = a2;

          }

          else if (fabs(a1->val) < 0.001 * fabs(a2->val))

          {  /* only second column can be eliminated, because first one

                has too small constraint coefficient */

             apq = a2, apr = a1;

          }

          else

          {  /* both columns are appropriate; choose that one which is

                shorter to minimize fill-in */

             if (npp_col_nnz(npp, a1->col) <= npp_col_nnz(npp, a2->col))

             {  /* first column is shorter */

                apq = a1, apr = a2;

             }

             else

             {  /* second column is shorter */

                apq = a2, apr = a1;

             }

          }

       }

       /* now columns q and r have been chosen */

       q = apq->col, r = apr->col;

       /* create transformation stack entry */

       info = npp_push_tse(npp,

          rcv_eq_doublet, sizeof(struct eq_doublet));

       info->p = p->i;

       info->apq = apq->val;

       info->ptr = NULL;

       /* transform each row i (i != p), where a[i,q] != 0, to eliminate

          column q */

       for (aiq = q->ptr; aiq != NULL; aiq = next)

       {  next = aiq->c_next;

          if (aiq == apq) continue; /* skip row p */

          i = aiq->row; /* row i to be transformed */

          /* save constraint coefficient a[i,q] */

          if (npp->sol != GLP_MIP)

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

             lfe->ref = i->i;

             lfe->val = aiq->val;

             lfe->next = info->ptr;

             info->ptr = lfe;

          }

          /* find coefficient a[i,r] in row i */

          for (air = i->ptr; air != NULL; air = air->r_next)

             if (air->col == r) break;

          /* if a[i,r] does not exist, create a[i,r] = 0 */

          if (air == NULL)

             air = npp_add_aij(npp, i, r, 0.0);

          /* compute gamma[i] = a[i,q] / a[p,q] */

          gamma = aiq->val / apq->val;

          /* (row i) := (row i) - gamma[i] * (row p); see (3)-(6) */

          /* new a[i,q] is exact zero due to elimnation; remove it from

             row i */

          npp_del_aij(npp, aiq);

          /* compute new a[i,r] */

          air->val -= gamma * apr->val;

          /* if new a[i,r] is close to zero due to numeric cancelation,

             remove it from row i */

          if (fabs(air->val) <= 1e-10)

             npp_del_aij(npp, air);

          /* compute new lower and upper bounds of row i */

          if (i->lb == i->ub)

             i->lb = i->ub = (i->lb - gamma * p->lb);

          else

          {  if (i->lb != -DBL_MAX)

                i->lb -= gamma * p->lb;

             if (i->ub != +DBL_MAX)

                i->ub -= gamma * p->lb;

          }

       }

       return q;

 }

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

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

       struct eq_doublet *info = _info;

       NPPLFE *lfe;

       double gamma, temp;

       /* we assume that processing row p is followed by processing

          column q as singleton of type "implied slack variable", in

          which case row p must always be active equality constraint */

       if (npp->sol == GLP_SOL)

       {  if (npp->r_stat[info->p] != GLP_NS)

          {  npp_error();

             return 1;

          }

       }

       if (npp->sol != GLP_MIP)

       {  /* compute value of multiplier for row p; see (14) */

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

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

          {  gamma = lfe->val / info->apq; /* a[i,q] / a[p,q] */

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

          }

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

       }

       return 0;

 }

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

 *  NAME

 *

 *  npp_forcing_row - process forcing row

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  int npp_forcing_row(NPP *npp, NPPROW *p, int at);

 *

 *  DESCRIPTION

 *

 *  The routine npp_forcing row processes row p of general format:

 *

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

 *              j

 *

 *     l[j] <= x[j] <= u[j],                                          (2)

 *

 *  where L[p] <= U[p] and l[j] < u[j] for all a[p,j] != 0. It is also

 *  assumed that:

 *

 *  1) if at = 0 then |L[p] - U'[p]| <= eps, where U'[p] is implied

 *     row upper bound (see below), eps is an absolute tolerance for row

 *     value;

 *

 *  2) if at = 1 then |U[p] - L'[p]| <= eps, where L'[p] is implied

 *     row lower bound (see below).

 *

 *  RETURNS

 *

 *  0 - success;

 *

 *  1 - cannot fix columns due to too small constraint coefficients.

 *

 *  PROBLEM TRANSFORMATION

 *

 *  Implied lower and upper bounds of row (1) are determined by bounds

 *  of corresponding columns (variables) as follows:

 *

 *     L'[p] = inf sum a[p,j] x[j] =

 *                  j

 *                                                                    (3)

 *           =  sum  a[p,j] l[j] +  sum  a[p,j] u[j],

 *            j in Jp             j in Jn

 *

 *     U'[p] = sup sum a[p,j] x[j] =

 *                                                                    (4)

 *           =  sum  a[p,j] u[j] +  sum  a[p,j] l[j],

 *            j in Jp             j in Jn

 *

 *     Jp = {j: a[p,j] > 0},  Jn = {j: a[p,j] < 0}.                   (5)

 *

 *  If L[p] =~ U'[p] (at = 0), solution can be primal feasible only when

 *  all variables take their boundary values as defined by (4):

 *

 *            ( u[j], if j in Jp

 *     x[j] = <                                                       (6)

 *            ( l[j], if j in Jn

 *

 *  Similarly, if U[p] =~ L'[p] (at = 1), solution can be primal feasible

 *  only when all variables take their boundary values as defined by (3):

 *

 *            ( l[j], if j in Jp

 *     x[j] = <                                                       (7)

 *            ( u[j], if j in Jn

 *

 *  Condition (6) or (7) allows fixing all columns (variables x[j])

 *  in row (1) on their bounds and then removing them from the problem

 *  (see the routine npp_fixed_col). Due to this row p becomes redundant,

 *  so it can be replaced by equivalent free (unbounded) row and also

 *  removed from the problem (see the routine npp_free_row).

 *

 *  1. To apply this transformation row (1) should not have coefficients

 *     whose magnitude is too small, i.e. all a[p,j] should satisfy to

 *     the following condition:

 *

 *        |a[p,j]| >= eps * max(1, |a[p,k]|),                         (8)

 *                           k

 *     where eps is a relative tolerance for constraint coefficients.

 *     Otherwise, fixing columns may be numerically unreliable and may

 *     lead to wrong solution.

 *

 *  2. The routine fixes columns and remove bounds of row p, however,

 *     it does not remove the row and columns from the problem.

 *

 *  RECOVERING BASIC SOLUTION

 *

 *  In the transformed problem row p being inactive constraint is

 *  assigned status GLP_BS (as the result of transformation of free

 *  row), and all columns in this row are assigned status GLP_NS (as the

 *  result of transformation of fixed columns).

 *

 *  Note that in the dual system of the transformed (as well as original)

 *  problem every column j in row p corresponds to the following row:

 *

 *     sum  a[i,j] pi[i] + a[p,j] pi[p] + lambda[j] = c[j],           (9)

 *     i!=p

 *

 *  from which it follows that:

 *

 *     lambda[j] = c[j] - sum a[i,j] pi[i] - a[p,j] pi[p].           (10)

 *                        i!=p

 *

 *  In the transformed problem values of all multipliers pi[i] are known

 *  (including pi[i], whose value is zero, since row p is inactive).

 *  Thus, using formula (10) it is possible to compute values of

 *  multipliers lambda[j] for all columns in row p.

 *

 *  Note also that in the original problem all columns in row p are

 *  bounded, not fixed. So status GLP_NS assigned to every such column

 *  must be changed to GLP_NL or GLP_NU depending on which bound the

 *  corresponding column has been fixed. This status change may lead to

 *  dual feasibility violation for solution of the original problem,

 *  because now column multipliers must satisfy to the following

 *  condition:

 *

 *               ( >= 0, if status of column j is GLP_NL,

 *     lambda[j] <                                                   (11)

 *               ( <= 0, if status of column j is GLP_NU.

 *

 *  If this condition holds, solution to the original problem is the

 *  same as to the transformed problem. Otherwise, we have to perform

 *  one degenerate pivoting step of the primal simplex method to obtain

 *  dual feasible (hence, optimal) solution to the original problem as

 *  follows. If, on problem transformation, row p was made active on its

 *  lower bound (case at = 0), we change its status to GLP_NL (or GLP_NS)

 *  and start increasing its multiplier pi[p]. Otherwise, if row p was

 *  made active on its upper bound (case at = 1), we change its status

 *  to GLP_NU (or GLP_NS) and start decreasing pi[p]. From (10) it

 *  follows that:

 *

 *     delta lambda[j] = - a[p,j] * delta pi[p] = - a[p,j] pi[p].    (12)

 *

 *  Simple analysis of formulae (3)-(5) shows that changing pi[p] in the

 *  specified direction causes increasing lambda[j] for every column j

 *  assigned status GLP_NL (delta lambda[j] > 0) and decreasing lambda[j]

 *  for every column j assigned status GLP_NU (delta lambda[j] < 0). It

 *  is understood that once the last lambda[q], which violates condition

 *  (11), has reached zero, multipliers lambda[j] for all columns get

 *  valid signs. Such column q can be determined as follows. Let d[j] be

 *  initial value of lambda[j] (i.e. reduced cost of column j) in the

 *  transformed problem computed with formula (10) when pi[p] = 0. Then

 *  lambda[j] = d[j] + delta lambda[j], and from (12) it follows that

 *  lambda[j] becomes zero if:

 *

 *     delta lambda[j] = - a[p,j] pi[p] = - d[j]  ==>

 *                                                                   (13)

 *     pi[p] = d[j] / a[p,j].

 *

 *  Therefore, the last column q, for which lambda[q] becomes zero, can

 *  be determined from the following condition:

 *

 *     |d[q] / a[p,q]| = max  |pi[p]| = max  |d[j] / a[p,j]|,        (14)

 *                      j in D         j in D

 *

 *  where D is a set of columns j whose, reduced costs d[j] have invalid

 *  signs, i.e. violate condition (11). (Thus, if D is empty, solution

 *  to the original problem is the same as solution to the transformed

 *  problem, and no correction is needed as was noticed above.) In

 *  solution to the original problem column q is assigned status GLP_BS,

 *  since it replaces column of auxiliary variable of row p (becoming

 *  active) in the basis, and multiplier for row p is assigned its new

 *  value, which is pi[p] = d[q] / a[p,q]. Note that due to primal

 *  degeneracy values of all columns having non-zero coefficients in row

 *  p remain unchanged.

 *

 *  RECOVERING INTERIOR-POINT SOLUTION

 *

 *  Value of multiplier pi[p] in solution to the original problem is

 *  corrected in the same way as for basic solution. Values of all

 *  columns having non-zero coefficients in row p remain unchanged.

 *

 *  RECOVERING MIP SOLUTION

 *

 *  None needed. */

 struct forcing_col

 {     /* column fixed on its bound by forcing row */

       int j;

       /* column reference number */

       char stat;

       /* original column status:

          GLP_NL - fixed on lower bound

          GLP_NU - fixed on upper bound */

       double a;

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

       double c;

       /* objective coefficient c[j] */

       NPPLFE *ptr;

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

       struct forcing_col *next;

       /* pointer to another column fixed by forcing row */

 };

 struct forcing_row

 {     /* forcing row */

       int p;

       /* row reference number */

       char stat;

       /* status assigned to the row if it becomes active:

          GLP_NS - active equality constraint

          GLP_NL - inequality constraint with lower bound active

          GLP_NU - inequality constraint with upper bound active */

       struct forcing_col *ptr;

       /* list of all columns having non-zero constraint coefficient

          a[p,j] in the forcing row */

 };

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

 int npp_forcing_row(NPP *npp, NPPROW *p, int at)

 {     /* process forcing row */

       struct forcing_row *info;

       struct forcing_col *col = NULL;

       NPPCOL *j;

       NPPAIJ *apj, *aij;

       NPPLFE *lfe;

       double big;

       xassert(at == 0 || at == 1);

       /* determine maximal magnitude of the row coefficients */

       big = 1.0;

       for (apj = p->ptr; apj != NULL; apj = apj->r_next)

          if (big < fabs(apj->val)) big = fabs(apj->val);

       /* if there are too small coefficients in the row, transformation

          should not be applied */

       for (apj = p->ptr; apj != NULL; apj = apj->r_next)

          if (fabs(apj->val) < 1e-7 * big) return 1;

       /* create transformation stack entry */

       info = npp_push_tse(npp,

          rcv_forcing_row, sizeof(struct forcing_row));

       info->p = p->i;

       if (p->lb == p->ub)

       {  /* equality constraint */

          info->stat = GLP_NS;

       }

       else if (at == 0)

       {  /* inequality constraint; case L[p] = U'[p] */

          info->stat = GLP_NL;

          xassert(p->lb != -DBL_MAX);

       }

       else /* at == 1 */

       {  /* inequality constraint; case U[p] = L'[p] */

          info->stat = GLP_NU;

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

       }

       info->ptr = NULL;

       /* scan the forcing row, fix columns at corresponding bounds, and

          save column information (the latter is not needed for MIP) */

       for (apj = p->ptr; apj != NULL; apj = apj->r_next)

       {  /* column j has non-zero coefficient in the forcing row */

          j = apj->col;

          /* it must be non-fixed */

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

          /* allocate stack entry to save column information */

          if (npp->sol != GLP_MIP)

          {  col = dmp_get_atom(npp->stack, sizeof(struct forcing_col));

             col->j = j->j;

             col->stat = -1; /* will be set below */

             col->a = apj->val;

             col->c = j->coef;

             col->ptr = NULL;

             col->next = info->ptr;

             info->ptr = col;

          }

          /* fix column j */

          if (at == 0 && apj->val < 0.0 || at != 0 && apj->val > 0.0)

          {  /* at its lower bound */

             if (npp->sol != GLP_MIP)

                col->stat = GLP_NL;

             xassert(j->lb != -DBL_MAX);

             j->ub = j->lb;

          }

          else

          {  /* at its upper bound */

             if (npp->sol != GLP_MIP)

                col->stat = GLP_NU;

             xassert(j->ub != +DBL_MAX);

             j->lb = j->ub;

          }

          /* save column coefficients a[i,j], i != p */

          if (npp->sol != GLP_MIP)

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

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

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

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

                lfe->val = aij->val;

                lfe->next = col->ptr;

                col->ptr = lfe;

             }

          }

       }

       /* make the row free (unbounded) */

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

       return 0;

 }

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

 {     /* recover forcing row */

       struct forcing_row *info = _info;

       struct forcing_col *col, *piv;

       NPPLFE *lfe;

       double d, big, temp;

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

       /* initially solution to the original problem is the same as

          to the transformed problem, where row p is inactive constraint

          with pi[p] = 0, and all columns are non-basic */

       if (npp->sol == GLP_SOL)

       {  if (npp->r_stat[info->p] != GLP_BS)

          {  npp_error();

             return 1;

          }

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

          {  if (npp->c_stat[col->j] != GLP_NS)

             {  npp_error();

                return 1;

             }

             npp->c_stat[col->j] = col->stat; /* original status */

          }

       }

       /* compute reduced costs d[j] for all columns with formula (10)

          and store them in col.c instead objective coefficients */

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

       {  d = col->c;

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

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

          col->c = d;

       }

       /* consider columns j, whose multipliers lambda[j] has wrong

          sign in solution to the transformed problem (where lambda[j] =

          d[j]), and choose column q, whose multipler lambda[q] reaches

          zero last on changing row multiplier pi[p]; see (14) */

       piv = NULL, big = 0.0;

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

       {  d = col->c; /* d[j] */

          temp = fabs(d / col->a);

          if (col->stat == GLP_NL)

          {  /* column j has active lower bound */

             if (d < 0.0 && big < temp)

                piv = col, big = temp;

          }

          else if (col->stat == GLP_NU)

          {  /* column j has active upper bound */

             if (d > 0.0 && big < temp)

                piv = col, big = temp;

          }

          else

          {  npp_error();

             return 1;

          }

       }

       /* if column q does not exist, no correction is needed */

       if (piv != NULL)

       {  /* correct solution; row p becomes active constraint while

             column q becomes basic */

          if (npp->sol == GLP_SOL)

          {  npp->r_stat[info->p] = info->stat;

             npp->c_stat[piv->j] = GLP_BS;

          }

          /* assign new value to row multiplier pi[p] = d[p] / a[p,q] */

          npp->r_pi[info->p] = piv->c / piv->a;

       }

 done: return 0;

 }

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

 *  NAME

 *

 *  npp_analyze_row - perform general row analysis

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

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

 *

 *  DESCRIPTION

 *

 *  The routine npp_analyze_row performs analysis of row p of general

 *  format:

 *

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

 *              j

 *

 *     l[j] <= x[j] <= u[j],                                          (2)

 *

 *  where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0.

 *

 *  RETURNS

 *

 *  0x?0 - row lower bound does not exist or is redundant;

 *

 *  0x?1 - row lower bound can be active;

 *

 *  0x?2 - row lower bound is a forcing bound;

 *

 *  0x0? - row upper bound does not exist or is redundant;

 *

 *  0x1? - row upper bound can be active;

 *

 *  0x2? - row upper bound is a forcing bound;

 *

 *  0x33 - row bounds are inconsistent with column bounds.

 *

 *  ALGORITHM

 *

 *  Analysis of row (1) is based on analysis of its implied lower and

 *  upper bounds, which are determined by bounds of corresponding columns

 *  (variables) as follows:

 *

 *     L'[p] = inf sum a[p,j] x[j] =

 *                  j

 *                                                                    (3)

 *           =  sum  a[p,j] l[j] +  sum  a[p,j] u[j],

 *            j in Jp             j in Jn

 *

 *     U'[p] = sup sum a[p,j] x[j] =

 *                                                                    (4)

 *           =  sum  a[p,j] u[j] +  sum  a[p,j] l[j],

 *            j in Jp             j in Jn

 *

 *     Jp = {j: a[p,j] > 0},  Jn = {j: a[p,j] < 0}.                   (5)

 *

 *  (Note that bounds of all columns in row p are assumed to be correct,

 *  so L'[p] <= U'[p].)

 *

 *  Analysis of row lower bound L[p] includes the following cases:

 *

 *  1) if L[p] > U'[p] + eps, where eps is an absolute tolerance for row

 *     value, row lower bound L[p] and implied row upper bound U'[p] are

 *     inconsistent, ergo, the problem has no primal feasible solution;

 *

 *  2) if U'[p] - eps <= L[p] <= U'[p] + eps, i.e. if L[p] =~ U'[p],

 *     the row is a forcing row on its lower bound (see description of

 *     the routine npp_forcing_row);

 *

 *  3) if L[p] > L'[p] + eps, row lower bound L[p] can be active (this

 *     conclusion does not account other rows in the problem);

 *

 *  4) if L[p] <= L'[p] + eps, row lower bound L[p] cannot be active, so

 *     it is redundant and can be removed (replaced by -oo).

 *

 *  Analysis of row upper bound U[p] is performed in a similar way and

 *  includes the following cases:

 *

 *  1) if U[p] < L'[p] - eps, row upper bound U[p] and implied row lower

 *     bound L'[p] are inconsistent, ergo the problem has no primal

 *     feasible solution;

 *

 *  2) if L'[p] - eps <= U[p] <= L'[p] + eps, i.e. if U[p] =~ L'[p],

 *     the row is a forcing row on its upper bound (see description of

 *     the routine npp_forcing_row);

 *

 *  3) if U[p] < U'[p] - eps, row upper bound U[p] can be active (this

 *     conclusion does not account other rows in the problem);

 *

 *  4) if U[p] >= U'[p] - eps, row upper bound U[p] cannot be active, so

 *     it is redundant and can be removed (replaced by +oo). */

 int npp_analyze_row(NPP *npp, NPPROW *p)

 {     /* perform general row analysis */

       NPPAIJ *aij;

       int ret = 0x00;

       double l, u, eps;

       xassert(npp == npp);

       /* compute implied lower bound L'[p]; see (3) */

       l = 0.0;

       for (aij = p->ptr; aij != NULL; aij = aij->r_next)

       {  if (aij->val > 0.0)

          {  if (aij->col->lb == -DBL_MAX)

             {  l = -DBL_MAX;

                break;

             }

             l += aij->val * aij->col->lb;

          }

          else /* aij->val < 0.0 */

          {  if (aij->col->ub == +DBL_MAX)

             {  l = -DBL_MAX;

                break;

             }

             l += aij->val * aij->col->ub;

          }

       }

       /* compute implied upper bound U'[p]; see (4) */

       u = 0.0;

       for (aij = p->ptr; aij != NULL; aij = aij->r_next)

       {  if (aij->val > 0.0)

          {  if (aij->col->ub == +DBL_MAX)

             {  u = +DBL_MAX;

                break;

             }

             u += aij->val * aij->col->ub;

          }

          else /* aij->val < 0.0 */

          {  if (aij->col->lb == -DBL_MAX)

             {  u = +DBL_MAX;

                break;

             }

             u += aij->val * aij->col->lb;

          }

       }

       /* column bounds are assumed correct, so L'[p] <= U'[p] */

       /* check if row lower bound is consistent */

       if (p->lb != -DBL_MAX)

       {  eps = 1e-3 + 1e-6 * fabs(p->lb);

          if (p->lb - eps > u)

          {  ret = 0x33;

             goto done;

          }

       }

       /* check if row upper bound is consistent */

       if (p->ub != +DBL_MAX)

       {  eps = 1e-3 + 1e-6 * fabs(p->ub);

          if (p->ub + eps < l)

          {  ret = 0x33;

             goto done;

          }

       }

       /* check if row lower bound can be active/forcing */

       if (p->lb != -DBL_MAX)

       {  eps = 1e-9 + 1e-12 * fabs(p->lb);

          if (p->lb - eps > l)

          {  if (p->lb + eps <= u)

                ret |= 0x01;

             else

                ret |= 0x02;

          }

       }

       /* check if row upper bound can be active/forcing */

       if (p->ub != +DBL_MAX)

       {  eps = 1e-9 + 1e-12 * fabs(p->ub);

          if (p->ub + eps < u)

          {  /* check if the upper bound is forcing */

             if (p->ub - eps >= l)

                ret |= 0x10;

             else

                ret |= 0x20;

          }

       }

 done: return ret;

 }

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

 *  NAME

 *

 *  npp_inactive_bound - remove row lower/upper inactive bound

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  void npp_inactive_bound(NPP *npp, NPPROW *p, int which);

 *

 *  DESCRIPTION

 *

 *  The routine npp_inactive_bound removes lower (if which = 0) or upper

 *  (if which = 1) bound of row p:

 *

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

 *

 *  which (bound) is assumed to be redundant.

 *

 *  PROBLEM TRANSFORMATION

 *

 *  If which = 0, current lower bound L[p] of row p is assigned -oo.

 *  If which = 1, current upper bound U[p] of row p is assigned +oo.

 *

 *  RECOVERING BASIC SOLUTION

 *

 *  If in solution to the transformed problem row p is inactive

 *  constraint (GLP_BS), its status is not changed in solution to the

 *  original problem. Otherwise, status of row p in solution to the

 *  original problem is defined by its type before transformation and

 *  its status in solution to the transformed problem as follows:

 *

 *     +---------------------+-------+---------------+---------------+

 *     |        Row          | Flag  | Row status in | Row status in |

 *     |        type         | which | transfmd soln | original soln |

 *     +---------------------+-------+---------------+---------------+

 *     |     sum >= L[p]     |   0   |    GLP_NF     |    GLP_NL     |

 *     |     sum <= U[p]     |   1   |    GLP_NF     |    GLP_NU     |

 *     | L[p] <= sum <= U[p] |   0   |    GLP_NU     |    GLP_NU     |

 *     | L[p] <= sum <= U[p] |   1   |    GLP_NL     |    GLP_NL     |

 *     |  sum = L[p] = U[p]  |   0   |    GLP_NU     |    GLP_NS     |

 *     |  sum = L[p] = U[p]  |   1   |    GLP_NL     |    GLP_NS     |

 *     +---------------------+-------+---------------+---------------+

 *

 *  RECOVERING INTERIOR-POINT SOLUTION

 *

 *  None needed.

 *

 *  RECOVERING MIP SOLUTION

 *

 *  None needed. */

 struct inactive_bound

 {     /* row inactive bound */

       int p;

       /* row reference number */

       char stat;

       /* row status (if active constraint) */

 };

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

 void npp_inactive_bound(NPP *npp, NPPROW *p, int which)

 {     /* remove row lower/upper inactive bound */

       struct inactive_bound *info;

       if (npp->sol == GLP_SOL)

       {  /* create transformation stack entry */

          info = npp_push_tse(npp,

             rcv_inactive_bound, sizeof(struct inactive_bound));

          info->p = p->i;

          if (p->ub == +DBL_MAX)

             info->stat = GLP_NL;

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

             info->stat = GLP_NU;

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

             info->stat = (char)(which == 0 ? GLP_NU : GLP_NL);

          else

             info->stat = GLP_NS;

       }

       /* remove row inactive bound */

       if (which == 0)

       {  xassert(p->lb != -DBL_MAX);

          p->lb = -DBL_MAX;

       }

       else if (which == 1)

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

          p->ub = +DBL_MAX;

       }

       else

          xassert(which != which);

       return;

 }

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

 {     /* recover row status */

       struct inactive_bound *info = _info;

       if (npp->sol != GLP_SOL)

       {  npp_error();

          return 1;

       }

       if (npp->r_stat[info->p] == GLP_BS)

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

       else

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

       return 0;

 }

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

 *  NAME

 *

 *  npp_implied_bounds - determine implied column bounds

 *

 *  SYNOPSIS

 *

 *  #include "glpnpp.h"

 *  void npp_implied_bounds(NPP *npp, NPPROW *p);

 *

 *  DESCRIPTION

 *

 *  The routine npp_implied_bounds inspects general row (constraint) p:

 *

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

 *

 *     l[j] <= x[j] <= u[j],                                          (2)

 *

 *  where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0, to compute

 *  implied bounds of columns (variables x[j]) in this row.

 *

 *  The routine stores implied column bounds l'[j] and u'[j] in column

 *  descriptors (NPPCOL); it does not change current column bounds l[j]

 *  and u[j]. (Implied column bounds can be then used to strengthen the

 *  current column bounds; see the routines npp_implied_lower and

 *  npp_implied_upper).

 *

 *  ALGORITHM

 *

 *  Current column bounds (2) define implied lower and upper bounds of

 *  row (1) as follows:

 *

 *     L'[p] = inf sum a[p,j] x[j] =

 *                  j

 *                                                                    (3)

 *           =  sum  a[p,j] l[j] +  sum  a[p,j] u[j],

 *            j in Jp             j in Jn

 *

 *     U'[p] = sup sum a[p,j] x[j] =

 *                                                                    (4)

 *           =  sum  a[p,j] u[j] +  sum  a[p,j] l[j],

 *            j in Jp             j in Jn

 *

 *     Jp = {j: a[p,j] > 0},  Jn = {j: a[p,j] < 0}.                   (5)

 *

 *  (Note that bounds of all columns in row p are assumed to be correct,

 *  so L'[p] <= U'[p].)

 *

 *  If L[p] > L'[p] and/or U[p] < U'[p], the lower and/or upper bound of

 *  row (1) can be active, in which case such row defines implied bounds

 *  of its variables.

 *

 *  Let x[k] be some variable having in row (1) coefficient a[p,k] != 0.

 *  Consider a case when row lower bound can be active (L[p] > L'[p]):

 *

 *     sum a[p,j] x[j] >= L[p]  ==>

 *      j

 *

 *     sum a[p,j] x[j] + a[p,k] x[k] >= L[p]  ==>

 *     j!=k

 *                                                                    (6)

 *     a[p,k] x[k] >= L[p] - sum a[p,j] x[j]  ==>

 *                           j!=k

 *

 *     a[p,k] x[k] >= L[p,k],

 *

 *  where

 *

 *     L[p,k] = inf(L[p] - sum a[p,j] x[j]) =

 *                         j!=k

 *

 *            = L[p] - sup sum a[p,j] x[j] =                          (7)

 *                         j!=k

 *

 *            = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j].

 *                    j in Jp\{k}       j in Jn\{k}

 *

 *  Thus:

 *

 *     x[k] >= l'[k] = L[p,k] / a[p,k],  if a[p,k] > 0,               (8)

 *

 *     x[k] <= u'[k] = L[p,k] / a[p,k],  if a[p,k] < 0.               (9)

 *

 *  where l'[k] and u'[k] are implied lower and upper bounds of variable

 *  x[k], resp.

 *

 *  Now consider a similar case when row upper bound can be active

 *  (U[p] < U'[p]):

 *

 *     sum a[p,j] x[j] <= U[p]  ==>

 *      j

 *

 *     sum a[p,j] x[j] + a[p,k] x[k] <= U[p]  ==>

 *     j!=k

 *                                                                   (10)

 *     a[p,k] x[k] <= U[p] - sum a[p,j] x[j]  ==>

 *                           j!=k

 *

 *     a[p,k] x[k] <= U[p,k],

 *

 *  where:

 *

 *     U[p,k] = sup(U[p] - sum a[p,j] x[j]) =

 *                         j!=k

 *

 *            = U[p] - inf sum a[p,j] x[j] =                         (11)

 *                         j!=k

 *

 *            = U[p] - sum a[p,j] l[j] - sum a[p,j] u[j].

 *                    j in Jp\{k}       j in Jn\{k}

 *

 *  Thus:

 *

 *     x[k] <= u'[k] = U[p,k] / a[p,k],  if a[p,k] > 0,              (12)

 *

 *     x[k] >= l'[k] = U[p,k] / a[p,k],  if a[p,k] < 0.              (13)

 *

 *  Note that in formulae (8), (9), (12), and (13) coefficient a[p,k]

 *  must not be too small in magnitude relatively to other non-zero

 *  coefficients in row (1), i.e. the following condition must hold:

 *

 *     |a[p,k]| >= eps * max(1, |a[p,j]|),                           (14)

 *                        j

 *

 *  where eps is a relative tolerance for constraint coefficients.

 *  Otherwise the implied column bounds can be numerical inreliable. For