/* glpnpp02.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_free_row - process free (unbounded) row

*

*  SYNOPSIS

*

*  #include "glpnpp.h"

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

*

*  DESCRIPTION

*

*  The routine npp_free_row processes row p, which is free (i.e. has

*  no finite bounds):

*

*     -inf < sum a[p,j] x[j] < +inf.                                 (1)

*             j

*

*  PROBLEM TRANSFORMATION

*

*  Constraint (1) cannot be active, so it is redundant and can be

*  removed from the original problem.

*

*  Removing row p leads to removing a column of multiplier pi[p] for

*  this row in the dual system. Since row p has no bounds, pi[p] = 0,

*  so removing the column does not affect the dual solution.

*

*  RECOVERING BASIC SOLUTION

*

*  In solution to the original problem row p is inactive constraint,

*  so it is assigned status GLP_BS, and multiplier pi[p] is assigned

*  zero value.

*

*  RECOVERING INTERIOR-POINT SOLUTION

*

*  In solution to the original problem row p is inactive constraint,

*  so its multiplier pi[p] is assigned zero value.

*

*  RECOVERING MIP SOLUTION

*

*  None needed. */

struct free_row

{     /* free (unbounded) row */

      int p;

      /* row reference number */

};

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

void npp_free_row(NPP *npp, NPPROW *p)

{     /* process free (unbounded) row */

      struct free_row *info;

      /* the row must be free */

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

      /* create transformation stack entry */

      info = npp_push_tse(npp,

         rcv_free_row, sizeof(struct free_row));

      info->p = p->i;

      /* remove the row from the problem */

      npp_del_row(npp, p);

      return;

}

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

{     /* recover free (unbounded) row */

      struct free_row *info = _info;

      if (npp->sol == GLP_SOL)

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

      if (npp->sol != GLP_MIP)

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

      return 0;

}

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

*  NAME

*

*  npp_geq_row - process row of 'not less than' type

*

*  SYNOPSIS

*

*  #include "glpnpp.h"

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

*

*  DESCRIPTION

*

*  The routine npp_geq_row processes row p, which is 'not less than'

*  inequality constraint:

*

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

*              j

*

*  where L[p] < U[p], and upper bound may not exist (U[p] = +oo).

*

*  PROBLEM TRANSFORMATION

*

*  Constraint (1) can be replaced by equality constraint:

*

*     sum a[p,j] x[j] - s = L[p],                                    (2)

*      j

*

*  where

*

*     0 <= s (<= U[p] - L[p])                                        (3)

*

*  is a non-negative surplus variable.

*

*  Since in the primal system there appears column s having the only

*  non-zero coefficient in row p, in the dual system there appears a

*  new row:

*

*     (-1) pi[p] + lambda = 0,                                       (4)

*

*  where (-1) is coefficient of column s in row p, pi[p] is multiplier

*  of row p, lambda is multiplier of column q, 0 is coefficient of

*  column s in the objective row.

*

*  RECOVERING BASIC SOLUTION

*

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

*  by its status and status of column q in solution to the transformed

*  problem as follows:

*

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

*     |         Transformed problem          | Original problem |

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

*     | Status of row p | Status of column s | Status of row p  |

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

*     |     GLP_BS      |       GLP_BS       |       N/A        |

*     |     GLP_BS      |       GLP_NL       |      GLP_BS      |

*     |     GLP_BS      |       GLP_NU       |      GLP_BS      |

*     |     GLP_NS      |       GLP_BS       |      GLP_BS      |

*     |     GLP_NS      |       GLP_NL       |      GLP_NL      |

*     |     GLP_NS      |       GLP_NU       |      GLP_NU      |

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

*

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

*  is the same as in solution to the transformed problem.

*

*  1. In solution to the transformed problem row p and column q cannot

*     be basic at the same time; otherwise the basis matrix would have

*     two linear dependent columns: unity column of auxiliary variable

*     of row p and unity column of variable s.

*

*  2. Though in the transformed problem row p is equality constraint,

*     it may be basic due to primal degenerate solution.

*

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

