You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
388 lines
14 KiB
388 lines
14 KiB
/* spxchuzr.c */
|
|
|
|
/***********************************************************************
|
|
* This code is part of GLPK (GNU Linear Programming Kit).
|
|
*
|
|
* Copyright (C) 2015 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;
|
|
}
|
|
|
|
/* eof */
|