1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpssx02.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,479 @@
1.4 +/* glpssx02.c */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpenv.h"
1.29 +#include "glpssx.h"
1.30 +
1.31 +static void show_progress(SSX *ssx, int phase)
1.32 +{ /* this auxiliary routine displays information about progress of
1.33 + the search */
1.34 + int i, def = 0;
1.35 + for (i = 1; i <= ssx->m; i++)
1.36 + if (ssx->type[ssx->Q_col[i]] == SSX_FX) def++;
1.37 + xprintf("%s%6d: %s = %22.15g (%d)\n", phase == 1 ? " " : "*",
1.38 + ssx->it_cnt, phase == 1 ? "infsum" : "objval",
1.39 + mpq_get_d(ssx->bbar[0]), def);
1.40 +#if 0
1.41 + ssx->tm_lag = utime();
1.42 +#else
1.43 + ssx->tm_lag = xtime();
1.44 +#endif
1.45 + return;
1.46 +}
1.47 +
1.48 +/*----------------------------------------------------------------------
1.49 +// ssx_phase_I - find primal feasible solution.
1.50 +//
1.51 +// This routine implements phase I of the primal simplex method.
1.52 +//
1.53 +// On exit the routine returns one of the following codes:
1.54 +//
1.55 +// 0 - feasible solution found;
1.56 +// 1 - problem has no feasible solution;
1.57 +// 2 - iterations limit exceeded;
1.58 +// 3 - time limit exceeded.
1.59 +----------------------------------------------------------------------*/
1.60 +
1.61 +int ssx_phase_I(SSX *ssx)
1.62 +{ int m = ssx->m;
1.63 + int n = ssx->n;
1.64 + int *type = ssx->type;
1.65 + mpq_t *lb = ssx->lb;
1.66 + mpq_t *ub = ssx->ub;
1.67 + mpq_t *coef = ssx->coef;
1.68 + int *A_ptr = ssx->A_ptr;
1.69 + int *A_ind = ssx->A_ind;
1.70 + mpq_t *A_val = ssx->A_val;
1.71 + int *Q_col = ssx->Q_col;
1.72 + mpq_t *bbar = ssx->bbar;
1.73 + mpq_t *pi = ssx->pi;
1.74 + mpq_t *cbar = ssx->cbar;
1.75 + int *orig_type, orig_dir;
1.76 + mpq_t *orig_lb, *orig_ub, *orig_coef;
1.77 + int i, k, ret;
1.78 + /* save components of the original LP problem, which are changed
1.79 + by the routine */
1.80 + orig_type = xcalloc(1+m+n, sizeof(int));
1.81 + orig_lb = xcalloc(1+m+n, sizeof(mpq_t));
1.82 + orig_ub = xcalloc(1+m+n, sizeof(mpq_t));
1.83 + orig_coef = xcalloc(1+m+n, sizeof(mpq_t));
1.84 + for (k = 1; k <= m+n; k++)
1.85 + { orig_type[k] = type[k];
1.86 + mpq_init(orig_lb[k]);
1.87 + mpq_set(orig_lb[k], lb[k]);
1.88 + mpq_init(orig_ub[k]);
1.89 + mpq_set(orig_ub[k], ub[k]);
1.90 + }
1.91 + orig_dir = ssx->dir;
1.92 + for (k = 0; k <= m+n; k++)
1.93 + { mpq_init(orig_coef[k]);
1.94 + mpq_set(orig_coef[k], coef[k]);
1.95 + }
1.96 + /* build an artificial basic solution, which is primal feasible,
1.97 + and also build an auxiliary objective function to minimize the
1.98 + sum of infeasibilities for the original problem */
1.99 + ssx->dir = SSX_MIN;
1.100 + for (k = 0; k <= m+n; k++) mpq_set_si(coef[k], 0, 1);
1.101 + mpq_set_si(bbar[0], 0, 1);
1.102 + for (i = 1; i <= m; i++)
1.103 + { int t;
1.104 + k = Q_col[i]; /* x[k] = xB[i] */
1.105 + t = type[k];
1.106 + if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
1.107 + { /* in the original problem x[k] has lower bound */
1.108 + if (mpq_cmp(bbar[i], lb[k]) < 0)
1.109 + { /* which is violated */
1.110 + type[k] = SSX_UP;
1.111 + mpq_set(ub[k], lb[k]);
1.112 + mpq_set_si(lb[k], 0, 1);
1.113 + mpq_set_si(coef[k], -1, 1);
1.114 + mpq_add(bbar[0], bbar[0], ub[k]);
1.115 + mpq_sub(bbar[0], bbar[0], bbar[i]);
1.116 + }
1.117 + }
1.118 + if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
1.119 + { /* in the original problem x[k] has upper bound */
1.120 + if (mpq_cmp(bbar[i], ub[k]) > 0)
1.121 + { /* which is violated */
1.122 + type[k] = SSX_LO;
1.123 + mpq_set(lb[k], ub[k]);
1.124 + mpq_set_si(ub[k], 0, 1);
1.125 + mpq_set_si(coef[k], +1, 1);
1.126 + mpq_add(bbar[0], bbar[0], bbar[i]);
1.127 + mpq_sub(bbar[0], bbar[0], lb[k]);
1.128 + }
1.129 + }
1.130 + }
1.131 + /* now the initial basic solution should be primal feasible due
1.132 + to changes of bounds of some basic variables, which turned to
1.133 + implicit artifical variables */
1.134 + /* compute simplex multipliers and reduced costs */
1.135 + ssx_eval_pi(ssx);
1.136 + ssx_eval_cbar(ssx);
1.137 + /* display initial progress of the search */
1.138 + show_progress(ssx, 1);
1.139 + /* main loop starts here */
1.140 + for (;;)
1.141 + { /* display current progress of the search */
1.142 +#if 0
1.143 + if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
1.144 +#else
1.145 + if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
1.146 +#endif
1.147 + show_progress(ssx, 1);
1.148 + /* we do not need to wait until all artificial variables have
1.149 + left the basis */
1.150 + if (mpq_sgn(bbar[0]) == 0)
1.151 + { /* the sum of infeasibilities is zero, therefore the current
1.152 + solution is primal feasible for the original problem */
1.153 + ret = 0;
1.154 + break;
1.155 + }
1.156 + /* check if the iterations limit has been exhausted */
1.157 + if (ssx->it_lim == 0)
1.158 + { ret = 2;
1.159 + break;
1.160 + }
1.161 + /* check if the time limit has been exhausted */
1.162 +#if 0
1.163 + if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
1.164 +#else
1.165 + if (ssx->tm_lim >= 0.0 &&
1.166 + ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
1.167 +#endif
1.168 + { ret = 3;
1.169 + break;
1.170 + }
1.171 + /* choose non-basic variable xN[q] */
1.172 + ssx_chuzc(ssx);
1.173 + /* if xN[q] cannot be chosen, the sum of infeasibilities is
1.174 + minimal but non-zero; therefore the original problem has no
1.175 + primal feasible solution */
1.176 + if (ssx->q == 0)
1.177 + { ret = 1;
1.178 + break;
1.179 + }
1.180 + /* compute q-th column of the simplex table */
1.181 + ssx_eval_col(ssx);
1.182 + /* choose basic variable xB[p] */
1.183 + ssx_chuzr(ssx);
1.184 + /* the sum of infeasibilities cannot be negative, therefore
1.185 + the auxiliary lp problem cannot have unbounded solution */
1.186 + xassert(ssx->p != 0);
1.187 + /* update values of basic variables */
1.188 + ssx_update_bbar(ssx);
1.189 + if (ssx->p > 0)
1.190 + { /* compute p-th row of the inverse inv(B) */
1.191 + ssx_eval_rho(ssx);
1.192 + /* compute p-th row of the simplex table */
1.193 + ssx_eval_row(ssx);
1.194 + xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
1.195 + /* update simplex multipliers */
1.196 + ssx_update_pi(ssx);
1.197 + /* update reduced costs of non-basic variables */
1.198 + ssx_update_cbar(ssx);
1.199 + }
1.200 + /* xB[p] is leaving the basis; if it is implicit artificial
1.201 + variable, the corresponding residual vanishes; therefore
1.202 + bounds of this variable should be restored to the original
1.203 + values */
1.204 + if (ssx->p > 0)
1.205 + { k = Q_col[ssx->p]; /* x[k] = xB[p] */
1.206 + if (type[k] != orig_type[k])
1.207 + { /* x[k] is implicit artificial variable */
1.208 + type[k] = orig_type[k];
1.209 + mpq_set(lb[k], orig_lb[k]);
1.210 + mpq_set(ub[k], orig_ub[k]);
1.211 + xassert(ssx->p_stat == SSX_NL || ssx->p_stat == SSX_NU);
1.212 + ssx->p_stat = (ssx->p_stat == SSX_NL ? SSX_NU : SSX_NL);
1.213 + if (type[k] == SSX_FX) ssx->p_stat = SSX_NS;
1.214 + /* nullify the objective coefficient at x[k] */
1.215 + mpq_set_si(coef[k], 0, 1);
1.216 + /* since coef[k] has been changed, we need to compute
1.217 + new reduced cost of x[k], which it will have in the
1.218 + adjacent basis */
1.219 + /* the formula d[j] = cN[j] - pi' * N[j] is used (note
1.220 + that the vector pi is not changed, because it depends
1.221 + on objective coefficients at basic variables, but in
1.222 + the adjacent basis, for which the vector pi has been
1.223 + just recomputed, x[k] is non-basic) */
1.224 + if (k <= m)
1.225 + { /* x[k] is auxiliary variable */
1.226 + mpq_neg(cbar[ssx->q], pi[k]);
1.227 + }
1.228 + else
1.229 + { /* x[k] is structural variable */
1.230 + int ptr;
1.231 + mpq_t temp;
1.232 + mpq_init(temp);
1.233 + mpq_set_si(cbar[ssx->q], 0, 1);
1.234 + for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
1.235 + { mpq_mul(temp, pi[A_ind[ptr]], A_val[ptr]);
1.236 + mpq_add(cbar[ssx->q], cbar[ssx->q], temp);
1.237 + }
1.238 + mpq_clear(temp);
1.239 + }
1.240 + }
1.241 + }
1.242 + /* jump to the adjacent vertex of the polyhedron */
1.243 + ssx_change_basis(ssx);
1.244 + /* one simplex iteration has been performed */
1.245 + if (ssx->it_lim > 0) ssx->it_lim--;
1.246 + ssx->it_cnt++;
1.247 + }
1.248 + /* display final progress of the search */
1.249 + show_progress(ssx, 1);
1.250 + /* restore components of the original problem, which were changed
1.251 + by the routine */
1.252 + for (k = 1; k <= m+n; k++)
1.253 + { type[k] = orig_type[k];
1.254 + mpq_set(lb[k], orig_lb[k]);
1.255 + mpq_clear(orig_lb[k]);
1.256 + mpq_set(ub[k], orig_ub[k]);
1.257 + mpq_clear(orig_ub[k]);
1.258 + }
1.259 + ssx->dir = orig_dir;
1.260 + for (k = 0; k <= m+n; k++)
1.261 + { mpq_set(coef[k], orig_coef[k]);
1.262 + mpq_clear(orig_coef[k]);
1.263 + }
1.264 + xfree(orig_type);
1.265 + xfree(orig_lb);
1.266 + xfree(orig_ub);
1.267 + xfree(orig_coef);
1.268 + /* return to the calling program */
1.269 + return ret;
1.270 +}
1.271 +
1.272 +/*----------------------------------------------------------------------
1.273 +// ssx_phase_II - find optimal solution.
1.274 +//
1.275 +// This routine implements phase II of the primal simplex method.
1.276 +//
1.277 +// On exit the routine returns one of the following codes:
1.278 +//
1.279 +// 0 - optimal solution found;
1.280 +// 1 - problem has unbounded solution;
1.281 +// 2 - iterations limit exceeded;
1.282 +// 3 - time limit exceeded.
1.283 +----------------------------------------------------------------------*/
1.284 +
1.285 +int ssx_phase_II(SSX *ssx)
1.286 +{ int ret;
1.287 + /* display initial progress of the search */
1.288 + show_progress(ssx, 2);
1.289 + /* main loop starts here */
1.290 + for (;;)
1.291 + { /* display current progress of the search */
1.292 +#if 0
1.293 + if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
1.294 +#else
1.295 + if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
1.296 +#endif
1.297 + show_progress(ssx, 2);
1.298 + /* check if the iterations limit has been exhausted */
1.299 + if (ssx->it_lim == 0)
1.300 + { ret = 2;
1.301 + break;
1.302 + }
1.303 + /* check if the time limit has been exhausted */
1.304 +#if 0
1.305 + if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
1.306 +#else
1.307 + if (ssx->tm_lim >= 0.0 &&
1.308 + ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
1.309 +#endif
1.310 + { ret = 3;
1.311 + break;
1.312 + }
1.313 + /* choose non-basic variable xN[q] */
1.314 + ssx_chuzc(ssx);
1.315 + /* if xN[q] cannot be chosen, the current basic solution is
1.316 + dual feasible and therefore optimal */
1.317 + if (ssx->q == 0)
1.318 + { ret = 0;
1.319 + break;
1.320 + }
1.321 + /* compute q-th column of the simplex table */
1.322 + ssx_eval_col(ssx);
1.323 + /* choose basic variable xB[p] */
1.324 + ssx_chuzr(ssx);
1.325 + /* if xB[p] cannot be chosen, the problem has no dual feasible
1.326 + solution (i.e. unbounded) */
1.327 + if (ssx->p == 0)
1.328 + { ret = 1;
1.329 + break;
1.330 + }
1.331 + /* update values of basic variables */
1.332 + ssx_update_bbar(ssx);
1.333 + if (ssx->p > 0)
1.334 + { /* compute p-th row of the inverse inv(B) */
1.335 + ssx_eval_rho(ssx);
1.336 + /* compute p-th row of the simplex table */
1.337 + ssx_eval_row(ssx);
1.338 + xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
1.339 +#if 0
1.340 + /* update simplex multipliers */
1.341 + ssx_update_pi(ssx);
1.342 +#endif
1.343 + /* update reduced costs of non-basic variables */
1.344 + ssx_update_cbar(ssx);
1.345 + }
1.346 + /* jump to the adjacent vertex of the polyhedron */
1.347 + ssx_change_basis(ssx);
1.348 + /* one simplex iteration has been performed */
1.349 + if (ssx->it_lim > 0) ssx->it_lim--;
1.350 + ssx->it_cnt++;
1.351 + }
1.352 + /* display final progress of the search */
1.353 + show_progress(ssx, 2);
1.354 + /* return to the calling program */
1.355 + return ret;
1.356 +}
1.357 +
1.358 +/*----------------------------------------------------------------------
1.359 +// ssx_driver - base driver to exact simplex method.
1.360 +//
1.361 +// This routine is a base driver to a version of the primal simplex
1.362 +// method using exact (bignum) arithmetic.
1.363 +//
1.364 +// On exit the routine returns one of the following codes:
1.365 +//
1.366 +// 0 - optimal solution found;
1.367 +// 1 - problem has no feasible solution;
1.368 +// 2 - problem has unbounded solution;
1.369 +// 3 - iterations limit exceeded (phase I);
1.370 +// 4 - iterations limit exceeded (phase II);
1.371 +// 5 - time limit exceeded (phase I);
1.372 +// 6 - time limit exceeded (phase II);
1.373 +// 7 - initial basis matrix is exactly singular.
1.374 +----------------------------------------------------------------------*/
1.375 +
1.376 +int ssx_driver(SSX *ssx)
1.377 +{ int m = ssx->m;
1.378 + int *type = ssx->type;
1.379 + mpq_t *lb = ssx->lb;
1.380 + mpq_t *ub = ssx->ub;
1.381 + int *Q_col = ssx->Q_col;
1.382 + mpq_t *bbar = ssx->bbar;
1.383 + int i, k, ret;
1.384 + ssx->tm_beg = xtime();
1.385 + /* factorize the initial basis matrix */
1.386 + if (ssx_factorize(ssx))
1.387 + { xprintf("Initial basis matrix is singular\n");
1.388 + ret = 7;
1.389 + goto done;
1.390 + }
1.391 + /* compute values of basic variables */
1.392 + ssx_eval_bbar(ssx);
1.393 + /* check if the initial basic solution is primal feasible */
1.394 + for (i = 1; i <= m; i++)
1.395 + { int t;
1.396 + k = Q_col[i]; /* x[k] = xB[i] */
1.397 + t = type[k];
1.398 + if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
1.399 + { /* x[k] has lower bound */
1.400 + if (mpq_cmp(bbar[i], lb[k]) < 0)
1.401 + { /* which is violated */
1.402 + break;
1.403 + }
1.404 + }
1.405 + if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
1.406 + { /* x[k] has upper bound */
1.407 + if (mpq_cmp(bbar[i], ub[k]) > 0)
1.408 + { /* which is violated */
1.409 + break;
1.410 + }
1.411 + }
1.412 + }
1.413 + if (i > m)
1.414 + { /* no basic variable violates its bounds */
1.415 + ret = 0;
1.416 + goto skip;
1.417 + }
1.418 + /* phase I: find primal feasible solution */
1.419 + ret = ssx_phase_I(ssx);
1.420 + switch (ret)
1.421 + { case 0:
1.422 + ret = 0;
1.423 + break;
1.424 + case 1:
1.425 + xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
1.426 + ret = 1;
1.427 + break;
1.428 + case 2:
1.429 + xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
1.430 + ret = 3;
1.431 + break;
1.432 + case 3:
1.433 + xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
1.434 + ret = 5;
1.435 + break;
1.436 + default:
1.437 + xassert(ret != ret);
1.438 + }
1.439 + /* compute values of basic variables (actually only the objective
1.440 + value needs to be computed) */
1.441 + ssx_eval_bbar(ssx);
1.442 +skip: /* compute simplex multipliers */
1.443 + ssx_eval_pi(ssx);
1.444 + /* compute reduced costs of non-basic variables */
1.445 + ssx_eval_cbar(ssx);
1.446 + /* if phase I failed, do not start phase II */
1.447 + if (ret != 0) goto done;
1.448 + /* phase II: find optimal solution */
1.449 + ret = ssx_phase_II(ssx);
1.450 + switch (ret)
1.451 + { case 0:
1.452 + xprintf("OPTIMAL SOLUTION FOUND\n");
1.453 + ret = 0;
1.454 + break;
1.455 + case 1:
1.456 + xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
1.457 + ret = 2;
1.458 + break;
1.459 + case 2:
1.460 + xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
1.461 + ret = 4;
1.462 + break;
1.463 + case 3:
1.464 + xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
1.465 + ret = 6;
1.466 + break;
1.467 + default:
1.468 + xassert(ret != ret);
1.469 + }
1.470 +done: /* decrease the time limit by the spent amount of time */
1.471 + if (ssx->tm_lim >= 0.0)
1.472 +#if 0
1.473 + { ssx->tm_lim -= utime() - ssx->tm_beg;
1.474 +#else
1.475 + { ssx->tm_lim -= xdifftime(xtime(), ssx->tm_beg);
1.476 +#endif
1.477 + if (ssx->tm_lim < 0.0) ssx->tm_lim = 0.0;
1.478 + }
1.479 + return ret;
1.480 +}
1.481 +
1.482 +/* eof */