{     /* inequality constraint row */

      int p;

      /* row reference number */

      int s;

      /* column reference number for slack/surplus variable */

};

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

void npp_geq_row(NPP *npp, NPPROW *p)

{     /* process row of 'not less than' type */

      struct ineq_row *info;

      NPPCOL *s;

      /* the row must have lower bound */

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

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

      /* create column for surplus variable */

      s = npp_add_col(npp);

      s->lb = 0.0;

      s->ub = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub - p->lb);

      /* and add it to the transformed problem */

      npp_add_aij(npp, p, s, -1.0);

      /* create transformation stack entry */

      info = npp_push_tse(npp,

         rcv_geq_row, sizeof(struct ineq_row));

      info->p = p->i;

      info->s = s->j;

      /* replace the row by equality constraint */

      p->ub = p->lb;

      return;

}

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

{     /* recover row of 'not less than' type */

      struct ineq_row *info = _info;

      if (npp->sol == GLP_SOL)

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

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

            {  npp_error();

               return 1;

            }

            else if (npp->c_stat[info->s] == GLP_NL ||

                     npp->c_stat[info->s] == GLP_NU)

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

            else

            {  npp_error();

               return 1;

            }

         }

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

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

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

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

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

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

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

            else

            {  npp_error();

               return 1;

            }

         }

         else

         {  npp_error();

            return 1;

         }

      }

      return 0;

}

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

*  NAME

*

*  npp_leq_row - process row of 'not greater than' type

*

*  SYNOPSIS

*

*  #include "glpnpp.h"

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

*

*  DESCRIPTION

*

*  The routine npp_leq_row processes row p, which is 'not greater than'

*  inequality constraint:

*

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

*                j

*

*  where L[p] < U[p], and lower bound may not exist (L[p] = +oo).

*

*  PROBLEM TRANSFORMATION

*

*  Constraint (1) can be replaced by equality constraint:

*

*     sum a[p,j] x[j] + s = L[p],                                    (2)

*      j

*

*  where

*

*     0 <= s (<= U[p] - L[p])                                        (3)

*

*  is a non-negative slack variable.

*

*  Since in the primal system there appears column s having the only

*  non-zero coefficient in row p, in the dual system there appears a

*  new row:

*

*     (+1) pi[p] + lambda = 0,                                       (4)

*

*  where (+1) is coefficient of column s in row p, pi[p] is multiplier

*  of row p, lambda is multiplier of column q, 0 is coefficient of

*  column s in the objective row.

*

*  RECOVERING BASIC SOLUTION

*

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

*  by its status and status of column q in solution to the transformed

*  problem as follows:

*

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

*     |         Transformed problem          | Original problem |

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

*     | Status of row p | Status of column s | Status of row p  |

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

*     |     GLP_BS      |       GLP_BS       |       N/A        |

*     |     GLP_BS      |       GLP_NL       |      GLP_BS      |

*     |     GLP_BS      |       GLP_NU       |      GLP_BS      |

*     |     GLP_NS      |       GLP_BS       |      GLP_BS      |

*     |     GLP_NS      |       GLP_NL       |      GLP_NU      |

*     |     GLP_NS      |       GLP_NU       |      GLP_NL      |

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

*

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

*  is the same as in solution to the transformed problem.

*

*  1. In solution to the transformed problem row p and column q cannot

*     be basic at the same time; otherwise the basis matrix would have

*     two linear dependent columns: unity column of auxiliary variable

*     of row p and unity column of variable s.

*

*  2. Though in the transformed problem row p is equality constraint,

*     it may be basic due to primal degeneracy.

*

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

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

void npp_leq_row(NPP *npp, NPPROW *p)

{     /* process row of 'not greater than' type */

      struct ineq_row *info;

      NPPCOL *s;

      /* the row must have upper bound */

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

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

      /* create column for slack variable */

      s = npp_add_col(npp);

      s->lb = 0.0;

      s->ub = (p->lb == -DBL_MAX ? +DBL_MAX : p->ub - p->lb);

      /* and add it to the transformed problem */

      npp_add_aij(npp, p, s, +1.0);

      /* create transformation stack entry */

      info = npp_push_tse(npp,

         rcv_leq_row, sizeof(struct ineq_row));

      info->p = p->i;

      info->s = s->j;

      /* replace the row by equality constraint */

      p->lb = p->ub;

      return;

}

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

