1 | /* glpssx02.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 | |
---|
28 | static void show_progress(SSX *ssx, int phase) |
---|
29 | { /* this auxiliary routine displays information about progress of |
---|
30 | the search */ |
---|
31 | int i, def = 0; |
---|
32 | for (i = 1; i <= ssx->m; i++) |
---|
33 | if (ssx->type[ssx->Q_col[i]] == SSX_FX) def++; |
---|
34 | xprintf("%s%6d: %s = %22.15g (%d)\n", phase == 1 ? " " : "*", |
---|
35 | ssx->it_cnt, phase == 1 ? "infsum" : "objval", |
---|
36 | mpq_get_d(ssx->bbar[0]), def); |
---|
37 | #if 0 |
---|
38 | ssx->tm_lag = utime(); |
---|
39 | #else |
---|
40 | ssx->tm_lag = xtime(); |
---|
41 | #endif |
---|
42 | return; |
---|
43 | } |
---|
44 | |
---|
45 | /*---------------------------------------------------------------------- |
---|
46 | // ssx_phase_I - find primal feasible solution. |
---|
47 | // |
---|
48 | // This routine implements phase I of the primal simplex method. |
---|
49 | // |
---|
50 | // On exit the routine returns one of the following codes: |
---|
51 | // |
---|
52 | // 0 - feasible solution found; |
---|
53 | // 1 - problem has no feasible solution; |
---|
54 | // 2 - iterations limit exceeded; |
---|
55 | // 3 - time limit exceeded. |
---|
56 | ----------------------------------------------------------------------*/ |
---|
57 | |
---|
58 | int ssx_phase_I(SSX *ssx) |
---|
59 | { int m = ssx->m; |
---|
60 | int n = ssx->n; |
---|
61 | int *type = ssx->type; |
---|
62 | mpq_t *lb = ssx->lb; |
---|
63 | mpq_t *ub = ssx->ub; |
---|
64 | mpq_t *coef = ssx->coef; |
---|
65 | int *A_ptr = ssx->A_ptr; |
---|
66 | int *A_ind = ssx->A_ind; |
---|
67 | mpq_t *A_val = ssx->A_val; |
---|
68 | int *Q_col = ssx->Q_col; |
---|
69 | mpq_t *bbar = ssx->bbar; |
---|
70 | mpq_t *pi = ssx->pi; |
---|
71 | mpq_t *cbar = ssx->cbar; |
---|
72 | int *orig_type, orig_dir; |
---|
73 | mpq_t *orig_lb, *orig_ub, *orig_coef; |
---|
74 | int i, k, ret; |
---|
75 | /* save components of the original LP problem, which are changed |
---|
76 | by the routine */ |
---|
77 | orig_type = xcalloc(1+m+n, sizeof(int)); |
---|
78 | orig_lb = xcalloc(1+m+n, sizeof(mpq_t)); |
---|
79 | orig_ub = xcalloc(1+m+n, sizeof(mpq_t)); |
---|
80 | orig_coef = xcalloc(1+m+n, sizeof(mpq_t)); |
---|
81 | for (k = 1; k <= m+n; k++) |
---|
82 | { orig_type[k] = type[k]; |
---|
83 | mpq_init(orig_lb[k]); |
---|
84 | mpq_set(orig_lb[k], lb[k]); |
---|
85 | mpq_init(orig_ub[k]); |
---|
86 | mpq_set(orig_ub[k], ub[k]); |
---|
87 | } |
---|
88 | orig_dir = ssx->dir; |
---|
89 | for (k = 0; k <= m+n; k++) |
---|
90 | { mpq_init(orig_coef[k]); |
---|
91 | mpq_set(orig_coef[k], coef[k]); |
---|
92 | } |
---|
93 | /* build an artificial basic solution, which is primal feasible, |
---|
94 | and also build an auxiliary objective function to minimize the |
---|
95 | sum of infeasibilities for the original problem */ |
---|
96 | ssx->dir = SSX_MIN; |
---|
97 | for (k = 0; k <= m+n; k++) mpq_set_si(coef[k], 0, 1); |
---|
98 | mpq_set_si(bbar[0], 0, 1); |
---|
99 | for (i = 1; i <= m; i++) |
---|
100 | { int t; |
---|
101 | k = Q_col[i]; /* x[k] = xB[i] */ |
---|
102 | t = type[k]; |
---|
103 | if (t == SSX_LO || t == SSX_DB || t == SSX_FX) |
---|
104 | { /* in the original problem x[k] has lower bound */ |
---|
105 | if (mpq_cmp(bbar[i], lb[k]) < 0) |
---|
106 | { /* which is violated */ |
---|
107 | type[k] = SSX_UP; |
---|
108 | mpq_set(ub[k], lb[k]); |
---|
109 | mpq_set_si(lb[k], 0, 1); |
---|
110 | mpq_set_si(coef[k], -1, 1); |
---|
111 | mpq_add(bbar[0], bbar[0], ub[k]); |
---|
112 | mpq_sub(bbar[0], bbar[0], bbar[i]); |
---|
113 | } |
---|
114 | } |
---|
115 | if (t == SSX_UP || t == SSX_DB || t == SSX_FX) |
---|
116 | { /* in the original problem x[k] has upper bound */ |
---|
117 | if (mpq_cmp(bbar[i], ub[k]) > 0) |
---|
118 | { /* which is violated */ |
---|
119 | type[k] = SSX_LO; |
---|
120 | mpq_set(lb[k], ub[k]); |
---|
121 | mpq_set_si(ub[k], 0, 1); |
---|
122 | mpq_set_si(coef[k], +1, 1); |
---|
123 | mpq_add(bbar[0], bbar[0], bbar[i]); |
---|
124 | mpq_sub(bbar[0], bbar[0], lb[k]); |
---|
125 | } |
---|
126 | } |
---|
127 | } |
---|
128 | /* now the initial basic solution should be primal feasible due |
---|
129 | to changes of bounds of some basic variables, which turned to |
---|
130 | implicit artifical variables */ |
---|
131 | /* compute simplex multipliers and reduced costs */ |
---|
132 | ssx_eval_pi(ssx); |
---|
133 | ssx_eval_cbar(ssx); |
---|
134 | /* display initial progress of the search */ |
---|
135 | show_progress(ssx, 1); |
---|
136 | /* main loop starts here */ |
---|
137 | for (;;) |
---|
138 | { /* display current progress of the search */ |
---|
139 | #if 0 |
---|
140 | if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001) |
---|
141 | #else |
---|
142 | if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001) |
---|
143 | #endif |
---|
144 | show_progress(ssx, 1); |
---|
145 | /* we do not need to wait until all artificial variables have |
---|
146 | left the basis */ |
---|
147 | if (mpq_sgn(bbar[0]) == 0) |
---|
148 | { /* the sum of infeasibilities is zero, therefore the current |
---|
149 | solution is primal feasible for the original problem */ |
---|
150 | ret = 0; |
---|
151 | break; |
---|
152 | } |
---|
153 | /* check if the iterations limit has been exhausted */ |
---|
154 | if (ssx->it_lim == 0) |
---|
155 | { ret = 2; |
---|
156 | break; |
---|
157 | } |
---|
158 | /* check if the time limit has been exhausted */ |
---|
159 | #if 0 |
---|
160 | if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg) |
---|
161 | #else |
---|
162 | if (ssx->tm_lim >= 0.0 && |
---|
163 | ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg)) |
---|
164 | #endif |
---|
165 | { ret = 3; |
---|
166 | break; |
---|
167 | } |
---|
168 | /* choose non-basic variable xN[q] */ |
---|
169 | ssx_chuzc(ssx); |
---|
170 | /* if xN[q] cannot be chosen, the sum of infeasibilities is |
---|
171 | minimal but non-zero; therefore the original problem has no |
---|
172 | primal feasible solution */ |
---|
173 | if (ssx->q == 0) |
---|
174 | { ret = 1; |
---|
175 | break; |
---|
176 | } |
---|
177 | /* compute q-th column of the simplex table */ |
---|
178 | ssx_eval_col(ssx); |
---|
179 | /* choose basic variable xB[p] */ |
---|
180 | ssx_chuzr(ssx); |
---|
181 | /* the sum of infeasibilities cannot be negative, therefore |
---|
182 | the auxiliary lp problem cannot have unbounded solution */ |
---|
183 | xassert(ssx->p != 0); |
---|
184 | /* update values of basic variables */ |
---|
185 | ssx_update_bbar(ssx); |
---|
186 | if (ssx->p > 0) |
---|
187 | { /* compute p-th row of the inverse inv(B) */ |
---|
188 | ssx_eval_rho(ssx); |
---|
189 | /* compute p-th row of the simplex table */ |
---|
190 | ssx_eval_row(ssx); |
---|
191 | xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0); |
---|
192 | /* update simplex multipliers */ |
---|
193 | ssx_update_pi(ssx); |
---|
194 | /* update reduced costs of non-basic variables */ |
---|
195 | ssx_update_cbar(ssx); |
---|
196 | } |
---|
197 | /* xB[p] is leaving the basis; if it is implicit artificial |
---|
198 | variable, the corresponding residual vanishes; therefore |
---|
199 | bounds of this variable should be restored to the original |
---|
200 | values */ |
---|
201 | if (ssx->p > 0) |
---|
202 | { k = Q_col[ssx->p]; /* x[k] = xB[p] */ |
---|
203 | if (type[k] != orig_type[k]) |
---|
204 | { /* x[k] is implicit artificial variable */ |
---|
205 | type[k] = orig_type[k]; |
---|
206 | mpq_set(lb[k], orig_lb[k]); |
---|
207 | mpq_set(ub[k], orig_ub[k]); |
---|
208 | xassert(ssx->p_stat == SSX_NL || ssx->p_stat == SSX_NU); |
---|
209 | ssx->p_stat = (ssx->p_stat == SSX_NL ? SSX_NU : SSX_NL); |
---|
210 | if (type[k] == SSX_FX) ssx->p_stat = SSX_NS; |
---|
211 | /* nullify the objective coefficient at x[k] */ |
---|
212 | mpq_set_si(coef[k], 0, 1); |
---|
213 | /* since coef[k] has been changed, we need to compute |
---|
214 | new reduced cost of x[k], which it will have in the |
---|
215 | adjacent basis */ |
---|
216 | /* the formula d[j] = cN[j] - pi' * N[j] is used (note |
---|
217 | that the vector pi is not changed, because it depends |
---|
218 | on objective coefficients at basic variables, but in |
---|
219 | the adjacent basis, for which the vector pi has been |
---|
220 | just recomputed, x[k] is non-basic) */ |
---|
221 | if (k <= m) |
---|
222 | { /* x[k] is auxiliary variable */ |
---|
223 | mpq_neg(cbar[ssx->q], pi[k]); |
---|
224 | } |
---|
225 | else |
---|
226 | { /* x[k] is structural variable */ |
---|
227 | int ptr; |
---|
228 | mpq_t temp; |
---|
229 | mpq_init(temp); |
---|
230 | mpq_set_si(cbar[ssx->q], 0, 1); |
---|
231 | for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++) |
---|
232 | { mpq_mul(temp, pi[A_ind[ptr]], A_val[ptr]); |
---|
233 | mpq_add(cbar[ssx->q], cbar[ssx->q], temp); |
---|
234 | } |
---|
235 | mpq_clear(temp); |
---|
236 | } |
---|
237 | } |
---|
238 | } |
---|
239 | /* jump to the adjacent vertex of the polyhedron */ |
---|
240 | ssx_change_basis(ssx); |
---|
241 | /* one simplex iteration has been performed */ |
---|
242 | if (ssx->it_lim > 0) ssx->it_lim--; |
---|
243 | ssx->it_cnt++; |
---|
244 | } |
---|
245 | /* display final progress of the search */ |
---|
246 | show_progress(ssx, 1); |
---|
247 | /* restore components of the original problem, which were changed |
---|
248 | by the routine */ |
---|
249 | for (k = 1; k <= m+n; k++) |
---|
250 | { type[k] = orig_type[k]; |
---|
251 | mpq_set(lb[k], orig_lb[k]); |
---|
252 | mpq_clear(orig_lb[k]); |
---|
253 | mpq_set(ub[k], orig_ub[k]); |
---|
254 | mpq_clear(orig_ub[k]); |
---|
255 | } |
---|
256 | ssx->dir = orig_dir; |
---|
257 | for (k = 0; k <= m+n; k++) |
---|
258 | { mpq_set(coef[k], orig_coef[k]); |
---|
259 | mpq_clear(orig_coef[k]); |
---|
260 | } |
---|
261 | xfree(orig_type); |
---|
262 | xfree(orig_lb); |
---|
263 | xfree(orig_ub); |
---|
264 | xfree(orig_coef); |
---|
265 | /* return to the calling program */ |
---|
266 | return ret; |
---|
267 | } |
---|
268 | |
---|
269 | /*---------------------------------------------------------------------- |
---|
270 | // ssx_phase_II - find optimal solution. |
---|
271 | // |
---|
272 | // This routine implements phase II of the primal simplex method. |
---|
273 | // |
---|
274 | // On exit the routine returns one of the following codes: |
---|
275 | // |
---|
276 | // 0 - optimal solution found; |
---|
277 | // 1 - problem has unbounded solution; |
---|
278 | // 2 - iterations limit exceeded; |
---|
279 | // 3 - time limit exceeded. |
---|
280 | ----------------------------------------------------------------------*/ |
---|
281 | |
---|
282 | int ssx_phase_II(SSX *ssx) |
---|
283 | { int ret; |
---|
284 | /* display initial progress of the search */ |
---|
285 | show_progress(ssx, 2); |
---|
286 | /* main loop starts here */ |
---|
287 | for (;;) |
---|
288 | { /* display current progress of the search */ |
---|
289 | #if 0 |
---|
290 | if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001) |
---|
291 | #else |
---|
292 | if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001) |
---|
293 | #endif |
---|
294 | show_progress(ssx, 2); |
---|
295 | /* check if the iterations limit has been exhausted */ |
---|
296 | if (ssx->it_lim == 0) |
---|
297 | { ret = 2; |
---|
298 | break; |
---|
299 | } |
---|
300 | /* check if the time limit has been exhausted */ |
---|
301 | #if 0 |
---|
302 | if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg) |
---|
303 | #else |
---|
304 | if (ssx->tm_lim >= 0.0 && |
---|
305 | ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg)) |
---|
306 | #endif |
---|
307 | { ret = 3; |
---|
308 | break; |
---|
309 | } |
---|
310 | /* choose non-basic variable xN[q] */ |
---|
311 | ssx_chuzc(ssx); |
---|
312 | /* if xN[q] cannot be chosen, the current basic solution is |
---|
313 | dual feasible and therefore optimal */ |
---|
314 | if (ssx->q == 0) |
---|
315 | { ret = 0; |
---|
316 | break; |
---|
317 | } |
---|
318 | /* compute q-th column of the simplex table */ |
---|
319 | ssx_eval_col(ssx); |
---|
320 | /* choose basic variable xB[p] */ |
---|
321 | ssx_chuzr(ssx); |
---|
322 | /* if xB[p] cannot be chosen, the problem has no dual feasible |
---|
323 | solution (i.e. unbounded) */ |
---|
324 | if (ssx->p == 0) |
---|
325 | { ret = 1; |
---|
326 | break; |
---|
327 | } |
---|
328 | /* update values of basic variables */ |
---|
329 | ssx_update_bbar(ssx); |
---|
330 | if (ssx->p > 0) |
---|
331 | { /* compute p-th row of the inverse inv(B) */ |
---|
332 | ssx_eval_rho(ssx); |
---|
333 | /* compute p-th row of the simplex table */ |
---|
334 | ssx_eval_row(ssx); |
---|
335 | xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0); |
---|
336 | #if 0 |
---|
337 | /* update simplex multipliers */ |
---|
338 | ssx_update_pi(ssx); |
---|
339 | #endif |
---|
340 | /* update reduced costs of non-basic variables */ |
---|
341 | ssx_update_cbar(ssx); |
---|
342 | } |
---|
343 | /* jump to the adjacent vertex of the polyhedron */ |
---|
344 | ssx_change_basis(ssx); |
---|
345 | /* one simplex iteration has been performed */ |
---|
346 | if (ssx->it_lim > 0) ssx->it_lim--; |
---|
347 | ssx->it_cnt++; |
---|
348 | } |
---|
349 | /* display final progress of the search */ |
---|
350 | show_progress(ssx, 2); |
---|
351 | /* return to the calling program */ |
---|
352 | return ret; |
---|
353 | } |
---|
354 | |
---|
355 | /*---------------------------------------------------------------------- |
---|
356 | // ssx_driver - base driver to exact simplex method. |
---|
357 | // |
---|
358 | // This routine is a base driver to a version of the primal simplex |
---|
359 | // method using exact (bignum) arithmetic. |
---|
360 | // |
---|
361 | // On exit the routine returns one of the following codes: |
---|
362 | // |
---|
363 | // 0 - optimal solution found; |
---|
364 | // 1 - problem has no feasible solution; |
---|
365 | // 2 - problem has unbounded solution; |
---|
366 | // 3 - iterations limit exceeded (phase I); |
---|
367 | // 4 - iterations limit exceeded (phase II); |
---|
368 | // 5 - time limit exceeded (phase I); |
---|
369 | // 6 - time limit exceeded (phase II); |
---|
370 | // 7 - initial basis matrix is exactly singular. |
---|
371 | ----------------------------------------------------------------------*/ |
---|
372 | |
---|
373 | int ssx_driver(SSX *ssx) |
---|
374 | { int m = ssx->m; |
---|
375 | int *type = ssx->type; |
---|
376 | mpq_t *lb = ssx->lb; |
---|
377 | mpq_t *ub = ssx->ub; |
---|
378 | int *Q_col = ssx->Q_col; |
---|
379 | mpq_t *bbar = ssx->bbar; |
---|
380 | int i, k, ret; |
---|
381 | ssx->tm_beg = xtime(); |
---|
382 | /* factorize the initial basis matrix */ |
---|
383 | if (ssx_factorize(ssx)) |
---|
384 | { xprintf("Initial basis matrix is singular\n"); |
---|
385 | ret = 7; |
---|
386 | goto done; |
---|
387 | } |
---|
388 | /* compute values of basic variables */ |
---|
389 | ssx_eval_bbar(ssx); |
---|
390 | /* check if the initial basic solution is primal feasible */ |
---|
391 | for (i = 1; i <= m; i++) |
---|
392 | { int t; |
---|
393 | k = Q_col[i]; /* x[k] = xB[i] */ |
---|
394 | t = type[k]; |
---|
395 | if (t == SSX_LO || t == SSX_DB || t == SSX_FX) |
---|
396 | { /* x[k] has lower bound */ |
---|
397 | if (mpq_cmp(bbar[i], lb[k]) < 0) |
---|
398 | { /* which is violated */ |
---|
399 | break; |
---|
400 | } |
---|
401 | } |
---|
402 | if (t == SSX_UP || t == SSX_DB || t == SSX_FX) |
---|
403 | { /* x[k] has upper bound */ |
---|
404 | if (mpq_cmp(bbar[i], ub[k]) > 0) |
---|
405 | { /* which is violated */ |
---|
406 | break; |
---|
407 | } |
---|
408 | } |
---|
409 | } |
---|
410 | if (i > m) |
---|
411 | { /* no basic variable violates its bounds */ |
---|
412 | ret = 0; |
---|
413 | goto skip; |
---|
414 | } |
---|
415 | /* phase I: find primal feasible solution */ |
---|
416 | ret = ssx_phase_I(ssx); |
---|
417 | switch (ret) |
---|
418 | { case 0: |
---|
419 | ret = 0; |
---|
420 | break; |
---|
421 | case 1: |
---|
422 | xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n"); |
---|
423 | ret = 1; |
---|
424 | break; |
---|
425 | case 2: |
---|
426 | xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n"); |
---|
427 | ret = 3; |
---|
428 | break; |
---|
429 | case 3: |
---|
430 | xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n"); |
---|
431 | ret = 5; |
---|
432 | break; |
---|
433 | default: |
---|
434 | xassert(ret != ret); |
---|
435 | } |
---|
436 | /* compute values of basic variables (actually only the objective |
---|
437 | value needs to be computed) */ |
---|
438 | ssx_eval_bbar(ssx); |
---|
439 | skip: /* compute simplex multipliers */ |
---|
440 | ssx_eval_pi(ssx); |
---|
441 | /* compute reduced costs of non-basic variables */ |
---|
442 | ssx_eval_cbar(ssx); |
---|
443 | /* if phase I failed, do not start phase II */ |
---|
444 | if (ret != 0) goto done; |
---|
445 | /* phase II: find optimal solution */ |
---|
446 | ret = ssx_phase_II(ssx); |
---|
447 | switch (ret) |
---|
448 | { case 0: |
---|
449 | xprintf("OPTIMAL SOLUTION FOUND\n"); |
---|
450 | ret = 0; |
---|
451 | break; |
---|
452 | case 1: |
---|
453 | xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n"); |
---|
454 | ret = 2; |
---|
455 | break; |
---|
456 | case 2: |
---|
457 | xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n"); |
---|
458 | ret = 4; |
---|
459 | break; |
---|
460 | case 3: |
---|
461 | xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n"); |
---|
462 | ret = 6; |
---|
463 | break; |
---|
464 | default: |
---|
465 | xassert(ret != ret); |
---|
466 | } |
---|
467 | done: /* decrease the time limit by the spent amount of time */ |
---|
468 | if (ssx->tm_lim >= 0.0) |
---|
469 | #if 0 |
---|
470 | { ssx->tm_lim -= utime() - ssx->tm_beg; |
---|
471 | #else |
---|
472 | { ssx->tm_lim -= xdifftime(xtime(), ssx->tm_beg); |
---|
473 | #endif |
---|
474 | if (ssx->tm_lim < 0.0) ssx->tm_lim = 0.0; |
---|
475 | } |
---|
476 | return ret; |
---|
477 | } |
---|
478 | |
---|
479 | /* eof */ |
---|