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.
479 lines
16 KiB
479 lines
16 KiB
/* glpssx02.c (simplex method, rational arithmetic) */
|
|
|
|
/***********************************************************************
|
|
* This code is part of GLPK (GNU Linear Programming Kit).
|
|
*
|
|
* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
|
|
* 2009, 2010, 2011, 2013 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 "glpssx.h"
|
|
|
|
static void show_progress(SSX *ssx, int phase)
|
|
{ /* this auxiliary routine displays information about progress of
|
|
the search */
|
|
int i, def = 0;
|
|
for (i = 1; i <= ssx->m; i++)
|
|
if (ssx->type[ssx->Q_col[i]] == SSX_FX) def++;
|
|
xprintf("%s%6d: %s = %22.15g (%d)\n", phase == 1 ? " " : "*",
|
|
ssx->it_cnt, phase == 1 ? "infsum" : "objval",
|
|
mpq_get_d(ssx->bbar[0]), def);
|
|
#if 0
|
|
ssx->tm_lag = utime();
|
|
#else
|
|
ssx->tm_lag = xtime();
|
|
#endif
|
|
return;
|
|
}
|
|
|
|
/*----------------------------------------------------------------------
|
|
// ssx_phase_I - find primal feasible solution.
|
|
//
|
|
// This routine implements phase I of the primal simplex method.
|
|
//
|
|
// On exit the routine returns one of the following codes:
|
|
//
|
|
// 0 - feasible solution found;
|
|
// 1 - problem has no feasible solution;
|
|
// 2 - iterations limit exceeded;
|
|
// 3 - time limit exceeded.
|
|
----------------------------------------------------------------------*/
|
|
|
|
int ssx_phase_I(SSX *ssx)
|
|
{ int m = ssx->m;
|
|
int n = ssx->n;
|
|
int *type = ssx->type;
|
|
mpq_t *lb = ssx->lb;
|
|
mpq_t *ub = ssx->ub;
|
|
mpq_t *coef = ssx->coef;
|
|
int *A_ptr = ssx->A_ptr;
|
|
int *A_ind = ssx->A_ind;
|
|
mpq_t *A_val = ssx->A_val;
|
|
int *Q_col = ssx->Q_col;
|
|
mpq_t *bbar = ssx->bbar;
|
|
mpq_t *pi = ssx->pi;
|
|
mpq_t *cbar = ssx->cbar;
|
|
int *orig_type, orig_dir;
|
|
mpq_t *orig_lb, *orig_ub, *orig_coef;
|
|
int i, k, ret;
|
|
/* save components of the original LP problem, which are changed
|
|
by the routine */
|
|
orig_type = xcalloc(1+m+n, sizeof(int));
|
|
orig_lb = xcalloc(1+m+n, sizeof(mpq_t));
|
|
orig_ub = xcalloc(1+m+n, sizeof(mpq_t));
|
|
orig_coef = xcalloc(1+m+n, sizeof(mpq_t));
|
|
for (k = 1; k <= m+n; k++)
|
|
{ orig_type[k] = type[k];
|
|
mpq_init(orig_lb[k]);
|
|
mpq_set(orig_lb[k], lb[k]);
|
|
mpq_init(orig_ub[k]);
|
|
mpq_set(orig_ub[k], ub[k]);
|
|
}
|
|
orig_dir = ssx->dir;
|
|
for (k = 0; k <= m+n; k++)
|
|
{ mpq_init(orig_coef[k]);
|
|
mpq_set(orig_coef[k], coef[k]);
|
|
}
|
|
/* build an artificial basic solution, which is primal feasible,
|
|
and also build an auxiliary objective function to minimize the
|
|
sum of infeasibilities for the original problem */
|
|
ssx->dir = SSX_MIN;
|
|
for (k = 0; k <= m+n; k++) mpq_set_si(coef[k], 0, 1);
|
|
mpq_set_si(bbar[0], 0, 1);
|
|
for (i = 1; i <= m; i++)
|
|
{ int t;
|
|
k = Q_col[i]; /* x[k] = xB[i] */
|
|
t = type[k];
|
|
if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
|
|
{ /* in the original problem x[k] has lower bound */
|
|
if (mpq_cmp(bbar[i], lb[k]) < 0)
|
|
{ /* which is violated */
|
|
type[k] = SSX_UP;
|
|
mpq_set(ub[k], lb[k]);
|
|
mpq_set_si(lb[k], 0, 1);
|
|
mpq_set_si(coef[k], -1, 1);
|
|
mpq_add(bbar[0], bbar[0], ub[k]);
|
|
mpq_sub(bbar[0], bbar[0], bbar[i]);
|
|
}
|
|
}
|
|
if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
|
|
{ /* in the original problem x[k] has upper bound */
|
|
if (mpq_cmp(bbar[i], ub[k]) > 0)
|
|
{ /* which is violated */
|
|
type[k] = SSX_LO;
|
|
mpq_set(lb[k], ub[k]);
|
|
mpq_set_si(ub[k], 0, 1);
|
|
mpq_set_si(coef[k], +1, 1);
|
|
mpq_add(bbar[0], bbar[0], bbar[i]);
|
|
mpq_sub(bbar[0], bbar[0], lb[k]);
|
|
}
|
|
}
|
|
}
|
|
/* now the initial basic solution should be primal feasible due
|
|
to changes of bounds of some basic variables, which turned to
|
|
implicit artifical variables */
|
|
/* compute simplex multipliers and reduced costs */
|
|
ssx_eval_pi(ssx);
|
|
ssx_eval_cbar(ssx);
|
|
/* display initial progress of the search */
|
|
show_progress(ssx, 1);
|
|
/* main loop starts here */
|
|
for (;;)
|
|
{ /* display current progress of the search */
|
|
#if 0
|
|
if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
|
|
#else
|
|
if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
|
|
#endif
|
|
show_progress(ssx, 1);
|
|
/* we do not need to wait until all artificial variables have
|
|
left the basis */
|
|
if (mpq_sgn(bbar[0]) == 0)
|
|
{ /* the sum of infeasibilities is zero, therefore the current
|
|
solution is primal feasible for the original problem */
|
|
ret = 0;
|
|
break;
|
|
}
|
|
/* check if the iterations limit has been exhausted */
|
|
if (ssx->it_lim == 0)
|
|
{ ret = 2;
|
|
break;
|
|
}
|
|
/* check if the time limit has been exhausted */
|
|
#if 0
|
|
if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
|
|
#else
|
|
if (ssx->tm_lim >= 0.0 &&
|
|
ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
|
|
#endif
|
|
{ ret = 3;
|
|
break;
|
|
}
|
|
/* choose non-basic variable xN[q] */
|
|
ssx_chuzc(ssx);
|
|
/* if xN[q] cannot be chosen, the sum of infeasibilities is
|
|
minimal but non-zero; therefore the original problem has no
|
|
primal feasible solution */
|
|
if (ssx->q == 0)
|
|
{ ret = 1;
|
|
break;
|
|
}
|
|
/* compute q-th column of the simplex table */
|
|
ssx_eval_col(ssx);
|
|
/* choose basic variable xB[p] */
|
|
ssx_chuzr(ssx);
|
|
/* the sum of infeasibilities cannot be negative, therefore
|
|
the auxiliary lp problem cannot have unbounded solution */
|
|
xassert(ssx->p != 0);
|
|
/* update values of basic variables */
|
|
ssx_update_bbar(ssx);
|
|
if (ssx->p > 0)
|
|
{ /* compute p-th row of the inverse inv(B) */
|
|
ssx_eval_rho(ssx);
|
|
/* compute p-th row of the simplex table */
|
|
ssx_eval_row(ssx);
|
|
xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
|
|
/* update simplex multipliers */
|
|
ssx_update_pi(ssx);
|
|
/* update reduced costs of non-basic variables */
|
|
ssx_update_cbar(ssx);
|
|
}
|
|
/* xB[p] is leaving the basis; if it is implicit artificial
|
|
variable, the corresponding residual vanishes; therefore
|
|
bounds of this variable should be restored to the original
|
|
values */
|
|
if (ssx->p > 0)
|
|
{ k = Q_col[ssx->p]; /* x[k] = xB[p] */
|
|
if (type[k] != orig_type[k])
|
|
{ /* x[k] is implicit artificial variable */
|
|
type[k] = orig_type[k];
|
|
mpq_set(lb[k], orig_lb[k]);
|
|
mpq_set(ub[k], orig_ub[k]);
|
|
xassert(ssx->p_stat == SSX_NL || ssx->p_stat == SSX_NU);
|
|
ssx->p_stat = (ssx->p_stat == SSX_NL ? SSX_NU : SSX_NL);
|
|
if (type[k] == SSX_FX) ssx->p_stat = SSX_NS;
|
|
/* nullify the objective coefficient at x[k] */
|
|
mpq_set_si(coef[k], 0, 1);
|
|
/* since coef[k] has been changed, we need to compute
|
|
new reduced cost of x[k], which it will have in the
|
|
adjacent basis */
|
|
/* the formula d[j] = cN[j] - pi' * N[j] is used (note
|
|
that the vector pi is not changed, because it depends
|
|
on objective coefficients at basic variables, but in
|
|
the adjacent basis, for which the vector pi has been
|
|
just recomputed, x[k] is non-basic) */
|
|
if (k <= m)
|
|
{ /* x[k] is auxiliary variable */
|
|
mpq_neg(cbar[ssx->q], pi[k]);
|
|
}
|
|
else
|
|
{ /* x[k] is structural variable */
|
|
int ptr;
|
|
mpq_t temp;
|
|
mpq_init(temp);
|
|
mpq_set_si(cbar[ssx->q], 0, 1);
|
|
for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
|
|
{ mpq_mul(temp, pi[A_ind[ptr]], A_val[ptr]);
|
|
mpq_add(cbar[ssx->q], cbar[ssx->q], temp);
|
|
}
|
|
mpq_clear(temp);
|
|
}
|
|
}
|
|
}
|
|
/* jump to the adjacent vertex of the polyhedron */
|
|
ssx_change_basis(ssx);
|
|
/* one simplex iteration has been performed */
|
|
if (ssx->it_lim > 0) ssx->it_lim--;
|
|
ssx->it_cnt++;
|
|
}
|
|
/* display final progress of the search */
|
|
show_progress(ssx, 1);
|
|
/* restore components of the original problem, which were changed
|
|
by the routine */
|
|
for (k = 1; k <= m+n; k++)
|
|
{ type[k] = orig_type[k];
|
|
mpq_set(lb[k], orig_lb[k]);
|
|
mpq_clear(orig_lb[k]);
|
|
mpq_set(ub[k], orig_ub[k]);
|
|
mpq_clear(orig_ub[k]);
|
|
}
|
|
ssx->dir = orig_dir;
|
|
for (k = 0; k <= m+n; k++)
|
|
{ mpq_set(coef[k], orig_coef[k]);
|
|
mpq_clear(orig_coef[k]);
|
|
}
|
|
xfree(orig_type);
|
|
xfree(orig_lb);
|
|
xfree(orig_ub);
|
|
xfree(orig_coef);
|
|
/* return to the calling program */
|
|
return ret;
|
|
}
|
|
|
|
/*----------------------------------------------------------------------
|
|
// ssx_phase_II - find optimal solution.
|
|
//
|
|
// This routine implements phase II of the primal simplex method.
|
|
//
|
|
// On exit the routine returns one of the following codes:
|
|
//
|
|
// 0 - optimal solution found;
|
|
// 1 - problem has unbounded solution;
|
|
// 2 - iterations limit exceeded;
|
|
// 3 - time limit exceeded.
|
|
----------------------------------------------------------------------*/
|
|
|
|
int ssx_phase_II(SSX *ssx)
|
|
{ int ret;
|
|
/* display initial progress of the search */
|
|
show_progress(ssx, 2);
|
|
/* main loop starts here */
|
|
for (;;)
|
|
{ /* display current progress of the search */
|
|
#if 0
|
|
if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
|
|
#else
|
|
if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
|
|
#endif
|
|
show_progress(ssx, 2);
|
|
/* check if the iterations limit has been exhausted */
|
|
if (ssx->it_lim == 0)
|
|
{ ret = 2;
|
|
break;
|
|
}
|
|
/* check if the time limit has been exhausted */
|
|
#if 0
|
|
if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
|
|
#else
|
|
if (ssx->tm_lim >= 0.0 &&
|
|
ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
|
|
#endif
|
|
{ ret = 3;
|
|
break;
|
|
}
|
|
/* choose non-basic variable xN[q] */
|
|
ssx_chuzc(ssx);
|
|
/* if xN[q] cannot be chosen, the current basic solution is
|
|
dual feasible and therefore optimal */
|
|
if (ssx->q == 0)
|
|
{ ret = 0;
|
|
break;
|
|
}
|
|
/* compute q-th column of the simplex table */
|
|
ssx_eval_col(ssx);
|
|
/* choose basic variable xB[p] */
|
|
ssx_chuzr(ssx);
|
|
/* if xB[p] cannot be chosen, the problem has no dual feasible
|
|
solution (i.e. unbounded) */
|
|
if (ssx->p == 0)
|
|
{ ret = 1;
|
|
break;
|
|
}
|
|
/* update values of basic variables */
|
|
ssx_update_bbar(ssx);
|
|
if (ssx->p > 0)
|
|
{ /* compute p-th row of the inverse inv(B) */
|
|
ssx_eval_rho(ssx);
|
|
/* compute p-th row of the simplex table */
|
|
ssx_eval_row(ssx);
|
|
xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
|
|
#if 0
|
|
/* update simplex multipliers */
|
|
ssx_update_pi(ssx);
|
|
#endif
|
|
/* update reduced costs of non-basic variables */
|
|
ssx_update_cbar(ssx);
|
|
}
|
|
/* jump to the adjacent vertex of the polyhedron */
|
|
ssx_change_basis(ssx);
|
|
/* one simplex iteration has been performed */
|
|
if (ssx->it_lim > 0) ssx->it_lim--;
|
|
ssx->it_cnt++;
|
|
}
|
|
/* display final progress of the search */
|
|
show_progress(ssx, 2);
|
|
/* return to the calling program */
|
|
return ret;
|
|
}
|
|
|
|
/*----------------------------------------------------------------------
|
|
// ssx_driver - base driver to exact simplex method.
|
|
//
|
|
// This routine is a base driver to a version of the primal simplex
|
|
// method using exact (bignum) arithmetic.
|
|
//
|
|
// On exit the routine returns one of the following codes:
|
|
//
|
|
// 0 - optimal solution found;
|
|
// 1 - problem has no feasible solution;
|
|
// 2 - problem has unbounded solution;
|
|
// 3 - iterations limit exceeded (phase I);
|
|
// 4 - iterations limit exceeded (phase II);
|
|
// 5 - time limit exceeded (phase I);
|
|
// 6 - time limit exceeded (phase II);
|
|
// 7 - initial basis matrix is exactly singular.
|
|
----------------------------------------------------------------------*/
|
|
|
|
int ssx_driver(SSX *ssx)
|
|
{ int m = ssx->m;
|
|
int *type = ssx->type;
|
|
mpq_t *lb = ssx->lb;
|
|
mpq_t *ub = ssx->ub;
|
|
int *Q_col = ssx->Q_col;
|
|
mpq_t *bbar = ssx->bbar;
|
|
int i, k, ret;
|
|
ssx->tm_beg = xtime();
|
|
/* factorize the initial basis matrix */
|
|
if (ssx_factorize(ssx))
|
|
{ xprintf("Initial basis matrix is singular\n");
|
|
ret = 7;
|
|
goto done;
|
|
}
|
|
/* compute values of basic variables */
|
|
ssx_eval_bbar(ssx);
|
|
/* check if the initial basic solution is primal feasible */
|
|
for (i = 1; i <= m; i++)
|
|
{ int t;
|
|
k = Q_col[i]; /* x[k] = xB[i] */
|
|
t = type[k];
|
|
if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
|
|
{ /* x[k] has lower bound */
|
|
if (mpq_cmp(bbar[i], lb[k]) < 0)
|
|
{ /* which is violated */
|
|
break;
|
|
}
|
|
}
|
|
if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
|
|
{ /* x[k] has upper bound */
|
|
if (mpq_cmp(bbar[i], ub[k]) > 0)
|
|
{ /* which is violated */
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
if (i > m)
|
|
{ /* no basic variable violates its bounds */
|
|
ret = 0;
|
|
goto skip;
|
|
}
|
|
/* phase I: find primal feasible solution */
|
|
ret = ssx_phase_I(ssx);
|
|
switch (ret)
|
|
{ case 0:
|
|
ret = 0;
|
|
break;
|
|
case 1:
|
|
xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
|
|
ret = 1;
|
|
break;
|
|
case 2:
|
|
xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
|
|
ret = 3;
|
|
break;
|
|
case 3:
|
|
xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
|
|
ret = 5;
|
|
break;
|
|
default:
|
|
xassert(ret != ret);
|
|
}
|
|
/* compute values of basic variables (actually only the objective
|
|
value needs to be computed) */
|
|
ssx_eval_bbar(ssx);
|
|
skip: /* compute simplex multipliers */
|
|
ssx_eval_pi(ssx);
|
|
/* compute reduced costs of non-basic variables */
|
|
ssx_eval_cbar(ssx);
|
|
/* if phase I failed, do not start phase II */
|
|
if (ret != 0) goto done;
|
|
/* phase II: find optimal solution */
|
|
ret = ssx_phase_II(ssx);
|
|
switch (ret)
|
|
{ case 0:
|
|
xprintf("OPTIMAL SOLUTION FOUND\n");
|
|
ret = 0;
|
|
break;
|
|
case 1:
|
|
xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
|
|
ret = 2;
|
|
break;
|
|
case 2:
|
|
xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
|
|
ret = 4;
|
|
break;
|
|
case 3:
|
|
xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
|
|
ret = 6;
|
|
break;
|
|
default:
|
|
xassert(ret != ret);
|
|
}
|
|
done: /* decrease the time limit by the spent amount of time */
|
|
if (ssx->tm_lim >= 0.0)
|
|
#if 0
|
|
{ ssx->tm_lim -= utime() - ssx->tm_beg;
|
|
#else
|
|
{ ssx->tm_lim -= xdifftime(xtime(), ssx->tm_beg);
|
|
#endif
|
|
if (ssx->tm_lim < 0.0) ssx->tm_lim = 0.0;
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
/* eof */
|