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