lemon-project-template-glpk

comparison deps/glpk/src/glpssx01.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6bb69fc997ad
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, 2011 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 */