{     /* recover row of 'not greater than' type */

      struct ineq_row *info = _info;

      if (npp->sol == GLP_SOL)

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

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

            {  npp_error();

               return 1;

            }

            else if (npp->c_stat[info->s] == GLP_NL ||

                     npp->c_stat[info->s] == GLP_NU)

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

            else

            {  npp_error();

               return 1;

            }

         }

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

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

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

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

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

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

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

            else

            {  npp_error();

               return 1;

            }

         }

         else

         {  npp_error();

            return 1;

         }

      }

      return 0;

}

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

*  NAME

*

*  npp_free_col - process free (unbounded) column

*

*  SYNOPSIS

*

*  #include "glpnpp.h"

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

*

*  DESCRIPTION

*

*  The routine npp_free_col processes column q, which is free (i.e. has

*  no finite bounds):

*

*     -oo < x[q] < +oo.                                              (1)

*

*  PROBLEM TRANSFORMATION

*

*  Free (unbounded) variable can be replaced by the difference of two

*  non-negative variables:

*

*     x[q] = s' - s'',   s', s'' >= 0.                               (2)

*

*  Assuming that in the transformed problem x[q] becomes s',

*  transformation (2) causes new column s'' to appear, which differs

*  from column s' only in the sign of coefficients in constraint and

*  objective rows. Thus, if in the dual system the following row

*  corresponds to column s':

*

*     sum a[i,q] pi[i] + lambda' = c[q],                             (3)

*      i

*

*  the row which corresponds to column s'' is the following:

*

*     sum (-a[i,q]) pi[i] + lambda'' = -c[q].                        (4)

*      i

*

*  Then from (3) and (4) it follows that:

*

*     lambda' + lambda'' = 0   =>   lambda' = lmabda'' = 0,          (5)

*

*  where lambda' and lambda'' are multipliers for columns s' and s'',

*  resp.

*

*  RECOVERING BASIC SOLUTION

*

*  With respect to (5) status of column q in solution to the original

*  problem is determined by statuses of columns s' and s'' in solution

*  to the transformed problem as follows:

*

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

*     |         Transformed problem          | Original problem |

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

*     | Status of col s' | Status of col s'' | Status of col q  |

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

*     |      GLP_BS      |      GLP_BS       |       N/A        |

*     |      GLP_BS      |      GLP_NL       |      GLP_BS      |

*     |      GLP_NL      |      GLP_BS       |      GLP_BS      |

*     |      GLP_NL      |      GLP_NL       |      GLP_NF      |

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

*

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

*

*  1. In solution to the transformed problem columns s' and s'' cannot

*     be basic at the same time, because they differ only in the sign,

*     hence, are linear dependent.

*

*  2. Though column q is free, it can be non-basic due to dual

*     degeneracy.

*

*  3. If column q is integral, columns s' and s'' are also integral.

*

*  RECOVERING INTERIOR-POINT SOLUTION

*

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

*

*  RECOVERING MIP SOLUTION

*

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

struct free_col

{     /* free (unbounded) column */

      int q;

      /* column reference number for variables x[q] and s' */

      int s;

      /* column reference number for variable s'' */

};

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

void npp_free_col(NPP *npp, NPPCOL *q)

{     /* process free (unbounded) column */

      struct free_col *info;

      NPPCOL *s;

      NPPAIJ *aij;

      /* the column must be free */

      xassert(q->lb == -DBL_MAX && q->ub == +DBL_MAX);

      /* variable x[q] becomes s' */

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

      /* create variable s'' */

      s = npp_add_col(npp);

      s->is_int = q->is_int;

      s->lb = 0.0, s->ub = +DBL_MAX;

      /* duplicate objective coefficient */

      s->coef = -q->coef;

      /* duplicate column of the constraint matrix */

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

         npp_add_aij(npp, aij->row, s, -aij->val);

      /* create transformation stack entry */

      info = npp_push_tse(npp,

         rcv_free_col, sizeof(struct free_col));

      info->q = q->j;

      info->s = s->j;

      return;

}

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

