alpar@9: /* glpnpp02.c */ alpar@9: alpar@9: /*********************************************************************** alpar@9: * This code is part of GLPK (GNU Linear Programming Kit). alpar@9: * alpar@9: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@9: * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, alpar@9: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@9: * E-mail: . alpar@9: * alpar@9: * GLPK is free software: you can redistribute it and/or modify it alpar@9: * under the terms of the GNU General Public License as published by alpar@9: * the Free Software Foundation, either version 3 of the License, or alpar@9: * (at your option) any later version. alpar@9: * alpar@9: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@9: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@9: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@9: * License for more details. alpar@9: * alpar@9: * You should have received a copy of the GNU General Public License alpar@9: * along with GLPK. If not, see . alpar@9: ***********************************************************************/ alpar@9: alpar@9: #include "glpnpp.h" alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_free_row - process free (unbounded) row alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_free_row(NPP *npp, NPPROW *p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_free_row processes row p, which is free (i.e. has alpar@9: * no finite bounds): alpar@9: * alpar@9: * -inf < sum a[p,j] x[j] < +inf. (1) alpar@9: * j alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Constraint (1) cannot be active, so it is redundant and can be alpar@9: * removed from the original problem. alpar@9: * alpar@9: * Removing row p leads to removing a column of multiplier pi[p] for alpar@9: * this row in the dual system. Since row p has no bounds, pi[p] = 0, alpar@9: * so removing the column does not affect the dual solution. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * In solution to the original problem row p is inactive constraint, alpar@9: * so it is assigned status GLP_BS, and multiplier pi[p] is assigned alpar@9: * zero value. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * In solution to the original problem row p is inactive constraint, alpar@9: * so its multiplier pi[p] is assigned zero value. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: struct free_row alpar@9: { /* free (unbounded) row */ alpar@9: int p; alpar@9: /* row reference number */ alpar@9: }; alpar@9: alpar@9: static int rcv_free_row(NPP *npp, void *info); alpar@9: alpar@9: void npp_free_row(NPP *npp, NPPROW *p) alpar@9: { /* process free (unbounded) row */ alpar@9: struct free_row *info; alpar@9: /* the row must be free */ alpar@9: xassert(p->lb == -DBL_MAX && p->ub == +DBL_MAX); alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_free_row, sizeof(struct free_row)); alpar@9: info->p = p->i; alpar@9: /* remove the row from the problem */ alpar@9: npp_del_row(npp, p); alpar@9: return; alpar@9: } alpar@9: alpar@9: static int rcv_free_row(NPP *npp, void *_info) alpar@9: { /* recover free (unbounded) row */ alpar@9: struct free_row *info = _info; alpar@9: if (npp->sol == GLP_SOL) alpar@9: npp->r_stat[info->p] = GLP_BS; alpar@9: if (npp->sol != GLP_MIP) alpar@9: npp->r_pi[info->p] = 0.0; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_geq_row - process row of 'not less than' type alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_geq_row(NPP *npp, NPPROW *p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_geq_row processes row p, which is 'not less than' alpar@9: * inequality constraint: alpar@9: * alpar@9: * L[p] <= sum a[p,j] x[j] (<= U[p]), (1) alpar@9: * j alpar@9: * alpar@9: * where L[p] < U[p], and upper bound may not exist (U[p] = +oo). alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Constraint (1) can be replaced by equality constraint: alpar@9: * alpar@9: * sum a[p,j] x[j] - s = L[p], (2) alpar@9: * j alpar@9: * alpar@9: * where alpar@9: * alpar@9: * 0 <= s (<= U[p] - L[p]) (3) alpar@9: * alpar@9: * is a non-negative surplus variable. alpar@9: * alpar@9: * Since in the primal system there appears column s having the only alpar@9: * non-zero coefficient in row p, in the dual system there appears a alpar@9: * new row: alpar@9: * alpar@9: * (-1) pi[p] + lambda = 0, (4) alpar@9: * alpar@9: * where (-1) is coefficient of column s in row p, pi[p] is multiplier alpar@9: * of row p, lambda is multiplier of column q, 0 is coefficient of alpar@9: * column s in the objective row. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * Status of row p in solution to the original problem is determined alpar@9: * by its status and status of column q in solution to the transformed alpar@9: * problem as follows: alpar@9: * alpar@9: * +--------------------------------------+------------------+ alpar@9: * | Transformed problem | Original problem | alpar@9: * +-----------------+--------------------+------------------+ alpar@9: * | Status of row p | Status of column s | Status of row p | alpar@9: * +-----------------+--------------------+------------------+ alpar@9: * | GLP_BS | GLP_BS | N/A | alpar@9: * | GLP_BS | GLP_NL | GLP_BS | alpar@9: * | GLP_BS | GLP_NU | GLP_BS | alpar@9: * | GLP_NS | GLP_BS | GLP_BS | alpar@9: * | GLP_NS | GLP_NL | GLP_NL | alpar@9: * | GLP_NS | GLP_NU | GLP_NU | alpar@9: * +-----------------+--------------------+------------------+ alpar@9: * alpar@9: * Value of row multiplier pi[p] in solution to the original problem alpar@9: * is the same as in solution to the transformed problem. alpar@9: * alpar@9: * 1. In solution to the transformed problem row p and column q cannot alpar@9: * be basic at the same time; otherwise the basis matrix would have alpar@9: * two linear dependent columns: unity column of auxiliary variable alpar@9: * of row p and unity column of variable s. alpar@9: * alpar@9: * 2. Though in the transformed problem row p is equality constraint, alpar@9: * it may be basic due to primal degenerate solution. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Value of row multiplier pi[p] in solution to the original problem alpar@9: * is the same as in solution to the transformed problem. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: struct ineq_row alpar@9: { /* inequality constraint row */ alpar@9: int p; alpar@9: /* row reference number */ alpar@9: int s; alpar@9: /* column reference number for slack/surplus variable */ alpar@9: }; alpar@9: alpar@9: static int rcv_geq_row(NPP *npp, void *info); alpar@9: alpar@9: void npp_geq_row(NPP *npp, NPPROW *p) alpar@9: { /* process row of 'not less than' type */ alpar@9: struct ineq_row *info; alpar@9: NPPCOL *s; alpar@9: /* the row must have lower bound */ alpar@9: xassert(p->lb != -DBL_MAX); alpar@9: xassert(p->lb < p->ub); alpar@9: /* create column for surplus variable */ alpar@9: s = npp_add_col(npp); alpar@9: s->lb = 0.0; alpar@9: s->ub = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub - p->lb); alpar@9: /* and add it to the transformed problem */ alpar@9: npp_add_aij(npp, p, s, -1.0); alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_geq_row, sizeof(struct ineq_row)); alpar@9: info->p = p->i; alpar@9: info->s = s->j; alpar@9: /* replace the row by equality constraint */ alpar@9: p->ub = p->lb; alpar@9: return; alpar@9: } alpar@9: alpar@9: static int rcv_geq_row(NPP *npp, void *_info) alpar@9: { /* recover row of 'not less than' type */ alpar@9: struct ineq_row *info = _info; alpar@9: if (npp->sol == GLP_SOL) alpar@9: { if (npp->r_stat[info->p] == GLP_BS) alpar@9: { if (npp->c_stat[info->s] == GLP_BS) alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: else if (npp->c_stat[info->s] == GLP_NL || alpar@9: npp->c_stat[info->s] == GLP_NU) alpar@9: npp->r_stat[info->p] = GLP_BS; alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: else if (npp->r_stat[info->p] == GLP_NS) alpar@9: { if (npp->c_stat[info->s] == GLP_BS) alpar@9: npp->r_stat[info->p] = GLP_BS; alpar@9: else if (npp->c_stat[info->s] == GLP_NL) alpar@9: npp->r_stat[info->p] = GLP_NL; alpar@9: else if (npp->c_stat[info->s] == GLP_NU) alpar@9: npp->r_stat[info->p] = GLP_NU; alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_leq_row - process row of 'not greater than' type alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_leq_row(NPP *npp, NPPROW *p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_leq_row processes row p, which is 'not greater than' alpar@9: * inequality constraint: alpar@9: * alpar@9: * (L[p] <=) sum a[p,j] x[j] <= U[p], (1) alpar@9: * j alpar@9: * alpar@9: * where L[p] < U[p], and lower bound may not exist (L[p] = +oo). alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Constraint (1) can be replaced by equality constraint: alpar@9: * alpar@9: * sum a[p,j] x[j] + s = L[p], (2) alpar@9: * j alpar@9: * alpar@9: * where alpar@9: * alpar@9: * 0 <= s (<= U[p] - L[p]) (3) alpar@9: * alpar@9: * is a non-negative slack variable. alpar@9: * alpar@9: * Since in the primal system there appears column s having the only alpar@9: * non-zero coefficient in row p, in the dual system there appears a alpar@9: * new row: alpar@9: * alpar@9: * (+1) pi[p] + lambda = 0, (4) alpar@9: * alpar@9: * where (+1) is coefficient of column s in row p, pi[p] is multiplier alpar@9: * of row p, lambda is multiplier of column q, 0 is coefficient of alpar@9: * column s in the objective row. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * Status of row p in solution to the original problem is determined alpar@9: * by its status and status of column q in solution to the transformed alpar@9: * problem as follows: alpar@9: * alpar@9: * +--------------------------------------+------------------+ alpar@9: * | Transformed problem | Original problem | alpar@9: * +-----------------+--------------------+------------------+ alpar@9: * | Status of row p | Status of column s | Status of row p | alpar@9: * +-----------------+--------------------+------------------+ alpar@9: * | GLP_BS | GLP_BS | N/A | alpar@9: * | GLP_BS | GLP_NL | GLP_BS | alpar@9: * | GLP_BS | GLP_NU | GLP_BS | alpar@9: * | GLP_NS | GLP_BS | GLP_BS | alpar@9: * | GLP_NS | GLP_NL | GLP_NU | alpar@9: * | GLP_NS | GLP_NU | GLP_NL | alpar@9: * +-----------------+--------------------+------------------+ alpar@9: * alpar@9: * Value of row multiplier pi[p] in solution to the original problem alpar@9: * is the same as in solution to the transformed problem. alpar@9: * alpar@9: * 1. In solution to the transformed problem row p and column q cannot alpar@9: * be basic at the same time; otherwise the basis matrix would have alpar@9: * two linear dependent columns: unity column of auxiliary variable alpar@9: * of row p and unity column of variable s. alpar@9: * alpar@9: * 2. Though in the transformed problem row p is equality constraint, alpar@9: * it may be basic due to primal degeneracy. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Value of row multiplier pi[p] in solution to the original problem alpar@9: * is the same as in solution to the transformed problem. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: static int rcv_leq_row(NPP *npp, void *info); alpar@9: alpar@9: void npp_leq_row(NPP *npp, NPPROW *p) alpar@9: { /* process row of 'not greater than' type */ alpar@9: struct ineq_row *info; alpar@9: NPPCOL *s; alpar@9: /* the row must have upper bound */ alpar@9: xassert(p->ub != +DBL_MAX); alpar@9: xassert(p->lb < p->ub); alpar@9: /* create column for slack variable */ alpar@9: s = npp_add_col(npp); alpar@9: s->lb = 0.0; alpar@9: s->ub = (p->lb == -DBL_MAX ? +DBL_MAX : p->ub - p->lb); alpar@9: /* and add it to the transformed problem */ alpar@9: npp_add_aij(npp, p, s, +1.0); alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_leq_row, sizeof(struct ineq_row)); alpar@9: info->p = p->i; alpar@9: info->s = s->j; alpar@9: /* replace the row by equality constraint */ alpar@9: p->lb = p->ub; alpar@9: return; alpar@9: } alpar@9: alpar@9: static int rcv_leq_row(NPP *npp, void *_info) alpar@9: { /* recover row of 'not greater than' type */ alpar@9: struct ineq_row *info = _info; alpar@9: if (npp->sol == GLP_SOL) alpar@9: { if (npp->r_stat[info->p] == GLP_BS) alpar@9: { if (npp->c_stat[info->s] == GLP_BS) alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: else if (npp->c_stat[info->s] == GLP_NL || alpar@9: npp->c_stat[info->s] == GLP_NU) alpar@9: npp->r_stat[info->p] = GLP_BS; alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: else if (npp->r_stat[info->p] == GLP_NS) alpar@9: { if (npp->c_stat[info->s] == GLP_BS) alpar@9: npp->r_stat[info->p] = GLP_BS; alpar@9: else if (npp->c_stat[info->s] == GLP_NL) alpar@9: npp->r_stat[info->p] = GLP_NU; alpar@9: else if (npp->c_stat[info->s] == GLP_NU) alpar@9: npp->r_stat[info->p] = GLP_NL; alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_free_col - process free (unbounded) column alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_free_col(NPP *npp, NPPCOL *q); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_free_col processes column q, which is free (i.e. has alpar@9: * no finite bounds): alpar@9: * alpar@9: * -oo < x[q] < +oo. (1) alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Free (unbounded) variable can be replaced by the difference of two alpar@9: * non-negative variables: alpar@9: * alpar@9: * x[q] = s' - s'', s', s'' >= 0. (2) alpar@9: * alpar@9: * Assuming that in the transformed problem x[q] becomes s', alpar@9: * transformation (2) causes new column s'' to appear, which differs alpar@9: * from column s' only in the sign of coefficients in constraint and alpar@9: * objective rows. Thus, if in the dual system the following row alpar@9: * corresponds to column s': alpar@9: * alpar@9: * sum a[i,q] pi[i] + lambda' = c[q], (3) alpar@9: * i alpar@9: * alpar@9: * the row which corresponds to column s'' is the following: alpar@9: * alpar@9: * sum (-a[i,q]) pi[i] + lambda'' = -c[q]. (4) alpar@9: * i alpar@9: * alpar@9: * Then from (3) and (4) it follows that: alpar@9: * alpar@9: * lambda' + lambda'' = 0 => lambda' = lmabda'' = 0, (5) alpar@9: * alpar@9: * where lambda' and lambda'' are multipliers for columns s' and s'', alpar@9: * resp. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * With respect to (5) status of column q in solution to the original alpar@9: * problem is determined by statuses of columns s' and s'' in solution alpar@9: * to the transformed problem as follows: alpar@9: * alpar@9: * +--------------------------------------+------------------+ alpar@9: * | Transformed problem | Original problem | alpar@9: * +------------------+-------------------+------------------+ alpar@9: * | Status of col s' | Status of col s'' | Status of col q | alpar@9: * +------------------+-------------------+------------------+ alpar@9: * | GLP_BS | GLP_BS | N/A | alpar@9: * | GLP_BS | GLP_NL | GLP_BS | alpar@9: * | GLP_NL | GLP_BS | GLP_BS | alpar@9: * | GLP_NL | GLP_NL | GLP_NF | alpar@9: * +------------------+-------------------+------------------+ alpar@9: * alpar@9: * Value of column q is computed with formula (2). alpar@9: * alpar@9: * 1. In solution to the transformed problem columns s' and s'' cannot alpar@9: * be basic at the same time, because they differ only in the sign, alpar@9: * hence, are linear dependent. alpar@9: * alpar@9: * 2. Though column q is free, it can be non-basic due to dual alpar@9: * degeneracy. alpar@9: * alpar@9: * 3. If column q is integral, columns s' and s'' are also integral. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Value of column q is computed with formula (2). alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * Value of column q is computed with formula (2). */ alpar@9: alpar@9: struct free_col alpar@9: { /* free (unbounded) column */ alpar@9: int q; alpar@9: /* column reference number for variables x[q] and s' */ alpar@9: int s; alpar@9: /* column reference number for variable s'' */ alpar@9: }; alpar@9: alpar@9: static int rcv_free_col(NPP *npp, void *info); alpar@9: alpar@9: void npp_free_col(NPP *npp, NPPCOL *q) alpar@9: { /* process free (unbounded) column */ alpar@9: struct free_col *info; alpar@9: NPPCOL *s; alpar@9: NPPAIJ *aij; alpar@9: /* the column must be free */ alpar@9: xassert(q->lb == -DBL_MAX && q->ub == +DBL_MAX); alpar@9: /* variable x[q] becomes s' */ alpar@9: q->lb = 0.0, q->ub = +DBL_MAX; alpar@9: /* create variable s'' */ alpar@9: s = npp_add_col(npp); alpar@9: s->is_int = q->is_int; alpar@9: s->lb = 0.0, s->ub = +DBL_MAX; alpar@9: /* duplicate objective coefficient */ alpar@9: s->coef = -q->coef; alpar@9: /* duplicate column of the constraint matrix */ alpar@9: for (aij = q->ptr; aij != NULL; aij = aij->c_next) alpar@9: npp_add_aij(npp, aij->row, s, -aij->val); alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_free_col, sizeof(struct free_col)); alpar@9: info->q = q->j; alpar@9: info->s = s->j; alpar@9: return; alpar@9: } alpar@9: alpar@9: static int rcv_free_col(NPP *npp, void *_info) alpar@9: { /* recover free (unbounded) column */ alpar@9: struct free_col *info = _info; alpar@9: if (npp->sol == GLP_SOL) alpar@9: { if (npp->c_stat[info->q] == GLP_BS) alpar@9: { if (npp->c_stat[info->s] == GLP_BS) alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: else if (npp->c_stat[info->s] == GLP_NL) alpar@9: npp->c_stat[info->q] = GLP_BS; alpar@9: else alpar@9: { npp_error(); alpar@9: return -1; alpar@9: } alpar@9: } alpar@9: else if (npp->c_stat[info->q] == GLP_NL) alpar@9: { if (npp->c_stat[info->s] == GLP_BS) alpar@9: npp->c_stat[info->q] = GLP_BS; alpar@9: else if (npp->c_stat[info->s] == GLP_NL) alpar@9: npp->c_stat[info->q] = GLP_NF; alpar@9: else alpar@9: { npp_error(); alpar@9: return -1; alpar@9: } alpar@9: } alpar@9: else alpar@9: { npp_error(); alpar@9: return -1; alpar@9: } alpar@9: } alpar@9: /* compute value of x[q] with formula (2) */ alpar@9: npp->c_value[info->q] -= npp->c_value[info->s]; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_lbnd_col - process column with (non-zero) lower bound alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_lbnd_col(NPP *npp, NPPCOL *q); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_lbnd_col processes column q, which has (non-zero) alpar@9: * lower bound: alpar@9: * alpar@9: * l[q] <= x[q] (<= u[q]), (1) alpar@9: * alpar@9: * where l[q] < u[q], and upper bound may not exist (u[q] = +oo). alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Column q can be replaced as follows: alpar@9: * alpar@9: * x[q] = l[q] + s, (2) alpar@9: * alpar@9: * where alpar@9: * alpar@9: * 0 <= s (<= u[q] - l[q]) (3) alpar@9: * alpar@9: * is a non-negative variable. alpar@9: * alpar@9: * Substituting x[q] from (2) into the objective row, we have: alpar@9: * alpar@9: * z = sum c[j] x[j] + c0 = alpar@9: * j alpar@9: * alpar@9: * = sum c[j] x[j] + c[q] x[q] + c0 = alpar@9: * j!=q alpar@9: * alpar@9: * = sum c[j] x[j] + c[q] (l[q] + s) + c0 = alpar@9: * j!=q alpar@9: * alpar@9: * = sum c[j] x[j] + c[q] s + c~0, alpar@9: * alpar@9: * where alpar@9: * alpar@9: * c~0 = c0 + c[q] l[q] (4) alpar@9: * alpar@9: * is the constant term of the objective in the transformed problem. alpar@9: * Similarly, substituting x[q] into constraint row i, we have: alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] <= U[i] ==> alpar@9: * j alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==> alpar@9: * j!=q alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] + a[i,q] (l[q] + s) <= U[i] ==> alpar@9: * j!=q alpar@9: * alpar@9: * L~[i] <= sum a[i,j] x[j] + a[i,q] s <= U~[i], alpar@9: * j!=q alpar@9: * alpar@9: * where alpar@9: * alpar@9: * L~[i] = L[i] - a[i,q] l[q], U~[i] = U[i] - a[i,q] l[q] (5) alpar@9: * alpar@9: * are lower and upper bounds of row i in the transformed problem, alpar@9: * resp. alpar@9: * alpar@9: * Transformation (2) does not affect the dual system. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * Status of column q in solution to the original problem is the same alpar@9: * as in solution to the transformed problem (GLP_BS, GLP_NL or GLP_NU). alpar@9: * Value of column q is computed with formula (2). alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Value of column q is computed with formula (2). alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * Value of column q is computed with formula (2). */ alpar@9: alpar@9: struct bnd_col alpar@9: { /* bounded column */ alpar@9: int q; alpar@9: /* column reference number for variables x[q] and s */ alpar@9: double bnd; alpar@9: /* lower/upper bound l[q] or u[q] */ alpar@9: }; alpar@9: alpar@9: static int rcv_lbnd_col(NPP *npp, void *info); alpar@9: alpar@9: void npp_lbnd_col(NPP *npp, NPPCOL *q) alpar@9: { /* process column with (non-zero) lower bound */ alpar@9: struct bnd_col *info; alpar@9: NPPROW *i; alpar@9: NPPAIJ *aij; alpar@9: /* the column must have non-zero lower bound */ alpar@9: xassert(q->lb != 0.0); alpar@9: xassert(q->lb != -DBL_MAX); alpar@9: xassert(q->lb < q->ub); alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_lbnd_col, sizeof(struct bnd_col)); alpar@9: info->q = q->j; alpar@9: info->bnd = q->lb; alpar@9: /* substitute x[q] into objective row */ alpar@9: npp->c0 += q->coef * q->lb; alpar@9: /* substitute x[q] into constraint rows */ alpar@9: for (aij = q->ptr; aij != NULL; aij = aij->c_next) alpar@9: { i = aij->row; alpar@9: if (i->lb == i->ub) alpar@9: i->ub = (i->lb -= aij->val * q->lb); alpar@9: else alpar@9: { if (i->lb != -DBL_MAX) alpar@9: i->lb -= aij->val * q->lb; alpar@9: if (i->ub != +DBL_MAX) alpar@9: i->ub -= aij->val * q->lb; alpar@9: } alpar@9: } alpar@9: /* column x[q] becomes column s */ alpar@9: if (q->ub != +DBL_MAX) alpar@9: q->ub -= q->lb; alpar@9: q->lb = 0.0; alpar@9: return; alpar@9: } alpar@9: alpar@9: static int rcv_lbnd_col(NPP *npp, void *_info) alpar@9: { /* recover column with (non-zero) lower bound */ alpar@9: struct bnd_col *info = _info; alpar@9: if (npp->sol == GLP_SOL) alpar@9: { if (npp->c_stat[info->q] == GLP_BS || alpar@9: npp->c_stat[info->q] == GLP_NL || alpar@9: npp->c_stat[info->q] == GLP_NU) alpar@9: npp->c_stat[info->q] = npp->c_stat[info->q]; alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: /* compute value of x[q] with formula (2) */ alpar@9: npp->c_value[info->q] = info->bnd + npp->c_value[info->q]; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_ubnd_col - process column with upper bound alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_ubnd_col(NPP *npp, NPPCOL *q); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_ubnd_col processes column q, which has upper bound: alpar@9: * alpar@9: * (l[q] <=) x[q] <= u[q], (1) alpar@9: * alpar@9: * where l[q] < u[q], and lower bound may not exist (l[q] = -oo). alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Column q can be replaced as follows: alpar@9: * alpar@9: * x[q] = u[q] - s, (2) alpar@9: * alpar@9: * where alpar@9: * alpar@9: * 0 <= s (<= u[q] - l[q]) (3) alpar@9: * alpar@9: * is a non-negative variable. alpar@9: * alpar@9: * Substituting x[q] from (2) into the objective row, we have: alpar@9: * alpar@9: * z = sum c[j] x[j] + c0 = alpar@9: * j alpar@9: * alpar@9: * = sum c[j] x[j] + c[q] x[q] + c0 = alpar@9: * j!=q alpar@9: * alpar@9: * = sum c[j] x[j] + c[q] (u[q] - s) + c0 = alpar@9: * j!=q alpar@9: * alpar@9: * = sum c[j] x[j] - c[q] s + c~0, alpar@9: * alpar@9: * where alpar@9: * alpar@9: * c~0 = c0 + c[q] u[q] (4) alpar@9: * alpar@9: * is the constant term of the objective in the transformed problem. alpar@9: * Similarly, substituting x[q] into constraint row i, we have: alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] <= U[i] ==> alpar@9: * j alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==> alpar@9: * j!=q alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] + a[i,q] (u[q] - s) <= U[i] ==> alpar@9: * j!=q alpar@9: * alpar@9: * L~[i] <= sum a[i,j] x[j] - a[i,q] s <= U~[i], alpar@9: * j!=q alpar@9: * alpar@9: * where alpar@9: * alpar@9: * L~[i] = L[i] - a[i,q] u[q], U~[i] = U[i] - a[i,q] u[q] (5) alpar@9: * alpar@9: * are lower and upper bounds of row i in the transformed problem, alpar@9: * resp. alpar@9: * alpar@9: * Note that in the transformed problem coefficients c[q] and a[i,q] alpar@9: * change their sign. Thus, the row of the dual system corresponding to alpar@9: * column q: alpar@9: * alpar@9: * sum a[i,q] pi[i] + lambda[q] = c[q] (6) alpar@9: * i alpar@9: * alpar@9: * in the transformed problem becomes the following: alpar@9: * alpar@9: * sum (-a[i,q]) pi[i] + lambda[s] = -c[q]. (7) alpar@9: * i alpar@9: * alpar@9: * Therefore: alpar@9: * alpar@9: * lambda[q] = - lambda[s], (8) alpar@9: * alpar@9: * where lambda[q] is multiplier for column q, lambda[s] is multiplier alpar@9: * for column s. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * With respect to (8) status of column q in solution to the original alpar@9: * problem is determined by status of column s in solution to the alpar@9: * transformed problem as follows: alpar@9: * alpar@9: * +-----------------------+--------------------+ alpar@9: * | Status of column s | Status of column q | alpar@9: * | (transformed problem) | (original problem) | alpar@9: * +-----------------------+--------------------+ alpar@9: * | GLP_BS | GLP_BS | alpar@9: * | GLP_NL | GLP_NU | alpar@9: * | GLP_NU | GLP_NL | alpar@9: * +-----------------------+--------------------+ alpar@9: * alpar@9: * Value of column q is computed with formula (2). alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Value of column q is computed with formula (2). alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * Value of column q is computed with formula (2). */ alpar@9: alpar@9: static int rcv_ubnd_col(NPP *npp, void *info); alpar@9: alpar@9: void npp_ubnd_col(NPP *npp, NPPCOL *q) alpar@9: { /* process column with upper bound */ alpar@9: struct bnd_col *info; alpar@9: NPPROW *i; alpar@9: NPPAIJ *aij; alpar@9: /* the column must have upper bound */ alpar@9: xassert(q->ub != +DBL_MAX); alpar@9: xassert(q->lb < q->ub); alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_ubnd_col, sizeof(struct bnd_col)); alpar@9: info->q = q->j; alpar@9: info->bnd = q->ub; alpar@9: /* substitute x[q] into objective row */ alpar@9: npp->c0 += q->coef * q->ub; alpar@9: q->coef = -q->coef; alpar@9: /* substitute x[q] into constraint rows */ alpar@9: for (aij = q->ptr; aij != NULL; aij = aij->c_next) alpar@9: { i = aij->row; alpar@9: if (i->lb == i->ub) alpar@9: i->ub = (i->lb -= aij->val * q->ub); alpar@9: else alpar@9: { if (i->lb != -DBL_MAX) alpar@9: i->lb -= aij->val * q->ub; alpar@9: if (i->ub != +DBL_MAX) alpar@9: i->ub -= aij->val * q->ub; alpar@9: } alpar@9: aij->val = -aij->val; alpar@9: } alpar@9: /* column x[q] becomes column s */ alpar@9: if (q->lb != -DBL_MAX) alpar@9: q->ub -= q->lb; alpar@9: else alpar@9: q->ub = +DBL_MAX; alpar@9: q->lb = 0.0; alpar@9: return; alpar@9: } alpar@9: alpar@9: static int rcv_ubnd_col(NPP *npp, void *_info) alpar@9: { /* recover column with upper bound */ alpar@9: struct bnd_col *info = _info; alpar@9: if (npp->sol == GLP_BS) alpar@9: { if (npp->c_stat[info->q] == GLP_BS) alpar@9: npp->c_stat[info->q] = GLP_BS; alpar@9: else if (npp->c_stat[info->q] == GLP_NL) alpar@9: npp->c_stat[info->q] = GLP_NU; alpar@9: else if (npp->c_stat[info->q] == GLP_NU) alpar@9: npp->c_stat[info->q] = GLP_NL; alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: /* compute value of x[q] with formula (2) */ alpar@9: npp->c_value[info->q] = info->bnd - npp->c_value[info->q]; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_dbnd_col - process non-negative column with upper bound alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_dbnd_col(NPP *npp, NPPCOL *q); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_dbnd_col processes column q, which is non-negative alpar@9: * and has upper bound: alpar@9: * alpar@9: * 0 <= x[q] <= u[q], (1) alpar@9: * alpar@9: * where u[q] > 0. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Upper bound of column q can be replaced by the following equality alpar@9: * constraint: alpar@9: * alpar@9: * x[q] + s = u[q], (2) alpar@9: * alpar@9: * where s >= 0 is a non-negative complement variable. alpar@9: * alpar@9: * Since in the primal system along with new row (2) there appears a alpar@9: * new column s having the only non-zero coefficient in this row, in alpar@9: * the dual system there appears a new row: alpar@9: * alpar@9: * (+1)pi + lambda[s] = 0, (3) alpar@9: * alpar@9: * where (+1) is coefficient at column s in row (2), pi is multiplier alpar@9: * for row (2), lambda[s] is multiplier for column s, 0 is coefficient alpar@9: * at column s in the objective row. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * Status of column q in solution to the original problem is determined alpar@9: * by its status and status of column s in solution to the transformed alpar@9: * problem as follows: alpar@9: * alpar@9: * +-----------------------------------+------------------+ alpar@9: * | Transformed problem | Original problem | alpar@9: * +-----------------+-----------------+------------------+ alpar@9: * | Status of col q | Status of col s | Status of col q | alpar@9: * +-----------------+-----------------+------------------+ alpar@9: * | GLP_BS | GLP_BS | GLP_BS | alpar@9: * | GLP_BS | GLP_NL | GLP_NU | alpar@9: * | GLP_NL | GLP_BS | GLP_NL | alpar@9: * | GLP_NL | GLP_NL | GLP_NL (*) | alpar@9: * +-----------------+-----------------+------------------+ alpar@9: * alpar@9: * Value of column q in solution to the original problem is the same as alpar@9: * in solution to the transformed problem. alpar@9: * alpar@9: * 1. Formally, in solution to the transformed problem columns q and s alpar@9: * cannot be non-basic at the same time, since the constraint (2) alpar@9: * would be violated. However, if u[q] is close to zero, violation alpar@9: * may be less than a working precision even if both columns q and s alpar@9: * are non-basic. In this degenerate case row (2) can be only basic, alpar@9: * i.e. non-active constraint (otherwise corresponding row of the alpar@9: * basis matrix would be zero). This allows to pivot out auxiliary alpar@9: * variable and pivot in column s, in which case the row becomes alpar@9: * active while column s becomes basic. alpar@9: * alpar@9: * 2. If column q is integral, column s is also integral. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Value of column q in solution to the original problem is the same as alpar@9: * in solution to the transformed problem. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * Value of column q in solution to the original problem is the same as alpar@9: * in solution to the transformed problem. */ alpar@9: alpar@9: struct dbnd_col alpar@9: { /* double-bounded column */ alpar@9: int q; alpar@9: /* column reference number for variable x[q] */ alpar@9: int s; alpar@9: /* column reference number for complement variable s */ alpar@9: }; alpar@9: alpar@9: static int rcv_dbnd_col(NPP *npp, void *info); alpar@9: alpar@9: void npp_dbnd_col(NPP *npp, NPPCOL *q) alpar@9: { /* process non-negative column with upper bound */ alpar@9: struct dbnd_col *info; alpar@9: NPPROW *p; alpar@9: NPPCOL *s; alpar@9: /* the column must be non-negative with upper bound */ alpar@9: xassert(q->lb == 0.0); alpar@9: xassert(q->ub > 0.0); alpar@9: xassert(q->ub != +DBL_MAX); alpar@9: /* create variable s */ alpar@9: s = npp_add_col(npp); alpar@9: s->is_int = q->is_int; alpar@9: s->lb = 0.0, s->ub = +DBL_MAX; alpar@9: /* create equality constraint (2) */ alpar@9: p = npp_add_row(npp); alpar@9: p->lb = p->ub = q->ub; alpar@9: npp_add_aij(npp, p, q, +1.0); alpar@9: npp_add_aij(npp, p, s, +1.0); alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_dbnd_col, sizeof(struct dbnd_col)); alpar@9: info->q = q->j; alpar@9: info->s = s->j; alpar@9: /* remove upper bound of x[q] */ alpar@9: q->ub = +DBL_MAX; alpar@9: return; alpar@9: } alpar@9: alpar@9: static int rcv_dbnd_col(NPP *npp, void *_info) alpar@9: { /* recover non-negative column with upper bound */ alpar@9: struct dbnd_col *info = _info; alpar@9: if (npp->sol == GLP_BS) alpar@9: { if (npp->c_stat[info->q] == GLP_BS) alpar@9: { if (npp->c_stat[info->s] == GLP_BS) alpar@9: npp->c_stat[info->q] = GLP_BS; alpar@9: else if (npp->c_stat[info->s] == GLP_NL) alpar@9: npp->c_stat[info->q] = GLP_NU; alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: else if (npp->c_stat[info->q] == GLP_NL) alpar@9: { if (npp->c_stat[info->s] == GLP_BS || alpar@9: npp->c_stat[info->s] == GLP_NL) alpar@9: npp->c_stat[info->q] = GLP_NL; alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_fixed_col - process fixed column alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_fixed_col(NPP *npp, NPPCOL *q); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_fixed_col processes column q, which is fixed: alpar@9: * alpar@9: * x[q] = s[q], (1) alpar@9: * alpar@9: * where s[q] is a fixed column value. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * The value of a fixed column can be substituted into the objective alpar@9: * and constraint rows that allows removing the column from the problem. alpar@9: * alpar@9: * Substituting x[q] = s[q] into the objective row, we have: alpar@9: * alpar@9: * z = sum c[j] x[j] + c0 = alpar@9: * j alpar@9: * alpar@9: * = sum c[j] x[j] + c[q] x[q] + c0 = alpar@9: * j!=q alpar@9: * alpar@9: * = sum c[j] x[j] + c[q] s[q] + c0 = alpar@9: * j!=q alpar@9: * alpar@9: * = sum c[j] x[j] + c~0, alpar@9: * j!=q alpar@9: * alpar@9: * where alpar@9: * alpar@9: * c~0 = c0 + c[q] s[q] (2) alpar@9: * alpar@9: * is the constant term of the objective in the transformed problem. alpar@9: * Similarly, substituting x[q] = s[q] into constraint row i, we have: alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] <= U[i] ==> alpar@9: * j alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==> alpar@9: * j!=q alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] + a[i,q] s[q] <= U[i] ==> alpar@9: * j!=q alpar@9: * alpar@9: * L~[i] <= sum a[i,j] x[j] + a[i,q] s <= U~[i], alpar@9: * j!=q alpar@9: * alpar@9: * where alpar@9: * alpar@9: * L~[i] = L[i] - a[i,q] s[q], U~[i] = U[i] - a[i,q] s[q] (3) alpar@9: * alpar@9: * are lower and upper bounds of row i in the transformed problem, alpar@9: * resp. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * Column q is assigned status GLP_NS and its value is assigned s[q]. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Value of column q is assigned s[q]. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * Value of column q is assigned s[q]. */ alpar@9: alpar@9: struct fixed_col alpar@9: { /* fixed column */ alpar@9: int q; alpar@9: /* column reference number for variable x[q] */ alpar@9: double s; alpar@9: /* value, at which x[q] is fixed */ alpar@9: }; alpar@9: alpar@9: static int rcv_fixed_col(NPP *npp, void *info); alpar@9: alpar@9: void npp_fixed_col(NPP *npp, NPPCOL *q) alpar@9: { /* process fixed column */ alpar@9: struct fixed_col *info; alpar@9: NPPROW *i; alpar@9: NPPAIJ *aij; alpar@9: /* the column must be fixed */ alpar@9: xassert(q->lb == q->ub); alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_fixed_col, sizeof(struct fixed_col)); alpar@9: info->q = q->j; alpar@9: info->s = q->lb; alpar@9: /* substitute x[q] = s[q] into objective row */ alpar@9: npp->c0 += q->coef * q->lb; alpar@9: /* substitute x[q] = s[q] into constraint rows */ alpar@9: for (aij = q->ptr; aij != NULL; aij = aij->c_next) alpar@9: { i = aij->row; alpar@9: if (i->lb == i->ub) alpar@9: i->ub = (i->lb -= aij->val * q->lb); alpar@9: else alpar@9: { if (i->lb != -DBL_MAX) alpar@9: i->lb -= aij->val * q->lb; alpar@9: if (i->ub != +DBL_MAX) alpar@9: i->ub -= aij->val * q->lb; alpar@9: } alpar@9: } alpar@9: /* remove the column from the problem */ alpar@9: npp_del_col(npp, q); alpar@9: return; alpar@9: } alpar@9: alpar@9: static int rcv_fixed_col(NPP *npp, void *_info) alpar@9: { /* recover fixed column */ alpar@9: struct fixed_col *info = _info; alpar@9: if (npp->sol == GLP_SOL) alpar@9: npp->c_stat[info->q] = GLP_NS; alpar@9: npp->c_value[info->q] = info->s; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_make_equality - process row with almost identical bounds alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_make_equality(NPP *npp, NPPROW *p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_make_equality processes row p: alpar@9: * alpar@9: * L[p] <= sum a[p,j] x[j] <= U[p], (1) alpar@9: * j alpar@9: * alpar@9: * where -oo < L[p] < U[p] < +oo, i.e. which is double-sided inequality alpar@9: * constraint. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - row bounds have not been changed; alpar@9: * alpar@9: * 1 - row has been replaced by equality constraint. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * If bounds of row (1) are very close to each other: alpar@9: * alpar@9: * U[p] - L[p] <= eps, (2) alpar@9: * alpar@9: * where eps is an absolute tolerance for row value, the row can be alpar@9: * replaced by the following almost equivalent equiality constraint: alpar@9: * alpar@9: * sum a[p,j] x[j] = b, (3) alpar@9: * j alpar@9: * alpar@9: * where b = (L[p] + U[p]) / 2. If the right-hand side in (3) happens alpar@9: * to be very close to its nearest integer: alpar@9: * alpar@9: * |b - floor(b + 0.5)| <= eps, (4) alpar@9: * alpar@9: * it is reasonable to use this nearest integer as the right-hand side. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * Status of row p in solution to the original problem is determined alpar@9: * by its status and the sign of its multiplier pi[p] in solution to alpar@9: * the transformed problem as follows: alpar@9: * alpar@9: * +-----------------------+---------+--------------------+ alpar@9: * | Status of row p | Sign of | Status of row p | alpar@9: * | (transformed problem) | pi[p] | (original problem) | alpar@9: * +-----------------------+---------+--------------------+ alpar@9: * | GLP_BS | + / - | GLP_BS | alpar@9: * | GLP_NS | + | GLP_NL | alpar@9: * | GLP_NS | - | GLP_NU | alpar@9: * +-----------------------+---------+--------------------+ alpar@9: * alpar@9: * Value of row multiplier pi[p] in solution to the original problem is alpar@9: * the same as in solution to the transformed problem. alpar@9: * alpar@9: * RECOVERING INTERIOR POINT SOLUTION alpar@9: * alpar@9: * Value of row multiplier pi[p] in solution to the original problem is alpar@9: * the same as in solution to the transformed problem. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: struct make_equality alpar@9: { /* row with almost identical bounds */ alpar@9: int p; alpar@9: /* row reference number */ alpar@9: }; alpar@9: alpar@9: static int rcv_make_equality(NPP *npp, void *info); alpar@9: alpar@9: int npp_make_equality(NPP *npp, NPPROW *p) alpar@9: { /* process row with almost identical bounds */ alpar@9: struct make_equality *info; alpar@9: double b, eps, nint; alpar@9: /* the row must be double-sided inequality */ alpar@9: xassert(p->lb != -DBL_MAX); alpar@9: xassert(p->ub != +DBL_MAX); alpar@9: xassert(p->lb < p->ub); alpar@9: /* check row bounds */ alpar@9: eps = 1e-9 + 1e-12 * fabs(p->lb); alpar@9: if (p->ub - p->lb > eps) return 0; alpar@9: /* row bounds are very close to each other */ alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_make_equality, sizeof(struct make_equality)); alpar@9: info->p = p->i; alpar@9: /* compute right-hand side */ alpar@9: b = 0.5 * (p->ub + p->lb); alpar@9: nint = floor(b + 0.5); alpar@9: if (fabs(b - nint) <= eps) b = nint; alpar@9: /* replace row p by almost equivalent equality constraint */ alpar@9: p->lb = p->ub = b; alpar@9: return 1; alpar@9: } alpar@9: alpar@9: int rcv_make_equality(NPP *npp, void *_info) alpar@9: { /* recover row with almost identical bounds */ alpar@9: struct make_equality *info = _info; alpar@9: if (npp->sol == GLP_SOL) alpar@9: { if (npp->r_stat[info->p] == GLP_BS) alpar@9: npp->r_stat[info->p] = GLP_BS; alpar@9: else if (npp->r_stat[info->p] == GLP_NS) alpar@9: { if (npp->r_pi[info->p] >= 0.0) alpar@9: npp->r_stat[info->p] = GLP_NL; alpar@9: else alpar@9: npp->r_stat[info->p] = GLP_NU; alpar@9: } alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_make_fixed - process column with almost identical bounds alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_make_fixed(NPP *npp, NPPCOL *q); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_make_fixed processes column q: alpar@9: * alpar@9: * l[q] <= x[q] <= u[q], (1) alpar@9: * alpar@9: * where -oo < l[q] < u[q] < +oo, i.e. which has both lower and upper alpar@9: * bounds. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - column bounds have not been changed; alpar@9: * alpar@9: * 1 - column has been fixed. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * If bounds of column (1) are very close to each other: alpar@9: * alpar@9: * u[q] - l[q] <= eps, (2) alpar@9: * alpar@9: * where eps is an absolute tolerance for column value, the column can alpar@9: * be fixed: alpar@9: * alpar@9: * x[q] = s[q], (3) alpar@9: * alpar@9: * where s[q] = (l[q] + u[q]) / 2. And if the fixed column value s[q] alpar@9: * happens to be very close to its nearest integer: alpar@9: * alpar@9: * |s[q] - floor(s[q] + 0.5)| <= eps, (4) alpar@9: * alpar@9: * it is reasonable to use this nearest integer as the fixed value. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * In the dual system of the original (as well as transformed) problem alpar@9: * column q corresponds to the following row: alpar@9: * alpar@9: * sum a[i,q] pi[i] + lambda[q] = c[q]. (5) alpar@9: * i alpar@9: * alpar@9: * Since multipliers pi[i] are known for all rows from solution to the alpar@9: * transformed problem, formula (5) allows computing value of multiplier alpar@9: * (reduced cost) for column q: alpar@9: * alpar@9: * lambda[q] = c[q] - sum a[i,q] pi[i]. (6) alpar@9: * i alpar@9: * alpar@9: * Status of column q in solution to the original problem is determined alpar@9: * by its status and the sign of its multiplier lambda[q] in solution to alpar@9: * the transformed problem as follows: alpar@9: * alpar@9: * +-----------------------+-----------+--------------------+ alpar@9: * | Status of column q | Sign of | Status of column q | alpar@9: * | (transformed problem) | lambda[q] | (original problem) | alpar@9: * +-----------------------+-----------+--------------------+ alpar@9: * | GLP_BS | + / - | GLP_BS | alpar@9: * | GLP_NS | + | GLP_NL | alpar@9: * | GLP_NS | - | GLP_NU | alpar@9: * +-----------------------+-----------+--------------------+ alpar@9: * alpar@9: * Value of column q in solution to the original problem is the same as alpar@9: * in solution to the transformed problem. alpar@9: * alpar@9: * RECOVERING INTERIOR POINT SOLUTION alpar@9: * alpar@9: * Value of column q in solution to the original problem is the same as alpar@9: * in solution to the transformed problem. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: struct make_fixed alpar@9: { /* column with almost identical bounds */ alpar@9: int q; alpar@9: /* column reference number */ alpar@9: double c; alpar@9: /* objective coefficient at x[q] */ alpar@9: NPPLFE *ptr; alpar@9: /* list of non-zero coefficients a[i,q] */ alpar@9: }; alpar@9: alpar@9: static int rcv_make_fixed(NPP *npp, void *info); alpar@9: alpar@9: int npp_make_fixed(NPP *npp, NPPCOL *q) alpar@9: { /* process column with almost identical bounds */ alpar@9: struct make_fixed *info; alpar@9: NPPAIJ *aij; alpar@9: NPPLFE *lfe; alpar@9: double s, eps, nint; alpar@9: /* the column must be double-bounded */ alpar@9: xassert(q->lb != -DBL_MAX); alpar@9: xassert(q->ub != +DBL_MAX); alpar@9: xassert(q->lb < q->ub); alpar@9: /* check column bounds */ alpar@9: eps = 1e-9 + 1e-12 * fabs(q->lb); alpar@9: if (q->ub - q->lb > eps) return 0; alpar@9: /* column bounds are very close to each other */ alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_make_fixed, sizeof(struct make_fixed)); alpar@9: info->q = q->j; alpar@9: info->c = q->coef; alpar@9: info->ptr = NULL; alpar@9: /* save column coefficients a[i,q] (needed for basic solution alpar@9: only) */ alpar@9: if (npp->sol == GLP_SOL) alpar@9: { for (aij = q->ptr; aij != NULL; aij = aij->c_next) alpar@9: { lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE)); alpar@9: lfe->ref = aij->row->i; alpar@9: lfe->val = aij->val; alpar@9: lfe->next = info->ptr; alpar@9: info->ptr = lfe; alpar@9: } alpar@9: } alpar@9: /* compute column fixed value */ alpar@9: s = 0.5 * (q->ub + q->lb); alpar@9: nint = floor(s + 0.5); alpar@9: if (fabs(s - nint) <= eps) s = nint; alpar@9: /* make column q fixed */ alpar@9: q->lb = q->ub = s; alpar@9: return 1; alpar@9: } alpar@9: alpar@9: static int rcv_make_fixed(NPP *npp, void *_info) alpar@9: { /* recover column with almost identical bounds */ alpar@9: struct make_fixed *info = _info; alpar@9: NPPLFE *lfe; alpar@9: double lambda; alpar@9: if (npp->sol == GLP_SOL) alpar@9: { if (npp->c_stat[info->q] == GLP_BS) alpar@9: npp->c_stat[info->q] = GLP_BS; alpar@9: else if (npp->c_stat[info->q] == GLP_NS) alpar@9: { /* compute multiplier for column q with formula (6) */ alpar@9: lambda = info->c; alpar@9: for (lfe = info->ptr; lfe != NULL; lfe = lfe->next) alpar@9: lambda -= lfe->val * npp->r_pi[lfe->ref]; alpar@9: /* assign status to non-basic column */ alpar@9: if (lambda >= 0.0) alpar@9: npp->c_stat[info->q] = GLP_NL; alpar@9: else alpar@9: npp->c_stat[info->q] = GLP_NU; alpar@9: } alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /* eof */