lemon-project-template-glpk

annotate deps/glpk/src/glpssx01.c @ 11:4fc6ad2fb8a6

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