1 /* glpapi07.c (exact simplex solver) */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 * Copyright (C) 2007-2017 Free Software Foundation, Inc.
6 * Written by Andrew Makhorin <mao@gnu.org>.
7 *
8 * GLPK is free software: you can redistribute it and/or modify it
9 * under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * GLPK is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
20 ***********************************************************************/
21
22 #include "draft.h"
23 #include "glpssx.h"
24 #include "misc.h"
25 #include "prob.h"
26
27 /***********************************************************************
28 * NAME
29 *
30 * glp_exact - solve LP problem in exact arithmetic
31 *
32 * SYNOPSIS
33 *
34 * int glp_exact(glp_prob *lp, const glp_smcp *parm);
35 *
36 * DESCRIPTION
37 *
38 * The routine glp_exact is a tentative implementation of the primal
39 * two-phase simplex method based on exact (rational) arithmetic. It is
40 * similar to the routine glp_simplex, however, for all internal
41 * computations it uses arithmetic of rational numbers, which is exact
42 * in mathematical sense, i.e. free of round-off errors unlike floating
43 * point arithmetic.
44 *
45 * Note that the routine glp_exact uses inly two control parameters
46 * passed in the structure glp_smcp, namely, it_lim and tm_lim.
47 *
48 * RETURNS
49 *
50 * 0 The LP problem instance has been successfully solved. This code
51 * does not necessarily mean that the solver has found optimal
52 * solution. It only means that the solution process was successful.
53 *
54 * GLP_EBADB
55 * Unable to start the search, because the initial basis specified
56 * in the problem object is invalid--the number of basic (auxiliary
57 * and structural) variables is not the same as the number of rows in
58 * the problem object.
59 *
60 * GLP_ESING
61 * Unable to start the search, because the basis matrix correspodning
62 * to the initial basis is exactly singular.
63 *
64 * GLP_EBOUND
65 * Unable to start the search, because some double-bounded variables
66 * have incorrect bounds.
67 *
68 * GLP_EFAIL
69 * The problem has no rows/columns.
70 *
71 * GLP_EITLIM
72 * The search was prematurely terminated, because the simplex
73 * iteration limit has been exceeded.
74 *
75 * GLP_ETMLIM
76 * The search was prematurely terminated, because the time limit has
77 * been exceeded. */
78
set_d_eps(mpq_t x,double val)79 static void set_d_eps(mpq_t x, double val)
80 { /* convert double val to rational x obtaining a more adequate
81 fraction than provided by mpq_set_d due to allowing a small
82 approximation error specified by a given relative tolerance;
83 for example, mpq_set_d would give the following
84 1/3 ~= 0.333333333333333314829616256247391... ->
85 -> 6004799503160661/18014398509481984
86 while this routine gives exactly 1/3 */
87 int s, n, j;
88 double f, p, q, eps = 1e-9;
89 mpq_t temp;
90 xassert(-DBL_MAX <= val && val <= +DBL_MAX);
91 #if 1 /* 30/VII-2008 */
92 if (val == floor(val))
93 { /* if val is integral, do not approximate */
94 mpq_set_d(x, val);
95 goto done;
96 }
97 #endif
98 if (val > 0.0)
99 s = +1;
100 else if (val < 0.0)
101 s = -1;
102 else
103 { mpq_set_si(x, 0, 1);
104 goto done;
105 }
106 f = frexp(fabs(val), &n);
107 /* |val| = f * 2^n, where 0.5 <= f < 1.0 */
108 fp2rat(f, 0.1 * eps, &p, &q);
109 /* f ~= p / q, where p and q are integers */
110 mpq_init(temp);
111 mpq_set_d(x, p);
112 mpq_set_d(temp, q);
113 mpq_div(x, x, temp);
114 mpq_set_si(temp, 1, 1);
115 for (j = 1; j <= abs(n); j++)
116 mpq_add(temp, temp, temp);
117 if (n > 0)
118 mpq_mul(x, x, temp);
119 else if (n < 0)
120 mpq_div(x, x, temp);
121 mpq_clear(temp);
122 if (s < 0) mpq_neg(x, x);
123 /* check that the desired tolerance has been attained */
124 xassert(fabs(val - mpq_get_d(x)) <= eps * (1.0 + fabs(val)));
125 done: return;
126 }
127
load_data(SSX * ssx,glp_prob * lp)128 static void load_data(SSX *ssx, glp_prob *lp)
129 { /* load LP problem data into simplex solver workspace */
130 int m = ssx->m;
131 int n = ssx->n;
132 int nnz = ssx->A_ptr[n+1]-1;
133 int j, k, type, loc, len, *ind;
134 double lb, ub, coef, *val;
135 xassert(lp->m == m);
136 xassert(lp->n == n);
137 xassert(lp->nnz == nnz);
138 /* types and bounds of rows and columns */
139 for (k = 1; k <= m+n; k++)
140 { if (k <= m)
141 { type = lp->row[k]->type;
142 lb = lp->row[k]->lb;
143 ub = lp->row[k]->ub;
144 }
145 else
146 { type = lp->col[k-m]->type;
147 lb = lp->col[k-m]->lb;
148 ub = lp->col[k-m]->ub;
149 }
150 switch (type)
151 { case GLP_FR: type = SSX_FR; break;
152 case GLP_LO: type = SSX_LO; break;
153 case GLP_UP: type = SSX_UP; break;
154 case GLP_DB: type = SSX_DB; break;
155 case GLP_FX: type = SSX_FX; break;
156 default: xassert(type != type);
157 }
158 ssx->type[k] = type;
159 set_d_eps(ssx->lb[k], lb);
160 set_d_eps(ssx->ub[k], ub);
161 }
162 /* optimization direction */
163 switch (lp->dir)
164 { case GLP_MIN: ssx->dir = SSX_MIN; break;
165 case GLP_MAX: ssx->dir = SSX_MAX; break;
166 default: xassert(lp != lp);
167 }
168 /* objective coefficients */
169 for (k = 0; k <= m+n; k++)
170 { if (k == 0)
171 coef = lp->c0;
172 else if (k <= m)
173 coef = 0.0;
174 else
175 coef = lp->col[k-m]->coef;
176 set_d_eps(ssx->coef[k], coef);
177 }
178 /* constraint coefficients */
179 ind = xcalloc(1+m, sizeof(int));
180 val = xcalloc(1+m, sizeof(double));
181 loc = 0;
182 for (j = 1; j <= n; j++)
183 { ssx->A_ptr[j] = loc+1;
184 len = glp_get_mat_col(lp, j, ind, val);
185 for (k = 1; k <= len; k++)
186 { loc++;
187 ssx->A_ind[loc] = ind[k];
188 set_d_eps(ssx->A_val[loc], val[k]);
189 }
190 }
191 xassert(loc == nnz);
192 xfree(ind);
193 xfree(val);
194 return;
195 }
196
load_basis(SSX * ssx,glp_prob * lp)197 static int load_basis(SSX *ssx, glp_prob *lp)
198 { /* load current LP basis into simplex solver workspace */
199 int m = ssx->m;
200 int n = ssx->n;
201 int *type = ssx->type;
202 int *stat = ssx->stat;
203 int *Q_row = ssx->Q_row;
204 int *Q_col = ssx->Q_col;
205 int i, j, k;
206 xassert(lp->m == m);
207 xassert(lp->n == n);
208 /* statuses of rows and columns */
209 for (k = 1; k <= m+n; k++)
210 { if (k <= m)
211 stat[k] = lp->row[k]->stat;
212 else
213 stat[k] = lp->col[k-m]->stat;
214 switch (stat[k])
215 { case GLP_BS:
216 stat[k] = SSX_BS;
217 break;
218 case GLP_NL:
219 stat[k] = SSX_NL;
220 xassert(type[k] == SSX_LO || type[k] == SSX_DB);
221 break;
222 case GLP_NU:
223 stat[k] = SSX_NU;
224 xassert(type[k] == SSX_UP || type[k] == SSX_DB);
225 break;
226 case GLP_NF:
227 stat[k] = SSX_NF;
228 xassert(type[k] == SSX_FR);
229 break;
230 case GLP_NS:
231 stat[k] = SSX_NS;
232 xassert(type[k] == SSX_FX);
233 break;
234 default:
235 xassert(stat != stat);
236 }
237 }
238 /* build permutation matix Q */
239 i = j = 0;
240 for (k = 1; k <= m+n; k++)
241 { if (stat[k] == SSX_BS)
242 { i++;
243 if (i > m) return 1;
244 Q_row[k] = i, Q_col[i] = k;
245 }
246 else
247 { j++;
248 if (j > n) return 1;
249 Q_row[k] = m+j, Q_col[m+j] = k;
250 }
251 }
252 xassert(i == m && j == n);
253 return 0;
254 }
255
glp_exact(glp_prob * lp,const glp_smcp * parm)256 int glp_exact(glp_prob *lp, const glp_smcp *parm)
257 { glp_smcp _parm;
258 SSX *ssx;
259 int m = lp->m;
260 int n = lp->n;
261 int nnz = lp->nnz;
262 int i, j, k, type, pst, dst, ret, stat;
263 double lb, ub, prim, dual, sum;
264 if (parm == NULL)
265 parm = &_parm, glp_init_smcp((glp_smcp *)parm);
266 /* check control parameters */
267 #if 1 /* 25/XI-2017 */
268 switch (parm->msg_lev)
269 { case GLP_MSG_OFF:
270 case GLP_MSG_ERR:
271 case GLP_MSG_ON:
272 case GLP_MSG_ALL:
273 case GLP_MSG_DBG:
274 break;
275 default:
276 xerror("glp_exact: msg_lev = %d; invalid parameter\n",
277 parm->msg_lev);
278 }
279 #endif
280 if (parm->it_lim < 0)
281 xerror("glp_exact: it_lim = %d; invalid parameter\n",
282 parm->it_lim);
283 if (parm->tm_lim < 0)
284 xerror("glp_exact: tm_lim = %d; invalid parameter\n",
285 parm->tm_lim);
286 /* the problem must have at least one row and one column */
287 if (!(m > 0 && n > 0))
288 #if 0 /* 25/XI-2017 */
289 { xprintf("glp_exact: problem has no rows/columns\n");
290 #else
291 { if (parm->msg_lev >= GLP_MSG_ERR)
292 xprintf("glp_exact: problem has no rows/columns\n");
293 #endif
294 return GLP_EFAIL;
295 }
296 #if 1
297 /* basic solution is currently undefined */
298 lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
299 lp->obj_val = 0.0;
300 lp->some = 0;
301 #endif
302 /* check that all double-bounded variables have correct bounds */
303 for (k = 1; k <= m+n; k++)
304 { if (k <= m)
305 { type = lp->row[k]->type;
306 lb = lp->row[k]->lb;
307 ub = lp->row[k]->ub;
308 }
309 else
310 { type = lp->col[k-m]->type;
311 lb = lp->col[k-m]->lb;
312 ub = lp->col[k-m]->ub;
313 }
314 if (type == GLP_DB && lb >= ub)
315 #if 0 /* 25/XI-2017 */
316 { xprintf("glp_exact: %s %d has invalid bounds\n",
317 k <= m ? "row" : "column", k <= m ? k : k-m);
318 #else
319 { if (parm->msg_lev >= GLP_MSG_ERR)
320 xprintf("glp_exact: %s %d has invalid bounds\n",
321 k <= m ? "row" : "column", k <= m ? k : k-m);
322 #endif
323 return GLP_EBOUND;
324 }
325 }
326 /* create the simplex solver workspace */
327 #if 1 /* 25/XI-2017 */
328 if (parm->msg_lev >= GLP_MSG_ALL)
329 {
330 #endif
331 xprintf("glp_exact: %d rows, %d columns, %d non-zeros\n",
332 m, n, nnz);
333 #ifdef HAVE_GMP
334 xprintf("GNU MP bignum library is being used\n");
335 #else
336 xprintf("GLPK bignum module is being used\n");
337 xprintf("(Consider installing GNU MP to attain a much better perf"
338 "ormance.)\n");
339 #endif
340 #if 1 /* 25/XI-2017 */
341 }
342 #endif
343 ssx = ssx_create(m, n, nnz);
344 /* load LP problem data into the workspace */
345 load_data(ssx, lp);
346 /* load current LP basis into the workspace */
347 if (load_basis(ssx, lp))
348 #if 0 /* 25/XI-2017 */
349 { xprintf("glp_exact: initial LP basis is invalid\n");
350 #else
351 { if (parm->msg_lev >= GLP_MSG_ERR)
352 xprintf("glp_exact: initial LP basis is invalid\n");
353 #endif
354 ret = GLP_EBADB;
355 goto done;
356 }
357 #if 0
358 /* inherit some control parameters from the LP object */
359 ssx->it_lim = lpx_get_int_parm(lp, LPX_K_ITLIM);
360 ssx->it_cnt = lpx_get_int_parm(lp, LPX_K_ITCNT);
361 ssx->tm_lim = lpx_get_real_parm(lp, LPX_K_TMLIM);
362 #else
363 #if 1 /* 25/XI-2017 */
364 ssx->msg_lev = parm->msg_lev;
365 #endif
366 ssx->it_lim = parm->it_lim;
367 ssx->it_cnt = lp->it_cnt;
368 ssx->tm_lim = (double)parm->tm_lim / 1000.0;
369 #endif
370 ssx->out_frq = 5.0;
371 ssx->tm_beg = xtime();
372 #if 0 /* 10/VI-2013 */
373 ssx->tm_lag = xlset(0);
374 #else
375 ssx->tm_lag = 0.0;
376 #endif
377 /* solve LP */
378 ret = ssx_driver(ssx);
379 #if 0
380 /* copy back some statistics to the LP object */
381 lpx_set_int_parm(lp, LPX_K_ITLIM, ssx->it_lim);
382 lpx_set_int_parm(lp, LPX_K_ITCNT, ssx->it_cnt);
383 lpx_set_real_parm(lp, LPX_K_TMLIM, ssx->tm_lim);
384 #else
385 lp->it_cnt = ssx->it_cnt;
386 #endif
387 /* analyze the return code */
388 switch (ret)
389 { case 0:
390 /* optimal solution found */
391 ret = 0;
392 pst = dst = GLP_FEAS;
393 break;
394 case 1:
395 /* problem has no feasible solution */
396 ret = 0;
397 pst = GLP_NOFEAS, dst = GLP_INFEAS;
398 break;
399 case 2:
400 /* problem has unbounded solution */
401 ret = 0;
402 pst = GLP_FEAS, dst = GLP_NOFEAS;
403 #if 1
404 xassert(1 <= ssx->q && ssx->q <= n);
405 lp->some = ssx->Q_col[m + ssx->q];
406 xassert(1 <= lp->some && lp->some <= m+n);
407 #endif
408 break;
409 case 3:
410 /* iteration limit exceeded (phase I) */
411 ret = GLP_EITLIM;
412 pst = dst = GLP_INFEAS;
413 break;
414 case 4:
415 /* iteration limit exceeded (phase II) */
416 ret = GLP_EITLIM;
417 pst = GLP_FEAS, dst = GLP_INFEAS;
418 break;
419 case 5:
420 /* time limit exceeded (phase I) */
421 ret = GLP_ETMLIM;
422 pst = dst = GLP_INFEAS;
423 break;
424 case 6:
425 /* time limit exceeded (phase II) */
426 ret = GLP_ETMLIM;
427 pst = GLP_FEAS, dst = GLP_INFEAS;
428 break;
429 case 7:
430 /* initial basis matrix is singular */
431 ret = GLP_ESING;
432 goto done;
433 default:
434 xassert(ret != ret);
435 }
436 /* store final basic solution components into LP object */
437 lp->pbs_stat = pst;
438 lp->dbs_stat = dst;
439 sum = lp->c0;
440 for (k = 1; k <= m+n; k++)
441 { if (ssx->stat[k] == SSX_BS)
442 { i = ssx->Q_row[k]; /* x[k] = xB[i] */
443 xassert(1 <= i && i <= m);
444 stat = GLP_BS;
445 prim = mpq_get_d(ssx->bbar[i]);
446 dual = 0.0;
447 }
448 else
449 { j = ssx->Q_row[k] - m; /* x[k] = xN[j] */
450 xassert(1 <= j && j <= n);
451 switch (ssx->stat[k])
452 { case SSX_NF:
453 stat = GLP_NF;
454 prim = 0.0;
455 break;
456 case SSX_NL:
457 stat = GLP_NL;
458 prim = mpq_get_d(ssx->lb[k]);
459 break;
460 case SSX_NU:
461 stat = GLP_NU;
462 prim = mpq_get_d(ssx->ub[k]);
463 break;
464 case SSX_NS:
465 stat = GLP_NS;
466 prim = mpq_get_d(ssx->lb[k]);
467 break;
468 default:
469 xassert(ssx != ssx);
470 }
471 dual = mpq_get_d(ssx->cbar[j]);
472 }
473 if (k <= m)
474 { glp_set_row_stat(lp, k, stat);
475 lp->row[k]->prim = prim;
476 lp->row[k]->dual = dual;
477 }
478 else
479 { glp_set_col_stat(lp, k-m, stat);
480 lp->col[k-m]->prim = prim;
481 lp->col[k-m]->dual = dual;
482 sum += lp->col[k-m]->coef * prim;
483 }
484 }
485 lp->obj_val = sum;
486 done: /* delete the simplex solver workspace */
487 ssx_delete(ssx);
488 #if 1 /* 23/XI-2015 */
489 xassert(gmp_pool_count() == 0);
490 gmp_free_mem();
491 #endif
492 /* return to the application program */
493 return ret;
494 }
495
496 /* eof */
497