lemon-project-template-glpk

comparison deps/glpk/src/glpapi07.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:4d97f2fabf5b
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, 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 "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 */