COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpapi07.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 14 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 14.5 KB
Line 
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
80static 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)));
126done: return;
127}
128
129static 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
198static 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
257int 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);
445done: /* delete the simplex solver workspace */
446      ssx_delete(ssx);
447      /* return to the application program */
448      return ret;
449}
450
451/* eof */
Note: See TracBrowser for help on using the repository browser.