diff options
Diffstat (limited to 'test/monniaux/glpk-4.65/src/simplex/spxchuzr.c')
| -rw-r--r-- | test/monniaux/glpk-4.65/src/simplex/spxchuzr.c | 594 |
1 files changed, 594 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/simplex/spxchuzr.c b/test/monniaux/glpk-4.65/src/simplex/spxchuzr.c new file mode 100644 index 00000000..8bef77ba --- /dev/null +++ b/test/monniaux/glpk-4.65/src/simplex/spxchuzr.c @@ -0,0 +1,594 @@ +/* spxchuzr.c */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2015-2018 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 "spxchuzr.h" + +/*********************************************************************** +* spx_chuzr_std - choose basic variable (textbook ratio test) +* +* This routine implements an improved textbook ratio test to choose +* basic variable xB[p]. +* +* The parameter phase specifies the search phase: +* +* 1 - searching for feasible basic solution. In this case the routine +* uses artificial bounds of basic variables that correspond to +* breakpoints of the penalty function: +* +* ( lB[i], if cB[i] = 0 +* ( +* lB'[i] = { uB[i], if cB[i] > 0 +* ( +* ( -inf, if cB[i] < 0 +* +* ( uB[i], if cB[i] = 0 +* ( +* uB'[i] = { +inf, if cB[i] > 0 +* ( +* ( lB[i], if cB[i] < 0 +* +* where lB[i] and uB[i] are original bounds of variable xB[i], +* cB[i] is the penalty (objective) coefficient of that variable. +* +* 2 - searching for optimal basic solution. In this case the routine +* uses original bounds of basic variables. +* +* Current values of basic variables should be placed in the array +* locations beta[1], ..., beta[m]. +* +* The parameter 1 <= q <= n-m specifies the index of non-basic +* variable xN[q] chosen. +* +* The parameter s specifies the direction in which xN[q] changes: +* s = +1.0 means xN[q] increases, and s = -1.0 means xN[q] decreases. +* (Thus, the corresponding ray parameter is theta = s (xN[q] - f[q]), +* where f[q] is the active bound of xN[q] in the current basis.) +* +* Elements of q-th simplex table column T[q] = (t[i,q]) corresponding +* to non-basic variable xN[q] should be placed in the array locations +* tcol[1], ..., tcol[m]. +* +* The parameter tol_piv specifies a tolerance for elements of the +* simplex table column T[q]. If |t[i,q]| < tol_piv, basic variable +* xB[i] is skipped, i.e. it is assumed that it does not depend on the +* ray parameter theta. +* +* The parameters tol and tol1 specify tolerances used to increase the +* choice freedom by simulating an artificial degeneracy as follows. +* If beta[i] <= lB[i] + delta[i], where delta[i] = tol + tol1 |lB[i]|, +* it is assumed that beta[i] is exactly the same as lB[i]. Similarly, +* if beta[i] >= uB[i] - delta[i], where delta[i] = tol + tol1 |uB[i]|, +* it is assumed that beta[i] is exactly the same as uB[i]. +* +* The routine determines the index 1 <= p <= m of basic variable xB[p] +* that reaches its (lower or upper) bound first on increasing the ray +* parameter theta, stores the bound flag (0 - lower bound or fixed +* value, 1 - upper bound) to the location pointed to by the pointer +* p_flag, and returns the index p. If non-basic variable xN[q] is +* double-bounded and reaches its opposite bound first, the routine +* returns (-1). And if the ray parameter may increase unlimitedly, the +* routine returns zero. +* +* Should note that the bound flag stored to the location pointed to by +* p_flag corresponds to the original (not artficial) bound of variable +* xB[p] and defines the active bound flag lp->flag[q] to be set in the +* adjacent basis for that basic variable. */ + +int spx_chuzr_std(SPXLP *lp, int phase, const double beta[/*1+m*/], + int q, double s, const double tcol[/*1+m*/], int *p_flag, + double tol_piv, double tol, double tol1) +{ int m = lp->m; + int n = lp->n; + double *c = lp->c; + double *l = lp->l; + double *u = lp->u; + int *head = lp->head; + int i, i_flag, k, p; + double alfa, biga, delta, lk, uk, teta, teta_min; + xassert(phase == 1 || phase == 2); + xassert(1 <= q && q <= n-m); + xassert(s == +1.0 || s == -1.0); + /* determine initial teta_min */ + k = head[m+q]; /* x[k] = xN[q] */ + if (l[k] == -DBL_MAX || u[k] == +DBL_MAX) + { /* xN[q] has no opposite bound */ + p = 0, *p_flag = 0, teta_min = DBL_MAX, biga = 0.0; + } + else + { /* xN[q] have both lower and upper bounds */ + p = -1, *p_flag = 0, teta_min = fabs(l[k] - u[k]), biga = 1.0; + } + /* walk thru the list of basic variables */ + for (i = 1; i <= m; i++) + { k = head[i]; /* x[k] = xB[i] */ + /* determine alfa such that delta xB[i] = alfa * teta */ + alfa = s * tcol[i]; + if (alfa <= -tol_piv) + { /* xB[i] decreases */ + /* determine actual lower bound of xB[i] */ + if (phase == 1 && c[k] < 0.0) + { /* xB[i] has no actual lower bound */ + continue; + } + else if (phase == 1 && c[k] > 0.0) + { /* actual lower bound of xB[i] is its upper bound */ + lk = u[k]; + xassert(lk != +DBL_MAX); + i_flag = 1; + } + else + { /* actual lower bound of xB[i] is its original bound */ + lk = l[k]; + if (lk == -DBL_MAX) + continue; + i_flag = 0; + } + /* determine teta on which xB[i] reaches its lower bound */ + delta = tol + tol1 * (lk >= 0.0 ? +lk : -lk); + if (beta[i] <= lk + delta) + teta = 0.0; + else + teta = (lk - beta[i]) / alfa; + } + else if (alfa >= +tol_piv) + { /* xB[i] increases */ + /* determine actual upper bound of xB[i] */ + if (phase == 1 && c[k] < 0.0) + { /* actual upper bound of xB[i] is its lower bound */ + uk = l[k]; + xassert(uk != -DBL_MAX); + i_flag = 0; + } + else if (phase == 1 && c[k] > 0.0) + { /* xB[i] has no actual upper bound */ + continue; + } + else + { /* actual upper bound of xB[i] is its original bound */ + uk = u[k]; + if (uk == +DBL_MAX) + continue; + i_flag = 1; + } + /* determine teta on which xB[i] reaches its upper bound */ + delta = tol + tol1 * (uk >= 0.0 ? +uk : -uk); + if (beta[i] >= uk - delta) + teta = 0.0; + else + teta = (uk - beta[i]) / alfa; + } + else + { /* xB[i] does not depend on teta */ + continue; + } + /* choose basic variable xB[p] for which teta is minimal */ + xassert(teta >= 0.0); + alfa = (alfa >= 0.0 ? +alfa : -alfa); + if (teta_min > teta || (teta_min == teta && biga < alfa)) + p = i, *p_flag = i_flag, teta_min = teta, biga = alfa; + } + /* if xB[p] is fixed variable, adjust its bound flag */ + if (p > 0) + { k = head[p]; + if (l[k] == u[k]) + *p_flag = 0; + } + return p; +} + +/*********************************************************************** +* spx_chuzr_harris - choose basic variable (Harris' ratio test) +* +* This routine implements Harris' ratio test to choose basic variable +* xB[p]. +* +* All the parameters, except tol and tol1, as well as the returned +* value have the same meaning as for the routine spx_chuzr_std (see +* above). +* +* The parameters tol and tol1 specify tolerances on bound violations +* for basic variables. For the lower bound of basic variable xB[i] the +* tolerance is delta[i] = tol + tol1 |lB[i]|, and for the upper bound +* the tolerance is delta[i] = tol + tol1 |uB[i]|. */ + +int spx_chuzr_harris(SPXLP *lp, int phase, const double beta[/*1+m*/], + int q, double s, const double tcol[/*1+m*/], int *p_flag, + double tol_piv, double tol, double tol1) +{ int m = lp->m; + int n = lp->n; + double *c = lp->c; + double *l = lp->l; + double *u = lp->u; + int *head = lp->head; + int i, i_flag, k, p; + double alfa, biga, delta, lk, uk, teta, teta_min; + xassert(phase == 1 || phase == 2); + xassert(1 <= q && q <= n-m); + xassert(s == +1.0 || s == -1.0); + /*--------------------------------------------------------------*/ + /* first pass: determine teta_min for relaxed bounds */ + /*--------------------------------------------------------------*/ + teta_min = DBL_MAX; + /* walk thru the list of basic variables */ + for (i = 1; i <= m; i++) + { k = head[i]; /* x[k] = xB[i] */ + /* determine alfa such that delta xB[i] = alfa * teta */ + alfa = s * tcol[i]; + if (alfa <= -tol_piv) + { /* xB[i] decreases */ + /* determine actual lower bound of xB[i] */ + if (phase == 1 && c[k] < 0.0) + { /* xB[i] has no actual lower bound */ + continue; + } + else if (phase == 1 && c[k] > 0.0) + { /* actual lower bound of xB[i] is its upper bound */ + lk = u[k]; + xassert(lk != +DBL_MAX); + } + else + { /* actual lower bound of xB[i] is its original bound */ + lk = l[k]; + if (lk == -DBL_MAX) + continue; + } + /* determine teta on which xB[i] reaches its relaxed lower + * bound */ + delta = tol + tol1 * (lk >= 0.0 ? +lk : -lk); + if (beta[i] < lk) + teta = - delta / alfa; + else + teta = ((lk - delta) - beta[i]) / alfa; + } + else if (alfa >= +tol_piv) + { /* xB[i] increases */ + /* determine actual upper bound of xB[i] */ + if (phase == 1 && c[k] < 0.0) + { /* actual upper bound of xB[i] is its lower bound */ + uk = l[k]; + xassert(uk != -DBL_MAX); + } + else if (phase == 1 && c[k] > 0.0) + { /* xB[i] has no actual upper bound */ + continue; + } + else + { /* actual upper bound of xB[i] is its original bound */ + uk = u[k]; + if (uk == +DBL_MAX) + continue; + } + /* determine teta on which xB[i] reaches its relaxed upper + * bound */ + delta = tol + tol1 * (uk >= 0.0 ? +uk : -uk); + if (beta[i] > uk) + teta = + delta / alfa; + else + teta = ((uk + delta) - beta[i]) / alfa; + } + else + { /* xB[i] does not depend on teta */ + continue; + } + xassert(teta >= 0.0); + if (teta_min > teta) + teta_min = teta; + } + /*--------------------------------------------------------------*/ + /* second pass: choose basic variable xB[p] */ + /*--------------------------------------------------------------*/ + k = head[m+q]; /* x[k] = xN[q] */ + if (l[k] != -DBL_MAX && u[k] != +DBL_MAX) + { /* xN[q] has both lower and upper bounds */ + if (fabs(l[k] - u[k]) <= teta_min) + { /* and reaches its opposite bound */ + p = -1, *p_flag = 0; + goto done; + } + } + if (teta_min == DBL_MAX) + { /* teta may increase unlimitedly */ + p = 0, *p_flag = 0; + goto done; + } + /* nothing is chosen so far */ + p = 0, *p_flag = 0, biga = 0.0; + /* walk thru the list of basic variables */ + for (i = 1; i <= m; i++) + { k = head[i]; /* x[k] = xB[i] */ + /* determine alfa such that delta xB[i] = alfa * teta */ + alfa = s * tcol[i]; + if (alfa <= -tol_piv) + { /* xB[i] decreases */ + /* determine actual lower bound of xB[i] */ + if (phase == 1 && c[k] < 0.0) + { /* xB[i] has no actual lower bound */ + continue; + } + else if (phase == 1 && c[k] > 0.0) + { /* actual lower bound of xB[i] is its upper bound */ + lk = u[k]; + xassert(lk != +DBL_MAX); + i_flag = 1; + } + else + { /* actual lower bound of xB[i] is its original bound */ + lk = l[k]; + if (lk == -DBL_MAX) + continue; + i_flag = 0; + } + /* determine teta on which xB[i] reaches its lower bound */ + teta = (lk - beta[i]) / alfa; + } + else if (alfa >= +tol_piv) + { /* xB[i] increases */ + /* determine actual upper bound of xB[i] */ + if (phase == 1 && c[k] < 0.0) + { /* actual upper bound of xB[i] is its lower bound */ + uk = l[k]; + xassert(uk != -DBL_MAX); + i_flag = 0; + } + else if (phase == 1 && c[k] > 0.0) + { /* xB[i] has no actual upper bound */ + continue; + } + else + { /* actual upper bound of xB[i] is its original bound */ + uk = u[k]; + if (uk == +DBL_MAX) + continue; + i_flag = 1; + } + /* determine teta on which xB[i] reaches its upper bound */ + teta = (uk - beta[i]) / alfa; + } + else + { /* xB[i] does not depend on teta */ + continue; + } + /* choose basic variable for which teta is not greater than + * teta_min determined for relaxed bounds and which has best + * (largest in magnitude) pivot */ + alfa = (alfa >= 0.0 ? +alfa : -alfa); + if (teta <= teta_min && biga < alfa) + p = i, *p_flag = i_flag, biga = alfa; + } + /* something must be chosen */ + xassert(1 <= p && p <= m); + /* if xB[p] is fixed variable, adjust its bound flag */ + k = head[p]; + if (l[k] == u[k]) + *p_flag = 0; +done: return p; +} + +#if 1 /* 22/VI-2017 */ +/*********************************************************************** +* spx_ls_eval_bp - determine penalty function break points +* +* This routine determines break points of the penalty function (which +* is the sum of primal infeasibilities). +* +* The parameters lp, beta, q, dq, tcol, and tol_piv have the same +* meaning as for the routine spx_chuzr_std (see above). +* +* The routine stores the break-points determined to the array elements +* bp[1], ..., bp[nbp] in *arbitrary* order, where 0 <= nbp <= 2*m+1 is +* the number of break-points returned by the routine on exit. */ + +int spx_ls_eval_bp(SPXLP *lp, const double beta[/*1+m*/], + int q, double dq, const double tcol[/*1+m*/], double tol_piv, + SPXBP bp[/*1+2*m+1*/]) +{ int m = lp->m; + int n = lp->n; + double *c = lp->c; + double *l = lp->l; + double *u = lp->u; + int *head = lp->head; + int i, k, nbp; + double s, alfa; + xassert(1 <= q && q <= n-m); + xassert(dq != 0.0); + s = (dq < 0.0 ? +1.0 : -1.0); + nbp = 0; + /* if chosen non-basic variable xN[q] is double-bounded, include + * it in the list, because it can cross its opposite bound */ + k = head[m+q]; /* x[k] = xN[q] */ + if (l[k] != -DBL_MAX && u[k] != +DBL_MAX) + { nbp++; + bp[nbp].i = 0; + xassert(l[k] < u[k]); /* xN[q] cannot be fixed */ + bp[nbp].teta = u[k] - l[k]; + bp[nbp].dc = s; + } + /* build the list of all basic variables xB[i] that can cross + * their bound(s) for the ray parameter 0 <= teta < teta_max */ + for (i = 1; i <= m; i++) + { k = head[i]; /* x[k] = xB[i] */ + xassert(l[k] <= u[k]); + /* determine alfa such that (delta xB[i]) = alfa * teta */ + alfa = s * tcol[i]; + if (alfa >= +tol_piv) + { /* xB[i] increases on increasing teta */ + if (l[k] == u[k]) + { /* xB[i] is fixed at lB[i] = uB[i] */ + if (c[k] <= 0.0) + { /* increasing xB[i] can cross its fixed value lB[i], + * because currently xB[i] <= lB[i] */ + nbp++; + bp[nbp].i = +i; + bp[nbp].teta = (l[k] - beta[i]) / alfa; + /* if xB[i] > lB[i] then cB[i] = +1 */ + bp[nbp].dc = +1.0 - c[k]; + } + } + else + { if (l[k] != -DBL_MAX && c[k] < 0.0) + { /* increasing xB[i] can cross its lower bound lB[i], + * because currently xB[i] < lB[i] */ + nbp++; + bp[nbp].i = +i; + bp[nbp].teta = (l[k] - beta[i]) / alfa; + bp[nbp].dc = +1.0; + } + if (u[k] != +DBL_MAX && c[k] <= 0.0) + { /* increasing xB[i] can cross its upper bound uB[i], + * because currently xB[i] does not violate it */ + nbp++; + bp[nbp].i = -i; + bp[nbp].teta = (u[k] - beta[i]) / alfa; + bp[nbp].dc = +1.0; + } + } + } + else if (alfa <= -tol_piv) + { /* xB[i] decreases on increasing teta */ + if (l[k] == u[k]) + { /* xB[i] is fixed at lB[i] = uB[i] */ + if (c[k] >= 0.0) + { /* decreasing xB[i] can cross its fixed value lB[i], + * because currently xB[i] >= lB[i] */ + nbp++; + bp[nbp].i = +i; + bp[nbp].teta = (l[k] - beta[i]) / alfa; + /* if xB[i] < lB[i] then cB[i] = -1 */ + bp[nbp].dc = -1.0 - c[k]; + } + } + else + { if (l[k] != -DBL_MAX && c[k] >= 0.0) + { /* decreasing xB[i] can cross its lower bound lB[i], + * because currently xB[i] does not violate it */ + nbp++; + bp[nbp].i = +i; + bp[nbp].teta = (l[k] - beta[i]) / alfa; + bp[nbp].dc = -1.0; + } + if (u[k] != +DBL_MAX && c[k] > 0.0) + { /* decreasing xB[i] can cross its upper bound uB[i], + * because currently xB[i] > uB[i] */ + nbp++; + bp[nbp].i = -i; + bp[nbp].teta = (u[k] - beta[i]) / alfa; + bp[nbp].dc = -1.0; + } + } + } + else + { /* xB[i] does not depend on teta within a tolerance */ + continue; + } + /* teta < 0 may happen only due to round-off errors when the + * current value of xB[i] is *close* to its (lower or upper) + * bound; in this case we replace teta by exact zero */ + if (bp[nbp].teta < 0.0) + bp[nbp].teta = 0.0; + } + xassert(nbp <= 2*m+1); + return nbp; +} +#endif + +#if 1 /* 22/VI-2017 */ +/*********************************************************************** +* spx_ls_select_bp - select and process penalty function break points +* +* This routine selects a next portion of the penalty function break +* points and processes them. +* +* On entry to the routine it is assumed that break points bp[1], ..., +* bp[num] are already processed, and slope is the penalty function +* slope to the right of the last processed break point bp[num]. +* (Initially, when num = 0, slope should be specified as -fabs(d[q]), +* where d[q] is the reduced cost of chosen non-basic variable xN[q].) +* +* The routine selects break points among bp[num+1], ..., bp[nbp], for +* which teta <= teta_lim, and moves these break points to the array +* elements bp[num+1], ..., bp[num1], where num <= num1 <= 2*m+1 is the +* new number of processed break points returned by the routine on +* exit. Then the routine sorts the break points by ascending teta and +* computes the change of the penalty function relative to its value at +* teta = 0. +* +* On exit the routine also replaces the parameter slope with a new +* value that corresponds to the new last break-point bp[num1]. */ + +static int CDECL fcmp(const void *v1, const void *v2) +{ const SPXBP *p1 = v1, *p2 = v2; + if (p1->teta < p2->teta) + return -1; + else if (p1->teta > p2->teta) + return +1; + else + return 0; +} + +int spx_ls_select_bp(SPXLP *lp, const double tcol[/*1+m*/], + int nbp, SPXBP bp[/*1+m+m+1*/], int num, double *slope, double + teta_lim) +{ int m = lp->m; + int i, t, num1; + double teta, dz; + xassert(0 <= num && num <= nbp && nbp <= m+m+1); + /* select a new portion of break points */ + num1 = num; + for (t = num+1; t <= nbp; t++) + { if (bp[t].teta <= teta_lim) + { /* move break point to the beginning of the new portion */ + num1++; + i = bp[num1].i, teta = bp[num1].teta, dz = bp[num1].dc; + bp[num1].i = bp[t].i, bp[num1].teta = bp[t].teta, + bp[num1].dc = bp[t].dc; + bp[t].i = i, bp[t].teta = teta, bp[t].dc = dz; + } + } + /* sort new break points bp[num+1], ..., bp[num1] by ascending + * the ray parameter teta */ + if (num1 - num > 1) + qsort(&bp[num+1], num1 - num, sizeof(SPXBP), fcmp); + /* calculate the penalty function change at the new break points + * selected */ + for (t = num+1; t <= num1; t++) + { /* calculate the penalty function change relative to its value + * at break point bp[t-1] */ + dz = (*slope) * (bp[t].teta - (t == 1 ? 0.0 : bp[t-1].teta)); + /* calculate the penalty function change relative to its value + * at teta = 0 */ + bp[t].dz = (t == 1 ? 0.0 : bp[t-1].dz) + dz; + /* calculate a new slope of the penalty function to the right + * of the current break point bp[t] */ + i = (bp[t].i >= 0 ? bp[t].i : -bp[t].i); + xassert(0 <= i && i <= m); + if (i == 0) + *slope += fabs(1.0 * bp[t].dc); + else + *slope += fabs(tcol[i] * bp[t].dc); + } + return num1; +} +#endif + +/* eof */ |