{     /* recover free (unbounded) column */

      struct free_col *info = _info;

      if (npp->sol == GLP_SOL)

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

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

            {  npp_error();

               return 1;

            }

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

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

            else

            {  npp_error();

               return -1;

            }

         }

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

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

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

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

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

            else

            {  npp_error();

               return -1;

            }

         }

         else

         {  npp_error();

            return -1;

         }

      }

      /* compute value of x[q] with formula (2) */

      npp->c_value[info->q] -= npp->c_value[info->s];

      return 0;

}

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

*  NAME

*

*  npp_lbnd_col - process column with (non-zero) lower bound

*

*  SYNOPSIS

*

*  #include "glpnpp.h"

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

*

*  DESCRIPTION

*

*  The routine npp_lbnd_col processes column q, which has (non-zero)

*  lower bound:

*

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

*

*  where l[q] < u[q], and upper bound may not exist (u[q] = +oo).

*

*  PROBLEM TRANSFORMATION

*

*  Column q can be replaced as follows:

*

*     x[q] = l[q] + s,                                               (2)

*

*  where

*

*     0 <= s (<= u[q] - l[q])                                        (3)

*

*  is a non-negative variable.

*

*  Substituting x[q] from (2) into the objective row, we have:

*

*     z = sum c[j] x[j] + c0 =

*          j

*

*       = sum c[j] x[j] + c[q] x[q] + c0 =

*         j!=q

*

*       = sum c[j] x[j] + c[q] (l[q] + s) + c0 =

*         j!=q

*

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

*

*  where

*

*     c~0 = c0 + c[q] l[q]                                           (4)

*

*  is the constant term of the objective in the transformed problem.

*  Similarly, substituting x[q] into constraint row i, we have:

*

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

*              j

*

*     L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i]  ==>

*             j!=q

*

*     L[i] <= sum a[i,j] x[j] + a[i,q] (l[q] + s) <= U[i]  ==>

*             j!=q

*

*     L~[i] <= sum a[i,j] x[j] + a[i,q] s <= U~[i],

*              j!=q

*

*  where

*

*     L~[i] = L[i] - a[i,q] l[q],  U~[i] = U[i] - a[i,q] l[q]        (5)

*

*  are lower and upper bounds of row i in the transformed problem,

*  resp.

*

*  Transformation (2) does not affect the dual system.

*

*  RECOVERING BASIC SOLUTION

*

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

*  as in solution to the transformed problem (GLP_BS, GLP_NL or GLP_NU).

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

*

*  RECOVERING INTERIOR-POINT SOLUTION

*

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

*

*  RECOVERING MIP SOLUTION

*

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

struct bnd_col

{     /* bounded column */

      int q;

      /* column reference number for variables x[q] and s */

      double bnd;

      /* lower/upper bound l[q] or u[q] */

};

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

void npp_lbnd_col(NPP *npp, NPPCOL *q)

{     /* process column with (non-zero) lower bound */

      struct bnd_col *info;

      NPPROW *i;

      NPPAIJ *aij;

      /* the column must have non-zero lower bound */

      xassert(q->lb != 0.0);

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

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

      /* create transformation stack entry */

      info = npp_push_tse(npp,

         rcv_lbnd_col, sizeof(struct bnd_col));

      info->q = q->j;

      info->bnd = q->lb;

      /* substitute x[q] into objective row */

      npp->c0 += q->coef * q->lb;

      /* substitute x[q] into constraint rows */

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

      {  i = aij->row;

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

            i->ub = (i->lb -= aij->val * q->lb);

         else

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

               i->lb -= aij->val * q->lb;

            if (i->ub != +DBL_MAX)

               i->ub -= aij->val * q->lb;

         }

      }

      /* column x[q] becomes column s */

      if (q->ub != +DBL_MAX)

         q->ub -= q->lb;

      q->lb = 0.0;

      return;

}

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

