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