1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpssx01.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,839 @@
1.4 +/* glpssx01.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 +#define xfault xerror
1.31 +
1.32 +/*----------------------------------------------------------------------
1.33 +// ssx_create - create simplex solver workspace.
1.34 +//
1.35 +// This routine creates the workspace used by simplex solver routines,
1.36 +// and returns a pointer to it.
1.37 +//
1.38 +// Parameters m, n, and nnz specify, respectively, the number of rows,
1.39 +// columns, and non-zero constraint coefficients.
1.40 +//
1.41 +// This routine only allocates the memory for the workspace components,
1.42 +// so the workspace needs to be saturated by data. */
1.43 +
1.44 +SSX *ssx_create(int m, int n, int nnz)
1.45 +{ SSX *ssx;
1.46 + int i, j, k;
1.47 + if (m < 1)
1.48 + xfault("ssx_create: m = %d; invalid number of rows\n", m);
1.49 + if (n < 1)
1.50 + xfault("ssx_create: n = %d; invalid number of columns\n", n);
1.51 + if (nnz < 0)
1.52 + xfault("ssx_create: nnz = %d; invalid number of non-zero const"
1.53 + "raint coefficients\n", nnz);
1.54 + ssx = xmalloc(sizeof(SSX));
1.55 + ssx->m = m;
1.56 + ssx->n = n;
1.57 + ssx->type = xcalloc(1+m+n, sizeof(int));
1.58 + ssx->lb = xcalloc(1+m+n, sizeof(mpq_t));
1.59 + for (k = 1; k <= m+n; k++) mpq_init(ssx->lb[k]);
1.60 + ssx->ub = xcalloc(1+m+n, sizeof(mpq_t));
1.61 + for (k = 1; k <= m+n; k++) mpq_init(ssx->ub[k]);
1.62 + ssx->coef = xcalloc(1+m+n, sizeof(mpq_t));
1.63 + for (k = 0; k <= m+n; k++) mpq_init(ssx->coef[k]);
1.64 + ssx->A_ptr = xcalloc(1+n+1, sizeof(int));
1.65 + ssx->A_ptr[n+1] = nnz+1;
1.66 + ssx->A_ind = xcalloc(1+nnz, sizeof(int));
1.67 + ssx->A_val = xcalloc(1+nnz, sizeof(mpq_t));
1.68 + for (k = 1; k <= nnz; k++) mpq_init(ssx->A_val[k]);
1.69 + ssx->stat = xcalloc(1+m+n, sizeof(int));
1.70 + ssx->Q_row = xcalloc(1+m+n, sizeof(int));
1.71 + ssx->Q_col = xcalloc(1+m+n, sizeof(int));
1.72 + ssx->binv = bfx_create_binv();
1.73 + ssx->bbar = xcalloc(1+m, sizeof(mpq_t));
1.74 + for (i = 0; i <= m; i++) mpq_init(ssx->bbar[i]);
1.75 + ssx->pi = xcalloc(1+m, sizeof(mpq_t));
1.76 + for (i = 1; i <= m; i++) mpq_init(ssx->pi[i]);
1.77 + ssx->cbar = xcalloc(1+n, sizeof(mpq_t));
1.78 + for (j = 1; j <= n; j++) mpq_init(ssx->cbar[j]);
1.79 + ssx->rho = xcalloc(1+m, sizeof(mpq_t));
1.80 + for (i = 1; i <= m; i++) mpq_init(ssx->rho[i]);
1.81 + ssx->ap = xcalloc(1+n, sizeof(mpq_t));
1.82 + for (j = 1; j <= n; j++) mpq_init(ssx->ap[j]);
1.83 + ssx->aq = xcalloc(1+m, sizeof(mpq_t));
1.84 + for (i = 1; i <= m; i++) mpq_init(ssx->aq[i]);
1.85 + mpq_init(ssx->delta);
1.86 + return ssx;
1.87 +}
1.88 +
1.89 +/*----------------------------------------------------------------------
1.90 +// ssx_factorize - factorize the current basis matrix.
1.91 +//
1.92 +// This routine computes factorization of the current basis matrix B
1.93 +// and returns the singularity flag. If the matrix B is non-singular,
1.94 +// the flag is zero, otherwise non-zero. */
1.95 +
1.96 +static int basis_col(void *info, int j, int ind[], mpq_t val[])
1.97 +{ /* this auxiliary routine provides row indices and numeric values
1.98 + of non-zero elements in j-th column of the matrix B */
1.99 + SSX *ssx = info;
1.100 + int m = ssx->m;
1.101 + int n = ssx->n;
1.102 + int *A_ptr = ssx->A_ptr;
1.103 + int *A_ind = ssx->A_ind;
1.104 + mpq_t *A_val = ssx->A_val;
1.105 + int *Q_col = ssx->Q_col;
1.106 + int k, len, ptr;
1.107 + xassert(1 <= j && j <= m);
1.108 + k = Q_col[j]; /* x[k] = xB[j] */
1.109 + xassert(1 <= k && k <= m+n);
1.110 + /* j-th column of the matrix B is k-th column of the augmented
1.111 + constraint matrix (I | -A) */
1.112 + if (k <= m)
1.113 + { /* it is a column of the unity matrix I */
1.114 + len = 1, ind[1] = k, mpq_set_si(val[1], 1, 1);
1.115 + }
1.116 + else
1.117 + { /* it is a column of the original constraint matrix -A */
1.118 + len = 0;
1.119 + for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
1.120 + { len++;
1.121 + ind[len] = A_ind[ptr];
1.122 + mpq_neg(val[len], A_val[ptr]);
1.123 + }
1.124 + }
1.125 + return len;
1.126 +}
1.127 +
1.128 +int ssx_factorize(SSX *ssx)
1.129 +{ int ret;
1.130 + ret = bfx_factorize(ssx->binv, ssx->m, basis_col, ssx);
1.131 + return ret;
1.132 +}
1.133 +
1.134 +/*----------------------------------------------------------------------
1.135 +// ssx_get_xNj - determine value of non-basic variable.
1.136 +//
1.137 +// This routine determines the value of non-basic variable xN[j] in the
1.138 +// current basic solution defined as follows:
1.139 +//
1.140 +// 0, if xN[j] is free variable
1.141 +// lN[j], if xN[j] is on its lower bound
1.142 +// uN[j], if xN[j] is on its upper bound
1.143 +// lN[j] = uN[j], if xN[j] is fixed variable
1.144 +//
1.145 +// where lN[j] and uN[j] are lower and upper bounds of xN[j]. */
1.146 +
1.147 +void ssx_get_xNj(SSX *ssx, int j, mpq_t x)
1.148 +{ int m = ssx->m;
1.149 + int n = ssx->n;
1.150 + mpq_t *lb = ssx->lb;
1.151 + mpq_t *ub = ssx->ub;
1.152 + int *stat = ssx->stat;
1.153 + int *Q_col = ssx->Q_col;
1.154 + int k;
1.155 + xassert(1 <= j && j <= n);
1.156 + k = Q_col[m+j]; /* x[k] = xN[j] */
1.157 + xassert(1 <= k && k <= m+n);
1.158 + switch (stat[k])
1.159 + { case SSX_NL:
1.160 + /* xN[j] is on its lower bound */
1.161 + mpq_set(x, lb[k]); break;
1.162 + case SSX_NU:
1.163 + /* xN[j] is on its upper bound */
1.164 + mpq_set(x, ub[k]); break;
1.165 + case SSX_NF:
1.166 + /* xN[j] is free variable */
1.167 + mpq_set_si(x, 0, 1); break;
1.168 + case SSX_NS:
1.169 + /* xN[j] is fixed variable */
1.170 + mpq_set(x, lb[k]); break;
1.171 + default:
1.172 + xassert(stat != stat);
1.173 + }
1.174 + return;
1.175 +}
1.176 +
1.177 +/*----------------------------------------------------------------------
1.178 +// ssx_eval_bbar - compute values of basic variables.
1.179 +//
1.180 +// This routine computes values of basic variables xB in the current
1.181 +// basic solution as follows:
1.182 +//
1.183 +// beta = - inv(B) * N * xN,
1.184 +//
1.185 +// where B is the basis matrix, N is the matrix of non-basic columns,
1.186 +// xN is a vector of current values of non-basic variables. */
1.187 +
1.188 +void ssx_eval_bbar(SSX *ssx)
1.189 +{ int m = ssx->m;
1.190 + int n = ssx->n;
1.191 + mpq_t *coef = ssx->coef;
1.192 + int *A_ptr = ssx->A_ptr;
1.193 + int *A_ind = ssx->A_ind;
1.194 + mpq_t *A_val = ssx->A_val;
1.195 + int *Q_col = ssx->Q_col;
1.196 + mpq_t *bbar = ssx->bbar;
1.197 + int i, j, k, ptr;
1.198 + mpq_t x, temp;
1.199 + mpq_init(x);
1.200 + mpq_init(temp);
1.201 + /* bbar := 0 */
1.202 + for (i = 1; i <= m; i++)
1.203 + mpq_set_si(bbar[i], 0, 1);
1.204 + /* bbar := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n] */
1.205 + for (j = 1; j <= n; j++)
1.206 + { ssx_get_xNj(ssx, j, x);
1.207 + if (mpq_sgn(x) == 0) continue;
1.208 + k = Q_col[m+j]; /* x[k] = xN[j] */
1.209 + if (k <= m)
1.210 + { /* N[j] is a column of the unity matrix I */
1.211 + mpq_sub(bbar[k], bbar[k], x);
1.212 + }
1.213 + else
1.214 + { /* N[j] is a column of the original constraint matrix -A */
1.215 + for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
1.216 + { mpq_mul(temp, A_val[ptr], x);
1.217 + mpq_add(bbar[A_ind[ptr]], bbar[A_ind[ptr]], temp);
1.218 + }
1.219 + }
1.220 + }
1.221 + /* bbar := inv(B) * bbar */
1.222 + bfx_ftran(ssx->binv, bbar, 0);
1.223 +#if 1
1.224 + /* compute value of the objective function */
1.225 + /* bbar[0] := c[0] */
1.226 + mpq_set(bbar[0], coef[0]);
1.227 + /* bbar[0] := bbar[0] + sum{i in B} cB[i] * xB[i] */
1.228 + for (i = 1; i <= m; i++)
1.229 + { k = Q_col[i]; /* x[k] = xB[i] */
1.230 + if (mpq_sgn(coef[k]) == 0) continue;
1.231 + mpq_mul(temp, coef[k], bbar[i]);
1.232 + mpq_add(bbar[0], bbar[0], temp);
1.233 + }
1.234 + /* bbar[0] := bbar[0] + sum{j in N} cN[j] * xN[j] */
1.235 + for (j = 1; j <= n; j++)
1.236 + { k = Q_col[m+j]; /* x[k] = xN[j] */
1.237 + if (mpq_sgn(coef[k]) == 0) continue;
1.238 + ssx_get_xNj(ssx, j, x);
1.239 + mpq_mul(temp, coef[k], x);
1.240 + mpq_add(bbar[0], bbar[0], temp);
1.241 + }
1.242 +#endif
1.243 + mpq_clear(x);
1.244 + mpq_clear(temp);
1.245 + return;
1.246 +}
1.247 +
1.248 +/*----------------------------------------------------------------------
1.249 +// ssx_eval_pi - compute values of simplex multipliers.
1.250 +//
1.251 +// This routine computes values of simplex multipliers (shadow prices)
1.252 +// pi in the current basic solution as follows:
1.253 +//
1.254 +// pi = inv(B') * cB,
1.255 +//
1.256 +// where B' is a matrix transposed to the basis matrix B, cB is a vector
1.257 +// of objective coefficients at basic variables xB. */
1.258 +
1.259 +void ssx_eval_pi(SSX *ssx)
1.260 +{ int m = ssx->m;
1.261 + mpq_t *coef = ssx->coef;
1.262 + int *Q_col = ssx->Q_col;
1.263 + mpq_t *pi = ssx->pi;
1.264 + int i;
1.265 + /* pi := cB */
1.266 + for (i = 1; i <= m; i++) mpq_set(pi[i], coef[Q_col[i]]);
1.267 + /* pi := inv(B') * cB */
1.268 + bfx_btran(ssx->binv, pi);
1.269 + return;
1.270 +}
1.271 +
1.272 +/*----------------------------------------------------------------------
1.273 +// ssx_eval_dj - compute reduced cost of non-basic variable.
1.274 +//
1.275 +// This routine computes reduced cost d[j] of non-basic variable xN[j]
1.276 +// in the current basic solution as follows:
1.277 +//
1.278 +// d[j] = cN[j] - N[j] * pi,
1.279 +//
1.280 +// where cN[j] is an objective coefficient at xN[j], N[j] is a column
1.281 +// of the augmented constraint matrix (I | -A) corresponding to xN[j],
1.282 +// pi is the vector of simplex multipliers (shadow prices). */
1.283 +
1.284 +void ssx_eval_dj(SSX *ssx, int j, mpq_t dj)
1.285 +{ int m = ssx->m;
1.286 + int n = ssx->n;
1.287 + mpq_t *coef = ssx->coef;
1.288 + int *A_ptr = ssx->A_ptr;
1.289 + int *A_ind = ssx->A_ind;
1.290 + mpq_t *A_val = ssx->A_val;
1.291 + int *Q_col = ssx->Q_col;
1.292 + mpq_t *pi = ssx->pi;
1.293 + int k, ptr, end;
1.294 + mpq_t temp;
1.295 + mpq_init(temp);
1.296 + xassert(1 <= j && j <= n);
1.297 + k = Q_col[m+j]; /* x[k] = xN[j] */
1.298 + xassert(1 <= k && k <= m+n);
1.299 + /* j-th column of the matrix N is k-th column of the augmented
1.300 + constraint matrix (I | -A) */
1.301 + if (k <= m)
1.302 + { /* it is a column of the unity matrix I */
1.303 + mpq_sub(dj, coef[k], pi[k]);
1.304 + }
1.305 + else
1.306 + { /* it is a column of the original constraint matrix -A */
1.307 + mpq_set(dj, coef[k]);
1.308 + for (ptr = A_ptr[k-m], end = A_ptr[k-m+1]; ptr < end; ptr++)
1.309 + { mpq_mul(temp, A_val[ptr], pi[A_ind[ptr]]);
1.310 + mpq_add(dj, dj, temp);
1.311 + }
1.312 + }
1.313 + mpq_clear(temp);
1.314 + return;
1.315 +}
1.316 +
1.317 +/*----------------------------------------------------------------------
1.318 +// ssx_eval_cbar - compute reduced costs of all non-basic variables.
1.319 +//
1.320 +// This routine computes the vector of reduced costs pi in the current
1.321 +// basic solution for all non-basic variables, including fixed ones. */
1.322 +
1.323 +void ssx_eval_cbar(SSX *ssx)
1.324 +{ int n = ssx->n;
1.325 + mpq_t *cbar = ssx->cbar;
1.326 + int j;
1.327 + for (j = 1; j <= n; j++)
1.328 + ssx_eval_dj(ssx, j, cbar[j]);
1.329 + return;
1.330 +}
1.331 +
1.332 +/*----------------------------------------------------------------------
1.333 +// ssx_eval_rho - compute p-th row of the inverse.
1.334 +//
1.335 +// This routine computes p-th row of the matrix inv(B), where B is the
1.336 +// current basis matrix.
1.337 +//
1.338 +// p-th row of the inverse is computed using the following formula:
1.339 +//
1.340 +// rho = inv(B') * e[p],
1.341 +//
1.342 +// where B' is a matrix transposed to B, e[p] is a unity vector, which
1.343 +// contains one in p-th position. */
1.344 +
1.345 +void ssx_eval_rho(SSX *ssx)
1.346 +{ int m = ssx->m;
1.347 + int p = ssx->p;
1.348 + mpq_t *rho = ssx->rho;
1.349 + int i;
1.350 + xassert(1 <= p && p <= m);
1.351 + /* rho := 0 */
1.352 + for (i = 1; i <= m; i++) mpq_set_si(rho[i], 0, 1);
1.353 + /* rho := e[p] */
1.354 + mpq_set_si(rho[p], 1, 1);
1.355 + /* rho := inv(B') * rho */
1.356 + bfx_btran(ssx->binv, rho);
1.357 + return;
1.358 +}
1.359 +
1.360 +/*----------------------------------------------------------------------
1.361 +// ssx_eval_row - compute pivot row of the simplex table.
1.362 +//
1.363 +// This routine computes p-th (pivot) row of the current simplex table
1.364 +// A~ = - inv(B) * N using the following formula:
1.365 +//
1.366 +// A~[p] = - N' * inv(B') * e[p] = - N' * rho[p],
1.367 +//
1.368 +// where N' is a matrix transposed to the matrix N, rho[p] is p-th row
1.369 +// of the inverse inv(B). */
1.370 +
1.371 +void ssx_eval_row(SSX *ssx)
1.372 +{ int m = ssx->m;
1.373 + int n = ssx->n;
1.374 + int *A_ptr = ssx->A_ptr;
1.375 + int *A_ind = ssx->A_ind;
1.376 + mpq_t *A_val = ssx->A_val;
1.377 + int *Q_col = ssx->Q_col;
1.378 + mpq_t *rho = ssx->rho;
1.379 + mpq_t *ap = ssx->ap;
1.380 + int j, k, ptr;
1.381 + mpq_t temp;
1.382 + mpq_init(temp);
1.383 + for (j = 1; j <= n; j++)
1.384 + { /* ap[j] := - N'[j] * rho (inner product) */
1.385 + k = Q_col[m+j]; /* x[k] = xN[j] */
1.386 + if (k <= m)
1.387 + mpq_neg(ap[j], rho[k]);
1.388 + else
1.389 + { mpq_set_si(ap[j], 0, 1);
1.390 + for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
1.391 + { mpq_mul(temp, A_val[ptr], rho[A_ind[ptr]]);
1.392 + mpq_add(ap[j], ap[j], temp);
1.393 + }
1.394 + }
1.395 + }
1.396 + mpq_clear(temp);
1.397 + return;
1.398 +}
1.399 +
1.400 +/*----------------------------------------------------------------------
1.401 +// ssx_eval_col - compute pivot column of the simplex table.
1.402 +//
1.403 +// This routine computes q-th (pivot) column of the current simplex
1.404 +// table A~ = - inv(B) * N using the following formula:
1.405 +//
1.406 +// A~[q] = - inv(B) * N[q],
1.407 +//
1.408 +// where N[q] is q-th column of the matrix N corresponding to chosen
1.409 +// non-basic variable xN[q]. */
1.410 +
1.411 +void ssx_eval_col(SSX *ssx)
1.412 +{ int m = ssx->m;
1.413 + int n = ssx->n;
1.414 + int *A_ptr = ssx->A_ptr;
1.415 + int *A_ind = ssx->A_ind;
1.416 + mpq_t *A_val = ssx->A_val;
1.417 + int *Q_col = ssx->Q_col;
1.418 + int q = ssx->q;
1.419 + mpq_t *aq = ssx->aq;
1.420 + int i, k, ptr;
1.421 + xassert(1 <= q && q <= n);
1.422 + /* aq := 0 */
1.423 + for (i = 1; i <= m; i++) mpq_set_si(aq[i], 0, 1);
1.424 + /* aq := N[q] */
1.425 + k = Q_col[m+q]; /* x[k] = xN[q] */
1.426 + if (k <= m)
1.427 + { /* N[q] is a column of the unity matrix I */
1.428 + mpq_set_si(aq[k], 1, 1);
1.429 + }
1.430 + else
1.431 + { /* N[q] is a column of the original constraint matrix -A */
1.432 + for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
1.433 + mpq_neg(aq[A_ind[ptr]], A_val[ptr]);
1.434 + }
1.435 + /* aq := inv(B) * aq */
1.436 + bfx_ftran(ssx->binv, aq, 1);
1.437 + /* aq := - aq */
1.438 + for (i = 1; i <= m; i++) mpq_neg(aq[i], aq[i]);
1.439 + return;
1.440 +}
1.441 +
1.442 +/*----------------------------------------------------------------------
1.443 +// ssx_chuzc - choose pivot column.
1.444 +//
1.445 +// This routine chooses non-basic variable xN[q] whose reduced cost
1.446 +// indicates possible improving of the objective function to enter it
1.447 +// in the basis.
1.448 +//
1.449 +// Currently the standard (textbook) pricing is used, i.e. that
1.450 +// non-basic variable is preferred which has greatest reduced cost (in
1.451 +// magnitude).
1.452 +//
1.453 +// If xN[q] has been chosen, the routine stores its number q and also
1.454 +// sets the flag q_dir that indicates direction in which xN[q] has to
1.455 +// change (+1 means increasing, -1 means decreasing).
1.456 +//
1.457 +// If the choice cannot be made, because the current basic solution is
1.458 +// dual feasible, the routine sets the number q to 0. */
1.459 +
1.460 +void ssx_chuzc(SSX *ssx)
1.461 +{ int m = ssx->m;
1.462 + int n = ssx->n;
1.463 + int dir = (ssx->dir == SSX_MIN ? +1 : -1);
1.464 + int *Q_col = ssx->Q_col;
1.465 + int *stat = ssx->stat;
1.466 + mpq_t *cbar = ssx->cbar;
1.467 + int j, k, s, q, q_dir;
1.468 + double best, temp;
1.469 + /* nothing is chosen so far */
1.470 + q = 0, q_dir = 0, best = 0.0;
1.471 + /* look through the list of non-basic variables */
1.472 + for (j = 1; j <= n; j++)
1.473 + { k = Q_col[m+j]; /* x[k] = xN[j] */
1.474 + s = dir * mpq_sgn(cbar[j]);
1.475 + if ((stat[k] == SSX_NF || stat[k] == SSX_NL) && s < 0 ||
1.476 + (stat[k] == SSX_NF || stat[k] == SSX_NU) && s > 0)
1.477 + { /* reduced cost of xN[j] indicates possible improving of
1.478 + the objective function */
1.479 + temp = fabs(mpq_get_d(cbar[j]));
1.480 + xassert(temp != 0.0);
1.481 + if (q == 0 || best < temp)
1.482 + q = j, q_dir = - s, best = temp;
1.483 + }
1.484 + }
1.485 + ssx->q = q, ssx->q_dir = q_dir;
1.486 + return;
1.487 +}
1.488 +
1.489 +/*----------------------------------------------------------------------
1.490 +// ssx_chuzr - choose pivot row.
1.491 +//
1.492 +// This routine looks through elements of q-th column of the simplex
1.493 +// table and chooses basic variable xB[p] which should leave the basis.
1.494 +//
1.495 +// The choice is based on the standard (textbook) ratio test.
1.496 +//
1.497 +// If xB[p] has been chosen, the routine stores its number p and also
1.498 +// sets its non-basic status p_stat which should be assigned to xB[p]
1.499 +// when it has left the basis and become xN[q].
1.500 +//
1.501 +// Special case p < 0 means that xN[q] is double-bounded variable and
1.502 +// it reaches its opposite bound before any basic variable does that,
1.503 +// so the current basis remains unchanged.
1.504 +//
1.505 +// If the choice cannot be made, because xN[q] can infinitely change in
1.506 +// the feasible direction, the routine sets the number p to 0. */
1.507 +
1.508 +void ssx_chuzr(SSX *ssx)
1.509 +{ int m = ssx->m;
1.510 + int n = ssx->n;
1.511 + int *type = ssx->type;
1.512 + mpq_t *lb = ssx->lb;
1.513 + mpq_t *ub = ssx->ub;
1.514 + int *Q_col = ssx->Q_col;
1.515 + mpq_t *bbar = ssx->bbar;
1.516 + int q = ssx->q;
1.517 + mpq_t *aq = ssx->aq;
1.518 + int q_dir = ssx->q_dir;
1.519 + int i, k, s, t, p, p_stat;
1.520 + mpq_t teta, temp;
1.521 + mpq_init(teta);
1.522 + mpq_init(temp);
1.523 + xassert(1 <= q && q <= n);
1.524 + xassert(q_dir == +1 || q_dir == -1);
1.525 + /* nothing is chosen so far */
1.526 + p = 0, p_stat = 0;
1.527 + /* look through the list of basic variables */
1.528 + for (i = 1; i <= m; i++)
1.529 + { s = q_dir * mpq_sgn(aq[i]);
1.530 + if (s < 0)
1.531 + { /* xB[i] decreases */
1.532 + k = Q_col[i]; /* x[k] = xB[i] */
1.533 + t = type[k];
1.534 + if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
1.535 + { /* xB[i] has finite lower bound */
1.536 + mpq_sub(temp, bbar[i], lb[k]);
1.537 + mpq_div(temp, temp, aq[i]);
1.538 + mpq_abs(temp, temp);
1.539 + if (p == 0 || mpq_cmp(teta, temp) > 0)
1.540 + { p = i;
1.541 + p_stat = (t == SSX_FX ? SSX_NS : SSX_NL);
1.542 + mpq_set(teta, temp);
1.543 + }
1.544 + }
1.545 + }
1.546 + else if (s > 0)
1.547 + { /* xB[i] increases */
1.548 + k = Q_col[i]; /* x[k] = xB[i] */
1.549 + t = type[k];
1.550 + if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
1.551 + { /* xB[i] has finite upper bound */
1.552 + mpq_sub(temp, bbar[i], ub[k]);
1.553 + mpq_div(temp, temp, aq[i]);
1.554 + mpq_abs(temp, temp);
1.555 + if (p == 0 || mpq_cmp(teta, temp) > 0)
1.556 + { p = i;
1.557 + p_stat = (t == SSX_FX ? SSX_NS : SSX_NU);
1.558 + mpq_set(teta, temp);
1.559 + }
1.560 + }
1.561 + }
1.562 + /* if something has been chosen and the ratio test indicates
1.563 + exact degeneracy, the search can be finished */
1.564 + if (p != 0 && mpq_sgn(teta) == 0) break;
1.565 + }
1.566 + /* if xN[q] is double-bounded, check if it can reach its opposite
1.567 + bound before any basic variable */
1.568 + k = Q_col[m+q]; /* x[k] = xN[q] */
1.569 + if (type[k] == SSX_DB)
1.570 + { mpq_sub(temp, ub[k], lb[k]);
1.571 + if (p == 0 || mpq_cmp(teta, temp) > 0)
1.572 + { p = -1;
1.573 + p_stat = -1;
1.574 + mpq_set(teta, temp);
1.575 + }
1.576 + }
1.577 + ssx->p = p;
1.578 + ssx->p_stat = p_stat;
1.579 + /* if xB[p] has been chosen, determine its actual change in the
1.580 + adjacent basis (it has the same sign as q_dir) */
1.581 + if (p != 0)
1.582 + { xassert(mpq_sgn(teta) >= 0);
1.583 + if (q_dir > 0)
1.584 + mpq_set(ssx->delta, teta);
1.585 + else
1.586 + mpq_neg(ssx->delta, teta);
1.587 + }
1.588 + mpq_clear(teta);
1.589 + mpq_clear(temp);
1.590 + return;
1.591 +}
1.592 +
1.593 +/*----------------------------------------------------------------------
1.594 +// ssx_update_bbar - update values of basic variables.
1.595 +//
1.596 +// This routine recomputes the current values of basic variables for
1.597 +// the adjacent basis.
1.598 +//
1.599 +// The simplex table for the current basis is the following:
1.600 +//
1.601 +// xB[i] = sum{j in 1..n} alfa[i,j] * xN[q], i = 1,...,m
1.602 +//
1.603 +// therefore
1.604 +//
1.605 +// delta xB[i] = alfa[i,q] * delta xN[q], i = 1,...,m
1.606 +//
1.607 +// where delta xN[q] = xN.new[q] - xN[q] is the change of xN[q] in the
1.608 +// adjacent basis, and delta xB[i] = xB.new[i] - xB[i] is the change of
1.609 +// xB[i]. This gives formulae for recomputing values of xB[i]:
1.610 +//
1.611 +// xB.new[p] = xN[q] + delta xN[q]
1.612 +//
1.613 +// (because xN[q] becomes xB[p] in the adjacent basis), and
1.614 +//
1.615 +// xB.new[i] = xB[i] + alfa[i,q] * delta xN[q], i != p
1.616 +//
1.617 +// for other basic variables. */
1.618 +
1.619 +void ssx_update_bbar(SSX *ssx)
1.620 +{ int m = ssx->m;
1.621 + int n = ssx->n;
1.622 + mpq_t *bbar = ssx->bbar;
1.623 + mpq_t *cbar = ssx->cbar;
1.624 + int p = ssx->p;
1.625 + int q = ssx->q;
1.626 + mpq_t *aq = ssx->aq;
1.627 + int i;
1.628 + mpq_t temp;
1.629 + mpq_init(temp);
1.630 + xassert(1 <= q && q <= n);
1.631 + if (p < 0)
1.632 + { /* xN[q] is double-bounded and goes to its opposite bound */
1.633 + /* nop */;
1.634 + }
1.635 + else
1.636 + { /* xN[q] becomes xB[p] in the adjacent basis */
1.637 + /* xB.new[p] = xN[q] + delta xN[q] */
1.638 + xassert(1 <= p && p <= m);
1.639 + ssx_get_xNj(ssx, q, temp);
1.640 + mpq_add(bbar[p], temp, ssx->delta);
1.641 + }
1.642 + /* update values of other basic variables depending on xN[q] */
1.643 + for (i = 1; i <= m; i++)
1.644 + { if (i == p) continue;
1.645 + /* xB.new[i] = xB[i] + alfa[i,q] * delta xN[q] */
1.646 + if (mpq_sgn(aq[i]) == 0) continue;
1.647 + mpq_mul(temp, aq[i], ssx->delta);
1.648 + mpq_add(bbar[i], bbar[i], temp);
1.649 + }
1.650 +#if 1
1.651 + /* update value of the objective function */
1.652 + /* z.new = z + d[q] * delta xN[q] */
1.653 + mpq_mul(temp, cbar[q], ssx->delta);
1.654 + mpq_add(bbar[0], bbar[0], temp);
1.655 +#endif
1.656 + mpq_clear(temp);
1.657 + return;
1.658 +}
1.659 +
1.660 +/*----------------------------------------------------------------------
1.661 +-- ssx_update_pi - update simplex multipliers.
1.662 +--
1.663 +-- This routine recomputes the vector of simplex multipliers for the
1.664 +-- adjacent basis. */
1.665 +
1.666 +void ssx_update_pi(SSX *ssx)
1.667 +{ int m = ssx->m;
1.668 + int n = ssx->n;
1.669 + mpq_t *pi = ssx->pi;
1.670 + mpq_t *cbar = ssx->cbar;
1.671 + int p = ssx->p;
1.672 + int q = ssx->q;
1.673 + mpq_t *aq = ssx->aq;
1.674 + mpq_t *rho = ssx->rho;
1.675 + int i;
1.676 + mpq_t new_dq, temp;
1.677 + mpq_init(new_dq);
1.678 + mpq_init(temp);
1.679 + xassert(1 <= p && p <= m);
1.680 + xassert(1 <= q && q <= n);
1.681 + /* compute d[q] in the adjacent basis */
1.682 + mpq_div(new_dq, cbar[q], aq[p]);
1.683 + /* update the vector of simplex multipliers */
1.684 + for (i = 1; i <= m; i++)
1.685 + { if (mpq_sgn(rho[i]) == 0) continue;
1.686 + mpq_mul(temp, new_dq, rho[i]);
1.687 + mpq_sub(pi[i], pi[i], temp);
1.688 + }
1.689 + mpq_clear(new_dq);
1.690 + mpq_clear(temp);
1.691 + return;
1.692 +}
1.693 +
1.694 +/*----------------------------------------------------------------------
1.695 +// ssx_update_cbar - update reduced costs of non-basic variables.
1.696 +//
1.697 +// This routine recomputes the vector of reduced costs of non-basic
1.698 +// variables for the adjacent basis. */
1.699 +
1.700 +void ssx_update_cbar(SSX *ssx)
1.701 +{ int m = ssx->m;
1.702 + int n = ssx->n;
1.703 + mpq_t *cbar = ssx->cbar;
1.704 + int p = ssx->p;
1.705 + int q = ssx->q;
1.706 + mpq_t *ap = ssx->ap;
1.707 + int j;
1.708 + mpq_t temp;
1.709 + mpq_init(temp);
1.710 + xassert(1 <= p && p <= m);
1.711 + xassert(1 <= q && q <= n);
1.712 + /* compute d[q] in the adjacent basis */
1.713 + /* d.new[q] = d[q] / alfa[p,q] */
1.714 + mpq_div(cbar[q], cbar[q], ap[q]);
1.715 + /* update reduced costs of other non-basic variables */
1.716 + for (j = 1; j <= n; j++)
1.717 + { if (j == q) continue;
1.718 + /* d.new[j] = d[j] - (alfa[p,j] / alfa[p,q]) * d[q] */
1.719 + if (mpq_sgn(ap[j]) == 0) continue;
1.720 + mpq_mul(temp, ap[j], cbar[q]);
1.721 + mpq_sub(cbar[j], cbar[j], temp);
1.722 + }
1.723 + mpq_clear(temp);
1.724 + return;
1.725 +}
1.726 +
1.727 +/*----------------------------------------------------------------------
1.728 +// ssx_change_basis - change current basis to adjacent one.
1.729 +//
1.730 +// This routine changes the current basis to the adjacent one swapping
1.731 +// basic variable xB[p] and non-basic variable xN[q]. */
1.732 +
1.733 +void ssx_change_basis(SSX *ssx)
1.734 +{ int m = ssx->m;
1.735 + int n = ssx->n;
1.736 + int *type = ssx->type;
1.737 + int *stat = ssx->stat;
1.738 + int *Q_row = ssx->Q_row;
1.739 + int *Q_col = ssx->Q_col;
1.740 + int p = ssx->p;
1.741 + int q = ssx->q;
1.742 + int p_stat = ssx->p_stat;
1.743 + int k, kp, kq;
1.744 + if (p < 0)
1.745 + { /* special case: xN[q] goes to its opposite bound */
1.746 + xassert(1 <= q && q <= n);
1.747 + k = Q_col[m+q]; /* x[k] = xN[q] */
1.748 + xassert(type[k] == SSX_DB);
1.749 + switch (stat[k])
1.750 + { case SSX_NL:
1.751 + stat[k] = SSX_NU;
1.752 + break;
1.753 + case SSX_NU:
1.754 + stat[k] = SSX_NL;
1.755 + break;
1.756 + default:
1.757 + xassert(stat != stat);
1.758 + }
1.759 + }
1.760 + else
1.761 + { /* xB[p] leaves the basis, xN[q] enters the basis */
1.762 + xassert(1 <= p && p <= m);
1.763 + xassert(1 <= q && q <= n);
1.764 + kp = Q_col[p]; /* x[kp] = xB[p] */
1.765 + kq = Q_col[m+q]; /* x[kq] = xN[q] */
1.766 + /* check non-basic status of xB[p] which becomes xN[q] */
1.767 + switch (type[kp])
1.768 + { case SSX_FR:
1.769 + xassert(p_stat == SSX_NF);
1.770 + break;
1.771 + case SSX_LO:
1.772 + xassert(p_stat == SSX_NL);
1.773 + break;
1.774 + case SSX_UP:
1.775 + xassert(p_stat == SSX_NU);
1.776 + break;
1.777 + case SSX_DB:
1.778 + xassert(p_stat == SSX_NL || p_stat == SSX_NU);
1.779 + break;
1.780 + case SSX_FX:
1.781 + xassert(p_stat == SSX_NS);
1.782 + break;
1.783 + default:
1.784 + xassert(type != type);
1.785 + }
1.786 + /* swap xB[p] and xN[q] */
1.787 + stat[kp] = (char)p_stat, stat[kq] = SSX_BS;
1.788 + Q_row[kp] = m+q, Q_row[kq] = p;
1.789 + Q_col[p] = kq, Q_col[m+q] = kp;
1.790 + /* update factorization of the basis matrix */
1.791 + if (bfx_update(ssx->binv, p))
1.792 + { if (ssx_factorize(ssx))
1.793 + xassert(("Internal error: basis matrix is singular", 0));
1.794 + }
1.795 + }
1.796 + return;
1.797 +}
1.798 +
1.799 +/*----------------------------------------------------------------------
1.800 +// ssx_delete - delete simplex solver workspace.
1.801 +//
1.802 +// This routine deletes the simplex solver workspace freeing all the
1.803 +// memory allocated to this object. */
1.804 +
1.805 +void ssx_delete(SSX *ssx)
1.806 +{ int m = ssx->m;
1.807 + int n = ssx->n;
1.808 + int nnz = ssx->A_ptr[n+1]-1;
1.809 + int i, j, k;
1.810 + xfree(ssx->type);
1.811 + for (k = 1; k <= m+n; k++) mpq_clear(ssx->lb[k]);
1.812 + xfree(ssx->lb);
1.813 + for (k = 1; k <= m+n; k++) mpq_clear(ssx->ub[k]);
1.814 + xfree(ssx->ub);
1.815 + for (k = 0; k <= m+n; k++) mpq_clear(ssx->coef[k]);
1.816 + xfree(ssx->coef);
1.817 + xfree(ssx->A_ptr);
1.818 + xfree(ssx->A_ind);
1.819 + for (k = 1; k <= nnz; k++) mpq_clear(ssx->A_val[k]);
1.820 + xfree(ssx->A_val);
1.821 + xfree(ssx->stat);
1.822 + xfree(ssx->Q_row);
1.823 + xfree(ssx->Q_col);
1.824 + bfx_delete_binv(ssx->binv);
1.825 + for (i = 0; i <= m; i++) mpq_clear(ssx->bbar[i]);
1.826 + xfree(ssx->bbar);
1.827 + for (i = 1; i <= m; i++) mpq_clear(ssx->pi[i]);
1.828 + xfree(ssx->pi);
1.829 + for (j = 1; j <= n; j++) mpq_clear(ssx->cbar[j]);
1.830 + xfree(ssx->cbar);
1.831 + for (i = 1; i <= m; i++) mpq_clear(ssx->rho[i]);
1.832 + xfree(ssx->rho);
1.833 + for (j = 1; j <= n; j++) mpq_clear(ssx->ap[j]);
1.834 + xfree(ssx->ap);
1.835 + for (i = 1; i <= m; i++) mpq_clear(ssx->aq[i]);
1.836 + xfree(ssx->aq);
1.837 + mpq_clear(ssx->delta);
1.838 + xfree(ssx);
1.839 + return;
1.840 +}
1.841 +
1.842 +/* eof */