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