{     /* recover column with (non-zero) lower bound */

      struct bnd_col *info = _info;

      if (npp->sol == GLP_SOL)

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

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

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

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

         else

         {  npp_error();

            return 1;

         }

      }

      /* compute value of x[q] with formula (2) */

      npp->c_value[info->q] = info->bnd + npp->c_value[info->q];

      return 0;

}

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

*  NAME

*

*  npp_ubnd_col - process column with upper bound

*

*  SYNOPSIS

*

*  #include "glpnpp.h"

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

*

*  DESCRIPTION

*

*  The routine npp_ubnd_col processes column q, which has upper bound:

*

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

*

*  where l[q] < u[q], and lower bound may not exist (l[q] = -oo).

*

*  PROBLEM TRANSFORMATION

*

*  Column q can be replaced as follows:

*

*     x[q] = u[q] - s,                                               (2)

*

*  where

*

*     0 <= s (<= u[q] - l[q])                                        (3)

*

*  is a non-negative variable.

*

*  Substituting x[q] from (2) into the objective row, we have:

*

*     z = sum c[j] x[j] + c0 =

*          j

*

*       = sum c[j] x[j] + c[q] x[q] + c0 =

*         j!=q

*

*       = sum c[j] x[j] + c[q] (u[q] - s) + c0 =

*         j!=q

*

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

*

*  where

*

*     c~0 = c0 + c[q] u[q]                                           (4)

*

*  is the constant term of the objective in the transformed problem.

*  Similarly, substituting x[q] into constraint row i, we have:

*

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

*              j

*

*     L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i]  ==>

*             j!=q

*

*     L[i] <= sum a[i,j] x[j] + a[i,q] (u[q] - s) <= U[i]  ==>

*             j!=q

*

*     L~[i] <= sum a[i,j] x[j] - a[i,q] s <= U~[i],

*              j!=q

*

*  where

*

*     L~[i] = L[i] - a[i,q] u[q],  U~[i] = U[i] - a[i,q] u[q]        (5)

*

*  are lower and upper bounds of row i in the transformed problem,

*  resp.

*

*  Note that in the transformed problem coefficients c[q] and a[i,q]

*  change their sign. Thus, the row of the dual system corresponding to

*  column q:

*

*     sum a[i,q] pi[i] + lambda[q] = c[q]                            (6)

*      i

*

*  in the transformed problem becomes the following:

*

*     sum (-a[i,q]) pi[i] + lambda[s] = -c[q].                       (7)

*      i

*

*  Therefore:

*

*     lambda[q] = - lambda[s],                                       (8)

*

*  where lambda[q] is multiplier for column q, lambda[s] is multiplier

*  for column s.

*

*  RECOVERING BASIC SOLUTION

*

*  With respect to (8) status of column q in solution to the original

*  problem is determined by status of column s in solution to the

*  transformed problem as follows:

*

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

*     |  Status of column s   | Status of column q |

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

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

*     |        GLP_BS         |       GLP_BS       |

*     |        GLP_NL         |       GLP_NU       |

*     |        GLP_NU         |       GLP_NL       |

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

*

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

*

*  RECOVERING INTERIOR-POINT SOLUTION

*

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

*

*  RECOVERING MIP SOLUTION

*

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

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

void npp_ubnd_col(NPP *npp, NPPCOL *q)

{     /* process column with upper bound */

      struct bnd_col *info;

      NPPROW *i;

      NPPAIJ *aij;

      /* the column must have upper bound */

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

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

      /* create transformation stack entry */

      info = npp_push_tse(npp,

         rcv_ubnd_col, sizeof(struct bnd_col));

      info->q = q->j;

      info->bnd = q->ub;

      /* substitute x[q] into objective row */

      npp->c0 += q->coef * q->ub;

      q->coef = -q->coef;

      /* substitute x[q] into constraint rows */

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

      {  i = aij->row;

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

            i->ub = (i->lb -= aij->val * q->ub);

         else

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

               i->lb -= aij->val * q->ub;

            if (i->ub != +DBL_MAX)

               i->ub -= aij->val * q->ub;

         }

         aij->val = -aij->val;

      }

      /* column x[q] becomes column s */

      if (q->lb != -DBL_MAX)

         q->ub -= q->lb;

      else

         q->ub = +DBL_MAX;

      q->lb = 0.0;

      return;

}

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

