lemon-project-template-glpk

comparison deps/glpk/src/glpssx02.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:d9b77317608a
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, 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
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 */