1 | /* glpapi07.c (exact simplex solver) */ |
---|
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 "glpapi.h" |
---|
26 | #include "glpssx.h" |
---|
27 | |
---|
28 | /*********************************************************************** |
---|
29 | * NAME |
---|
30 | * |
---|
31 | * glp_exact - solve LP problem in exact arithmetic |
---|
32 | * |
---|
33 | * SYNOPSIS |
---|
34 | * |
---|
35 | * int glp_exact(glp_prob *lp, const glp_smcp *parm); |
---|
36 | * |
---|
37 | * DESCRIPTION |
---|
38 | * |
---|
39 | * The routine glp_exact is a tentative implementation of the primal |
---|
40 | * two-phase simplex method based on exact (rational) arithmetic. It is |
---|
41 | * similar to the routine glp_simplex, however, for all internal |
---|
42 | * computations it uses arithmetic of rational numbers, which is exact |
---|
43 | * in mathematical sense, i.e. free of round-off errors unlike floating |
---|
44 | * point arithmetic. |
---|
45 | * |
---|
46 | * Note that the routine glp_exact uses inly two control parameters |
---|
47 | * passed in the structure glp_smcp, namely, it_lim and tm_lim. |
---|
48 | * |
---|
49 | * RETURNS |
---|
50 | * |
---|
51 | * 0 The LP problem instance has been successfully solved. This code |
---|
52 | * does not necessarily mean that the solver has found optimal |
---|
53 | * solution. It only means that the solution process was successful. |
---|
54 | * |
---|
55 | * GLP_EBADB |
---|
56 | * Unable to start the search, because the initial basis specified |
---|
57 | * in the problem object is invalid--the number of basic (auxiliary |
---|
58 | * and structural) variables is not the same as the number of rows in |
---|
59 | * the problem object. |
---|
60 | * |
---|
61 | * GLP_ESING |
---|
62 | * Unable to start the search, because the basis matrix correspodning |
---|
63 | * to the initial basis is exactly singular. |
---|
64 | * |
---|
65 | * GLP_EBOUND |
---|
66 | * Unable to start the search, because some double-bounded variables |
---|
67 | * have incorrect bounds. |
---|
68 | * |
---|
69 | * GLP_EFAIL |
---|
70 | * The problem has no rows/columns. |
---|
71 | * |
---|
72 | * GLP_EITLIM |
---|
73 | * The search was prematurely terminated, because the simplex |
---|
74 | * iteration limit has been exceeded. |
---|
75 | * |
---|
76 | * GLP_ETMLIM |
---|
77 | * The search was prematurely terminated, because the time limit has |
---|
78 | * been exceeded. */ |
---|
79 | |
---|
80 | static void set_d_eps(mpq_t x, double val) |
---|
81 | { /* convert double val to rational x obtaining a more adequate |
---|
82 | fraction than provided by mpq_set_d due to allowing a small |
---|
83 | approximation error specified by a given relative tolerance; |
---|
84 | for example, mpq_set_d would give the following |
---|
85 | 1/3 ~= 0.333333333333333314829616256247391... -> |
---|
86 | -> 6004799503160661/18014398509481984 |
---|
87 | while this routine gives exactly 1/3 */ |
---|
88 | int s, n, j; |
---|
89 | double f, p, q, eps = 1e-9; |
---|
90 | mpq_t temp; |
---|
91 | xassert(-DBL_MAX <= val && val <= +DBL_MAX); |
---|
92 | #if 1 /* 30/VII-2008 */ |
---|
93 | if (val == floor(val)) |
---|
94 | { /* if val is integral, do not approximate */ |
---|
95 | mpq_set_d(x, val); |
---|
96 | goto done; |
---|
97 | } |
---|
98 | #endif |
---|
99 | if (val > 0.0) |
---|
100 | s = +1; |
---|
101 | else if (val < 0.0) |
---|
102 | s = -1; |
---|
103 | else |
---|
104 | { mpq_set_si(x, 0, 1); |
---|
105 | goto done; |
---|
106 | } |
---|
107 | f = frexp(fabs(val), &n); |
---|
108 | /* |val| = f * 2^n, where 0.5 <= f < 1.0 */ |
---|
109 | fp2rat(f, 0.1 * eps, &p, &q); |
---|
110 | /* f ~= p / q, where p and q are integers */ |
---|
111 | mpq_init(temp); |
---|
112 | mpq_set_d(x, p); |
---|
113 | mpq_set_d(temp, q); |
---|
114 | mpq_div(x, x, temp); |
---|
115 | mpq_set_si(temp, 1, 1); |
---|
116 | for (j = 1; j <= abs(n); j++) |
---|
117 | mpq_add(temp, temp, temp); |
---|
118 | if (n > 0) |
---|
119 | mpq_mul(x, x, temp); |
---|
120 | else if (n < 0) |
---|
121 | mpq_div(x, x, temp); |
---|
122 | mpq_clear(temp); |
---|
123 | if (s < 0) mpq_neg(x, x); |
---|
124 | /* check that the desired tolerance has been attained */ |
---|
125 | xassert(fabs(val - mpq_get_d(x)) <= eps * (1.0 + fabs(val))); |
---|
126 | done: return; |
---|
127 | } |
---|
128 | |
---|
129 | static void load_data(SSX *ssx, LPX *lp) |
---|
130 | { /* load LP problem data into simplex solver workspace */ |
---|
131 | int m = ssx->m; |
---|
132 | int n = ssx->n; |
---|
133 | int nnz = ssx->A_ptr[n+1]-1; |
---|
134 | int j, k, type, loc, len, *ind; |
---|
135 | double lb, ub, coef, *val; |
---|
136 | xassert(lpx_get_num_rows(lp) == m); |
---|
137 | xassert(lpx_get_num_cols(lp) == n); |
---|
138 | xassert(lpx_get_num_nz(lp) == nnz); |
---|
139 | /* types and bounds of rows and columns */ |
---|
140 | for (k = 1; k <= m+n; k++) |
---|
141 | { if (k <= m) |
---|
142 | { type = lpx_get_row_type(lp, k); |
---|
143 | lb = lpx_get_row_lb(lp, k); |
---|
144 | ub = lpx_get_row_ub(lp, k); |
---|
145 | } |
---|
146 | else |
---|
147 | { type = lpx_get_col_type(lp, k-m); |
---|
148 | lb = lpx_get_col_lb(lp, k-m); |
---|
149 | ub = lpx_get_col_ub(lp, k-m); |
---|
150 | } |
---|
151 | switch (type) |
---|
152 | { case LPX_FR: type = SSX_FR; break; |
---|
153 | case LPX_LO: type = SSX_LO; break; |
---|
154 | case LPX_UP: type = SSX_UP; break; |
---|
155 | case LPX_DB: type = SSX_DB; break; |
---|
156 | case LPX_FX: type = SSX_FX; break; |
---|
157 | default: xassert(type != type); |
---|
158 | } |
---|
159 | ssx->type[k] = type; |
---|
160 | set_d_eps(ssx->lb[k], lb); |
---|
161 | set_d_eps(ssx->ub[k], ub); |
---|
162 | } |
---|
163 | /* optimization direction */ |
---|
164 | switch (lpx_get_obj_dir(lp)) |
---|
165 | { case LPX_MIN: ssx->dir = SSX_MIN; break; |
---|
166 | case LPX_MAX: ssx->dir = SSX_MAX; break; |
---|
167 | default: xassert(lp != lp); |
---|
168 | } |
---|
169 | /* objective coefficients */ |
---|
170 | for (k = 0; k <= m+n; k++) |
---|
171 | { if (k == 0) |
---|
172 | coef = lpx_get_obj_coef(lp, 0); |
---|
173 | else if (k <= m) |
---|
174 | coef = 0.0; |
---|
175 | else |
---|
176 | coef = lpx_get_obj_coef(lp, k-m); |
---|
177 | set_d_eps(ssx->coef[k], coef); |
---|
178 | } |
---|
179 | /* constraint coefficients */ |
---|
180 | ind = xcalloc(1+m, sizeof(int)); |
---|
181 | val = xcalloc(1+m, sizeof(double)); |
---|
182 | loc = 0; |
---|
183 | for (j = 1; j <= n; j++) |
---|
184 | { ssx->A_ptr[j] = loc+1; |
---|
185 | len = lpx_get_mat_col(lp, j, ind, val); |
---|
186 | for (k = 1; k <= len; k++) |
---|
187 | { loc++; |
---|
188 | ssx->A_ind[loc] = ind[k]; |
---|
189 | set_d_eps(ssx->A_val[loc], val[k]); |
---|
190 | } |
---|
191 | } |
---|
192 | xassert(loc == nnz); |
---|
193 | xfree(ind); |
---|
194 | xfree(val); |
---|
195 | return; |
---|
196 | } |
---|
197 | |
---|
198 | static int load_basis(SSX *ssx, LPX *lp) |
---|
199 | { /* load current LP basis into simplex solver workspace */ |
---|
200 | int m = ssx->m; |
---|
201 | int n = ssx->n; |
---|
202 | int *type = ssx->type; |
---|
203 | int *stat = ssx->stat; |
---|
204 | int *Q_row = ssx->Q_row; |
---|
205 | int *Q_col = ssx->Q_col; |
---|
206 | int i, j, k; |
---|
207 | xassert(lpx_get_num_rows(lp) == m); |
---|
208 | xassert(lpx_get_num_cols(lp) == n); |
---|
209 | /* statuses of rows and columns */ |
---|
210 | for (k = 1; k <= m+n; k++) |
---|
211 | { if (k <= m) |
---|
212 | stat[k] = lpx_get_row_stat(lp, k); |
---|
213 | else |
---|
214 | stat[k] = lpx_get_col_stat(lp, k-m); |
---|
215 | switch (stat[k]) |
---|
216 | { case LPX_BS: |
---|
217 | stat[k] = SSX_BS; |
---|
218 | break; |
---|
219 | case LPX_NL: |
---|
220 | stat[k] = SSX_NL; |
---|
221 | xassert(type[k] == SSX_LO || type[k] == SSX_DB); |
---|
222 | break; |
---|
223 | case LPX_NU: |
---|
224 | stat[k] = SSX_NU; |
---|
225 | xassert(type[k] == SSX_UP || type[k] == SSX_DB); |
---|
226 | break; |
---|
227 | case LPX_NF: |
---|
228 | stat[k] = SSX_NF; |
---|
229 | xassert(type[k] == SSX_FR); |
---|
230 | break; |
---|
231 | case LPX_NS: |
---|
232 | stat[k] = SSX_NS; |
---|
233 | xassert(type[k] == SSX_FX); |
---|
234 | break; |
---|
235 | default: |
---|
236 | xassert(stat != stat); |
---|
237 | } |
---|
238 | } |
---|
239 | /* build permutation matix Q */ |
---|
240 | i = j = 0; |
---|
241 | for (k = 1; k <= m+n; k++) |
---|
242 | { if (stat[k] == SSX_BS) |
---|
243 | { i++; |
---|
244 | if (i > m) return 1; |
---|
245 | Q_row[k] = i, Q_col[i] = k; |
---|
246 | } |
---|
247 | else |
---|
248 | { j++; |
---|
249 | if (j > n) return 1; |
---|
250 | Q_row[k] = m+j, Q_col[m+j] = k; |
---|
251 | } |
---|
252 | } |
---|
253 | xassert(i == m && j == n); |
---|
254 | return 0; |
---|
255 | } |
---|
256 | |
---|
257 | int glp_exact(glp_prob *lp, const glp_smcp *parm) |
---|
258 | { glp_smcp _parm; |
---|
259 | SSX *ssx; |
---|
260 | int m = lpx_get_num_rows(lp); |
---|
261 | int n = lpx_get_num_cols(lp); |
---|
262 | int nnz = lpx_get_num_nz(lp); |
---|
263 | int i, j, k, type, pst, dst, ret, *stat; |
---|
264 | double lb, ub, *prim, *dual, sum; |
---|
265 | if (parm == NULL) |
---|
266 | parm = &_parm, glp_init_smcp((glp_smcp *)parm); |
---|
267 | /* check control parameters */ |
---|
268 | if (parm->it_lim < 0) |
---|
269 | xerror("glp_exact: it_lim = %d; invalid parameter\n", |
---|
270 | parm->it_lim); |
---|
271 | if (parm->tm_lim < 0) |
---|
272 | xerror("glp_exact: tm_lim = %d; invalid parameter\n", |
---|
273 | parm->tm_lim); |
---|
274 | /* the problem must have at least one row and one column */ |
---|
275 | if (!(m > 0 && n > 0)) |
---|
276 | { xprintf("glp_exact: problem has no rows/columns\n"); |
---|
277 | return GLP_EFAIL; |
---|
278 | } |
---|
279 | #if 1 |
---|
280 | /* basic solution is currently undefined */ |
---|
281 | lp->pbs_stat = lp->dbs_stat = GLP_UNDEF; |
---|
282 | lp->obj_val = 0.0; |
---|
283 | lp->some = 0; |
---|
284 | #endif |
---|
285 | /* check that all double-bounded variables have correct bounds */ |
---|
286 | for (k = 1; k <= m+n; k++) |
---|
287 | { if (k <= m) |
---|
288 | { type = lpx_get_row_type(lp, k); |
---|
289 | lb = lpx_get_row_lb(lp, k); |
---|
290 | ub = lpx_get_row_ub(lp, k); |
---|
291 | } |
---|
292 | else |
---|
293 | { type = lpx_get_col_type(lp, k-m); |
---|
294 | lb = lpx_get_col_lb(lp, k-m); |
---|
295 | ub = lpx_get_col_ub(lp, k-m); |
---|
296 | } |
---|
297 | if (type == LPX_DB && lb >= ub) |
---|
298 | { xprintf("glp_exact: %s %d has invalid bounds\n", |
---|
299 | k <= m ? "row" : "column", k <= m ? k : k-m); |
---|
300 | return GLP_EBOUND; |
---|
301 | } |
---|
302 | } |
---|
303 | /* create the simplex solver workspace */ |
---|
304 | xprintf("glp_exact: %d rows, %d columns, %d non-zeros\n", |
---|
305 | m, n, nnz); |
---|
306 | #ifdef HAVE_GMP |
---|
307 | xprintf("GNU MP bignum library is being used\n"); |
---|
308 | #else |
---|
309 | xprintf("GLPK bignum module is being used\n"); |
---|
310 | xprintf("(Consider installing GNU MP to attain a much better perf" |
---|
311 | "ormance.)\n"); |
---|
312 | #endif |
---|
313 | ssx = ssx_create(m, n, nnz); |
---|
314 | /* load LP problem data into the workspace */ |
---|
315 | load_data(ssx, lp); |
---|
316 | /* load current LP basis into the workspace */ |
---|
317 | if (load_basis(ssx, lp)) |
---|
318 | { xprintf("glp_exact: initial LP basis is invalid\n"); |
---|
319 | ret = GLP_EBADB; |
---|
320 | goto done; |
---|
321 | } |
---|
322 | /* inherit some control parameters from the LP object */ |
---|
323 | #if 0 |
---|
324 | ssx->it_lim = lpx_get_int_parm(lp, LPX_K_ITLIM); |
---|
325 | ssx->it_cnt = lpx_get_int_parm(lp, LPX_K_ITCNT); |
---|
326 | ssx->tm_lim = lpx_get_real_parm(lp, LPX_K_TMLIM); |
---|
327 | #else |
---|
328 | ssx->it_lim = parm->it_lim; |
---|
329 | ssx->it_cnt = lp->it_cnt; |
---|
330 | ssx->tm_lim = (double)parm->tm_lim / 1000.0; |
---|
331 | #endif |
---|
332 | ssx->out_frq = 5.0; |
---|
333 | ssx->tm_beg = xtime(); |
---|
334 | ssx->tm_lag = xlset(0); |
---|
335 | /* solve LP */ |
---|
336 | ret = ssx_driver(ssx); |
---|
337 | /* copy back some statistics to the LP object */ |
---|
338 | #if 0 |
---|
339 | lpx_set_int_parm(lp, LPX_K_ITLIM, ssx->it_lim); |
---|
340 | lpx_set_int_parm(lp, LPX_K_ITCNT, ssx->it_cnt); |
---|
341 | lpx_set_real_parm(lp, LPX_K_TMLIM, ssx->tm_lim); |
---|
342 | #else |
---|
343 | lp->it_cnt = ssx->it_cnt; |
---|
344 | #endif |
---|
345 | /* analyze the return code */ |
---|
346 | switch (ret) |
---|
347 | { case 0: |
---|
348 | /* optimal solution found */ |
---|
349 | ret = 0; |
---|
350 | pst = LPX_P_FEAS, dst = LPX_D_FEAS; |
---|
351 | break; |
---|
352 | case 1: |
---|
353 | /* problem has no feasible solution */ |
---|
354 | ret = 0; |
---|
355 | pst = LPX_P_NOFEAS, dst = LPX_D_INFEAS; |
---|
356 | break; |
---|
357 | case 2: |
---|
358 | /* problem has unbounded solution */ |
---|
359 | ret = 0; |
---|
360 | pst = LPX_P_FEAS, dst = LPX_D_NOFEAS; |
---|
361 | #if 1 |
---|
362 | xassert(1 <= ssx->q && ssx->q <= n); |
---|
363 | lp->some = ssx->Q_col[m + ssx->q]; |
---|
364 | xassert(1 <= lp->some && lp->some <= m+n); |
---|
365 | #endif |
---|
366 | break; |
---|
367 | case 3: |
---|
368 | /* iteration limit exceeded (phase I) */ |
---|
369 | ret = GLP_EITLIM; |
---|
370 | pst = LPX_P_INFEAS, dst = LPX_D_INFEAS; |
---|
371 | break; |
---|
372 | case 4: |
---|
373 | /* iteration limit exceeded (phase II) */ |
---|
374 | ret = GLP_EITLIM; |
---|
375 | pst = LPX_P_FEAS, dst = LPX_D_INFEAS; |
---|
376 | break; |
---|
377 | case 5: |
---|
378 | /* time limit exceeded (phase I) */ |
---|
379 | ret = GLP_ETMLIM; |
---|
380 | pst = LPX_P_INFEAS, dst = LPX_D_INFEAS; |
---|
381 | break; |
---|
382 | case 6: |
---|
383 | /* time limit exceeded (phase II) */ |
---|
384 | ret = GLP_ETMLIM; |
---|
385 | pst = LPX_P_FEAS, dst = LPX_D_INFEAS; |
---|
386 | break; |
---|
387 | case 7: |
---|
388 | /* initial basis matrix is singular */ |
---|
389 | ret = GLP_ESING; |
---|
390 | goto done; |
---|
391 | default: |
---|
392 | xassert(ret != ret); |
---|
393 | } |
---|
394 | /* obtain final basic solution components */ |
---|
395 | stat = xcalloc(1+m+n, sizeof(int)); |
---|
396 | prim = xcalloc(1+m+n, sizeof(double)); |
---|
397 | dual = xcalloc(1+m+n, sizeof(double)); |
---|
398 | for (k = 1; k <= m+n; k++) |
---|
399 | { if (ssx->stat[k] == SSX_BS) |
---|
400 | { i = ssx->Q_row[k]; /* x[k] = xB[i] */ |
---|
401 | xassert(1 <= i && i <= m); |
---|
402 | stat[k] = LPX_BS; |
---|
403 | prim[k] = mpq_get_d(ssx->bbar[i]); |
---|
404 | dual[k] = 0.0; |
---|
405 | } |
---|
406 | else |
---|
407 | { j = ssx->Q_row[k] - m; /* x[k] = xN[j] */ |
---|
408 | xassert(1 <= j && j <= n); |
---|
409 | switch (ssx->stat[k]) |
---|
410 | { case SSX_NF: |
---|
411 | stat[k] = LPX_NF; |
---|
412 | prim[k] = 0.0; |
---|
413 | break; |
---|
414 | case SSX_NL: |
---|
415 | stat[k] = LPX_NL; |
---|
416 | prim[k] = mpq_get_d(ssx->lb[k]); |
---|
417 | break; |
---|
418 | case SSX_NU: |
---|
419 | stat[k] = LPX_NU; |
---|
420 | prim[k] = mpq_get_d(ssx->ub[k]); |
---|
421 | break; |
---|
422 | case SSX_NS: |
---|
423 | stat[k] = LPX_NS; |
---|
424 | prim[k] = mpq_get_d(ssx->lb[k]); |
---|
425 | break; |
---|
426 | default: |
---|
427 | xassert(ssx != ssx); |
---|
428 | } |
---|
429 | dual[k] = mpq_get_d(ssx->cbar[j]); |
---|
430 | } |
---|
431 | } |
---|
432 | /* and store them into the LP object */ |
---|
433 | pst = pst - LPX_P_UNDEF + GLP_UNDEF; |
---|
434 | dst = dst - LPX_D_UNDEF + GLP_UNDEF; |
---|
435 | for (k = 1; k <= m+n; k++) |
---|
436 | stat[k] = stat[k] - LPX_BS + GLP_BS; |
---|
437 | sum = lpx_get_obj_coef(lp, 0); |
---|
438 | for (j = 1; j <= n; j++) |
---|
439 | sum += lpx_get_obj_coef(lp, j) * prim[m+j]; |
---|
440 | lpx_put_solution(lp, 1, &pst, &dst, &sum, |
---|
441 | &stat[0], &prim[0], &dual[0], &stat[m], &prim[m], &dual[m]); |
---|
442 | xfree(stat); |
---|
443 | xfree(prim); |
---|
444 | xfree(dual); |
---|
445 | done: /* delete the simplex solver workspace */ |
---|
446 | ssx_delete(ssx); |
---|
447 | /* return to the application program */ |
---|
448 | return ret; |
---|
449 | } |
---|
450 | |
---|
451 | /* eof */ |
---|