alpar@9: /* glpnpp03.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_empty_row - process empty row alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_empty_row(NPP *npp, NPPROW *p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_empty_row processes row p, which is empty, i.e. alpar@9: * coefficients at all columns in this row are zero: alpar@9: * alpar@9: * L[p] <= sum 0 x[j] <= U[p], (1) alpar@9: * alpar@9: * where L[p] <= U[p]. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - success; alpar@9: * alpar@9: * 1 - problem has no primal feasible solution. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * If the following conditions hold: alpar@9: * alpar@9: * L[p] <= +eps, U[p] >= -eps, (2) alpar@9: * alpar@9: * where eps is an absolute tolerance for row value, the row p is alpar@9: * redundant. In this case it can be replaced by equivalent redundant alpar@9: * row, which is free (unbounded), and then removed from the problem. alpar@9: * Otherwise, the row p is infeasible and, thus, the problem has no alpar@9: * primal feasible solution. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * See the routine npp_free_row. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * See the routine npp_free_row. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: int npp_empty_row(NPP *npp, NPPROW *p) alpar@9: { /* process empty row */ alpar@9: double eps = 1e-3; alpar@9: /* the row must be empty */ alpar@9: xassert(p->ptr == NULL); alpar@9: /* check primal feasibility */ alpar@9: if (p->lb > +eps || p->ub < -eps) alpar@9: return 1; alpar@9: /* replace the row by equivalent free (unbounded) row */ alpar@9: p->lb = -DBL_MAX, p->ub = +DBL_MAX; alpar@9: /* and process it */ alpar@9: npp_free_row(npp, p); alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_empty_col - process empty column alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_empty_col(NPP *npp, NPPCOL *q); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_empty_col processes column q: alpar@9: * alpar@9: * l[q] <= x[q] <= u[q], (1) alpar@9: * alpar@9: * where l[q] <= u[q], which is empty, i.e. has zero coefficients in alpar@9: * all constraint rows. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - success; alpar@9: * alpar@9: * 1 - problem has no dual feasible solution. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * The row of the dual system corresponding to the empty column is the alpar@9: * following: alpar@9: * alpar@9: * sum 0 pi[i] + lambda[q] = c[q], (2) alpar@9: * i alpar@9: * alpar@9: * from which it follows that: alpar@9: * alpar@9: * lambda[q] = c[q]. (3) alpar@9: * alpar@9: * If the following condition holds: alpar@9: * alpar@9: * c[q] < - eps, (4) alpar@9: * alpar@9: * where eps is an absolute tolerance for column multiplier, the lower alpar@9: * column bound l[q] must be active to provide dual feasibility (note alpar@9: * that being preprocessed the problem is always minimization). In this alpar@9: * case the column can be fixed on its lower bound and removed from the alpar@9: * problem (if the column is integral, its bounds are also assumed to alpar@9: * be integral). And if the column has no lower bound (l[q] = -oo), the alpar@9: * problem has no dual feasible solution. alpar@9: * alpar@9: * If the following condition holds: alpar@9: * alpar@9: * c[q] > + eps, (5) alpar@9: * alpar@9: * the upper column bound u[q] must be active to provide dual alpar@9: * feasibility. In this case the column can be fixed on its upper bound alpar@9: * and removed from the problem. And if the column has no upper bound alpar@9: * (u[q] = +oo), the problem has no dual feasible solution. alpar@9: * alpar@9: * Finally, if the following condition holds: alpar@9: * alpar@9: * - eps <= c[q] <= +eps, (6) alpar@9: * alpar@9: * dual feasibility does not depend on a particular value of column q. alpar@9: * In this case the column can be fixed either on its lower bound (if alpar@9: * l[q] > -oo) or on its upper bound (if u[q] < +oo) or at zero (if the alpar@9: * column is unbounded) and then removed from the problem. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * See the routine npp_fixed_col. Having been recovered the column alpar@9: * is assigned status GLP_NS. However, if actually it is not fixed alpar@9: * (l[q] < u[q]), its status should be changed to GLP_NL, GLP_NU, or alpar@9: * GLP_NF depending on which bound it was fixed on transformation stage. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * See the routine npp_fixed_col. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * See the routine npp_fixed_col. */ alpar@9: alpar@9: struct empty_col alpar@9: { /* empty column */ alpar@9: int q; alpar@9: /* column reference number */ alpar@9: char stat; alpar@9: /* status in basic solution */ alpar@9: }; alpar@9: alpar@9: static int rcv_empty_col(NPP *npp, void *info); alpar@9: alpar@9: int npp_empty_col(NPP *npp, NPPCOL *q) alpar@9: { /* process empty column */ alpar@9: struct empty_col *info; alpar@9: double eps = 1e-3; alpar@9: /* the column must be empty */ alpar@9: xassert(q->ptr == NULL); alpar@9: /* check dual feasibility */ alpar@9: if (q->coef > +eps && q->lb == -DBL_MAX) alpar@9: return 1; alpar@9: if (q->coef < -eps && q->ub == +DBL_MAX) alpar@9: return 1; alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_empty_col, sizeof(struct empty_col)); alpar@9: info->q = q->j; alpar@9: /* fix the column */ alpar@9: if (q->lb == -DBL_MAX && q->ub == +DBL_MAX) alpar@9: { /* free column */ alpar@9: info->stat = GLP_NF; alpar@9: q->lb = q->ub = 0.0; alpar@9: } alpar@9: else if (q->ub == +DBL_MAX) alpar@9: lo: { /* column with lower bound */ alpar@9: info->stat = GLP_NL; alpar@9: q->ub = q->lb; alpar@9: } alpar@9: else if (q->lb == -DBL_MAX) alpar@9: up: { /* column with upper bound */ alpar@9: info->stat = GLP_NU; alpar@9: q->lb = q->ub; alpar@9: } alpar@9: else if (q->lb != q->ub) alpar@9: { /* double-bounded column */ alpar@9: if (q->coef >= +DBL_EPSILON) goto lo; alpar@9: if (q->coef <= -DBL_EPSILON) goto up; alpar@9: if (fabs(q->lb) <= fabs(q->ub)) goto lo; else goto up; alpar@9: } alpar@9: else alpar@9: { /* fixed column */ alpar@9: info->stat = GLP_NS; alpar@9: } alpar@9: /* process fixed column */ alpar@9: npp_fixed_col(npp, q); alpar@9: return 0; alpar@9: } alpar@9: alpar@9: static int rcv_empty_col(NPP *npp, void *_info) alpar@9: { /* recover empty column */ alpar@9: struct empty_col *info = _info; alpar@9: if (npp->sol == GLP_SOL) alpar@9: npp->c_stat[info->q] = info->stat; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_implied_value - process implied column value alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_implied_value(NPP *npp, NPPCOL *q, double s); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * For column q: alpar@9: * alpar@9: * l[q] <= x[q] <= u[q], (1) alpar@9: * alpar@9: * where l[q] < u[q], the routine npp_implied_value processes its alpar@9: * implied value s[q]. If this implied value satisfies to the current alpar@9: * column bounds and integrality condition, the routine fixes column q alpar@9: * at the given point. Note that the column is kept in the problem in alpar@9: * any case. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - column has been fixed; alpar@9: * alpar@9: * 1 - implied value violates to current column bounds; alpar@9: * alpar@9: * 2 - implied value violates integrality condition. alpar@9: * alpar@9: * ALGORITHM alpar@9: * alpar@9: * Implied column value s[q] satisfies to the current column bounds if alpar@9: * the following condition holds: alpar@9: * alpar@9: * l[q] - eps <= s[q] <= u[q] + eps, (2) alpar@9: * alpar@9: * where eps is an absolute tolerance for column value. If the column alpar@9: * is integral, the following condition also must hold: alpar@9: * alpar@9: * |s[q] - floor(s[q]+0.5)| <= eps, (3) alpar@9: * alpar@9: * where floor(s[q]+0.5) is the nearest integer to s[q]. alpar@9: * alpar@9: * If both condition (2) and (3) are satisfied, the column can be fixed alpar@9: * at the value s[q], or, if it is integral, at floor(s[q]+0.5). alpar@9: * Otherwise, if s[q] violates (2) or (3), the problem has no feasible alpar@9: * solution. alpar@9: * alpar@9: * Note: If s[q] is close to l[q] or u[q], it seems to be reasonable to alpar@9: * fix the column at its lower or upper bound, resp. rather than at the alpar@9: * implied value. */ alpar@9: alpar@9: int npp_implied_value(NPP *npp, NPPCOL *q, double s) alpar@9: { /* process implied column value */ alpar@9: double eps, nint; alpar@9: xassert(npp == npp); alpar@9: /* column must not be fixed */ alpar@9: xassert(q->lb < q->ub); alpar@9: /* check integrality */ alpar@9: if (q->is_int) alpar@9: { nint = floor(s + 0.5); alpar@9: if (fabs(s - nint) <= 1e-5) alpar@9: s = nint; alpar@9: else alpar@9: return 2; alpar@9: } alpar@9: /* check current column lower bound */ alpar@9: if (q->lb != -DBL_MAX) alpar@9: { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb)); alpar@9: if (s < q->lb - eps) return 1; alpar@9: /* if s[q] is close to l[q], fix column at its lower bound alpar@9: rather than at the implied value */ alpar@9: if (s < q->lb + 1e-3 * eps) alpar@9: { q->ub = q->lb; alpar@9: return 0; alpar@9: } alpar@9: } alpar@9: /* check current column upper bound */ alpar@9: if (q->ub != +DBL_MAX) alpar@9: { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub)); alpar@9: if (s > q->ub + eps) return 1; alpar@9: /* if s[q] is close to u[q], fix column at its upper bound alpar@9: rather than at the implied value */ alpar@9: if (s > q->ub - 1e-3 * eps) alpar@9: { q->lb = q->ub; alpar@9: return 0; alpar@9: } alpar@9: } alpar@9: /* fix column at the implied value */ alpar@9: q->lb = q->ub = s; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_eq_singlet - process row singleton (equality constraint) alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_eq_singlet(NPP *npp, NPPROW *p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_eq_singlet processes row p, which is equiality alpar@9: * constraint having the only non-zero coefficient: alpar@9: * alpar@9: * a[p,q] x[q] = b. (1) alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - success; alpar@9: * alpar@9: * 1 - problem has no primal feasible solution; alpar@9: * alpar@9: * 2 - problem has no integer feasible solution. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * The equality constraint defines implied value of column q: alpar@9: * alpar@9: * x[q] = s[q] = b / a[p,q]. (2) alpar@9: * alpar@9: * If the implied value s[q] satisfies to the column bounds (see the alpar@9: * routine npp_implied_value), the column can be fixed at s[q] and alpar@9: * removed from the problem. In this case row p becomes redundant, so alpar@9: * it can be replaced by equivalent free row and also removed from the alpar@9: * problem. alpar@9: * alpar@9: * Note that the routine removes from the problem only row p. Column q alpar@9: * becomes fixed, however, it is kept in the problem. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * In solution to the original problem row p is assigned status GLP_NS alpar@9: * (active equality constraint), and column q is assigned status GLP_BS alpar@9: * (basic column). alpar@9: * alpar@9: * Multiplier for row p can be computed as follows. In the dual system alpar@9: * of the original problem column q corresponds to the following row: alpar@9: * alpar@9: * sum a[i,q] pi[i] + lambda[q] = c[q] ==> alpar@9: * i alpar@9: * alpar@9: * sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q]. alpar@9: * i!=p alpar@9: * alpar@9: * Therefore: alpar@9: * alpar@9: * 1 alpar@9: * pi[p] = ------ (c[q] - lambda[q] - sum a[i,q] pi[i]), (3) alpar@9: * a[p,q] i!=q alpar@9: * alpar@9: * where lambda[q] = 0 (since column[q] is basic), and pi[i] for all alpar@9: * i != p are known in solution to the transformed problem. alpar@9: * alpar@9: * Value of column q in solution to the original problem is assigned alpar@9: * its implied value s[q]. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Multiplier for row p is computed with formula (3). Value of column alpar@9: * q is assigned its implied value s[q]. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * Value of column q is assigned its implied value s[q]. */ alpar@9: alpar@9: struct eq_singlet alpar@9: { /* row singleton (equality constraint) */ alpar@9: int p; alpar@9: /* row reference number */ alpar@9: int q; alpar@9: /* column reference number */ alpar@9: double apq; alpar@9: /* constraint coefficient a[p,q] */ 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], i != p */ alpar@9: }; alpar@9: alpar@9: static int rcv_eq_singlet(NPP *npp, void *info); alpar@9: alpar@9: int npp_eq_singlet(NPP *npp, NPPROW *p) alpar@9: { /* process row singleton (equality constraint) */ alpar@9: struct eq_singlet *info; alpar@9: NPPCOL *q; alpar@9: NPPAIJ *aij; alpar@9: NPPLFE *lfe; alpar@9: int ret; alpar@9: double s; alpar@9: /* the row must be singleton equality constraint */ alpar@9: xassert(p->lb == p->ub); alpar@9: xassert(p->ptr != NULL && p->ptr->r_next == NULL); alpar@9: /* compute and process implied column value */ alpar@9: aij = p->ptr; alpar@9: q = aij->col; alpar@9: s = p->lb / aij->val; alpar@9: ret = npp_implied_value(npp, q, s); alpar@9: xassert(0 <= ret && ret <= 2); alpar@9: if (ret != 0) return ret; alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_eq_singlet, sizeof(struct eq_singlet)); alpar@9: info->p = p->i; alpar@9: info->q = q->j; alpar@9: info->apq = aij->val; alpar@9: info->c = q->coef; alpar@9: info->ptr = NULL; alpar@9: /* save column coefficients a[i,q], i != p (not needed for MIP alpar@9: solution) */ alpar@9: if (npp->sol != GLP_MIP) alpar@9: { for (aij = q->ptr; aij != NULL; aij = aij->c_next) alpar@9: { if (aij->row == p) continue; /* skip a[p,q] */ 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: /* remove the row from the problem */ alpar@9: npp_del_row(npp, p); alpar@9: return 0; alpar@9: } alpar@9: alpar@9: static int rcv_eq_singlet(NPP *npp, void *_info) alpar@9: { /* recover row singleton (equality constraint) */ alpar@9: struct eq_singlet *info = _info; alpar@9: NPPLFE *lfe; alpar@9: double temp; alpar@9: if (npp->sol == GLP_SOL) alpar@9: { /* column q must be already recovered as GLP_NS */ alpar@9: if (npp->c_stat[info->q] != GLP_NS) alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: npp->r_stat[info->p] = GLP_NS; alpar@9: npp->c_stat[info->q] = GLP_BS; alpar@9: } alpar@9: if (npp->sol != GLP_MIP) alpar@9: { /* compute multiplier for row p with formula (3) */ alpar@9: temp = info->c; alpar@9: for (lfe = info->ptr; lfe != NULL; lfe = lfe->next) alpar@9: temp -= lfe->val * npp->r_pi[lfe->ref]; alpar@9: npp->r_pi[info->p] = temp / info->apq; alpar@9: } alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_implied_lower - process implied column lower bound alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_implied_lower(NPP *npp, NPPCOL *q, double l); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * For column q: alpar@9: * alpar@9: * l[q] <= x[q] <= u[q], (1) alpar@9: * alpar@9: * where l[q] < u[q], the routine npp_implied_lower processes its alpar@9: * implied lower bound l'[q]. As the result the current column lower alpar@9: * bound may increase. Note that the column is kept in the problem in alpar@9: * any case. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - current column lower bound has not changed; alpar@9: * alpar@9: * 1 - current column lower bound has changed, but not significantly; alpar@9: * alpar@9: * 2 - current column lower bound has significantly changed; alpar@9: * alpar@9: * 3 - column has been fixed on its upper bound; alpar@9: * alpar@9: * 4 - implied lower bound violates current column upper bound. alpar@9: * alpar@9: * ALGORITHM alpar@9: * alpar@9: * If column q is integral, before processing its implied lower bound alpar@9: * should be rounded up: alpar@9: * alpar@9: * ( floor(l'[q]+0.5), if |l'[q] - floor(l'[q]+0.5)| <= eps alpar@9: * l'[q] := < (2) alpar@9: * ( ceil(l'[q]), otherwise alpar@9: * alpar@9: * where floor(l'[q]+0.5) is the nearest integer to l'[q], ceil(l'[q]) alpar@9: * is smallest integer not less than l'[q], and eps is an absolute alpar@9: * tolerance for column value. alpar@9: * alpar@9: * Processing implied column lower bound l'[q] includes the following alpar@9: * cases: alpar@9: * alpar@9: * 1) if l'[q] < l[q] + eps, implied lower bound is redundant; alpar@9: * alpar@9: * 2) if l[q] + eps <= l[q] <= u[q] + eps, current column lower bound alpar@9: * l[q] can be strengthened by replacing it with l'[q]. If in this alpar@9: * case new column lower bound becomes close to current column upper alpar@9: * bound u[q], the column can be fixed on its upper bound; alpar@9: * alpar@9: * 3) if l'[q] > u[q] + eps, implied lower bound violates current alpar@9: * column upper bound u[q], in which case the problem has no primal alpar@9: * feasible solution. */ alpar@9: alpar@9: int npp_implied_lower(NPP *npp, NPPCOL *q, double l) alpar@9: { /* process implied column lower bound */ alpar@9: int ret; alpar@9: double eps, nint; alpar@9: xassert(npp == npp); alpar@9: /* column must not be fixed */ alpar@9: xassert(q->lb < q->ub); alpar@9: /* implied lower bound must be finite */ alpar@9: xassert(l != -DBL_MAX); alpar@9: /* if column is integral, round up l'[q] */ alpar@9: if (q->is_int) alpar@9: { nint = floor(l + 0.5); alpar@9: if (fabs(l - nint) <= 1e-5) alpar@9: l = nint; alpar@9: else alpar@9: l = ceil(l); alpar@9: } alpar@9: /* check current column lower bound */ alpar@9: if (q->lb != -DBL_MAX) alpar@9: { eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->lb)); alpar@9: if (l < q->lb + eps) alpar@9: { ret = 0; /* redundant */ alpar@9: goto done; alpar@9: } alpar@9: } alpar@9: /* check current column upper bound */ alpar@9: if (q->ub != +DBL_MAX) alpar@9: { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub)); alpar@9: if (l > q->ub + eps) alpar@9: { ret = 4; /* infeasible */ alpar@9: goto done; alpar@9: } alpar@9: /* if l'[q] is close to u[q], fix column at its upper bound */ alpar@9: if (l > q->ub - 1e-3 * eps) alpar@9: { q->lb = q->ub; alpar@9: ret = 3; /* fixed */ alpar@9: goto done; alpar@9: } alpar@9: } alpar@9: /* check if column lower bound changes significantly */ alpar@9: if (q->lb == -DBL_MAX) alpar@9: ret = 2; /* significantly */ alpar@9: else if (q->is_int && l > q->lb + 0.5) alpar@9: ret = 2; /* significantly */ alpar@9: else if (l > q->lb + 0.30 * (1.0 + fabs(q->lb))) alpar@9: ret = 2; /* significantly */ alpar@9: else alpar@9: ret = 1; /* not significantly */ alpar@9: /* set new column lower bound */ alpar@9: q->lb = l; alpar@9: done: return ret; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_implied_upper - process implied column upper bound alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_implied_upper(NPP *npp, NPPCOL *q, double u); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * For column q: alpar@9: * alpar@9: * l[q] <= x[q] <= u[q], (1) alpar@9: * alpar@9: * where l[q] < u[q], the routine npp_implied_upper processes its alpar@9: * implied upper bound u'[q]. As the result the current column upper alpar@9: * bound may decrease. Note that the column is kept in the problem in alpar@9: * any case. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - current column upper bound has not changed; alpar@9: * alpar@9: * 1 - current column upper bound has changed, but not significantly; alpar@9: * alpar@9: * 2 - current column upper bound has significantly changed; alpar@9: * alpar@9: * 3 - column has been fixed on its lower bound; alpar@9: * alpar@9: * 4 - implied upper bound violates current column lower bound. alpar@9: * alpar@9: * ALGORITHM alpar@9: * alpar@9: * If column q is integral, before processing its implied upper bound alpar@9: * should be rounded down: alpar@9: * alpar@9: * ( floor(u'[q]+0.5), if |u'[q] - floor(l'[q]+0.5)| <= eps alpar@9: * u'[q] := < (2) alpar@9: * ( floor(l'[q]), otherwise alpar@9: * alpar@9: * where floor(u'[q]+0.5) is the nearest integer to u'[q], alpar@9: * floor(u'[q]) is largest integer not greater than u'[q], and eps is alpar@9: * an absolute tolerance for column value. alpar@9: * alpar@9: * Processing implied column upper bound u'[q] includes the following alpar@9: * cases: alpar@9: * alpar@9: * 1) if u'[q] > u[q] - eps, implied upper bound is redundant; alpar@9: * alpar@9: * 2) if l[q] - eps <= u[q] <= u[q] - eps, current column upper bound alpar@9: * u[q] can be strengthened by replacing it with u'[q]. If in this alpar@9: * case new column upper bound becomes close to current column lower alpar@9: * bound, the column can be fixed on its lower bound; alpar@9: * alpar@9: * 3) if u'[q] < l[q] - eps, implied upper bound violates current alpar@9: * column lower bound l[q], in which case the problem has no primal alpar@9: * feasible solution. */ alpar@9: alpar@9: int npp_implied_upper(NPP *npp, NPPCOL *q, double u) alpar@9: { int ret; alpar@9: double eps, nint; alpar@9: xassert(npp == npp); alpar@9: /* column must not be fixed */ alpar@9: xassert(q->lb < q->ub); alpar@9: /* implied upper bound must be finite */ alpar@9: xassert(u != +DBL_MAX); alpar@9: /* if column is integral, round down u'[q] */ alpar@9: if (q->is_int) alpar@9: { nint = floor(u + 0.5); alpar@9: if (fabs(u - nint) <= 1e-5) alpar@9: u = nint; alpar@9: else alpar@9: u = floor(u); alpar@9: } alpar@9: /* check current column upper bound */ alpar@9: if (q->ub != +DBL_MAX) alpar@9: { eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->ub)); alpar@9: if (u > q->ub - eps) alpar@9: { ret = 0; /* redundant */ alpar@9: goto done; alpar@9: } alpar@9: } alpar@9: /* check current column lower bound */ alpar@9: if (q->lb != -DBL_MAX) alpar@9: { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb)); alpar@9: if (u < q->lb - eps) alpar@9: { ret = 4; /* infeasible */ alpar@9: goto done; alpar@9: } alpar@9: /* if u'[q] is close to l[q], fix column at its lower bound */ alpar@9: if (u < q->lb + 1e-3 * eps) alpar@9: { q->ub = q->lb; alpar@9: ret = 3; /* fixed */ alpar@9: goto done; alpar@9: } alpar@9: } alpar@9: /* check if column upper bound changes significantly */ alpar@9: if (q->ub == +DBL_MAX) alpar@9: ret = 2; /* significantly */ alpar@9: else if (q->is_int && u < q->ub - 0.5) alpar@9: ret = 2; /* significantly */ alpar@9: else if (u < q->ub - 0.30 * (1.0 + fabs(q->ub))) alpar@9: ret = 2; /* significantly */ alpar@9: else alpar@9: ret = 1; /* not significantly */ alpar@9: /* set new column upper bound */ alpar@9: q->ub = u; alpar@9: done: return ret; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_ineq_singlet - process row singleton (inequality constraint) alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_ineq_singlet(NPP *npp, NPPROW *p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_ineq_singlet processes row p, which is inequality alpar@9: * constraint having the only non-zero coefficient: alpar@9: * alpar@9: * L[p] <= a[p,q] * x[q] <= U[p], (1) alpar@9: * alpar@9: * where L[p] < U[p], L[p] > -oo and/or U[p] < +oo. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - current column bounds have not changed; alpar@9: * alpar@9: * 1 - current column bounds have changed, but not significantly; alpar@9: * alpar@9: * 2 - current column bounds have significantly changed; alpar@9: * alpar@9: * 3 - column has been fixed on its lower or upper bound; alpar@9: * alpar@9: * 4 - problem has no primal feasible solution. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Inequality constraint (1) defines implied bounds of column q: alpar@9: * alpar@9: * ( L[p] / a[p,q], if a[p,q] > 0 alpar@9: * l'[q] = < (2) alpar@9: * ( U[p] / a[p,q], if a[p,q] < 0 alpar@9: * alpar@9: * ( U[p] / a[p,q], if a[p,q] > 0 alpar@9: * u'[q] = < (3) alpar@9: * ( L[p] / a[p,q], if a[p,q] < 0 alpar@9: * alpar@9: * If these implied bounds do not violate current bounds of column q: alpar@9: * alpar@9: * l[q] <= x[q] <= u[q], (4) alpar@9: * alpar@9: * they can be used to strengthen the current column bounds: alpar@9: * alpar@9: * l[q] := max(l[q], l'[q]), (5) alpar@9: * alpar@9: * u[q] := min(u[q], u'[q]). (6) alpar@9: * alpar@9: * (See the routines npp_implied_lower and npp_implied_upper.) alpar@9: * alpar@9: * Once bounds of row p (1) have been carried over column q, the row alpar@9: * becomes redundant, so it can be replaced by equivalent free row and alpar@9: * removed from the problem. alpar@9: * alpar@9: * Note that the routine removes from the problem only row p. Column q, alpar@9: * even it has been fixed, is kept in the problem. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * Note that the row in the dual system corresponding to column q is alpar@9: * the following: alpar@9: * alpar@9: * sum a[i,q] pi[i] + lambda[q] = c[q] ==> alpar@9: * i alpar@9: * (7) alpar@9: * sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q], alpar@9: * i!=p alpar@9: * alpar@9: * where pi[i] for all i != p are known in solution to the transformed alpar@9: * problem. Row p does not exist in the transformed problem, so it has alpar@9: * zero multiplier there. This allows computing multiplier for column q alpar@9: * in solution to the transformed problem: alpar@9: * alpar@9: * lambda~[q] = c[q] - sum a[i,q] pi[i]. (8) alpar@9: * i!=p alpar@9: * alpar@9: * Let in solution to the transformed problem column q be non-basic alpar@9: * with lower bound active (GLP_NL, lambda~[q] >= 0), and this lower alpar@9: * bound be implied one l'[q]. From the original problem's standpoint alpar@9: * this then means that actually the original column lower bound l[q] alpar@9: * is inactive, and active is that row bound L[p] or U[p] that defines alpar@9: * the implied bound l'[q] (2). In this case in solution to the alpar@9: * original problem column q is assigned status GLP_BS while row p is alpar@9: * assigned status GLP_NL (if a[p,q] > 0) or GLP_NU (if a[p,q] < 0). alpar@9: * Since now column q is basic, its multiplier lambda[q] is zero. This alpar@9: * allows using (7) and (8) to find multiplier for row p in solution to alpar@9: * the original problem: alpar@9: * alpar@9: * 1 alpar@9: * pi[p] = ------ (c[q] - sum a[i,q] pi[i]) = lambda~[q] / a[p,q] (9) alpar@9: * a[p,q] i!=p alpar@9: * alpar@9: * Now let in solution to the transformed problem column q be non-basic alpar@9: * with upper bound active (GLP_NU, lambda~[q] <= 0), and this upper alpar@9: * bound be implied one u'[q]. As in the previous case this then means alpar@9: * that from the original problem's standpoint actually the original alpar@9: * column upper bound u[q] is inactive, and active is that row bound alpar@9: * L[p] or U[p] that defines the implied bound u'[q] (3). In this case alpar@9: * in solution to the original problem column q is assigned status alpar@9: * GLP_BS, row p is assigned status GLP_NU (if a[p,q] > 0) or GLP_NL alpar@9: * (if a[p,q] < 0), and its multiplier is computed with formula (9). alpar@9: * alpar@9: * Strengthening bounds of column q according to (5) and (6) may make alpar@9: * it fixed. Thus, if in solution to the transformed problem column q is alpar@9: * non-basic and fixed (GLP_NS), we can suppose that if lambda~[q] > 0, alpar@9: * column q has active lower bound (GLP_NL), and if lambda~[q] < 0, alpar@9: * column q has active upper bound (GLP_NU), reducing this case to two alpar@9: * previous ones. If, however, lambda~[q] is close to zero or alpar@9: * corresponding bound of row p does not exist (this may happen if alpar@9: * lambda~[q] has wrong sign due to round-off errors, in which case it alpar@9: * is expected to be close to zero, since solution is assumed to be dual alpar@9: * feasible), column q can be assigned status GLP_BS (basic), and row p alpar@9: * can be made active on its existing bound. In the latter case row alpar@9: * multiplier pi[p] computed with formula (9) will be also close to alpar@9: * zero, and dual feasibility will be kept. alpar@9: * alpar@9: * In all other cases, namely, if in solution to the transformed alpar@9: * problem column q is basic (GLP_BS), or non-basic with original lower alpar@9: * bound l[q] active (GLP_NL), or non-basic with original upper bound alpar@9: * u[q] active (GLP_NU), constraint (1) is inactive. So in solution to alpar@9: * the original problem status of column q remains unchanged, row p is alpar@9: * assigned status GLP_BS, and its multiplier pi[p] is assigned zero alpar@9: * value. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * First, value of multiplier for column q in solution to the original alpar@9: * problem is computed with formula (8). If lambda~[q] > 0 and column q alpar@9: * has implied lower bound, or if lambda~[q] < 0 and column q has alpar@9: * implied upper bound, this means that from the original problem's alpar@9: * standpoint actually row p has corresponding active bound, in which alpar@9: * case its multiplier pi[p] is computed with formula (9). In other alpar@9: * cases, when the sign of lambda~[q] corresponds to original bound of alpar@9: * column q, or when lambda~[q] =~ 0, value of row multiplier pi[p] is alpar@9: * assigned zero value. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: struct ineq_singlet alpar@9: { /* row singleton (inequality constraint) */ alpar@9: int p; alpar@9: /* row reference number */ alpar@9: int q; alpar@9: /* column reference number */ alpar@9: double apq; alpar@9: /* constraint coefficient a[p,q] */ alpar@9: double c; alpar@9: /* objective coefficient at x[q] */ alpar@9: double lb; alpar@9: /* row lower bound */ alpar@9: double ub; alpar@9: /* row upper bound */ alpar@9: char lb_changed; alpar@9: /* this flag is set if column lower bound was changed */ alpar@9: char ub_changed; alpar@9: /* this flag is set if column upper bound was changed */ alpar@9: NPPLFE *ptr; alpar@9: /* list of non-zero coefficients a[i,q], i != p */ alpar@9: }; alpar@9: alpar@9: static int rcv_ineq_singlet(NPP *npp, void *info); alpar@9: alpar@9: int npp_ineq_singlet(NPP *npp, NPPROW *p) alpar@9: { /* process row singleton (inequality constraint) */ alpar@9: struct ineq_singlet *info; alpar@9: NPPCOL *q; alpar@9: NPPAIJ *apq, *aij; alpar@9: NPPLFE *lfe; alpar@9: int lb_changed, ub_changed; alpar@9: double ll, uu; alpar@9: /* the row must be singleton inequality constraint */ alpar@9: xassert(p->lb != -DBL_MAX || p->ub != +DBL_MAX); alpar@9: xassert(p->lb < p->ub); alpar@9: xassert(p->ptr != NULL && p->ptr->r_next == NULL); alpar@9: /* compute implied column bounds */ alpar@9: apq = p->ptr; alpar@9: q = apq->col; alpar@9: xassert(q->lb < q->ub); alpar@9: if (apq->val > 0.0) alpar@9: { ll = (p->lb == -DBL_MAX ? -DBL_MAX : p->lb / apq->val); alpar@9: uu = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub / apq->val); alpar@9: } alpar@9: else alpar@9: { ll = (p->ub == +DBL_MAX ? -DBL_MAX : p->ub / apq->val); alpar@9: uu = (p->lb == -DBL_MAX ? +DBL_MAX : p->lb / apq->val); alpar@9: } alpar@9: /* process implied column lower bound */ alpar@9: if (ll == -DBL_MAX) alpar@9: lb_changed = 0; alpar@9: else alpar@9: { lb_changed = npp_implied_lower(npp, q, ll); alpar@9: xassert(0 <= lb_changed && lb_changed <= 4); alpar@9: if (lb_changed == 4) return 4; /* infeasible */ alpar@9: } alpar@9: /* process implied column upper bound */ alpar@9: if (uu == +DBL_MAX) alpar@9: ub_changed = 0; alpar@9: else if (lb_changed == 3) alpar@9: { /* column was fixed on its upper bound due to l'[q] = u[q] */ alpar@9: /* note that L[p] < U[p], so l'[q] = u[q] < u'[q] */ alpar@9: ub_changed = 0; alpar@9: } alpar@9: else alpar@9: { ub_changed = npp_implied_upper(npp, q, uu); alpar@9: xassert(0 <= ub_changed && ub_changed <= 4); alpar@9: if (ub_changed == 4) return 4; /* infeasible */ alpar@9: } alpar@9: /* if neither lower nor upper column bound was changed, the row alpar@9: is originally redundant and can be replaced by free row */ alpar@9: if (!lb_changed && !ub_changed) alpar@9: { p->lb = -DBL_MAX, p->ub = +DBL_MAX; alpar@9: npp_free_row(npp, p); alpar@9: return 0; alpar@9: } alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_ineq_singlet, sizeof(struct ineq_singlet)); alpar@9: info->p = p->i; alpar@9: info->q = q->j; alpar@9: info->apq = apq->val; alpar@9: info->c = q->coef; alpar@9: info->lb = p->lb; alpar@9: info->ub = p->ub; alpar@9: info->lb_changed = (char)lb_changed; alpar@9: info->ub_changed = (char)ub_changed; alpar@9: info->ptr = NULL; alpar@9: /* save column coefficients a[i,q], i != p (not needed for MIP alpar@9: solution) */ alpar@9: if (npp->sol != GLP_MIP) alpar@9: { for (aij = q->ptr; aij != NULL; aij = aij->c_next) alpar@9: { if (aij == apq) continue; /* skip a[p,q] */ 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: /* remove the row from the problem */ alpar@9: npp_del_row(npp, p); alpar@9: return lb_changed >= ub_changed ? lb_changed : ub_changed; alpar@9: } alpar@9: alpar@9: static int rcv_ineq_singlet(NPP *npp, void *_info) alpar@9: { /* recover row singleton (inequality constraint) */ alpar@9: struct ineq_singlet *info = _info; alpar@9: NPPLFE *lfe; alpar@9: double lambda; alpar@9: if (npp->sol == GLP_MIP) goto done; alpar@9: /* compute lambda~[q] in solution to the transformed problem alpar@9: with formula (8) */ 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: if (npp->sol == GLP_SOL) alpar@9: { /* recover basic solution */ alpar@9: if (npp->c_stat[info->q] == GLP_BS) alpar@9: { /* column q is basic, so row p is inactive */ alpar@9: npp->r_stat[info->p] = GLP_BS; alpar@9: npp->r_pi[info->p] = 0.0; alpar@9: } alpar@9: else if (npp->c_stat[info->q] == GLP_NL) alpar@9: nl: { /* column q is non-basic with lower bound active */ alpar@9: if (info->lb_changed) alpar@9: { /* it is implied bound, so actually row p is active alpar@9: while column q is basic */ alpar@9: npp->r_stat[info->p] = alpar@9: (char)(info->apq > 0.0 ? GLP_NL : GLP_NU); alpar@9: npp->c_stat[info->q] = GLP_BS; alpar@9: npp->r_pi[info->p] = lambda / info->apq; alpar@9: } alpar@9: else alpar@9: { /* it is original bound, so row p is inactive */ alpar@9: npp->r_stat[info->p] = GLP_BS; alpar@9: npp->r_pi[info->p] = 0.0; alpar@9: } alpar@9: } alpar@9: else if (npp->c_stat[info->q] == GLP_NU) alpar@9: nu: { /* column q is non-basic with upper bound active */ alpar@9: if (info->ub_changed) alpar@9: { /* it is implied bound, so actually row p is active alpar@9: while column q is basic */ alpar@9: npp->r_stat[info->p] = alpar@9: (char)(info->apq > 0.0 ? GLP_NU : GLP_NL); alpar@9: npp->c_stat[info->q] = GLP_BS; alpar@9: npp->r_pi[info->p] = lambda / info->apq; alpar@9: } alpar@9: else alpar@9: { /* it is original bound, so row p is inactive */ alpar@9: npp->r_stat[info->p] = GLP_BS; alpar@9: npp->r_pi[info->p] = 0.0; alpar@9: } alpar@9: } alpar@9: else if (npp->c_stat[info->q] == GLP_NS) alpar@9: { /* column q is non-basic and fixed; note, however, that in alpar@9: in the original problem it is non-fixed */ alpar@9: if (lambda > +1e-7) alpar@9: { if (info->apq > 0.0 && info->lb != -DBL_MAX || alpar@9: info->apq < 0.0 && info->ub != +DBL_MAX || alpar@9: !info->lb_changed) alpar@9: { /* either corresponding bound of row p exists or alpar@9: column q remains non-basic with its original lower alpar@9: bound active */ alpar@9: npp->c_stat[info->q] = GLP_NL; alpar@9: goto nl; alpar@9: } alpar@9: } alpar@9: if (lambda < -1e-7) alpar@9: { if (info->apq > 0.0 && info->ub != +DBL_MAX || alpar@9: info->apq < 0.0 && info->lb != -DBL_MAX || alpar@9: !info->ub_changed) alpar@9: { /* either corresponding bound of row p exists or alpar@9: column q remains non-basic with its original upper alpar@9: bound active */ alpar@9: npp->c_stat[info->q] = GLP_NU; alpar@9: goto nu; alpar@9: } alpar@9: } alpar@9: /* either lambda~[q] is close to zero, or corresponding alpar@9: bound of row p does not exist, because lambda~[q] has alpar@9: wrong sign due to round-off errors; in the latter case alpar@9: lambda~[q] is also assumed to be close to zero; so, we alpar@9: can make row p active on its existing bound and column q alpar@9: basic; pi[p] will have wrong sign, but it also will be alpar@9: close to zero (rarus casus of dual degeneracy) */ alpar@9: if (info->lb != -DBL_MAX && info->ub == +DBL_MAX) alpar@9: { /* row lower bound exists, but upper bound doesn't */ alpar@9: npp->r_stat[info->p] = GLP_NL; alpar@9: } alpar@9: else if (info->lb == -DBL_MAX && info->ub != +DBL_MAX) alpar@9: { /* row upper bound exists, but lower bound doesn't */ alpar@9: npp->r_stat[info->p] = GLP_NU; alpar@9: } alpar@9: else if (info->lb != -DBL_MAX && info->ub != +DBL_MAX) alpar@9: { /* both row lower and upper bounds exist */ alpar@9: /* to choose proper active row bound we should not use alpar@9: lambda~[q], because its value being close to zero is alpar@9: unreliable; so we choose that bound which provides alpar@9: primal feasibility for original constraint (1) */ alpar@9: if (info->apq * npp->c_value[info->q] <= alpar@9: 0.5 * (info->lb + info->ub)) 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: npp->c_stat[info->q] = GLP_BS; alpar@9: npp->r_pi[info->p] = lambda / info->apq; alpar@9: } alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: if (npp->sol == GLP_IPT) alpar@9: { /* recover interior-point solution */ alpar@9: if (lambda > +DBL_EPSILON && info->lb_changed || alpar@9: lambda < -DBL_EPSILON && info->ub_changed) alpar@9: { /* actually row p has corresponding active bound */ alpar@9: npp->r_pi[info->p] = lambda / info->apq; alpar@9: } alpar@9: else alpar@9: { /* either bounds of column q are both inactive or its alpar@9: original bound is active */ alpar@9: npp->r_pi[info->p] = 0.0; alpar@9: } alpar@9: } alpar@9: done: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_implied_slack - process column singleton (implied slack variable) alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_implied_slack(NPP *npp, NPPCOL *q); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_implied_slack processes column q: alpar@9: * alpar@9: * l[q] <= x[q] <= u[q], (1) alpar@9: * alpar@9: * where l[q] < u[q], having the only non-zero coefficient in row p, alpar@9: * which is equality constraint: alpar@9: * alpar@9: * sum a[p,j] x[j] + a[p,q] x[q] = b. (2) alpar@9: * j!=q alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * (If x[q] is integral, this transformation must not be used.) alpar@9: * alpar@9: * The term a[p,q] x[q] in constraint (2) can be considered as a slack alpar@9: * variable that allows to carry bounds of column q over row p and then alpar@9: * remove column q from the problem. alpar@9: * alpar@9: * Constraint (2) can be written as follows: alpar@9: * alpar@9: * sum a[p,j] x[j] = b - a[p,q] x[q]. (3) alpar@9: * j!=q alpar@9: * alpar@9: * According to (1) constraint (3) is equivalent to the following alpar@9: * inequality constraint: alpar@9: * alpar@9: * L[p] <= sum a[p,j] x[j] <= U[p], (4) alpar@9: * j!=q alpar@9: * alpar@9: * where alpar@9: * alpar@9: * ( b - a[p,q] u[q], if a[p,q] > 0 alpar@9: * L[p] = < (5) alpar@9: * ( b - a[p,q] l[q], if a[p,q] < 0 alpar@9: * alpar@9: * ( b - a[p,q] l[q], if a[p,q] > 0 alpar@9: * U[p] = < (6) alpar@9: * ( b - a[p,q] u[q], if a[p,q] < 0 alpar@9: * alpar@9: * From (2) it follows that: alpar@9: * alpar@9: * 1 alpar@9: * x[q] = ------ (b - sum a[p,j] x[j]). (7) alpar@9: * a[p,q] j!=q alpar@9: * alpar@9: * In order to eliminate x[q] from the objective row we substitute it alpar@9: * from (6) to that row: alpar@9: * alpar@9: * z = sum c[j] x[j] + c[q] x[q] + c[0] = alpar@9: * j!=q alpar@9: * 1 alpar@9: * = sum c[j] x[j] + c[q] [------ (b - sum a[p,j] x[j])] + c0 = alpar@9: * j!=q a[p,q] j!=q alpar@9: * alpar@9: * = sum c~[j] x[j] + c~[0], alpar@9: * j!=q alpar@9: * a[p,j] b alpar@9: * c~[j] = c[j] - c[q] ------, c~0 = c0 - c[q] ------ (8) alpar@9: * a[p,q] a[p,q] alpar@9: * alpar@9: * are values of objective coefficients and constant term, resp., in alpar@9: * the transformed problem. alpar@9: * alpar@9: * Note that column q is column singleton, so in the dual system of the alpar@9: * original problem it corresponds to the following row singleton: alpar@9: * alpar@9: * a[p,q] pi[p] + lambda[q] = c[q]. (9) alpar@9: * alpar@9: * In the transformed problem row (9) would be the following: alpar@9: * alpar@9: * a[p,q] pi~[p] + lambda[q] = c~[q] = 0. (10) alpar@9: * alpar@9: * Subtracting (10) from (9) we have: alpar@9: * alpar@9: * a[p,q] (pi[p] - pi~[p]) = c[q] alpar@9: * alpar@9: * that gives the following formula to compute multiplier for row p in alpar@9: * solution to the original problem using its value in solution to the alpar@9: * transformed problem: alpar@9: * alpar@9: * pi[p] = pi~[p] + c[q] / a[p,q]. (11) alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * Status of column q in solution to the original problem is defined alpar@9: * by status of row p in solution to the transformed problem and the alpar@9: * sign of coefficient a[p,q] in the original inequality constraint (2) alpar@9: * as follows: alpar@9: * alpar@9: * +-----------------------+---------+--------------------+ alpar@9: * | Status of row p | Sign of | Status of column q | alpar@9: * | (transformed problem) | a[p,q] | (original problem) | alpar@9: * +-----------------------+---------+--------------------+ alpar@9: * | GLP_BS | + / - | GLP_BS | alpar@9: * | GLP_NL | + | GLP_NU | alpar@9: * | GLP_NL | - | GLP_NL | alpar@9: * | GLP_NU | + | GLP_NL | alpar@9: * | GLP_NU | - | GLP_NU | alpar@9: * | GLP_NF | + / - | GLP_NF | alpar@9: * +-----------------------+---------+--------------------+ alpar@9: * alpar@9: * Value of column q is computed with formula (7). Since originally row alpar@9: * p is equality constraint, its status is assigned GLP_NS, and value of alpar@9: * its multiplier pi[p] is computed with formula (11). alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Value of column q is computed with formula (7). Row multiplier value alpar@9: * pi[p] is computed with formula (11). alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * Value of column q is computed with formula (7). */ alpar@9: alpar@9: struct implied_slack alpar@9: { /* column singleton (implied slack variable) */ alpar@9: int p; alpar@9: /* row reference number */ alpar@9: int q; alpar@9: /* column reference number */ alpar@9: double apq; alpar@9: /* constraint coefficient a[p,q] */ alpar@9: double b; alpar@9: /* right-hand side of original equality constraint */ alpar@9: double c; alpar@9: /* original objective coefficient at x[q] */ alpar@9: NPPLFE *ptr; alpar@9: /* list of non-zero coefficients a[p,j], j != q */ alpar@9: }; alpar@9: alpar@9: static int rcv_implied_slack(NPP *npp, void *info); alpar@9: alpar@9: void npp_implied_slack(NPP *npp, NPPCOL *q) alpar@9: { /* process column singleton (implied slack variable) */ alpar@9: struct implied_slack *info; alpar@9: NPPROW *p; alpar@9: NPPAIJ *aij; alpar@9: NPPLFE *lfe; alpar@9: /* the column must be non-integral non-fixed singleton */ alpar@9: xassert(!q->is_int); alpar@9: xassert(q->lb < q->ub); alpar@9: xassert(q->ptr != NULL && q->ptr->c_next == NULL); alpar@9: /* corresponding row must be equality constraint */ alpar@9: aij = q->ptr; alpar@9: p = aij->row; alpar@9: xassert(p->lb == p->ub); alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_implied_slack, sizeof(struct implied_slack)); alpar@9: info->p = p->i; alpar@9: info->q = q->j; alpar@9: info->apq = aij->val; alpar@9: info->b = p->lb; alpar@9: info->c = q->coef; alpar@9: info->ptr = NULL; alpar@9: /* save row coefficients a[p,j], j != q, and substitute x[q] alpar@9: into the objective row */ alpar@9: for (aij = p->ptr; aij != NULL; aij = aij->r_next) alpar@9: { if (aij->col == q) continue; /* skip a[p,q] */ alpar@9: lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE)); alpar@9: lfe->ref = aij->col->j; alpar@9: lfe->val = aij->val; alpar@9: lfe->next = info->ptr; alpar@9: info->ptr = lfe; alpar@9: aij->col->coef -= info->c * (aij->val / info->apq); alpar@9: } alpar@9: npp->c0 += info->c * (info->b / info->apq); alpar@9: /* compute new row bounds */ alpar@9: if (info->apq > 0.0) alpar@9: { p->lb = (q->ub == +DBL_MAX ? alpar@9: -DBL_MAX : info->b - info->apq * q->ub); alpar@9: p->ub = (q->lb == -DBL_MAX ? alpar@9: +DBL_MAX : info->b - info->apq * q->lb); alpar@9: } alpar@9: else alpar@9: { p->lb = (q->lb == -DBL_MAX ? alpar@9: -DBL_MAX : info->b - info->apq * q->lb); alpar@9: p->ub = (q->ub == +DBL_MAX ? alpar@9: +DBL_MAX : info->b - info->apq * q->ub); 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_implied_slack(NPP *npp, void *_info) alpar@9: { /* recover column singleton (implied slack variable) */ alpar@9: struct implied_slack *info = _info; alpar@9: NPPLFE *lfe; alpar@9: double temp; alpar@9: if (npp->sol == GLP_SOL) alpar@9: { /* assign statuses to row p and column q */ alpar@9: if (npp->r_stat[info->p] == GLP_BS || alpar@9: npp->r_stat[info->p] == GLP_NF) alpar@9: npp->c_stat[info->q] = npp->r_stat[info->p]; alpar@9: else if (npp->r_stat[info->p] == GLP_NL) alpar@9: npp->c_stat[info->q] = alpar@9: (char)(info->apq > 0.0 ? GLP_NU : GLP_NL); alpar@9: else if (npp->r_stat[info->p] == GLP_NU) alpar@9: npp->c_stat[info->q] = alpar@9: (char)(info->apq > 0.0 ? GLP_NL : GLP_NU); alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: npp->r_stat[info->p] = GLP_NS; alpar@9: } alpar@9: if (npp->sol != GLP_MIP) alpar@9: { /* compute multiplier for row p */ alpar@9: npp->r_pi[info->p] += info->c / info->apq; alpar@9: } alpar@9: /* compute value of column q */ alpar@9: temp = info->b; alpar@9: for (lfe = info->ptr; lfe != NULL; lfe = lfe->next) alpar@9: temp -= lfe->val * npp->c_value[lfe->ref]; alpar@9: npp->c_value[info->q] = temp / info->apq; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_implied_free - process column singleton (implied free variable) alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_implied_free(NPP *npp, NPPCOL *q); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_implied_free processes column q: alpar@9: * alpar@9: * l[q] <= x[q] <= u[q], (1) alpar@9: * alpar@9: * having non-zero coefficient in the only row p, which is inequality alpar@9: * constraint: alpar@9: * alpar@9: * L[p] <= sum a[p,j] x[j] + a[p,q] x[q] <= U[p], (2) alpar@9: * j!=q alpar@9: * alpar@9: * where l[q] < u[q], L[p] < U[p], L[p] > -oo and/or U[p] < +oo. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - success; alpar@9: * alpar@9: * 1 - column lower and/or upper bound(s) can be active; alpar@9: * alpar@9: * 2 - problem has no dual feasible solution. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Constraint (2) can be written as follows: alpar@9: * alpar@9: * L[p] - sum a[p,j] x[j] <= a[p,q] x[q] <= U[p] - sum a[p,j] x[j], alpar@9: * j!=q j!=q alpar@9: * alpar@9: * from which it follows that: alpar@9: * alpar@9: * alfa <= a[p,q] x[q] <= beta, (3) alpar@9: * alpar@9: * where alpar@9: * alpar@9: * alfa = inf(L[p] - sum a[p,j] x[j]) = alpar@9: * j!=q alpar@9: * alpar@9: * = L[p] - sup sum a[p,j] x[j] = (4) alpar@9: * j!=q alpar@9: * alpar@9: * = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j], alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * beta = sup(L[p] - sum a[p,j] x[j]) = alpar@9: * j!=q alpar@9: * alpar@9: * = L[p] - inf sum a[p,j] x[j] = (5) alpar@9: * j!=q alpar@9: * alpar@9: * = L[p] - sum a[p,j] l[j] - sum a[p,j] u[j], alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * Jp = {j != q: a[p,j] > 0}, Jn = {j != q: a[p,j] < 0}. (6) alpar@9: * alpar@9: * Inequality (3) defines implied bounds of variable x[q]: alpar@9: * alpar@9: * l'[q] <= x[q] <= u'[q], (7) alpar@9: * alpar@9: * where alpar@9: * alpar@9: * ( alfa / a[p,q], if a[p,q] > 0 alpar@9: * l'[q] = < (8a) alpar@9: * ( beta / a[p,q], if a[p,q] < 0 alpar@9: * alpar@9: * ( beta / a[p,q], if a[p,q] > 0 alpar@9: * u'[q] = < (8b) alpar@9: * ( alfa / a[p,q], if a[p,q] < 0 alpar@9: * alpar@9: * Thus, if l'[q] > l[q] - eps and u'[q] < u[q] + eps, where eps is alpar@9: * an absolute tolerance for column value, column bounds (1) cannot be alpar@9: * active, in which case column q can be replaced by equivalent free alpar@9: * (unbounded) column. alpar@9: * alpar@9: * Note that column q is column singleton, so in the dual system of the alpar@9: * original problem it corresponds to the following row singleton: alpar@9: * alpar@9: * a[p,q] pi[p] + lambda[q] = c[q], (9) alpar@9: * alpar@9: * from which it follows that: alpar@9: * alpar@9: * pi[p] = (c[q] - lambda[q]) / a[p,q]. (10) alpar@9: * alpar@9: * Let x[q] be implied free (unbounded) variable. Then column q can be alpar@9: * only basic, so its multiplier lambda[q] is equal to zero, and from alpar@9: * (10) we have: alpar@9: * alpar@9: * pi[p] = c[q] / a[p,q]. (11) alpar@9: * alpar@9: * There are possible three cases: alpar@9: * alpar@9: * 1) pi[p] < -eps, where eps is an absolute tolerance for row alpar@9: * multiplier. In this case, to provide dual feasibility of the alpar@9: * original problem, row p must be active on its lower bound, and alpar@9: * if its lower bound does not exist (L[p] = -oo), the problem has alpar@9: * no dual feasible solution; alpar@9: * alpar@9: * 2) pi[p] > +eps. In this case row p must be active on its upper alpar@9: * bound, and if its upper bound does not exist (U[p] = +oo), the alpar@9: * problem has no dual feasible solution; alpar@9: * alpar@9: * 3) -eps <= pi[p] <= +eps. In this case any (either lower or upper) alpar@9: * bound of row p can be active, because this does not affect dual alpar@9: * feasibility. alpar@9: * alpar@9: * Thus, in all three cases original inequality constraint (2) can be alpar@9: * replaced by equality constraint, where the right-hand side is either alpar@9: * lower or upper bound of row p, and bounds of column q can be removed alpar@9: * that makes it free (unbounded). (May note that this transformation alpar@9: * can be followed by transformation "Column singleton (implied slack alpar@9: * variable)" performed by the routine npp_implied_slack.) 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 in solution to the transformed problem and its bound, alpar@9: * which was choosen to be active: alpar@9: * alpar@9: * +-----------------------+--------+--------------------+ alpar@9: * | Status of row p | Active | Status of row p | alpar@9: * | (transformed problem) | bound | (original problem) | alpar@9: * +-----------------------+--------+--------------------+ alpar@9: * | GLP_BS | L[p] | GLP_BS | alpar@9: * | GLP_BS | U[p] | GLP_BS | alpar@9: * | GLP_NS | L[p] | GLP_NL | alpar@9: * | GLP_NS | U[p] | GLP_NU | alpar@9: * +-----------------------+--------+--------------------+ alpar@9: * alpar@9: * Value of row multiplier pi[p] (as well as value of column q) in alpar@9: * solution to the original problem is the same as in solution to the alpar@9: * 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 implied_free alpar@9: { /* column singleton (implied free variable) */ alpar@9: int p; alpar@9: /* row reference number */ alpar@9: char stat; alpar@9: /* row status: alpar@9: GLP_NL - active constraint on lower bound alpar@9: GLP_NU - active constraint on upper bound */ alpar@9: }; alpar@9: alpar@9: static int rcv_implied_free(NPP *npp, void *info); alpar@9: alpar@9: int npp_implied_free(NPP *npp, NPPCOL *q) alpar@9: { /* process column singleton (implied free variable) */ alpar@9: struct implied_free *info; alpar@9: NPPROW *p; alpar@9: NPPAIJ *apq, *aij; alpar@9: double alfa, beta, l, u, pi, eps; alpar@9: /* the column must be non-fixed singleton */ alpar@9: xassert(q->lb < q->ub); alpar@9: xassert(q->ptr != NULL && q->ptr->c_next == NULL); alpar@9: /* corresponding row must be inequality constraint */ alpar@9: apq = q->ptr; alpar@9: p = apq->row; alpar@9: xassert(p->lb != -DBL_MAX || p->ub != +DBL_MAX); alpar@9: xassert(p->lb < p->ub); alpar@9: /* compute alfa */ alpar@9: alfa = p->lb; alpar@9: if (alfa != -DBL_MAX) alpar@9: { for (aij = p->ptr; aij != NULL; aij = aij->r_next) alpar@9: { if (aij == apq) continue; /* skip a[p,q] */ alpar@9: if (aij->val > 0.0) alpar@9: { if (aij->col->ub == +DBL_MAX) alpar@9: { alfa = -DBL_MAX; alpar@9: break; alpar@9: } alpar@9: alfa -= aij->val * aij->col->ub; alpar@9: } alpar@9: else /* < 0.0 */ alpar@9: { if (aij->col->lb == -DBL_MAX) alpar@9: { alfa = -DBL_MAX; alpar@9: break; alpar@9: } alpar@9: alfa -= aij->val * aij->col->lb; alpar@9: } alpar@9: } alpar@9: } alpar@9: /* compute beta */ alpar@9: beta = p->ub; alpar@9: if (beta != +DBL_MAX) alpar@9: { for (aij = p->ptr; aij != NULL; aij = aij->r_next) alpar@9: { if (aij == apq) continue; /* skip a[p,q] */ alpar@9: if (aij->val > 0.0) alpar@9: { if (aij->col->lb == -DBL_MAX) alpar@9: { beta = +DBL_MAX; alpar@9: break; alpar@9: } alpar@9: beta -= aij->val * aij->col->lb; alpar@9: } alpar@9: else /* < 0.0 */ alpar@9: { if (aij->col->ub == +DBL_MAX) alpar@9: { beta = +DBL_MAX; alpar@9: break; alpar@9: } alpar@9: beta -= aij->val * aij->col->ub; alpar@9: } alpar@9: } alpar@9: } alpar@9: /* compute implied column lower bound l'[q] */ alpar@9: if (apq->val > 0.0) alpar@9: l = (alfa == -DBL_MAX ? -DBL_MAX : alfa / apq->val); alpar@9: else /* < 0.0 */ alpar@9: l = (beta == +DBL_MAX ? -DBL_MAX : beta / apq->val); alpar@9: /* compute implied column upper bound u'[q] */ alpar@9: if (apq->val > 0.0) alpar@9: u = (beta == +DBL_MAX ? +DBL_MAX : beta / apq->val); alpar@9: else alpar@9: u = (alfa == -DBL_MAX ? +DBL_MAX : alfa / apq->val); alpar@9: /* check if column lower bound l[q] can be active */ alpar@9: if (q->lb != -DBL_MAX) alpar@9: { eps = 1e-9 + 1e-12 * fabs(q->lb); alpar@9: if (l < q->lb - eps) return 1; /* yes, it can */ alpar@9: } alpar@9: /* check if column upper bound u[q] can be active */ alpar@9: if (q->ub != +DBL_MAX) alpar@9: { eps = 1e-9 + 1e-12 * fabs(q->ub); alpar@9: if (u > q->ub + eps) return 1; /* yes, it can */ alpar@9: } alpar@9: /* okay; make column q free (unbounded) */ alpar@9: q->lb = -DBL_MAX, q->ub = +DBL_MAX; alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_implied_free, sizeof(struct implied_free)); alpar@9: info->p = p->i; alpar@9: info->stat = -1; alpar@9: /* compute row multiplier pi[p] */ alpar@9: pi = q->coef / apq->val; alpar@9: /* check dual feasibility for row p */ alpar@9: if (pi > +DBL_EPSILON) alpar@9: { /* lower bound L[p] must be active */ alpar@9: if (p->lb != -DBL_MAX) alpar@9: nl: { info->stat = GLP_NL; alpar@9: p->ub = p->lb; alpar@9: } alpar@9: else alpar@9: { if (pi > +1e-5) return 2; /* dual infeasibility */ alpar@9: /* take a chance on U[p] */ alpar@9: xassert(p->ub != +DBL_MAX); alpar@9: goto nu; alpar@9: } alpar@9: } alpar@9: else if (pi < -DBL_EPSILON) alpar@9: { /* upper bound U[p] must be active */ alpar@9: if (p->ub != +DBL_MAX) alpar@9: nu: { info->stat = GLP_NU; alpar@9: p->lb = p->ub; alpar@9: } alpar@9: else alpar@9: { if (pi < -1e-5) return 2; /* dual infeasibility */ alpar@9: /* take a chance on L[p] */ alpar@9: xassert(p->lb != -DBL_MAX); alpar@9: goto nl; alpar@9: } alpar@9: } alpar@9: else alpar@9: { /* any bound (either L[p] or U[p]) can be made active */ alpar@9: if (p->ub == +DBL_MAX) alpar@9: { xassert(p->lb != -DBL_MAX); alpar@9: goto nl; alpar@9: } alpar@9: if (p->lb == -DBL_MAX) alpar@9: { xassert(p->ub != +DBL_MAX); alpar@9: goto nu; alpar@9: } alpar@9: if (fabs(p->lb) <= fabs(p->ub)) goto nl; else goto nu; alpar@9: } alpar@9: return 0; alpar@9: } alpar@9: alpar@9: static int rcv_implied_free(NPP *npp, void *_info) alpar@9: { /* recover column singleton (implied free variable) */ alpar@9: struct implied_free *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: { xassert(info->stat == GLP_NL || info->stat == GLP_NU); alpar@9: npp->r_stat[info->p] = info->stat; 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_eq_doublet - process row doubleton (equality constraint) alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * NPPCOL *npp_eq_doublet(NPP *npp, NPPROW *p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_eq_doublet processes row p, which is equality alpar@9: * constraint having exactly two non-zero coefficients: alpar@9: * alpar@9: * a[p,q] x[q] + a[p,r] x[r] = b. (1) alpar@9: * alpar@9: * As the result of processing one of columns q or r is eliminated from alpar@9: * all other rows and, thus, becomes column singleton of type "implied alpar@9: * slack variable". Row p is not changed and along with column q and r alpar@9: * remains in the problem. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * The routine npp_eq_doublet returns pointer to the descriptor of that alpar@9: * column q or r which has been eliminated. If, due to some reason, the alpar@9: * elimination was not performed, the routine returns NULL. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * First, we decide which column q or r will be eliminated. Let it be alpar@9: * column q. Consider i-th constraint row, where column q has non-zero alpar@9: * coefficient a[i,q] != 0: alpar@9: * alpar@9: * L[i] <= sum a[i,j] x[j] <= U[i]. (2) alpar@9: * j alpar@9: * alpar@9: * In order to eliminate column q from row (2) we subtract from it row alpar@9: * (1) multiplied by gamma[i] = a[i,q] / a[p,q], i.e. we replace in the alpar@9: * transformed problem row (2) by its linear combination with row (1). alpar@9: * This transformation changes only coefficients in columns q and r, alpar@9: * and bounds of row i as follows: alpar@9: * alpar@9: * a~[i,q] = a[i,q] - gamma[i] a[p,q] = 0, (3) alpar@9: * alpar@9: * a~[i,r] = a[i,r] - gamma[i] a[p,r], (4) alpar@9: * alpar@9: * L~[i] = L[i] - gamma[i] b, (5) alpar@9: * alpar@9: * U~[i] = U[i] - gamma[i] b. (6) alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * The transformation of the primal system of the original problem: alpar@9: * alpar@9: * L <= A x <= U (7) alpar@9: * alpar@9: * is equivalent to multiplying from the left a transformation matrix F alpar@9: * by components of this primal system, which in the transformed problem alpar@9: * becomes the following: alpar@9: * alpar@9: * F L <= F A x <= F U ==> L~ <= A~x <= U~. (8) alpar@9: * alpar@9: * The matrix F has the following structure: alpar@9: * alpar@9: * ( 1 -gamma[1] ) alpar@9: * ( ) alpar@9: * ( 1 -gamma[2] ) alpar@9: * ( ) alpar@9: * ( ... ... ) alpar@9: * ( ) alpar@9: * F = ( 1 -gamma[p-1] ) (9) alpar@9: * ( ) alpar@9: * ( 1 ) alpar@9: * ( ) alpar@9: * ( -gamma[p+1] 1 ) alpar@9: * ( ) alpar@9: * ( ... ... ) alpar@9: * alpar@9: * where its column containing elements -gamma[i] corresponds to row p alpar@9: * of the primal system. alpar@9: * alpar@9: * From (8) it follows that the dual system of the original problem: alpar@9: * alpar@9: * A'pi + lambda = c, (10) alpar@9: * alpar@9: * in the transformed problem becomes the following: alpar@9: * alpar@9: * A'F'inv(F')pi + lambda = c ==> (A~)'pi~ + lambda = c, (11) alpar@9: * alpar@9: * where: alpar@9: * alpar@9: * pi~ = inv(F')pi (12) alpar@9: * alpar@9: * is the vector of row multipliers in the transformed problem. Thus: alpar@9: * alpar@9: * pi = F'pi~. (13) alpar@9: * alpar@9: * Therefore, as it follows from (13), value of multiplier for row p in alpar@9: * solution to the original problem can be computed as follows: alpar@9: * alpar@9: * pi[p] = pi~[p] - sum gamma[i] pi~[i], (14) alpar@9: * i alpar@9: * alpar@9: * where pi~[i] = pi[i] is multiplier for row i (i != p). alpar@9: * alpar@9: * Note that the statuses of all rows and columns are not changed. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Multiplier for row p in solution to the original problem is computed alpar@9: * with formula (14). alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: struct eq_doublet alpar@9: { /* row doubleton (equality constraint) */ alpar@9: int p; alpar@9: /* row reference number */ alpar@9: double apq; alpar@9: /* constraint coefficient a[p,q] */ alpar@9: NPPLFE *ptr; alpar@9: /* list of non-zero coefficients a[i,q], i != p */ alpar@9: }; alpar@9: alpar@9: static int rcv_eq_doublet(NPP *npp, void *info); alpar@9: alpar@9: NPPCOL *npp_eq_doublet(NPP *npp, NPPROW *p) alpar@9: { /* process row doubleton (equality constraint) */ alpar@9: struct eq_doublet *info; alpar@9: NPPROW *i; alpar@9: NPPCOL *q, *r; alpar@9: NPPAIJ *apq, *apr, *aiq, *air, *next; alpar@9: NPPLFE *lfe; alpar@9: double gamma; alpar@9: /* the row must be doubleton equality constraint */ alpar@9: xassert(p->lb == p->ub); alpar@9: xassert(p->ptr != NULL && p->ptr->r_next != NULL && alpar@9: p->ptr->r_next->r_next == NULL); alpar@9: /* choose column to be eliminated */ alpar@9: { NPPAIJ *a1, *a2; alpar@9: a1 = p->ptr, a2 = a1->r_next; alpar@9: if (fabs(a2->val) < 0.001 * fabs(a1->val)) alpar@9: { /* only first column can be eliminated, because second one alpar@9: has too small constraint coefficient */ alpar@9: apq = a1, apr = a2; alpar@9: } alpar@9: else if (fabs(a1->val) < 0.001 * fabs(a2->val)) alpar@9: { /* only second column can be eliminated, because first one alpar@9: has too small constraint coefficient */ alpar@9: apq = a2, apr = a1; alpar@9: } alpar@9: else alpar@9: { /* both columns are appropriate; choose that one which is alpar@9: shorter to minimize fill-in */ alpar@9: if (npp_col_nnz(npp, a1->col) <= npp_col_nnz(npp, a2->col)) alpar@9: { /* first column is shorter */ alpar@9: apq = a1, apr = a2; alpar@9: } alpar@9: else alpar@9: { /* second column is shorter */ alpar@9: apq = a2, apr = a1; alpar@9: } alpar@9: } alpar@9: } alpar@9: /* now columns q and r have been chosen */ alpar@9: q = apq->col, r = apr->col; alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_eq_doublet, sizeof(struct eq_doublet)); alpar@9: info->p = p->i; alpar@9: info->apq = apq->val; alpar@9: info->ptr = NULL; alpar@9: /* transform each row i (i != p), where a[i,q] != 0, to eliminate alpar@9: column q */ alpar@9: for (aiq = q->ptr; aiq != NULL; aiq = next) alpar@9: { next = aiq->c_next; alpar@9: if (aiq == apq) continue; /* skip row p */ alpar@9: i = aiq->row; /* row i to be transformed */ alpar@9: /* save constraint coefficient a[i,q] */ alpar@9: if (npp->sol != GLP_MIP) alpar@9: { lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE)); alpar@9: lfe->ref = i->i; alpar@9: lfe->val = aiq->val; alpar@9: lfe->next = info->ptr; alpar@9: info->ptr = lfe; alpar@9: } alpar@9: /* find coefficient a[i,r] in row i */ alpar@9: for (air = i->ptr; air != NULL; air = air->r_next) alpar@9: if (air->col == r) break; alpar@9: /* if a[i,r] does not exist, create a[i,r] = 0 */ alpar@9: if (air == NULL) alpar@9: air = npp_add_aij(npp, i, r, 0.0); alpar@9: /* compute gamma[i] = a[i,q] / a[p,q] */ alpar@9: gamma = aiq->val / apq->val; alpar@9: /* (row i) := (row i) - gamma[i] * (row p); see (3)-(6) */ alpar@9: /* new a[i,q] is exact zero due to elimnation; remove it from alpar@9: row i */ alpar@9: npp_del_aij(npp, aiq); alpar@9: /* compute new a[i,r] */ alpar@9: air->val -= gamma * apr->val; alpar@9: /* if new a[i,r] is close to zero due to numeric cancelation, alpar@9: remove it from row i */ alpar@9: if (fabs(air->val) <= 1e-10) alpar@9: npp_del_aij(npp, air); alpar@9: /* compute new lower and upper bounds of row i */ alpar@9: if (i->lb == i->ub) alpar@9: i->lb = i->ub = (i->lb - gamma * p->lb); alpar@9: else alpar@9: { if (i->lb != -DBL_MAX) alpar@9: i->lb -= gamma * p->lb; alpar@9: if (i->ub != +DBL_MAX) alpar@9: i->ub -= gamma * p->lb; alpar@9: } alpar@9: } alpar@9: return q; alpar@9: } alpar@9: alpar@9: static int rcv_eq_doublet(NPP *npp, void *_info) alpar@9: { /* recover row doubleton (equality constraint) */ alpar@9: struct eq_doublet *info = _info; alpar@9: NPPLFE *lfe; alpar@9: double gamma, temp; alpar@9: /* we assume that processing row p is followed by processing alpar@9: column q as singleton of type "implied slack variable", in alpar@9: which case row p must always be active equality constraint */ alpar@9: if (npp->sol == GLP_SOL) alpar@9: { if (npp->r_stat[info->p] != GLP_NS) alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: if (npp->sol != GLP_MIP) alpar@9: { /* compute value of multiplier for row p; see (14) */ alpar@9: temp = npp->r_pi[info->p]; alpar@9: for (lfe = info->ptr; lfe != NULL; lfe = lfe->next) alpar@9: { gamma = lfe->val / info->apq; /* a[i,q] / a[p,q] */ alpar@9: temp -= gamma * npp->r_pi[lfe->ref]; alpar@9: } alpar@9: npp->r_pi[info->p] = temp; alpar@9: } alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_forcing_row - process forcing row alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_forcing_row(NPP *npp, NPPROW *p, int at); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_forcing row processes row p of general format: alpar@9: * alpar@9: * L[p] <= sum a[p,j] x[j] <= U[p], (1) alpar@9: * j alpar@9: * alpar@9: * l[j] <= x[j] <= u[j], (2) alpar@9: * alpar@9: * where L[p] <= U[p] and l[j] < u[j] for all a[p,j] != 0. It is also alpar@9: * assumed that: alpar@9: * alpar@9: * 1) if at = 0 then |L[p] - U'[p]| <= eps, where U'[p] is implied alpar@9: * row upper bound (see below), eps is an absolute tolerance for row alpar@9: * value; alpar@9: * alpar@9: * 2) if at = 1 then |U[p] - L'[p]| <= eps, where L'[p] is implied alpar@9: * row lower bound (see below). alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0 - success; alpar@9: * alpar@9: * 1 - cannot fix columns due to too small constraint coefficients. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * Implied lower and upper bounds of row (1) are determined by bounds alpar@9: * of corresponding columns (variables) as follows: alpar@9: * alpar@9: * L'[p] = inf sum a[p,j] x[j] = alpar@9: * j alpar@9: * (3) alpar@9: * = sum a[p,j] l[j] + sum a[p,j] u[j], alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * U'[p] = sup sum a[p,j] x[j] = alpar@9: * (4) alpar@9: * = sum a[p,j] u[j] + sum a[p,j] l[j], alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5) alpar@9: * alpar@9: * If L[p] =~ U'[p] (at = 0), solution can be primal feasible only when alpar@9: * all variables take their boundary values as defined by (4): alpar@9: * alpar@9: * ( u[j], if j in Jp alpar@9: * x[j] = < (6) alpar@9: * ( l[j], if j in Jn alpar@9: * alpar@9: * Similarly, if U[p] =~ L'[p] (at = 1), solution can be primal feasible alpar@9: * only when all variables take their boundary values as defined by (3): alpar@9: * alpar@9: * ( l[j], if j in Jp alpar@9: * x[j] = < (7) alpar@9: * ( u[j], if j in Jn alpar@9: * alpar@9: * Condition (6) or (7) allows fixing all columns (variables x[j]) alpar@9: * in row (1) on their bounds and then removing them from the problem alpar@9: * (see the routine npp_fixed_col). Due to this row p becomes redundant, alpar@9: * so it can be replaced by equivalent free (unbounded) row and also alpar@9: * removed from the problem (see the routine npp_free_row). alpar@9: * alpar@9: * 1. To apply this transformation row (1) should not have coefficients alpar@9: * whose magnitude is too small, i.e. all a[p,j] should satisfy to alpar@9: * the following condition: alpar@9: * alpar@9: * |a[p,j]| >= eps * max(1, |a[p,k]|), (8) alpar@9: * k alpar@9: * where eps is a relative tolerance for constraint coefficients. alpar@9: * Otherwise, fixing columns may be numerically unreliable and may alpar@9: * lead to wrong solution. alpar@9: * alpar@9: * 2. The routine fixes columns and remove bounds of row p, however, alpar@9: * it does not remove the row and columns from the problem. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * In the transformed problem row p being inactive constraint is alpar@9: * assigned status GLP_BS (as the result of transformation of free alpar@9: * row), and all columns in this row are assigned status GLP_NS (as the alpar@9: * result of transformation of fixed columns). alpar@9: * alpar@9: * Note that in the dual system of the transformed (as well as original) alpar@9: * problem every column j in row p corresponds to the following row: alpar@9: * alpar@9: * sum a[i,j] pi[i] + a[p,j] pi[p] + lambda[j] = c[j], (9) alpar@9: * i!=p alpar@9: * alpar@9: * from which it follows that: alpar@9: * alpar@9: * lambda[j] = c[j] - sum a[i,j] pi[i] - a[p,j] pi[p]. (10) alpar@9: * i!=p alpar@9: * alpar@9: * In the transformed problem values of all multipliers pi[i] are known alpar@9: * (including pi[i], whose value is zero, since row p is inactive). alpar@9: * Thus, using formula (10) it is possible to compute values of alpar@9: * multipliers lambda[j] for all columns in row p. alpar@9: * alpar@9: * Note also that in the original problem all columns in row p are alpar@9: * bounded, not fixed. So status GLP_NS assigned to every such column alpar@9: * must be changed to GLP_NL or GLP_NU depending on which bound the alpar@9: * corresponding column has been fixed. This status change may lead to alpar@9: * dual feasibility violation for solution of the original problem, alpar@9: * because now column multipliers must satisfy to the following alpar@9: * condition: alpar@9: * alpar@9: * ( >= 0, if status of column j is GLP_NL, alpar@9: * lambda[j] < (11) alpar@9: * ( <= 0, if status of column j is GLP_NU. alpar@9: * alpar@9: * If this condition holds, solution to the original problem is the alpar@9: * same as to the transformed problem. Otherwise, we have to perform alpar@9: * one degenerate pivoting step of the primal simplex method to obtain alpar@9: * dual feasible (hence, optimal) solution to the original problem as alpar@9: * follows. If, on problem transformation, row p was made active on its alpar@9: * lower bound (case at = 0), we change its status to GLP_NL (or GLP_NS) alpar@9: * and start increasing its multiplier pi[p]. Otherwise, if row p was alpar@9: * made active on its upper bound (case at = 1), we change its status alpar@9: * to GLP_NU (or GLP_NS) and start decreasing pi[p]. From (10) it alpar@9: * follows that: alpar@9: * alpar@9: * delta lambda[j] = - a[p,j] * delta pi[p] = - a[p,j] pi[p]. (12) alpar@9: * alpar@9: * Simple analysis of formulae (3)-(5) shows that changing pi[p] in the alpar@9: * specified direction causes increasing lambda[j] for every column j alpar@9: * assigned status GLP_NL (delta lambda[j] > 0) and decreasing lambda[j] alpar@9: * for every column j assigned status GLP_NU (delta lambda[j] < 0). It alpar@9: * is understood that once the last lambda[q], which violates condition alpar@9: * (11), has reached zero, multipliers lambda[j] for all columns get alpar@9: * valid signs. Such column q can be determined as follows. Let d[j] be alpar@9: * initial value of lambda[j] (i.e. reduced cost of column j) in the alpar@9: * transformed problem computed with formula (10) when pi[p] = 0. Then alpar@9: * lambda[j] = d[j] + delta lambda[j], and from (12) it follows that alpar@9: * lambda[j] becomes zero if: alpar@9: * alpar@9: * delta lambda[j] = - a[p,j] pi[p] = - d[j] ==> alpar@9: * (13) alpar@9: * pi[p] = d[j] / a[p,j]. alpar@9: * alpar@9: * Therefore, the last column q, for which lambda[q] becomes zero, can alpar@9: * be determined from the following condition: alpar@9: * alpar@9: * |d[q] / a[p,q]| = max |pi[p]| = max |d[j] / a[p,j]|, (14) alpar@9: * j in D j in D alpar@9: * alpar@9: * where D is a set of columns j whose, reduced costs d[j] have invalid alpar@9: * signs, i.e. violate condition (11). (Thus, if D is empty, solution alpar@9: * to the original problem is the same as solution to the transformed alpar@9: * problem, and no correction is needed as was noticed above.) In alpar@9: * solution to the original problem column q is assigned status GLP_BS, alpar@9: * since it replaces column of auxiliary variable of row p (becoming alpar@9: * active) in the basis, and multiplier for row p is assigned its new alpar@9: * value, which is pi[p] = d[q] / a[p,q]. Note that due to primal alpar@9: * degeneracy values of all columns having non-zero coefficients in row alpar@9: * p remain unchanged. alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * Value of multiplier pi[p] in solution to the original problem is alpar@9: * corrected in the same way as for basic solution. Values of all alpar@9: * columns having non-zero coefficients in row p remain unchanged. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: struct forcing_col alpar@9: { /* column fixed on its bound by forcing row */ alpar@9: int j; alpar@9: /* column reference number */ alpar@9: char stat; alpar@9: /* original column status: alpar@9: GLP_NL - fixed on lower bound alpar@9: GLP_NU - fixed on upper bound */ alpar@9: double a; alpar@9: /* constraint coefficient a[p,j] */ alpar@9: double c; alpar@9: /* objective coefficient c[j] */ alpar@9: NPPLFE *ptr; alpar@9: /* list of non-zero coefficients a[i,j], i != p */ alpar@9: struct forcing_col *next; alpar@9: /* pointer to another column fixed by forcing row */ alpar@9: }; alpar@9: alpar@9: struct forcing_row alpar@9: { /* forcing row */ alpar@9: int p; alpar@9: /* row reference number */ alpar@9: char stat; alpar@9: /* status assigned to the row if it becomes active: alpar@9: GLP_NS - active equality constraint alpar@9: GLP_NL - inequality constraint with lower bound active alpar@9: GLP_NU - inequality constraint with upper bound active */ alpar@9: struct forcing_col *ptr; alpar@9: /* list of all columns having non-zero constraint coefficient alpar@9: a[p,j] in the forcing row */ alpar@9: }; alpar@9: alpar@9: static int rcv_forcing_row(NPP *npp, void *info); alpar@9: alpar@9: int npp_forcing_row(NPP *npp, NPPROW *p, int at) alpar@9: { /* process forcing row */ alpar@9: struct forcing_row *info; alpar@9: struct forcing_col *col = NULL; alpar@9: NPPCOL *j; alpar@9: NPPAIJ *apj, *aij; alpar@9: NPPLFE *lfe; alpar@9: double big; alpar@9: xassert(at == 0 || at == 1); alpar@9: /* determine maximal magnitude of the row coefficients */ alpar@9: big = 1.0; alpar@9: for (apj = p->ptr; apj != NULL; apj = apj->r_next) alpar@9: if (big < fabs(apj->val)) big = fabs(apj->val); alpar@9: /* if there are too small coefficients in the row, transformation alpar@9: should not be applied */ alpar@9: for (apj = p->ptr; apj != NULL; apj = apj->r_next) alpar@9: if (fabs(apj->val) < 1e-7 * big) return 1; alpar@9: /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_forcing_row, sizeof(struct forcing_row)); alpar@9: info->p = p->i; alpar@9: if (p->lb == p->ub) alpar@9: { /* equality constraint */ alpar@9: info->stat = GLP_NS; alpar@9: } alpar@9: else if (at == 0) alpar@9: { /* inequality constraint; case L[p] = U'[p] */ alpar@9: info->stat = GLP_NL; alpar@9: xassert(p->lb != -DBL_MAX); alpar@9: } alpar@9: else /* at == 1 */ alpar@9: { /* inequality constraint; case U[p] = L'[p] */ alpar@9: info->stat = GLP_NU; alpar@9: xassert(p->ub != +DBL_MAX); alpar@9: } alpar@9: info->ptr = NULL; alpar@9: /* scan the forcing row, fix columns at corresponding bounds, and alpar@9: save column information (the latter is not needed for MIP) */ alpar@9: for (apj = p->ptr; apj != NULL; apj = apj->r_next) alpar@9: { /* column j has non-zero coefficient in the forcing row */ alpar@9: j = apj->col; alpar@9: /* it must be non-fixed */ alpar@9: xassert(j->lb < j->ub); alpar@9: /* allocate stack entry to save column information */ alpar@9: if (npp->sol != GLP_MIP) alpar@9: { col = dmp_get_atom(npp->stack, sizeof(struct forcing_col)); alpar@9: col->j = j->j; alpar@9: col->stat = -1; /* will be set below */ alpar@9: col->a = apj->val; alpar@9: col->c = j->coef; alpar@9: col->ptr = NULL; alpar@9: col->next = info->ptr; alpar@9: info->ptr = col; alpar@9: } alpar@9: /* fix column j */ alpar@9: if (at == 0 && apj->val < 0.0 || at != 0 && apj->val > 0.0) alpar@9: { /* at its lower bound */ alpar@9: if (npp->sol != GLP_MIP) alpar@9: col->stat = GLP_NL; alpar@9: xassert(j->lb != -DBL_MAX); alpar@9: j->ub = j->lb; alpar@9: } alpar@9: else alpar@9: { /* at its upper bound */ alpar@9: if (npp->sol != GLP_MIP) alpar@9: col->stat = GLP_NU; alpar@9: xassert(j->ub != +DBL_MAX); alpar@9: j->lb = j->ub; alpar@9: } alpar@9: /* save column coefficients a[i,j], i != p */ alpar@9: if (npp->sol != GLP_MIP) alpar@9: { for (aij = j->ptr; aij != NULL; aij = aij->c_next) alpar@9: { if (aij == apj) continue; /* skip a[p,j] */ 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 = col->ptr; alpar@9: col->ptr = lfe; alpar@9: } alpar@9: } alpar@9: } alpar@9: /* make the row free (unbounded) */ alpar@9: p->lb = -DBL_MAX, p->ub = +DBL_MAX; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: static int rcv_forcing_row(NPP *npp, void *_info) alpar@9: { /* recover forcing row */ alpar@9: struct forcing_row *info = _info; alpar@9: struct forcing_col *col, *piv; alpar@9: NPPLFE *lfe; alpar@9: double d, big, temp; alpar@9: if (npp->sol == GLP_MIP) goto done; alpar@9: /* initially solution to the original problem is the same as alpar@9: to the transformed problem, where row p is inactive constraint alpar@9: with pi[p] = 0, and all columns are non-basic */ alpar@9: if (npp->sol == GLP_SOL) alpar@9: { if (npp->r_stat[info->p] != GLP_BS) alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: for (col = info->ptr; col != NULL; col = col->next) alpar@9: { if (npp->c_stat[col->j] != GLP_NS) alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: npp->c_stat[col->j] = col->stat; /* original status */ alpar@9: } alpar@9: } alpar@9: /* compute reduced costs d[j] for all columns with formula (10) alpar@9: and store them in col.c instead objective coefficients */ alpar@9: for (col = info->ptr; col != NULL; col = col->next) alpar@9: { d = col->c; alpar@9: for (lfe = col->ptr; lfe != NULL; lfe = lfe->next) alpar@9: d -= lfe->val * npp->r_pi[lfe->ref]; alpar@9: col->c = d; alpar@9: } alpar@9: /* consider columns j, whose multipliers lambda[j] has wrong alpar@9: sign in solution to the transformed problem (where lambda[j] = alpar@9: d[j]), and choose column q, whose multipler lambda[q] reaches alpar@9: zero last on changing row multiplier pi[p]; see (14) */ alpar@9: piv = NULL, big = 0.0; alpar@9: for (col = info->ptr; col != NULL; col = col->next) alpar@9: { d = col->c; /* d[j] */ alpar@9: temp = fabs(d / col->a); alpar@9: if (col->stat == GLP_NL) alpar@9: { /* column j has active lower bound */ alpar@9: if (d < 0.0 && big < temp) alpar@9: piv = col, big = temp; alpar@9: } alpar@9: else if (col->stat == GLP_NU) alpar@9: { /* column j has active upper bound */ alpar@9: if (d > 0.0 && big < temp) alpar@9: piv = col, big = temp; alpar@9: } alpar@9: else alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: } alpar@9: /* if column q does not exist, no correction is needed */ alpar@9: if (piv != NULL) alpar@9: { /* correct solution; row p becomes active constraint while alpar@9: column q becomes basic */ alpar@9: if (npp->sol == GLP_SOL) alpar@9: { npp->r_stat[info->p] = info->stat; alpar@9: npp->c_stat[piv->j] = GLP_BS; alpar@9: } alpar@9: /* assign new value to row multiplier pi[p] = d[p] / a[p,q] */ alpar@9: npp->r_pi[info->p] = piv->c / piv->a; alpar@9: } alpar@9: done: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_analyze_row - perform general row analysis alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * int npp_analyze_row(NPP *npp, NPPROW *p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_analyze_row performs analysis of row p of general alpar@9: * format: alpar@9: * alpar@9: * L[p] <= sum a[p,j] x[j] <= U[p], (1) alpar@9: * j alpar@9: * alpar@9: * l[j] <= x[j] <= u[j], (2) alpar@9: * alpar@9: * where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0. alpar@9: * alpar@9: * RETURNS alpar@9: * alpar@9: * 0x?0 - row lower bound does not exist or is redundant; alpar@9: * alpar@9: * 0x?1 - row lower bound can be active; alpar@9: * alpar@9: * 0x?2 - row lower bound is a forcing bound; alpar@9: * alpar@9: * 0x0? - row upper bound does not exist or is redundant; alpar@9: * alpar@9: * 0x1? - row upper bound can be active; alpar@9: * alpar@9: * 0x2? - row upper bound is a forcing bound; alpar@9: * alpar@9: * 0x33 - row bounds are inconsistent with column bounds. alpar@9: * alpar@9: * ALGORITHM alpar@9: * alpar@9: * Analysis of row (1) is based on analysis of its implied lower and alpar@9: * upper bounds, which are determined by bounds of corresponding columns alpar@9: * (variables) as follows: alpar@9: * alpar@9: * L'[p] = inf sum a[p,j] x[j] = alpar@9: * j alpar@9: * (3) alpar@9: * = sum a[p,j] l[j] + sum a[p,j] u[j], alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * U'[p] = sup sum a[p,j] x[j] = alpar@9: * (4) alpar@9: * = sum a[p,j] u[j] + sum a[p,j] l[j], alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5) alpar@9: * alpar@9: * (Note that bounds of all columns in row p are assumed to be correct, alpar@9: * so L'[p] <= U'[p].) alpar@9: * alpar@9: * Analysis of row lower bound L[p] includes the following cases: alpar@9: * alpar@9: * 1) if L[p] > U'[p] + eps, where eps is an absolute tolerance for row alpar@9: * value, row lower bound L[p] and implied row upper bound U'[p] are alpar@9: * inconsistent, ergo, the problem has no primal feasible solution; alpar@9: * alpar@9: * 2) if U'[p] - eps <= L[p] <= U'[p] + eps, i.e. if L[p] =~ U'[p], alpar@9: * the row is a forcing row on its lower bound (see description of alpar@9: * the routine npp_forcing_row); alpar@9: * alpar@9: * 3) if L[p] > L'[p] + eps, row lower bound L[p] can be active (this alpar@9: * conclusion does not account other rows in the problem); alpar@9: * alpar@9: * 4) if L[p] <= L'[p] + eps, row lower bound L[p] cannot be active, so alpar@9: * it is redundant and can be removed (replaced by -oo). alpar@9: * alpar@9: * Analysis of row upper bound U[p] is performed in a similar way and alpar@9: * includes the following cases: alpar@9: * alpar@9: * 1) if U[p] < L'[p] - eps, row upper bound U[p] and implied row lower alpar@9: * bound L'[p] are inconsistent, ergo the problem has no primal alpar@9: * feasible solution; alpar@9: * alpar@9: * 2) if L'[p] - eps <= U[p] <= L'[p] + eps, i.e. if U[p] =~ L'[p], alpar@9: * the row is a forcing row on its upper bound (see description of alpar@9: * the routine npp_forcing_row); alpar@9: * alpar@9: * 3) if U[p] < U'[p] - eps, row upper bound U[p] can be active (this alpar@9: * conclusion does not account other rows in the problem); alpar@9: * alpar@9: * 4) if U[p] >= U'[p] - eps, row upper bound U[p] cannot be active, so alpar@9: * it is redundant and can be removed (replaced by +oo). */ alpar@9: alpar@9: int npp_analyze_row(NPP *npp, NPPROW *p) alpar@9: { /* perform general row analysis */ alpar@9: NPPAIJ *aij; alpar@9: int ret = 0x00; alpar@9: double l, u, eps; alpar@9: xassert(npp == npp); alpar@9: /* compute implied lower bound L'[p]; see (3) */ alpar@9: l = 0.0; alpar@9: for (aij = p->ptr; aij != NULL; aij = aij->r_next) alpar@9: { if (aij->val > 0.0) alpar@9: { if (aij->col->lb == -DBL_MAX) alpar@9: { l = -DBL_MAX; alpar@9: break; alpar@9: } alpar@9: l += aij->val * aij->col->lb; alpar@9: } alpar@9: else /* aij->val < 0.0 */ alpar@9: { if (aij->col->ub == +DBL_MAX) alpar@9: { l = -DBL_MAX; alpar@9: break; alpar@9: } alpar@9: l += aij->val * aij->col->ub; alpar@9: } alpar@9: } alpar@9: /* compute implied upper bound U'[p]; see (4) */ alpar@9: u = 0.0; alpar@9: for (aij = p->ptr; aij != NULL; aij = aij->r_next) alpar@9: { if (aij->val > 0.0) alpar@9: { if (aij->col->ub == +DBL_MAX) alpar@9: { u = +DBL_MAX; alpar@9: break; alpar@9: } alpar@9: u += aij->val * aij->col->ub; alpar@9: } alpar@9: else /* aij->val < 0.0 */ alpar@9: { if (aij->col->lb == -DBL_MAX) alpar@9: { u = +DBL_MAX; alpar@9: break; alpar@9: } alpar@9: u += aij->val * aij->col->lb; alpar@9: } alpar@9: } alpar@9: /* column bounds are assumed correct, so L'[p] <= U'[p] */ alpar@9: /* check if row lower bound is consistent */ alpar@9: if (p->lb != -DBL_MAX) alpar@9: { eps = 1e-3 + 1e-6 * fabs(p->lb); alpar@9: if (p->lb - eps > u) alpar@9: { ret = 0x33; alpar@9: goto done; alpar@9: } alpar@9: } alpar@9: /* check if row upper bound is consistent */ alpar@9: if (p->ub != +DBL_MAX) alpar@9: { eps = 1e-3 + 1e-6 * fabs(p->ub); alpar@9: if (p->ub + eps < l) alpar@9: { ret = 0x33; alpar@9: goto done; alpar@9: } alpar@9: } alpar@9: /* check if row lower bound can be active/forcing */ alpar@9: if (p->lb != -DBL_MAX) alpar@9: { eps = 1e-9 + 1e-12 * fabs(p->lb); alpar@9: if (p->lb - eps > l) alpar@9: { if (p->lb + eps <= u) alpar@9: ret |= 0x01; alpar@9: else alpar@9: ret |= 0x02; alpar@9: } alpar@9: } alpar@9: /* check if row upper bound can be active/forcing */ alpar@9: if (p->ub != +DBL_MAX) alpar@9: { eps = 1e-9 + 1e-12 * fabs(p->ub); alpar@9: if (p->ub + eps < u) alpar@9: { /* check if the upper bound is forcing */ alpar@9: if (p->ub - eps >= l) alpar@9: ret |= 0x10; alpar@9: else alpar@9: ret |= 0x20; alpar@9: } alpar@9: } alpar@9: done: return ret; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_inactive_bound - remove row lower/upper inactive bound alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_inactive_bound(NPP *npp, NPPROW *p, int which); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_inactive_bound removes lower (if which = 0) or upper alpar@9: * (if which = 1) bound of row p: alpar@9: * alpar@9: * L[p] <= sum a[p,j] x[j] <= U[p], alpar@9: * alpar@9: * which (bound) is assumed to be redundant. alpar@9: * alpar@9: * PROBLEM TRANSFORMATION alpar@9: * alpar@9: * If which = 0, current lower bound L[p] of row p is assigned -oo. alpar@9: * If which = 1, current upper bound U[p] of row p is assigned +oo. alpar@9: * alpar@9: * RECOVERING BASIC SOLUTION alpar@9: * alpar@9: * If in solution to the transformed problem row p is inactive alpar@9: * constraint (GLP_BS), its status is not changed in solution to the alpar@9: * original problem. Otherwise, status of row p in solution to the alpar@9: * original problem is defined by its type before transformation and alpar@9: * its status in solution to the transformed problem as follows: alpar@9: * alpar@9: * +---------------------+-------+---------------+---------------+ alpar@9: * | Row | Flag | Row status in | Row status in | alpar@9: * | type | which | transfmd soln | original soln | alpar@9: * +---------------------+-------+---------------+---------------+ alpar@9: * | sum >= L[p] | 0 | GLP_NF | GLP_NL | alpar@9: * | sum <= U[p] | 1 | GLP_NF | GLP_NU | alpar@9: * | L[p] <= sum <= U[p] | 0 | GLP_NU | GLP_NU | alpar@9: * | L[p] <= sum <= U[p] | 1 | GLP_NL | GLP_NL | alpar@9: * | sum = L[p] = U[p] | 0 | GLP_NU | GLP_NS | alpar@9: * | sum = L[p] = U[p] | 1 | GLP_NL | GLP_NS | alpar@9: * +---------------------+-------+---------------+---------------+ alpar@9: * alpar@9: * RECOVERING INTERIOR-POINT SOLUTION alpar@9: * alpar@9: * None needed. alpar@9: * alpar@9: * RECOVERING MIP SOLUTION alpar@9: * alpar@9: * None needed. */ alpar@9: alpar@9: struct inactive_bound alpar@9: { /* row inactive bound */ alpar@9: int p; alpar@9: /* row reference number */ alpar@9: char stat; alpar@9: /* row status (if active constraint) */ alpar@9: }; alpar@9: alpar@9: static int rcv_inactive_bound(NPP *npp, void *info); alpar@9: alpar@9: void npp_inactive_bound(NPP *npp, NPPROW *p, int which) alpar@9: { /* remove row lower/upper inactive bound */ alpar@9: struct inactive_bound *info; alpar@9: if (npp->sol == GLP_SOL) alpar@9: { /* create transformation stack entry */ alpar@9: info = npp_push_tse(npp, alpar@9: rcv_inactive_bound, sizeof(struct inactive_bound)); alpar@9: info->p = p->i; alpar@9: if (p->ub == +DBL_MAX) alpar@9: info->stat = GLP_NL; alpar@9: else if (p->lb == -DBL_MAX) alpar@9: info->stat = GLP_NU; alpar@9: else if (p->lb != p->ub) alpar@9: info->stat = (char)(which == 0 ? GLP_NU : GLP_NL); alpar@9: else alpar@9: info->stat = GLP_NS; alpar@9: } alpar@9: /* remove row inactive bound */ alpar@9: if (which == 0) alpar@9: { xassert(p->lb != -DBL_MAX); alpar@9: p->lb = -DBL_MAX; alpar@9: } alpar@9: else if (which == 1) alpar@9: { xassert(p->ub != +DBL_MAX); alpar@9: p->ub = +DBL_MAX; alpar@9: } alpar@9: else alpar@9: xassert(which != which); alpar@9: return; alpar@9: } alpar@9: alpar@9: static int rcv_inactive_bound(NPP *npp, void *_info) alpar@9: { /* recover row status */ alpar@9: struct inactive_bound *info = _info; alpar@9: if (npp->sol != GLP_SOL) alpar@9: { npp_error(); alpar@9: return 1; alpar@9: } alpar@9: if (npp->r_stat[info->p] == GLP_BS) alpar@9: npp->r_stat[info->p] = GLP_BS; alpar@9: else alpar@9: npp->r_stat[info->p] = info->stat; alpar@9: return 0; alpar@9: } alpar@9: alpar@9: /*********************************************************************** alpar@9: * NAME alpar@9: * alpar@9: * npp_implied_bounds - determine implied column bounds alpar@9: * alpar@9: * SYNOPSIS alpar@9: * alpar@9: * #include "glpnpp.h" alpar@9: * void npp_implied_bounds(NPP *npp, NPPROW *p); alpar@9: * alpar@9: * DESCRIPTION alpar@9: * alpar@9: * The routine npp_implied_bounds inspects general row (constraint) p: alpar@9: * alpar@9: * L[p] <= sum a[p,j] x[j] <= U[p], (1) alpar@9: * alpar@9: * l[j] <= x[j] <= u[j], (2) alpar@9: * alpar@9: * where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0, to compute alpar@9: * implied bounds of columns (variables x[j]) in this row. alpar@9: * alpar@9: * The routine stores implied column bounds l'[j] and u'[j] in column alpar@9: * descriptors (NPPCOL); it does not change current column bounds l[j] alpar@9: * and u[j]. (Implied column bounds can be then used to strengthen the alpar@9: * current column bounds; see the routines npp_implied_lower and alpar@9: * npp_implied_upper). alpar@9: * alpar@9: * ALGORITHM alpar@9: * alpar@9: * Current column bounds (2) define implied lower and upper bounds of alpar@9: * row (1) as follows: alpar@9: * alpar@9: * L'[p] = inf sum a[p,j] x[j] = alpar@9: * j alpar@9: * (3) alpar@9: * = sum a[p,j] l[j] + sum a[p,j] u[j], alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * U'[p] = sup sum a[p,j] x[j] = alpar@9: * (4) alpar@9: * = sum a[p,j] u[j] + sum a[p,j] l[j], alpar@9: * j in Jp j in Jn alpar@9: * alpar@9: * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5) alpar@9: * alpar@9: * (Note that bounds of all columns in row p are assumed to be correct, alpar@9: * so L'[p] <= U'[p].) alpar@9: * alpar@9: * If L[p] > L'[p] and/or U[p] < U'[p], the lower and/or upper bound of alpar@9: * row (1) can be active, in which case such row defines implied bounds alpar@9: * of its variables. alpar@9: * alpar@9: * Let x[k] be some variable having in row (1) coefficient a[p,k] != 0. alpar@9: * Consider a case when row lower bound can be active (L[p] > L'[p]): alpar@9: * alpar@9: * sum a[p,j] x[j] >= L[p] ==> alpar@9: * j alpar@9: * alpar@9: * sum a[p,j] x[j] + a[p,k] x[k] >= L[p] ==> alpar@9: * j!=k alpar@9: * (6) alpar@9: * a[p,k] x[k] >= L[p] - sum a[p,j] x[j] ==> alpar@9: * j!=k alpar@9: * alpar@9: * a[p,k] x[k] >= L[p,k], alpar@9: * alpar@9: * where alpar@9: * alpar@9: * L[p,k] = inf(L[p] - sum a[p,j] x[j]) = alpar@9: * j!=k alpar@9: * alpar@9: * = L[p] - sup sum a[p,j] x[j] = (7) alpar@9: * j!=k alpar@9: * alpar@9: * = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j]. alpar@9: * j in Jp\{k} j in Jn\{k} alpar@9: * alpar@9: * Thus: alpar@9: * alpar@9: * x[k] >= l'[k] = L[p,k] / a[p,k], if a[p,k] > 0, (8) alpar@9: * alpar@9: * x[k] <= u'[k] = L[p,k] / a[p,k], if a[p,k] < 0. (9) alpar@9: * alpar@9: * where l'[k] and u'[k] are implied lower and upper bounds of variable alpar@9: * x[k], resp. alpar@9: * alpar@9: * Now consider a similar case when row upper bound can be active alpar@9: * (U[p] < U'[p]): alpar@9: * alpar@9: * sum a[p,j] x[j] <= U[p] ==> alpar@9: * j alpar@9: * alpar@9: * sum a[p,j] x[j] + a[p,k] x[k] <= U[p] ==> alpar@9: * j!=k alpar@9: * (10) alpar@9: * a[p,k] x[k] <= U[p] - sum a[p,j] x[j] ==> alpar@9: * j!=k alpar@9: * alpar@9: * a[p,k] x[k] <= U[p,k], alpar@9: * alpar@9: * where: alpar@9: * alpar@9: * U[p,k] = sup(U[p] - sum a[p,j] x[j]) = alpar@9: * j!=k alpar@9: * alpar@9: * = U[p] - inf sum a[p,j] x[j] = (11) alpar@9: * j!=k alpar@9: * alpar@9: * = U[p] - sum a[p,j] l[j] - sum a[p,j] u[j]. alpar@9: * j in Jp\{k} j in Jn\{k} alpar@9: * alpar@9: * Thus: alpar@9: * alpar@9: * x[k] <= u'[k] = U[p,k] / a[p,k], if a[p,k] > 0, (12) alpar@9: * alpar@9: * x[k] >= l'[k] = U[p,k] / a[p,k], if a[p,k] < 0. (13) alpar@9: * alpar@9: * Note that in formulae (8), (9), (12), and (13) coefficient a[p,k] alpar@9: * must not be too small in magnitude relatively to other non-zero alpar@9: * coefficients in row (1), i.e. the following condition must hold: alpar@9: * alpar@9: * |a[p,k]| >= eps * max(1, |a[p,j]|), (14) alpar@9: * j alpar@9: * alpar@9: * where eps is a relative tolerance for constraint coefficients. alpar@9: * Otherwise the implied column bounds can be numerical inreliable. For alpar@9: * example, using formula (8) for the following inequality constraint: alpar@9: * alpar@9: * 1e-12 x1 - x2 - x3 >= 0, alpar@9: * alpar@9: * where x1 >= -1, x2, x3, >= 0, may lead to numerically unreliable alpar@9: * conclusion that x1 >= 0. alpar@9: * alpar@9: * Using formulae (8), (9), (12), and (13) to compute implied bounds alpar@9: * for one variable requires |J| operations, where J = {j: a[p,j] != 0}, alpar@9: * because this needs computing L[p,k] and U[p,k]. Thus, computing alpar@9: * implied bounds for all variables in row (1) would require |J|^2 alpar@9: * operations, that is not a good technique. However, the total number alpar@9: * of operations can be reduced to |J| as follows. alpar@9: * alpar@9: * Let a[p,k] > 0. Then from (7) and (11) we have: alpar@9: * alpar@9: * L[p,k] = L[p] - (U'[p] - a[p,k] u[k]) = alpar@9: * alpar@9: * = L[p] - U'[p] + a[p,k] u[k], alpar@9: * alpar@9: * U[p,k] = U[p] - (L'[p] - a[p,k] l[k]) = alpar@9: * alpar@9: * = U[p] - L'[p] + a[p,k] l[k], alpar@9: * alpar@9: * where L'[p] and U'[p] are implied row lower and upper bounds defined alpar@9: * by formulae (3) and (4). Substituting these expressions into (8) and alpar@9: * (12) gives: alpar@9: * alpar@9: * l'[k] = L[p,k] / a[p,k] = u[k] + (L[p] - U'[p]) / a[p,k], (15) alpar@9: * alpar@9: * u'[k] = U[p,k] / a[p,k] = l[k] + (U[p] - L'[p]) / a[p,k]. (16) alpar@9: * alpar@9: * Similarly, if a[p,k] < 0, according to (7) and (11) we have: alpar@9: * alpar@9: * L[p,k] = L[p] - (U'[p] - a[p,k] l[k]) = alpar@9: * alpar@9: * = L[p] - U'[p] + a[p,k] l[k], alpar@9: * alpar@9: * U[p,k] = U[p] - (L'[p] - a[p,k] u[k]) = alpar@9: * alpar@9: * = U[p] - L'[p] + a[p,k] u[k], alpar@9: * alpar@9: * and substituting these expressions into (8) and (12) gives: alpar@9: * alpar@9: * l'[k] = U[p,k] / a[p,k] = u[k] + (U[p] - L'[p]) / a[p,k], (17) alpar@9: * alpar@9: * u'[k] = L[p,k] / a[p,k] = l[k] + (L[p] - U'[p]) / a[p,k]. (18) alpar@9: * alpar@9: * Note that formulae (15)-(18) can be used only if L'[p] and U'[p] alpar@9: * exist. However, if for some variable x[j] it happens that l[j] = -oo alpar@9: * and/or u[j] = +oo, values of L'[p] (if a[p,j] > 0) and/or U'[p] (if alpar@9: * a[p,j] < 0) are undefined. Consider, therefore, the most general alpar@9: * situation, when some column bounds (2) may not exist. alpar@9: * alpar@9: * Let: alpar@9: * alpar@9: * J' = {j : (a[p,j] > 0 and l[j] = -oo) or alpar@9: * (19) alpar@9: * (a[p,j] < 0 and u[j] = +oo)}. alpar@9: * alpar@9: * Then (assuming that row upper bound U[p] can be active) the following alpar@9: * three cases are possible: alpar@9: * alpar@9: * 1) |J'| = 0. In this case L'[p] exists, thus, for all variables x[j] alpar@9: * in row (1) we can use formulae (16) and (17); alpar@9: * alpar@9: * 2) J' = {k}. In this case L'[p] = -oo, however, U[p,k] (11) exists, alpar@9: * so for variable x[k] we can use formulae (12) and (13). Note that alpar@9: * for all other variables x[j] (j != k) l'[j] = -oo (if a[p,j] < 0) alpar@9: * or u'[j] = +oo (if a[p,j] > 0); alpar@9: * alpar@9: * 3) |J'| > 1. In this case for all variables x[j] in row [1] we have alpar@9: * l'[j] = -oo (if a[p,j] < 0) or u'[j] = +oo (if a[p,j] > 0). alpar@9: * alpar@9: * Similarly, let: alpar@9: * alpar@9: * J'' = {j : (a[p,j] > 0 and u[j] = +oo) or alpar@9: * (20) alpar@9: * (a[p,j] < 0 and l[j] = -oo)}. alpar@9: * alpar@9: * Then (assuming that row lower bound L[p] can be active) the following alpar@9: * three cases are possible: alpar@9: * alpar@9: * 1) |J''| = 0. In this case U'[p] exists, thus, for all variables x[j] alpar@9: * in row (1) we can use formulae (15) and (18); alpar@9: * alpar@9: * 2) J'' = {k}. In this case U'[p] = +oo, however, L[p,k] (7) exists, alpar@9: * so for variable x[k] we can use formulae (8) and (9). Note that alpar@9: * for all other variables x[j] (j != k) l'[j] = -oo (if a[p,j] > 0) alpar@9: * or u'[j] = +oo (if a[p,j] < 0); alpar@9: * alpar@9: * 3) |J''| > 1. In this case for all variables x[j] in row (1) we have alpar@9: * l'[j] = -oo (if a[p,j] > 0) or u'[j] = +oo (if a[p,j] < 0). */ alpar@9: alpar@9: void npp_implied_bounds(NPP *npp, NPPROW *p) alpar@9: { NPPAIJ *apj, *apk; alpar@9: double big, eps, temp; alpar@9: xassert(npp == npp); alpar@9: /* initialize implied bounds for all variables and determine alpar@9: maximal magnitude of row coefficients a[p,j] */ alpar@9: big = 1.0; alpar@9: for (apj = p->ptr; apj != NULL; apj = apj->r_next) alpar@9: { apj->col->ll.ll = -DBL_MAX, apj->col->uu.uu = +DBL_MAX; alpar@9: if (big < fabs(apj->val)) big = fabs(apj->val); alpar@9: } alpar@9: eps = 1e-6 * big; alpar@9: /* process row lower bound (assuming that it can be active) */ alpar@9: if (p->lb != -DBL_MAX) alpar@9: { apk = NULL; alpar@9: for (apj = p->ptr; apj != NULL; apj = apj->r_next) alpar@9: { if (apj->val > 0.0 && apj->col->ub == +DBL_MAX || alpar@9: apj->val < 0.0 && apj->col->lb == -DBL_MAX) alpar@9: { if (apk == NULL) alpar@9: apk = apj; alpar@9: else alpar@9: goto skip1; alpar@9: } alpar@9: } alpar@9: /* if a[p,k] = NULL then |J'| = 0 else J' = { k } */ alpar@9: temp = p->lb; alpar@9: for (apj = p->ptr; apj != NULL; apj = apj->r_next) alpar@9: { if (apj == apk) alpar@9: /* skip a[p,k] */; alpar@9: else if (apj->val > 0.0) alpar@9: temp -= apj->val * apj->col->ub; alpar@9: else /* apj->val < 0.0 */ alpar@9: temp -= apj->val * apj->col->lb; alpar@9: } alpar@9: /* compute column implied bounds */ alpar@9: if (apk == NULL) alpar@9: { /* temp = L[p] - U'[p] */ alpar@9: for (apj = p->ptr; apj != NULL; apj = apj->r_next) alpar@9: { if (apj->val >= +eps) alpar@9: { /* l'[j] := u[j] + (L[p] - U'[p]) / a[p,j] */ alpar@9: apj->col->ll.ll = apj->col->ub + temp / apj->val; alpar@9: } alpar@9: else if (apj->val <= -eps) alpar@9: { /* u'[j] := l[j] + (L[p] - U'[p]) / a[p,j] */ alpar@9: apj->col->uu.uu = apj->col->lb + temp / apj->val; alpar@9: } alpar@9: } alpar@9: } alpar@9: else alpar@9: { /* temp = L[p,k] */ alpar@9: if (apk->val >= +eps) alpar@9: { /* l'[k] := L[p,k] / a[p,k] */ alpar@9: apk->col->ll.ll = temp / apk->val; alpar@9: } alpar@9: else if (apk->val <= -eps) alpar@9: { /* u'[k] := L[p,k] / a[p,k] */ alpar@9: apk->col->uu.uu = temp / apk->val; alpar@9: } alpar@9: } alpar@9: skip1: ; alpar@9: } alpar@9: /* process row upper bound (assuming that it can be active) */ alpar@9: if (p->ub != +DBL_MAX) alpar@9: { apk = NULL; alpar@9: for (apj = p->ptr; apj != NULL; apj = apj->r_next) alpar@9: { if (apj->val > 0.0 && apj->col->lb == -DBL_MAX || alpar@9: apj->val < 0.0 && apj->col->ub == +DBL_MAX) alpar@9: { if (apk == NULL) alpar@9: apk = apj; alpar@9: else alpar@9: goto skip2; alpar@9: } alpar@9: } alpar@9: /* if a[p,k] = NULL then |J''| = 0 else J'' = { k } */ alpar@9: temp = p->ub; alpar@9: for (apj = p->ptr; apj != NULL; apj = apj->r_next) alpar@9: { if (apj == apk) alpar@9: /* skip a[p,k] */; alpar@9: else if (apj->val > 0.0) alpar@9: temp -= apj->val * apj->col->lb; alpar@9: else /* apj->val < 0.0 */ alpar@9: temp -= apj->val * apj->col->ub; alpar@9: } alpar@9: /* compute column implied bounds */ alpar@9: if (apk == NULL) alpar@9: { /* temp = U[p] - L'[p] */ alpar@9: for (apj = p->ptr; apj != NULL; apj = apj->r_next) alpar@9: { if (apj->val >= +eps) alpar@9: { /* u'[j] := l[j] + (U[p] - L'[p]) / a[p,j] */ alpar@9: apj->col->uu.uu = apj->col->lb + temp / apj->val; alpar@9: } alpar@9: else if (apj->val <= -eps) alpar@9: { /* l'[j] := u[j] + (U[p] - L'[p]) / a[p,j] */ alpar@9: apj->col->ll.ll = apj->col->ub + temp / apj->val; alpar@9: } alpar@9: } alpar@9: } alpar@9: else alpar@9: { /* temp = U[p,k] */ alpar@9: if (apk->val >= +eps) alpar@9: { /* u'[k] := U[p,k] / a[p,k] */ alpar@9: apk->col->uu.uu = temp / apk->val; alpar@9: } alpar@9: else if (apk->val <= -eps) alpar@9: { /* l'[k] := U[p,k] / a[p,k] */ alpar@9: apk->col->ll.ll = temp / apk->val; alpar@9: } alpar@9: } alpar@9: skip2: ; alpar@9: } alpar@9: return; alpar@9: } alpar@9: alpar@9: /* eof */