{     /* recover column with upper bound */

      struct bnd_col *info = _info;

      if (npp->sol == GLP_BS)

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

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

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

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

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

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

         else

         {  npp_error();

            return 1;

         }

      }

      /* compute value of x[q] with formula (2) */

      npp->c_value[info->q] = info->bnd - npp->c_value[info->q];

      return 0;

}

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

*  NAME

*

*  npp_dbnd_col - process non-negative column with upper bound

*

*  SYNOPSIS

*

*  #include "glpnpp.h"

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

*

*  DESCRIPTION

*

*  The routine npp_dbnd_col processes column q, which is non-negative

*  and has upper bound:

*

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

*

*  where u[q] > 0.

*

*  PROBLEM TRANSFORMATION

*

*  Upper bound of column q can be replaced by the following equality

*  constraint:

*

*     x[q] + s = u[q],                                               (2)

*

*  where s >= 0 is a non-negative complement variable.

*

*  Since in the primal system along with new row (2) there appears a

*  new column s having the only non-zero coefficient in this row, in

*  the dual system there appears a new row:

*

*     (+1)pi + lambda[s] = 0,                                        (3)

*

*  where (+1) is coefficient at column s in row (2), pi is multiplier

*  for row (2), lambda[s] is multiplier for column s, 0 is coefficient

*  at column s in the objective row.

*

*  RECOVERING BASIC SOLUTION

*

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

*  by its status and status of column s in solution to the transformed

*  problem as follows:

*

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

*     |         Transformed problem       | Original problem |

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

*     | Status of col q | Status of col s | Status of col q  |

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

*     |     GLP_BS      |     GLP_BS      |      GLP_BS      |

*     |     GLP_BS      |     GLP_NL      |      GLP_NU      |

*     |     GLP_NL      |     GLP_BS      |      GLP_NL      |

*     |     GLP_NL      |     GLP_NL      |      GLP_NL (*)  |

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

*

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

*  in solution to the transformed problem.

*

*  1. Formally, in solution to the transformed problem columns q and s

*     cannot be non-basic at the same time, since the constraint (2)

*     would be violated. However, if u[q] is close to zero, violation

*     may be less than a working precision even if both columns q and s

*     are non-basic. In this degenerate case row (2) can be only basic,

*     i.e. non-active constraint (otherwise corresponding row of the

*     basis matrix would be zero). This allows to pivot out auxiliary

*     variable and pivot in column s, in which case the row becomes

*     active while column s becomes basic.

*

*  2. If column q is integral, column s is also integral.

*

*  RECOVERING INTERIOR-POINT SOLUTION

*

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

*  in solution to the transformed problem.

*

*  RECOVERING MIP SOLUTION

*

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

*  in solution to the transformed problem. */

struct dbnd_col

{     /* double-bounded column */

      int q;

      /* column reference number for variable x[q] */

      int s;

      /* column reference number for complement variable s */

};

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

void npp_dbnd_col(NPP *npp, NPPCOL *q)

{     /* process non-negative column with upper bound */

      struct dbnd_col *info;

      NPPROW *p;

      NPPCOL *s;

      /* the column must be non-negative with upper bound */

      xassert(q->lb == 0.0);

      xassert(q->ub > 0.0);

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

      /* create variable s */

      s = npp_add_col(npp);

      s->is_int = q->is_int;

      s->lb = 0.0, s->ub = +DBL_MAX;

      /* create equality constraint (2) */

      p = npp_add_row(npp);

      p->lb = p->ub = q->ub;

      npp_add_aij(npp, p, q, +1.0);

      npp_add_aij(npp, p, s, +1.0);

      /* create transformation stack entry */

      info = npp_push_tse(npp,

         rcv_dbnd_col, sizeof(struct dbnd_col));

      info->q = q->j;

      info->s = s->j;

      /* remove upper bound of x[q] */

      q->ub = +DBL_MAX;

      return;

}

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

{     /* recover non-negative column with upper bound */

      struct dbnd_col *info = _info;

      if (npp->sol == GLP_BS)

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