1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpssx01.c	Sun Nov 06 20:59:10 2011 +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, 2011 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 */