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 */ |
---|