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