[1] | 1 | /* glpssx01.c */ |
---|
| 2 | |
---|
| 3 | /*********************************************************************** |
---|
| 4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
| 5 | * |
---|
| 6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
| 7 | * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
---|
| 8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
| 9 | * E-mail: <mao@gnu.org>. |
---|
| 10 | * |
---|
| 11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
| 12 | * under the terms of the GNU General Public License as published by |
---|
| 13 | * the Free Software Foundation, either version 3 of the License, or |
---|
| 14 | * (at your option) any later version. |
---|
| 15 | * |
---|
| 16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
| 17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
| 18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
| 19 | * License for more details. |
---|
| 20 | * |
---|
| 21 | * You should have received a copy of the GNU General Public License |
---|
| 22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
| 23 | ***********************************************************************/ |
---|
| 24 | |
---|
| 25 | #include "glpenv.h" |
---|
| 26 | #include "glpssx.h" |
---|
| 27 | #define xfault xerror |
---|
| 28 | |
---|
| 29 | /*---------------------------------------------------------------------- |
---|
| 30 | // ssx_create - create simplex solver workspace. |
---|
| 31 | // |
---|
| 32 | // This routine creates the workspace used by simplex solver routines, |
---|
| 33 | // and returns a pointer to it. |
---|
| 34 | // |
---|
| 35 | // Parameters m, n, and nnz specify, respectively, the number of rows, |
---|
| 36 | // columns, and non-zero constraint coefficients. |
---|
| 37 | // |
---|
| 38 | // This routine only allocates the memory for the workspace components, |
---|
| 39 | // so the workspace needs to be saturated by data. */ |
---|
| 40 | |
---|
| 41 | SSX *ssx_create(int m, int n, int nnz) |
---|
| 42 | { SSX *ssx; |
---|
| 43 | int i, j, k; |
---|
| 44 | if (m < 1) |
---|
| 45 | xfault("ssx_create: m = %d; invalid number of rows\n", m); |
---|
| 46 | if (n < 1) |
---|
| 47 | xfault("ssx_create: n = %d; invalid number of columns\n", n); |
---|
| 48 | if (nnz < 0) |
---|
| 49 | xfault("ssx_create: nnz = %d; invalid number of non-zero const" |
---|
| 50 | "raint coefficients\n", nnz); |
---|
| 51 | ssx = xmalloc(sizeof(SSX)); |
---|
| 52 | ssx->m = m; |
---|
| 53 | ssx->n = n; |
---|
| 54 | ssx->type = xcalloc(1+m+n, sizeof(int)); |
---|
| 55 | ssx->lb = xcalloc(1+m+n, sizeof(mpq_t)); |
---|
| 56 | for (k = 1; k <= m+n; k++) mpq_init(ssx->lb[k]); |
---|
| 57 | ssx->ub = xcalloc(1+m+n, sizeof(mpq_t)); |
---|
| 58 | for (k = 1; k <= m+n; k++) mpq_init(ssx->ub[k]); |
---|
| 59 | ssx->coef = xcalloc(1+m+n, sizeof(mpq_t)); |
---|
| 60 | for (k = 0; k <= m+n; k++) mpq_init(ssx->coef[k]); |
---|
| 61 | ssx->A_ptr = xcalloc(1+n+1, sizeof(int)); |
---|
| 62 | ssx->A_ptr[n+1] = nnz+1; |
---|
| 63 | ssx->A_ind = xcalloc(1+nnz, sizeof(int)); |
---|
| 64 | ssx->A_val = xcalloc(1+nnz, sizeof(mpq_t)); |
---|
| 65 | for (k = 1; k <= nnz; k++) mpq_init(ssx->A_val[k]); |
---|
| 66 | ssx->stat = xcalloc(1+m+n, sizeof(int)); |
---|
| 67 | ssx->Q_row = xcalloc(1+m+n, sizeof(int)); |
---|
| 68 | ssx->Q_col = xcalloc(1+m+n, sizeof(int)); |
---|
| 69 | ssx->binv = bfx_create_binv(); |
---|
| 70 | ssx->bbar = xcalloc(1+m, sizeof(mpq_t)); |
---|
| 71 | for (i = 0; i <= m; i++) mpq_init(ssx->bbar[i]); |
---|
| 72 | ssx->pi = xcalloc(1+m, sizeof(mpq_t)); |
---|
| 73 | for (i = 1; i <= m; i++) mpq_init(ssx->pi[i]); |
---|
| 74 | ssx->cbar = xcalloc(1+n, sizeof(mpq_t)); |
---|
| 75 | for (j = 1; j <= n; j++) mpq_init(ssx->cbar[j]); |
---|
| 76 | ssx->rho = xcalloc(1+m, sizeof(mpq_t)); |
---|
| 77 | for (i = 1; i <= m; i++) mpq_init(ssx->rho[i]); |
---|
| 78 | ssx->ap = xcalloc(1+n, sizeof(mpq_t)); |
---|
| 79 | for (j = 1; j <= n; j++) mpq_init(ssx->ap[j]); |
---|
| 80 | ssx->aq = xcalloc(1+m, sizeof(mpq_t)); |
---|
| 81 | for (i = 1; i <= m; i++) mpq_init(ssx->aq[i]); |
---|
| 82 | mpq_init(ssx->delta); |
---|
| 83 | return ssx; |
---|
| 84 | } |
---|
| 85 | |
---|
| 86 | /*---------------------------------------------------------------------- |
---|
| 87 | // ssx_factorize - factorize the current basis matrix. |
---|
| 88 | // |
---|
| 89 | // This routine computes factorization of the current basis matrix B |
---|
| 90 | // and returns the singularity flag. If the matrix B is non-singular, |
---|
| 91 | // the flag is zero, otherwise non-zero. */ |
---|
| 92 | |
---|
| 93 | static int basis_col(void *info, int j, int ind[], mpq_t val[]) |
---|
| 94 | { /* this auxiliary routine provides row indices and numeric values |
---|
| 95 | of non-zero elements in j-th column of the matrix B */ |
---|
| 96 | SSX *ssx = info; |
---|
| 97 | int m = ssx->m; |
---|
| 98 | int n = ssx->n; |
---|
| 99 | int *A_ptr = ssx->A_ptr; |
---|
| 100 | int *A_ind = ssx->A_ind; |
---|
| 101 | mpq_t *A_val = ssx->A_val; |
---|
| 102 | int *Q_col = ssx->Q_col; |
---|
| 103 | int k, len, ptr; |
---|
| 104 | xassert(1 <= j && j <= m); |
---|
| 105 | k = Q_col[j]; /* x[k] = xB[j] */ |
---|
| 106 | xassert(1 <= k && k <= m+n); |
---|
| 107 | /* j-th column of the matrix B is k-th column of the augmented |
---|
| 108 | constraint matrix (I | -A) */ |
---|
| 109 | if (k <= m) |
---|
| 110 | { /* it is a column of the unity matrix I */ |
---|
| 111 | len = 1, ind[1] = k, mpq_set_si(val[1], 1, 1); |
---|
| 112 | } |
---|
| 113 | else |
---|
| 114 | { /* it is a column of the original constraint matrix -A */ |
---|
| 115 | len = 0; |
---|
| 116 | for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++) |
---|
| 117 | { len++; |
---|
| 118 | ind[len] = A_ind[ptr]; |
---|
| 119 | mpq_neg(val[len], A_val[ptr]); |
---|
| 120 | } |
---|
| 121 | } |
---|
| 122 | return len; |
---|
| 123 | } |
---|
| 124 | |
---|
| 125 | int ssx_factorize(SSX *ssx) |
---|
| 126 | { int ret; |
---|
| 127 | ret = bfx_factorize(ssx->binv, ssx->m, basis_col, ssx); |
---|
| 128 | return ret; |
---|
| 129 | } |
---|
| 130 | |
---|
| 131 | /*---------------------------------------------------------------------- |
---|
| 132 | // ssx_get_xNj - determine value of non-basic variable. |
---|
| 133 | // |
---|
| 134 | // This routine determines the value of non-basic variable xN[j] in the |
---|
| 135 | // current basic solution defined as follows: |
---|
| 136 | // |
---|
| 137 | // 0, if xN[j] is free variable |
---|
| 138 | // lN[j], if xN[j] is on its lower bound |
---|
| 139 | // uN[j], if xN[j] is on its upper bound |
---|
| 140 | // lN[j] = uN[j], if xN[j] is fixed variable |
---|
| 141 | // |
---|
| 142 | // where lN[j] and uN[j] are lower and upper bounds of xN[j]. */ |
---|
| 143 | |
---|
| 144 | void ssx_get_xNj(SSX *ssx, int j, mpq_t x) |
---|
| 145 | { int m = ssx->m; |
---|
| 146 | int n = ssx->n; |
---|
| 147 | mpq_t *lb = ssx->lb; |
---|
| 148 | mpq_t *ub = ssx->ub; |
---|
| 149 | int *stat = ssx->stat; |
---|
| 150 | int *Q_col = ssx->Q_col; |
---|
| 151 | int k; |
---|
| 152 | xassert(1 <= j && j <= n); |
---|
| 153 | k = Q_col[m+j]; /* x[k] = xN[j] */ |
---|
| 154 | xassert(1 <= k && k <= m+n); |
---|
| 155 | switch (stat[k]) |
---|
| 156 | { case SSX_NL: |
---|
| 157 | /* xN[j] is on its lower bound */ |
---|
| 158 | mpq_set(x, lb[k]); break; |
---|
| 159 | case SSX_NU: |
---|
| 160 | /* xN[j] is on its upper bound */ |
---|
| 161 | mpq_set(x, ub[k]); break; |
---|
| 162 | case SSX_NF: |
---|
| 163 | /* xN[j] is free variable */ |
---|
| 164 | mpq_set_si(x, 0, 1); break; |
---|
| 165 | case SSX_NS: |
---|
| 166 | /* xN[j] is fixed variable */ |
---|
| 167 | mpq_set(x, lb[k]); break; |
---|
| 168 | default: |
---|
| 169 | xassert(stat != stat); |
---|
| 170 | } |
---|
| 171 | return; |
---|
| 172 | } |
---|
| 173 | |
---|
| 174 | /*---------------------------------------------------------------------- |
---|
| 175 | // ssx_eval_bbar - compute values of basic variables. |
---|
| 176 | // |
---|
| 177 | // This routine computes values of basic variables xB in the current |
---|
| 178 | // basic solution as follows: |
---|
| 179 | // |
---|
| 180 | // beta = - inv(B) * N * xN, |
---|
| 181 | // |
---|
| 182 | // where B is the basis matrix, N is the matrix of non-basic columns, |
---|
| 183 | // xN is a vector of current values of non-basic variables. */ |
---|
| 184 | |
---|
| 185 | void ssx_eval_bbar(SSX *ssx) |
---|
| 186 | { int m = ssx->m; |
---|
| 187 | int n = ssx->n; |
---|
| 188 | mpq_t *coef = ssx->coef; |
---|
| 189 | int *A_ptr = ssx->A_ptr; |
---|
| 190 | int *A_ind = ssx->A_ind; |
---|
| 191 | mpq_t *A_val = ssx->A_val; |
---|
| 192 | int *Q_col = ssx->Q_col; |
---|
| 193 | mpq_t *bbar = ssx->bbar; |
---|
| 194 | int i, j, k, ptr; |
---|
| 195 | mpq_t x, temp; |
---|
| 196 | mpq_init(x); |
---|
| 197 | mpq_init(temp); |
---|
| 198 | /* bbar := 0 */ |
---|
| 199 | for (i = 1; i <= m; i++) |
---|
| 200 | mpq_set_si(bbar[i], 0, 1); |
---|
| 201 | /* bbar := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n] */ |
---|
| 202 | for (j = 1; j <= n; j++) |
---|
| 203 | { ssx_get_xNj(ssx, j, x); |
---|
| 204 | if (mpq_sgn(x) == 0) continue; |
---|
| 205 | k = Q_col[m+j]; /* x[k] = xN[j] */ |
---|
| 206 | if (k <= m) |
---|
| 207 | { /* N[j] is a column of the unity matrix I */ |
---|
| 208 | mpq_sub(bbar[k], bbar[k], x); |
---|
| 209 | } |
---|
| 210 | else |
---|
| 211 | { /* N[j] is a column of the original constraint matrix -A */ |
---|
| 212 | for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++) |
---|
| 213 | { mpq_mul(temp, A_val[ptr], x); |
---|
| 214 | mpq_add(bbar[A_ind[ptr]], bbar[A_ind[ptr]], temp); |
---|
| 215 | } |
---|
| 216 | } |
---|
| 217 | } |
---|
| 218 | /* bbar := inv(B) * bbar */ |
---|
| 219 | bfx_ftran(ssx->binv, bbar, 0); |
---|
| 220 | #if 1 |
---|
| 221 | /* compute value of the objective function */ |
---|
| 222 | /* bbar[0] := c[0] */ |
---|
| 223 | mpq_set(bbar[0], coef[0]); |
---|
| 224 | /* bbar[0] := bbar[0] + sum{i in B} cB[i] * xB[i] */ |
---|
| 225 | for (i = 1; i <= m; i++) |
---|
| 226 | { k = Q_col[i]; /* x[k] = xB[i] */ |
---|
| 227 | if (mpq_sgn(coef[k]) == 0) continue; |
---|
| 228 | mpq_mul(temp, coef[k], bbar[i]); |
---|
| 229 | mpq_add(bbar[0], bbar[0], temp); |
---|
| 230 | } |
---|
| 231 | /* bbar[0] := bbar[0] + sum{j in N} cN[j] * xN[j] */ |
---|
| 232 | for (j = 1; j <= n; j++) |
---|
| 233 | { k = Q_col[m+j]; /* x[k] = xN[j] */ |
---|
| 234 | if (mpq_sgn(coef[k]) == 0) continue; |
---|
| 235 | ssx_get_xNj(ssx, j, x); |
---|
| 236 | mpq_mul(temp, coef[k], x); |
---|
| 237 | mpq_add(bbar[0], bbar[0], temp); |
---|
| 238 | } |
---|
| 239 | #endif |
---|
| 240 | mpq_clear(x); |
---|
| 241 | mpq_clear(temp); |
---|
| 242 | return; |
---|
| 243 | } |
---|
| 244 | |
---|
| 245 | /*---------------------------------------------------------------------- |
---|
| 246 | // ssx_eval_pi - compute values of simplex multipliers. |
---|
| 247 | // |
---|
| 248 | // This routine computes values of simplex multipliers (shadow prices) |
---|
| 249 | // pi in the current basic solution as follows: |
---|
| 250 | // |
---|
| 251 | // pi = inv(B') * cB, |
---|
| 252 | // |
---|
| 253 | // where B' is a matrix transposed to the basis matrix B, cB is a vector |
---|
| 254 | // of objective coefficients at basic variables xB. */ |
---|
| 255 | |
---|
| 256 | void ssx_eval_pi(SSX *ssx) |
---|
| 257 | { int m = ssx->m; |
---|
| 258 | mpq_t *coef = ssx->coef; |
---|
| 259 | int *Q_col = ssx->Q_col; |
---|
| 260 | mpq_t *pi = ssx->pi; |
---|
| 261 | int i; |
---|
| 262 | /* pi := cB */ |
---|
| 263 | for (i = 1; i <= m; i++) mpq_set(pi[i], coef[Q_col[i]]); |
---|
| 264 | /* pi := inv(B') * cB */ |
---|
| 265 | bfx_btran(ssx->binv, pi); |
---|
| 266 | return; |
---|
| 267 | } |
---|
| 268 | |
---|
| 269 | /*---------------------------------------------------------------------- |
---|
| 270 | // ssx_eval_dj - compute reduced cost of non-basic variable. |
---|
| 271 | // |
---|
| 272 | // This routine computes reduced cost d[j] of non-basic variable xN[j] |
---|
| 273 | // in the current basic solution as follows: |
---|
| 274 | // |
---|
| 275 | // d[j] = cN[j] - N[j] * pi, |
---|
| 276 | // |
---|
| 277 | // where cN[j] is an objective coefficient at xN[j], N[j] is a column |
---|
| 278 | // of the augmented constraint matrix (I | -A) corresponding to xN[j], |
---|
| 279 | // pi is the vector of simplex multipliers (shadow prices). */ |
---|
| 280 | |
---|
| 281 | void ssx_eval_dj(SSX *ssx, int j, mpq_t dj) |
---|
| 282 | { int m = ssx->m; |
---|
| 283 | int n = ssx->n; |
---|
| 284 | mpq_t *coef = ssx->coef; |
---|
| 285 | int *A_ptr = ssx->A_ptr; |
---|
| 286 | int *A_ind = ssx->A_ind; |
---|
| 287 | mpq_t *A_val = ssx->A_val; |
---|
| 288 | int *Q_col = ssx->Q_col; |
---|
| 289 | mpq_t *pi = ssx->pi; |
---|
| 290 | int k, ptr, end; |
---|
| 291 | mpq_t temp; |
---|
| 292 | mpq_init(temp); |
---|
| 293 | xassert(1 <= j && j <= n); |
---|
| 294 | k = Q_col[m+j]; /* x[k] = xN[j] */ |
---|
| 295 | xassert(1 <= k && k <= m+n); |
---|
| 296 | /* j-th column of the matrix N is k-th column of the augmented |
---|
| 297 | constraint matrix (I | -A) */ |
---|
| 298 | if (k <= m) |
---|
| 299 | { /* it is a column of the unity matrix I */ |
---|
| 300 | mpq_sub(dj, coef[k], pi[k]); |
---|
| 301 | } |
---|
| 302 | else |
---|
| 303 | { /* it is a column of the original constraint matrix -A */ |
---|
| 304 | mpq_set(dj, coef[k]); |
---|
| 305 | for (ptr = A_ptr[k-m], end = A_ptr[k-m+1]; ptr < end; ptr++) |
---|
| 306 | { mpq_mul(temp, A_val[ptr], pi[A_ind[ptr]]); |
---|
| 307 | mpq_add(dj, dj, temp); |
---|
| 308 | } |
---|
| 309 | } |
---|
| 310 | mpq_clear(temp); |
---|
| 311 | return; |
---|
| 312 | } |
---|
| 313 | |
---|
| 314 | /*---------------------------------------------------------------------- |
---|
| 315 | // ssx_eval_cbar - compute reduced costs of all non-basic variables. |
---|
| 316 | // |
---|
| 317 | // This routine computes the vector of reduced costs pi in the current |
---|
| 318 | // basic solution for all non-basic variables, including fixed ones. */ |
---|
| 319 | |
---|
| 320 | void ssx_eval_cbar(SSX *ssx) |
---|
| 321 | { int n = ssx->n; |
---|
| 322 | mpq_t *cbar = ssx->cbar; |
---|
| 323 | int j; |
---|
| 324 | for (j = 1; j <= n; j++) |
---|
| 325 | ssx_eval_dj(ssx, j, cbar[j]); |
---|
| 326 | return; |
---|
| 327 | } |
---|
| 328 | |
---|
| 329 | /*---------------------------------------------------------------------- |
---|
| 330 | // ssx_eval_rho - compute p-th row of the inverse. |
---|
| 331 | // |
---|
| 332 | // This routine computes p-th row of the matrix inv(B), where B is the |
---|
| 333 | // current basis matrix. |
---|
| 334 | // |
---|
| 335 | // p-th row of the inverse is computed using the following formula: |
---|
| 336 | // |
---|
| 337 | // rho = inv(B') * e[p], |
---|
| 338 | // |
---|
| 339 | // where B' is a matrix transposed to B, e[p] is a unity vector, which |
---|
| 340 | // contains one in p-th position. */ |
---|
| 341 | |
---|
| 342 | void ssx_eval_rho(SSX *ssx) |
---|
| 343 | { int m = ssx->m; |
---|
| 344 | int p = ssx->p; |
---|
| 345 | mpq_t *rho = ssx->rho; |
---|
| 346 | int i; |
---|
| 347 | xassert(1 <= p && p <= m); |
---|
| 348 | /* rho := 0 */ |
---|
| 349 | for (i = 1; i <= m; i++) mpq_set_si(rho[i], 0, 1); |
---|
| 350 | /* rho := e[p] */ |
---|
| 351 | mpq_set_si(rho[p], 1, 1); |
---|
| 352 | /* rho := inv(B') * rho */ |
---|
| 353 | bfx_btran(ssx->binv, rho); |
---|
| 354 | return; |
---|
| 355 | } |
---|
| 356 | |
---|
| 357 | /*---------------------------------------------------------------------- |
---|
| 358 | // ssx_eval_row - compute pivot row of the simplex table. |
---|
| 359 | // |
---|
| 360 | // This routine computes p-th (pivot) row of the current simplex table |
---|
| 361 | // A~ = - inv(B) * N using the following formula: |
---|
| 362 | // |
---|
| 363 | // A~[p] = - N' * inv(B') * e[p] = - N' * rho[p], |
---|
| 364 | // |
---|
| 365 | // where N' is a matrix transposed to the matrix N, rho[p] is p-th row |
---|
| 366 | // of the inverse inv(B). */ |
---|
| 367 | |
---|
| 368 | void ssx_eval_row(SSX *ssx) |
---|
| 369 | { int m = ssx->m; |
---|
| 370 | int n = ssx->n; |
---|
| 371 | int *A_ptr = ssx->A_ptr; |
---|
| 372 | int *A_ind = ssx->A_ind; |
---|
| 373 | mpq_t *A_val = ssx->A_val; |
---|
| 374 | int *Q_col = ssx->Q_col; |
---|
| 375 | mpq_t *rho = ssx->rho; |
---|
| 376 | mpq_t *ap = ssx->ap; |
---|
| 377 | int j, k, ptr; |
---|
| 378 | mpq_t temp; |
---|
| 379 | mpq_init(temp); |
---|
| 380 | for (j = 1; j <= n; j++) |
---|
| 381 | { /* ap[j] := - N'[j] * rho (inner product) */ |
---|
| 382 | k = Q_col[m+j]; /* x[k] = xN[j] */ |
---|
| 383 | if (k <= m) |
---|
| 384 | mpq_neg(ap[j], rho[k]); |
---|
| 385 | else |
---|
| 386 | { mpq_set_si(ap[j], 0, 1); |
---|
| 387 | for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++) |
---|
| 388 | { mpq_mul(temp, A_val[ptr], rho[A_ind[ptr]]); |
---|
| 389 | mpq_add(ap[j], ap[j], temp); |
---|
| 390 | } |
---|
| 391 | } |
---|
| 392 | } |
---|
| 393 | mpq_clear(temp); |
---|
| 394 | return; |
---|
| 395 | } |
---|
| 396 | |
---|
| 397 | /*---------------------------------------------------------------------- |
---|
| 398 | // ssx_eval_col - compute pivot column of the simplex table. |
---|
| 399 | // |
---|
| 400 | // This routine computes q-th (pivot) column of the current simplex |
---|
| 401 | // table A~ = - inv(B) * N using the following formula: |
---|
| 402 | // |
---|
| 403 | // A~[q] = - inv(B) * N[q], |
---|
| 404 | // |
---|
| 405 | // where N[q] is q-th column of the matrix N corresponding to chosen |
---|
| 406 | // non-basic variable xN[q]. */ |
---|
| 407 | |
---|
| 408 | void ssx_eval_col(SSX *ssx) |
---|
| 409 | { int m = ssx->m; |
---|
| 410 | int n = ssx->n; |
---|
| 411 | int *A_ptr = ssx->A_ptr; |
---|
| 412 | int *A_ind = ssx->A_ind; |
---|
| 413 | mpq_t *A_val = ssx->A_val; |
---|
| 414 | int *Q_col = ssx->Q_col; |
---|
| 415 | int q = ssx->q; |
---|
| 416 | mpq_t *aq = ssx->aq; |
---|
| 417 | int i, k, ptr; |
---|
| 418 | xassert(1 <= q && q <= n); |
---|
| 419 | /* aq := 0 */ |
---|
| 420 | for (i = 1; i <= m; i++) mpq_set_si(aq[i], 0, 1); |
---|
| 421 | /* aq := N[q] */ |
---|
| 422 | k = Q_col[m+q]; /* x[k] = xN[q] */ |
---|
| 423 | if (k <= m) |
---|
| 424 | { /* N[q] is a column of the unity matrix I */ |
---|
| 425 | mpq_set_si(aq[k], 1, 1); |
---|
| 426 | } |
---|
| 427 | else |
---|
| 428 | { /* N[q] is a column of the original constraint matrix -A */ |
---|
| 429 | for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++) |
---|
| 430 | mpq_neg(aq[A_ind[ptr]], A_val[ptr]); |
---|
| 431 | } |
---|
| 432 | /* aq := inv(B) * aq */ |
---|
| 433 | bfx_ftran(ssx->binv, aq, 1); |
---|
| 434 | /* aq := - aq */ |
---|
| 435 | for (i = 1; i <= m; i++) mpq_neg(aq[i], aq[i]); |
---|
| 436 | return; |
---|
| 437 | } |
---|
| 438 | |
---|
| 439 | /*---------------------------------------------------------------------- |
---|
| 440 | // ssx_chuzc - choose pivot column. |
---|
| 441 | // |
---|
| 442 | // This routine chooses non-basic variable xN[q] whose reduced cost |
---|
| 443 | // indicates possible improving of the objective function to enter it |
---|
| 444 | // in the basis. |
---|
| 445 | // |
---|
| 446 | // Currently the standard (textbook) pricing is used, i.e. that |
---|
| 447 | // non-basic variable is preferred which has greatest reduced cost (in |
---|
| 448 | // magnitude). |
---|
| 449 | // |
---|
| 450 | // If xN[q] has been chosen, the routine stores its number q and also |
---|
| 451 | // sets the flag q_dir that indicates direction in which xN[q] has to |
---|
| 452 | // change (+1 means increasing, -1 means decreasing). |
---|
| 453 | // |
---|
| 454 | // If the choice cannot be made, because the current basic solution is |
---|
| 455 | // dual feasible, the routine sets the number q to 0. */ |
---|
| 456 | |
---|
| 457 | void ssx_chuzc(SSX *ssx) |
---|
| 458 | { int m = ssx->m; |
---|
| 459 | int n = ssx->n; |
---|
| 460 | int dir = (ssx->dir == SSX_MIN ? +1 : -1); |
---|
| 461 | int *Q_col = ssx->Q_col; |
---|
| 462 | int *stat = ssx->stat; |
---|
| 463 | mpq_t *cbar = ssx->cbar; |
---|
| 464 | int j, k, s, q, q_dir; |
---|
| 465 | double best, temp; |
---|
| 466 | /* nothing is chosen so far */ |
---|
| 467 | q = 0, q_dir = 0, best = 0.0; |
---|
| 468 | /* look through the list of non-basic variables */ |
---|
| 469 | for (j = 1; j <= n; j++) |
---|
| 470 | { k = Q_col[m+j]; /* x[k] = xN[j] */ |
---|
| 471 | s = dir * mpq_sgn(cbar[j]); |
---|
| 472 | if ((stat[k] == SSX_NF || stat[k] == SSX_NL) && s < 0 || |
---|
| 473 | (stat[k] == SSX_NF || stat[k] == SSX_NU) && s > 0) |
---|
| 474 | { /* reduced cost of xN[j] indicates possible improving of |
---|
| 475 | the objective function */ |
---|
| 476 | temp = fabs(mpq_get_d(cbar[j])); |
---|
| 477 | xassert(temp != 0.0); |
---|
| 478 | if (q == 0 || best < temp) |
---|
| 479 | q = j, q_dir = - s, best = temp; |
---|
| 480 | } |
---|
| 481 | } |
---|
| 482 | ssx->q = q, ssx->q_dir = q_dir; |
---|
| 483 | return; |
---|
| 484 | } |
---|
| 485 | |
---|
| 486 | /*---------------------------------------------------------------------- |
---|
| 487 | // ssx_chuzr - choose pivot row. |
---|
| 488 | // |
---|
| 489 | // This routine looks through elements of q-th column of the simplex |
---|
| 490 | // table and chooses basic variable xB[p] which should leave the basis. |
---|
| 491 | // |
---|
| 492 | // The choice is based on the standard (textbook) ratio test. |
---|
| 493 | // |
---|
| 494 | // If xB[p] has been chosen, the routine stores its number p and also |
---|
| 495 | // sets its non-basic status p_stat which should be assigned to xB[p] |
---|
| 496 | // when it has left the basis and become xN[q]. |
---|
| 497 | // |
---|
| 498 | // Special case p < 0 means that xN[q] is double-bounded variable and |
---|
| 499 | // it reaches its opposite bound before any basic variable does that, |
---|
| 500 | // so the current basis remains unchanged. |
---|
| 501 | // |
---|
| 502 | // If the choice cannot be made, because xN[q] can infinitely change in |
---|
| 503 | // the feasible direction, the routine sets the number p to 0. */ |
---|
| 504 | |
---|
| 505 | void ssx_chuzr(SSX *ssx) |
---|
| 506 | { int m = ssx->m; |
---|
| 507 | int n = ssx->n; |
---|
| 508 | int *type = ssx->type; |
---|
| 509 | mpq_t *lb = ssx->lb; |
---|
| 510 | mpq_t *ub = ssx->ub; |
---|
| 511 | int *Q_col = ssx->Q_col; |
---|
| 512 | mpq_t *bbar = ssx->bbar; |
---|
| 513 | int q = ssx->q; |
---|
| 514 | mpq_t *aq = ssx->aq; |
---|
| 515 | int q_dir = ssx->q_dir; |
---|
| 516 | int i, k, s, t, p, p_stat; |
---|
| 517 | mpq_t teta, temp; |
---|
| 518 | mpq_init(teta); |
---|
| 519 | mpq_init(temp); |
---|
| 520 | xassert(1 <= q && q <= n); |
---|
| 521 | xassert(q_dir == +1 || q_dir == -1); |
---|
| 522 | /* nothing is chosen so far */ |
---|
| 523 | p = 0, p_stat = 0; |
---|
| 524 | /* look through the list of basic variables */ |
---|
| 525 | for (i = 1; i <= m; i++) |
---|
| 526 | { s = q_dir * mpq_sgn(aq[i]); |
---|
| 527 | if (s < 0) |
---|
| 528 | { /* xB[i] decreases */ |
---|
| 529 | k = Q_col[i]; /* x[k] = xB[i] */ |
---|
| 530 | t = type[k]; |
---|
| 531 | if (t == SSX_LO || t == SSX_DB || t == SSX_FX) |
---|
| 532 | { /* xB[i] has finite lower bound */ |
---|
| 533 | mpq_sub(temp, bbar[i], lb[k]); |
---|
| 534 | mpq_div(temp, temp, aq[i]); |
---|
| 535 | mpq_abs(temp, temp); |
---|
| 536 | if (p == 0 || mpq_cmp(teta, temp) > 0) |
---|
| 537 | { p = i; |
---|
| 538 | p_stat = (t == SSX_FX ? SSX_NS : SSX_NL); |
---|
| 539 | mpq_set(teta, temp); |
---|
| 540 | } |
---|
| 541 | } |
---|
| 542 | } |
---|
| 543 | else if (s > 0) |
---|
| 544 | { /* xB[i] increases */ |
---|
| 545 | k = Q_col[i]; /* x[k] = xB[i] */ |
---|
| 546 | t = type[k]; |
---|
| 547 | if (t == SSX_UP || t == SSX_DB || t == SSX_FX) |
---|
| 548 | { /* xB[i] has finite upper bound */ |
---|
| 549 | mpq_sub(temp, bbar[i], ub[k]); |
---|
| 550 | mpq_div(temp, temp, aq[i]); |
---|
| 551 | mpq_abs(temp, temp); |
---|
| 552 | if (p == 0 || mpq_cmp(teta, temp) > 0) |
---|
| 553 | { p = i; |
---|
| 554 | p_stat = (t == SSX_FX ? SSX_NS : SSX_NU); |
---|
| 555 | mpq_set(teta, temp); |
---|
| 556 | } |
---|
| 557 | } |
---|
| 558 | } |
---|
| 559 | /* if something has been chosen and the ratio test indicates |
---|
| 560 | exact degeneracy, the search can be finished */ |
---|
| 561 | if (p != 0 && mpq_sgn(teta) == 0) break; |
---|
| 562 | } |
---|
| 563 | /* if xN[q] is double-bounded, check if it can reach its opposite |
---|
| 564 | bound before any basic variable */ |
---|
| 565 | k = Q_col[m+q]; /* x[k] = xN[q] */ |
---|
| 566 | if (type[k] == SSX_DB) |
---|
| 567 | { mpq_sub(temp, ub[k], lb[k]); |
---|
| 568 | if (p == 0 || mpq_cmp(teta, temp) > 0) |
---|
| 569 | { p = -1; |
---|
| 570 | p_stat = -1; |
---|
| 571 | mpq_set(teta, temp); |
---|
| 572 | } |
---|
| 573 | } |
---|
| 574 | ssx->p = p; |
---|
| 575 | ssx->p_stat = p_stat; |
---|
| 576 | /* if xB[p] has been chosen, determine its actual change in the |
---|
| 577 | adjacent basis (it has the same sign as q_dir) */ |
---|
| 578 | if (p != 0) |
---|
| 579 | { xassert(mpq_sgn(teta) >= 0); |
---|
| 580 | if (q_dir > 0) |
---|
| 581 | mpq_set(ssx->delta, teta); |
---|
| 582 | else |
---|
| 583 | mpq_neg(ssx->delta, teta); |
---|
| 584 | } |
---|
| 585 | mpq_clear(teta); |
---|
| 586 | mpq_clear(temp); |
---|
| 587 | return; |
---|
| 588 | } |
---|
| 589 | |
---|
| 590 | /*---------------------------------------------------------------------- |
---|
| 591 | // ssx_update_bbar - update values of basic variables. |
---|
| 592 | // |
---|
| 593 | // This routine recomputes the current values of basic variables for |
---|
| 594 | // the adjacent basis. |
---|
| 595 | // |
---|
| 596 | // The simplex table for the current basis is the following: |
---|
| 597 | // |
---|
| 598 | // xB[i] = sum{j in 1..n} alfa[i,j] * xN[q], i = 1,...,m |
---|
| 599 | // |
---|
| 600 | // therefore |
---|
| 601 | // |
---|
| 602 | // delta xB[i] = alfa[i,q] * delta xN[q], i = 1,...,m |
---|
| 603 | // |
---|
| 604 | // where delta xN[q] = xN.new[q] - xN[q] is the change of xN[q] in the |
---|
| 605 | // adjacent basis, and delta xB[i] = xB.new[i] - xB[i] is the change of |
---|
| 606 | // xB[i]. This gives formulae for recomputing values of xB[i]: |
---|
| 607 | // |
---|
| 608 | // xB.new[p] = xN[q] + delta xN[q] |
---|
| 609 | // |
---|
| 610 | // (because xN[q] becomes xB[p] in the adjacent basis), and |
---|
| 611 | // |
---|
| 612 | // xB.new[i] = xB[i] + alfa[i,q] * delta xN[q], i != p |
---|
| 613 | // |
---|
| 614 | // for other basic variables. */ |
---|
| 615 | |
---|
| 616 | void ssx_update_bbar(SSX *ssx) |
---|
| 617 | { int m = ssx->m; |
---|
| 618 | int n = ssx->n; |
---|
| 619 | mpq_t *bbar = ssx->bbar; |
---|
| 620 | mpq_t *cbar = ssx->cbar; |
---|
| 621 | int p = ssx->p; |
---|
| 622 | int q = ssx->q; |
---|
| 623 | mpq_t *aq = ssx->aq; |
---|
| 624 | int i; |
---|
| 625 | mpq_t temp; |
---|
| 626 | mpq_init(temp); |
---|
| 627 | xassert(1 <= q && q <= n); |
---|
| 628 | if (p < 0) |
---|
| 629 | { /* xN[q] is double-bounded and goes to its opposite bound */ |
---|
| 630 | /* nop */; |
---|
| 631 | } |
---|
| 632 | else |
---|
| 633 | { /* xN[q] becomes xB[p] in the adjacent basis */ |
---|
| 634 | /* xB.new[p] = xN[q] + delta xN[q] */ |
---|
| 635 | xassert(1 <= p && p <= m); |
---|
| 636 | ssx_get_xNj(ssx, q, temp); |
---|
| 637 | mpq_add(bbar[p], temp, ssx->delta); |
---|
| 638 | } |
---|
| 639 | /* update values of other basic variables depending on xN[q] */ |
---|
| 640 | for (i = 1; i <= m; i++) |
---|
| 641 | { if (i == p) continue; |
---|
| 642 | /* xB.new[i] = xB[i] + alfa[i,q] * delta xN[q] */ |
---|
| 643 | if (mpq_sgn(aq[i]) == 0) continue; |
---|
| 644 | mpq_mul(temp, aq[i], ssx->delta); |
---|
| 645 | mpq_add(bbar[i], bbar[i], temp); |
---|
| 646 | } |
---|
| 647 | #if 1 |
---|
| 648 | /* update value of the objective function */ |
---|
| 649 | /* z.new = z + d[q] * delta xN[q] */ |
---|
| 650 | mpq_mul(temp, cbar[q], ssx->delta); |
---|
| 651 | mpq_add(bbar[0], bbar[0], temp); |
---|
| 652 | #endif |
---|
| 653 | mpq_clear(temp); |
---|
| 654 | return; |
---|
| 655 | } |
---|
| 656 | |
---|
| 657 | /*---------------------------------------------------------------------- |
---|
| 658 | -- ssx_update_pi - update simplex multipliers. |
---|
| 659 | -- |
---|
| 660 | -- This routine recomputes the vector of simplex multipliers for the |
---|
| 661 | -- adjacent basis. */ |
---|
| 662 | |
---|
| 663 | void ssx_update_pi(SSX *ssx) |
---|
| 664 | { int m = ssx->m; |
---|
| 665 | int n = ssx->n; |
---|
| 666 | mpq_t *pi = ssx->pi; |
---|
| 667 | mpq_t *cbar = ssx->cbar; |
---|
| 668 | int p = ssx->p; |
---|
| 669 | int q = ssx->q; |
---|
| 670 | mpq_t *aq = ssx->aq; |
---|
| 671 | mpq_t *rho = ssx->rho; |
---|
| 672 | int i; |
---|
| 673 | mpq_t new_dq, temp; |
---|
| 674 | mpq_init(new_dq); |
---|
| 675 | mpq_init(temp); |
---|
| 676 | xassert(1 <= p && p <= m); |
---|
| 677 | xassert(1 <= q && q <= n); |
---|
| 678 | /* compute d[q] in the adjacent basis */ |
---|
| 679 | mpq_div(new_dq, cbar[q], aq[p]); |
---|
| 680 | /* update the vector of simplex multipliers */ |
---|
| 681 | for (i = 1; i <= m; i++) |
---|
| 682 | { if (mpq_sgn(rho[i]) == 0) continue; |
---|
| 683 | mpq_mul(temp, new_dq, rho[i]); |
---|
| 684 | mpq_sub(pi[i], pi[i], temp); |
---|
| 685 | } |
---|
| 686 | mpq_clear(new_dq); |
---|
| 687 | mpq_clear(temp); |
---|
| 688 | return; |
---|
| 689 | } |
---|
| 690 | |
---|
| 691 | /*---------------------------------------------------------------------- |
---|
| 692 | // ssx_update_cbar - update reduced costs of non-basic variables. |
---|
| 693 | // |
---|
| 694 | // This routine recomputes the vector of reduced costs of non-basic |
---|
| 695 | // variables for the adjacent basis. */ |
---|
| 696 | |
---|
| 697 | void ssx_update_cbar(SSX *ssx) |
---|
| 698 | { int m = ssx->m; |
---|
| 699 | int n = ssx->n; |
---|
| 700 | mpq_t *cbar = ssx->cbar; |
---|
| 701 | int p = ssx->p; |
---|
| 702 | int q = ssx->q; |
---|
| 703 | mpq_t *ap = ssx->ap; |
---|
| 704 | int j; |
---|
| 705 | mpq_t temp; |
---|
| 706 | mpq_init(temp); |
---|
| 707 | xassert(1 <= p && p <= m); |
---|
| 708 | xassert(1 <= q && q <= n); |
---|
| 709 | /* compute d[q] in the adjacent basis */ |
---|
| 710 | /* d.new[q] = d[q] / alfa[p,q] */ |
---|
| 711 | mpq_div(cbar[q], cbar[q], ap[q]); |
---|
| 712 | /* update reduced costs of other non-basic variables */ |
---|
| 713 | for (j = 1; j <= n; j++) |
---|
| 714 | { if (j == q) continue; |
---|
| 715 | /* d.new[j] = d[j] - (alfa[p,j] / alfa[p,q]) * d[q] */ |
---|
| 716 | if (mpq_sgn(ap[j]) == 0) continue; |
---|
| 717 | mpq_mul(temp, ap[j], cbar[q]); |
---|
| 718 | mpq_sub(cbar[j], cbar[j], temp); |
---|
| 719 | } |
---|
| 720 | mpq_clear(temp); |
---|
| 721 | return; |
---|
| 722 | } |
---|
| 723 | |
---|
| 724 | /*---------------------------------------------------------------------- |
---|
| 725 | // ssx_change_basis - change current basis to adjacent one. |
---|
| 726 | // |
---|
| 727 | // This routine changes the current basis to the adjacent one swapping |
---|
| 728 | // basic variable xB[p] and non-basic variable xN[q]. */ |
---|
| 729 | |
---|
| 730 | void ssx_change_basis(SSX *ssx) |
---|
| 731 | { int m = ssx->m; |
---|
| 732 | int n = ssx->n; |
---|
| 733 | int *type = ssx->type; |
---|
| 734 | int *stat = ssx->stat; |
---|
| 735 | int *Q_row = ssx->Q_row; |
---|
| 736 | int *Q_col = ssx->Q_col; |
---|
| 737 | int p = ssx->p; |
---|
| 738 | int q = ssx->q; |
---|
| 739 | int p_stat = ssx->p_stat; |
---|
| 740 | int k, kp, kq; |
---|
| 741 | if (p < 0) |
---|
| 742 | { /* special case: xN[q] goes to its opposite bound */ |
---|
| 743 | xassert(1 <= q && q <= n); |
---|
| 744 | k = Q_col[m+q]; /* x[k] = xN[q] */ |
---|
| 745 | xassert(type[k] == SSX_DB); |
---|
| 746 | switch (stat[k]) |
---|
| 747 | { case SSX_NL: |
---|
| 748 | stat[k] = SSX_NU; |
---|
| 749 | break; |
---|
| 750 | case SSX_NU: |
---|
| 751 | stat[k] = SSX_NL; |
---|
| 752 | break; |
---|
| 753 | default: |
---|
| 754 | xassert(stat != stat); |
---|
| 755 | } |
---|
| 756 | } |
---|
| 757 | else |
---|
| 758 | { /* xB[p] leaves the basis, xN[q] enters the basis */ |
---|
| 759 | xassert(1 <= p && p <= m); |
---|
| 760 | xassert(1 <= q && q <= n); |
---|
| 761 | kp = Q_col[p]; /* x[kp] = xB[p] */ |
---|
| 762 | kq = Q_col[m+q]; /* x[kq] = xN[q] */ |
---|
| 763 | /* check non-basic status of xB[p] which becomes xN[q] */ |
---|
| 764 | switch (type[kp]) |
---|
| 765 | { case SSX_FR: |
---|
| 766 | xassert(p_stat == SSX_NF); |
---|
| 767 | break; |
---|
| 768 | case SSX_LO: |
---|
| 769 | xassert(p_stat == SSX_NL); |
---|
| 770 | break; |
---|
| 771 | case SSX_UP: |
---|
| 772 | xassert(p_stat == SSX_NU); |
---|
| 773 | break; |
---|
| 774 | case SSX_DB: |
---|
| 775 | xassert(p_stat == SSX_NL || p_stat == SSX_NU); |
---|
| 776 | break; |
---|
| 777 | case SSX_FX: |
---|
| 778 | xassert(p_stat == SSX_NS); |
---|
| 779 | break; |
---|
| 780 | default: |
---|
| 781 | xassert(type != type); |
---|
| 782 | } |
---|
| 783 | /* swap xB[p] and xN[q] */ |
---|
| 784 | stat[kp] = (char)p_stat, stat[kq] = SSX_BS; |
---|
| 785 | Q_row[kp] = m+q, Q_row[kq] = p; |
---|
| 786 | Q_col[p] = kq, Q_col[m+q] = kp; |
---|
| 787 | /* update factorization of the basis matrix */ |
---|
| 788 | if (bfx_update(ssx->binv, p)) |
---|
| 789 | { if (ssx_factorize(ssx)) |
---|
| 790 | xassert(("Internal error: basis matrix is singular", 0)); |
---|
| 791 | } |
---|
| 792 | } |
---|
| 793 | return; |
---|
| 794 | } |
---|
| 795 | |
---|
| 796 | /*---------------------------------------------------------------------- |
---|
| 797 | // ssx_delete - delete simplex solver workspace. |
---|
| 798 | // |
---|
| 799 | // This routine deletes the simplex solver workspace freeing all the |
---|
| 800 | // memory allocated to this object. */ |
---|
| 801 | |
---|
| 802 | void ssx_delete(SSX *ssx) |
---|
| 803 | { int m = ssx->m; |
---|
| 804 | int n = ssx->n; |
---|
| 805 | int nnz = ssx->A_ptr[n+1]-1; |
---|
| 806 | int i, j, k; |
---|
| 807 | xfree(ssx->type); |
---|
| 808 | for (k = 1; k <= m+n; k++) mpq_clear(ssx->lb[k]); |
---|
| 809 | xfree(ssx->lb); |
---|
| 810 | for (k = 1; k <= m+n; k++) mpq_clear(ssx->ub[k]); |
---|
| 811 | xfree(ssx->ub); |
---|
| 812 | for (k = 0; k <= m+n; k++) mpq_clear(ssx->coef[k]); |
---|
| 813 | xfree(ssx->coef); |
---|
| 814 | xfree(ssx->A_ptr); |
---|
| 815 | xfree(ssx->A_ind); |
---|
| 816 | for (k = 1; k <= nnz; k++) mpq_clear(ssx->A_val[k]); |
---|
| 817 | xfree(ssx->A_val); |
---|
| 818 | xfree(ssx->stat); |
---|
| 819 | xfree(ssx->Q_row); |
---|
| 820 | xfree(ssx->Q_col); |
---|
| 821 | bfx_delete_binv(ssx->binv); |
---|
| 822 | for (i = 0; i <= m; i++) mpq_clear(ssx->bbar[i]); |
---|
| 823 | xfree(ssx->bbar); |
---|
| 824 | for (i = 1; i <= m; i++) mpq_clear(ssx->pi[i]); |
---|
| 825 | xfree(ssx->pi); |
---|
| 826 | for (j = 1; j <= n; j++) mpq_clear(ssx->cbar[j]); |
---|
| 827 | xfree(ssx->cbar); |
---|
| 828 | for (i = 1; i <= m; i++) mpq_clear(ssx->rho[i]); |
---|
| 829 | xfree(ssx->rho); |
---|
| 830 | for (j = 1; j <= n; j++) mpq_clear(ssx->ap[j]); |
---|
| 831 | xfree(ssx->ap); |
---|
| 832 | for (i = 1; i <= m; i++) mpq_clear(ssx->aq[i]); |
---|
| 833 | xfree(ssx->aq); |
---|
| 834 | mpq_clear(ssx->delta); |
---|
| 835 | xfree(ssx); |
---|
| 836 | return; |
---|
| 837 | } |
---|
| 838 | |
---|
| 839 | /* eof */ |
---|