COIN-OR::LEMON - Graph Library

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

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

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 16.4 KB
Line 
1/* glpssx02.c */
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 "glpenv.h"
26#include "glpssx.h"
27
28static void show_progress(SSX *ssx, int phase)
29{     /* this auxiliary routine displays information about progress of
30         the search */
31      int i, def = 0;
32      for (i = 1; i <= ssx->m; i++)
33         if (ssx->type[ssx->Q_col[i]] == SSX_FX) def++;
34      xprintf("%s%6d:   %s = %22.15g   (%d)\n", phase == 1 ? " " : "*",
35         ssx->it_cnt, phase == 1 ? "infsum" : "objval",
36         mpq_get_d(ssx->bbar[0]), def);
37#if 0
38      ssx->tm_lag = utime();
39#else
40      ssx->tm_lag = xtime();
41#endif
42      return;
43}
44
45/*----------------------------------------------------------------------
46// ssx_phase_I - find primal feasible solution.
47//
48// This routine implements phase I of the primal simplex method.
49//
50// On exit the routine returns one of the following codes:
51//
52// 0 - feasible solution found;
53// 1 - problem has no feasible solution;
54// 2 - iterations limit exceeded;
55// 3 - time limit exceeded.
56----------------------------------------------------------------------*/
57
58int ssx_phase_I(SSX *ssx)
59{     int m = ssx->m;
60      int n = ssx->n;
61      int *type = ssx->type;
62      mpq_t *lb = ssx->lb;
63      mpq_t *ub = ssx->ub;
64      mpq_t *coef = ssx->coef;
65      int *A_ptr = ssx->A_ptr;
66      int *A_ind = ssx->A_ind;
67      mpq_t *A_val = ssx->A_val;
68      int *Q_col = ssx->Q_col;
69      mpq_t *bbar = ssx->bbar;
70      mpq_t *pi = ssx->pi;
71      mpq_t *cbar = ssx->cbar;
72      int *orig_type, orig_dir;
73      mpq_t *orig_lb, *orig_ub, *orig_coef;
74      int i, k, ret;
75      /* save components of the original LP problem, which are changed
76         by the routine */
77      orig_type = xcalloc(1+m+n, sizeof(int));
78      orig_lb = xcalloc(1+m+n, sizeof(mpq_t));
79      orig_ub = xcalloc(1+m+n, sizeof(mpq_t));
80      orig_coef = xcalloc(1+m+n, sizeof(mpq_t));
81      for (k = 1; k <= m+n; k++)
82      {  orig_type[k] = type[k];
83         mpq_init(orig_lb[k]);
84         mpq_set(orig_lb[k], lb[k]);
85         mpq_init(orig_ub[k]);
86         mpq_set(orig_ub[k], ub[k]);
87      }
88      orig_dir = ssx->dir;
89      for (k = 0; k <= m+n; k++)
90      {  mpq_init(orig_coef[k]);
91         mpq_set(orig_coef[k], coef[k]);
92      }
93      /* build an artificial basic solution, which is primal feasible,
94         and also build an auxiliary objective function to minimize the
95         sum of infeasibilities for the original problem */
96      ssx->dir = SSX_MIN;
97      for (k = 0; k <= m+n; k++) mpq_set_si(coef[k], 0, 1);
98      mpq_set_si(bbar[0], 0, 1);
99      for (i = 1; i <= m; i++)
100      {  int t;
101         k = Q_col[i]; /* x[k] = xB[i] */
102         t = type[k];
103         if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
104         {  /* in the original problem x[k] has lower bound */
105            if (mpq_cmp(bbar[i], lb[k]) < 0)
106            {  /* which is violated */
107               type[k] = SSX_UP;
108               mpq_set(ub[k], lb[k]);
109               mpq_set_si(lb[k], 0, 1);
110               mpq_set_si(coef[k], -1, 1);
111               mpq_add(bbar[0], bbar[0], ub[k]);
112               mpq_sub(bbar[0], bbar[0], bbar[i]);
113            }
114         }
115         if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
116         {  /* in the original problem x[k] has upper bound */
117            if (mpq_cmp(bbar[i], ub[k]) > 0)
118            {  /* which is violated */
119               type[k] = SSX_LO;
120               mpq_set(lb[k], ub[k]);
121               mpq_set_si(ub[k], 0, 1);
122               mpq_set_si(coef[k], +1, 1);
123               mpq_add(bbar[0], bbar[0], bbar[i]);
124               mpq_sub(bbar[0], bbar[0], lb[k]);
125            }
126         }
127      }
128      /* now the initial basic solution should be primal feasible due
129         to changes of bounds of some basic variables, which turned to
130         implicit artifical variables */
131      /* compute simplex multipliers and reduced costs */
132      ssx_eval_pi(ssx);
133      ssx_eval_cbar(ssx);
134      /* display initial progress of the search */
135      show_progress(ssx, 1);
136      /* main loop starts here */
137      for (;;)
138      {  /* display current progress of the search */
139#if 0
140         if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
141#else
142         if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
143#endif
144            show_progress(ssx, 1);
145         /* we do not need to wait until all artificial variables have
146            left the basis */
147         if (mpq_sgn(bbar[0]) == 0)
148         {  /* the sum of infeasibilities is zero, therefore the current
149               solution is primal feasible for the original problem */
150            ret = 0;
151            break;
152         }
153         /* check if the iterations limit has been exhausted */
154         if (ssx->it_lim == 0)
155         {  ret = 2;
156            break;
157         }
158         /* check if the time limit has been exhausted */
159#if 0
160         if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
161#else
162         if (ssx->tm_lim >= 0.0 &&
163             ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
164#endif
165         {  ret = 3;
166            break;
167         }
168         /* choose non-basic variable xN[q] */
169         ssx_chuzc(ssx);
170         /* if xN[q] cannot be chosen, the sum of infeasibilities is
171            minimal but non-zero; therefore the original problem has no
172            primal feasible solution */
173         if (ssx->q == 0)
174         {  ret = 1;
175            break;
176         }
177         /* compute q-th column of the simplex table */
178         ssx_eval_col(ssx);
179         /* choose basic variable xB[p] */
180         ssx_chuzr(ssx);
181         /* the sum of infeasibilities cannot be negative, therefore
182            the auxiliary lp problem cannot have unbounded solution */
183         xassert(ssx->p != 0);
184         /* update values of basic variables */
185         ssx_update_bbar(ssx);
186         if (ssx->p > 0)
187         {  /* compute p-th row of the inverse inv(B) */
188            ssx_eval_rho(ssx);
189            /* compute p-th row of the simplex table */
190            ssx_eval_row(ssx);
191            xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
192            /* update simplex multipliers */
193            ssx_update_pi(ssx);
194            /* update reduced costs of non-basic variables */
195            ssx_update_cbar(ssx);
196         }
197         /* xB[p] is leaving the basis; if it is implicit artificial
198            variable, the corresponding residual vanishes; therefore
199            bounds of this variable should be restored to the original
200            values */
201         if (ssx->p > 0)
202         {  k = Q_col[ssx->p]; /* x[k] = xB[p] */
203            if (type[k] != orig_type[k])
204            {  /* x[k] is implicit artificial variable */
205               type[k] = orig_type[k];
206               mpq_set(lb[k], orig_lb[k]);
207               mpq_set(ub[k], orig_ub[k]);
208               xassert(ssx->p_stat == SSX_NL || ssx->p_stat == SSX_NU);
209               ssx->p_stat = (ssx->p_stat == SSX_NL ? SSX_NU : SSX_NL);
210               if (type[k] == SSX_FX) ssx->p_stat = SSX_NS;
211               /* nullify the objective coefficient at x[k] */
212               mpq_set_si(coef[k], 0, 1);
213               /* since coef[k] has been changed, we need to compute
214                  new reduced cost of x[k], which it will have in the
215                  adjacent basis */
216               /* the formula d[j] = cN[j] - pi' * N[j] is used (note
217                  that the vector pi is not changed, because it depends
218                  on objective coefficients at basic variables, but in
219                  the adjacent basis, for which the vector pi has been
220                  just recomputed, x[k] is non-basic) */
221               if (k <= m)
222               {  /* x[k] is auxiliary variable */
223                  mpq_neg(cbar[ssx->q], pi[k]);
224               }
225               else
226               {  /* x[k] is structural variable */
227                  int ptr;
228                  mpq_t temp;
229                  mpq_init(temp);
230                  mpq_set_si(cbar[ssx->q], 0, 1);
231                  for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
232                  {  mpq_mul(temp, pi[A_ind[ptr]], A_val[ptr]);
233                     mpq_add(cbar[ssx->q], cbar[ssx->q], temp);
234                  }
235                  mpq_clear(temp);
236               }
237            }
238         }
239         /* jump to the adjacent vertex of the polyhedron */
240         ssx_change_basis(ssx);
241         /* one simplex iteration has been performed */
242         if (ssx->it_lim > 0) ssx->it_lim--;
243         ssx->it_cnt++;
244      }
245      /* display final progress of the search */
246      show_progress(ssx, 1);
247      /* restore components of the original problem, which were changed
248         by the routine */
249      for (k = 1; k <= m+n; k++)
250      {  type[k] = orig_type[k];
251         mpq_set(lb[k], orig_lb[k]);
252         mpq_clear(orig_lb[k]);
253         mpq_set(ub[k], orig_ub[k]);
254         mpq_clear(orig_ub[k]);
255      }
256      ssx->dir = orig_dir;
257      for (k = 0; k <= m+n; k++)
258      {  mpq_set(coef[k], orig_coef[k]);
259         mpq_clear(orig_coef[k]);
260      }
261      xfree(orig_type);
262      xfree(orig_lb);
263      xfree(orig_ub);
264      xfree(orig_coef);
265      /* return to the calling program */
266      return ret;
267}
268
269/*----------------------------------------------------------------------
270// ssx_phase_II - find optimal solution.
271//
272// This routine implements phase II of the primal simplex method.
273//
274// On exit the routine returns one of the following codes:
275//
276// 0 - optimal solution found;
277// 1 - problem has unbounded solution;
278// 2 - iterations limit exceeded;
279// 3 - time limit exceeded.
280----------------------------------------------------------------------*/
281
282int ssx_phase_II(SSX *ssx)
283{     int ret;
284      /* display initial progress of the search */
285      show_progress(ssx, 2);
286      /* main loop starts here */
287      for (;;)
288      {  /* display current progress of the search */
289#if 0
290         if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
291#else
292         if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
293#endif
294            show_progress(ssx, 2);
295         /* check if the iterations limit has been exhausted */
296         if (ssx->it_lim == 0)
297         {  ret = 2;
298            break;
299         }
300         /* check if the time limit has been exhausted */
301#if 0
302         if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
303#else
304         if (ssx->tm_lim >= 0.0 &&
305             ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
306#endif
307         {  ret = 3;
308            break;
309         }
310         /* choose non-basic variable xN[q] */
311         ssx_chuzc(ssx);
312         /* if xN[q] cannot be chosen, the current basic solution is
313            dual feasible and therefore optimal */
314         if (ssx->q == 0)
315         {  ret = 0;
316            break;
317         }
318         /* compute q-th column of the simplex table */
319         ssx_eval_col(ssx);
320         /* choose basic variable xB[p] */
321         ssx_chuzr(ssx);
322         /* if xB[p] cannot be chosen, the problem has no dual feasible
323            solution (i.e. unbounded) */
324         if (ssx->p == 0)
325         {  ret = 1;
326            break;
327         }
328         /* update values of basic variables */
329         ssx_update_bbar(ssx);
330         if (ssx->p > 0)
331         {  /* compute p-th row of the inverse inv(B) */
332            ssx_eval_rho(ssx);
333            /* compute p-th row of the simplex table */
334            ssx_eval_row(ssx);
335            xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
336#if 0
337            /* update simplex multipliers */
338            ssx_update_pi(ssx);
339#endif
340            /* update reduced costs of non-basic variables */
341            ssx_update_cbar(ssx);
342         }
343         /* jump to the adjacent vertex of the polyhedron */
344         ssx_change_basis(ssx);
345         /* one simplex iteration has been performed */
346         if (ssx->it_lim > 0) ssx->it_lim--;
347         ssx->it_cnt++;
348      }
349      /* display final progress of the search */
350      show_progress(ssx, 2);
351      /* return to the calling program */
352      return ret;
353}
354
355/*----------------------------------------------------------------------
356// ssx_driver - base driver to exact simplex method.
357//
358// This routine is a base driver to a version of the primal simplex
359// method using exact (bignum) arithmetic.
360//
361// On exit the routine returns one of the following codes:
362//
363// 0 - optimal solution found;
364// 1 - problem has no feasible solution;
365// 2 - problem has unbounded solution;
366// 3 - iterations limit exceeded (phase I);
367// 4 - iterations limit exceeded (phase II);
368// 5 - time limit exceeded (phase I);
369// 6 - time limit exceeded (phase II);
370// 7 - initial basis matrix is exactly singular.
371----------------------------------------------------------------------*/
372
373int ssx_driver(SSX *ssx)
374{     int m = ssx->m;
375      int *type = ssx->type;
376      mpq_t *lb = ssx->lb;
377      mpq_t *ub = ssx->ub;
378      int *Q_col = ssx->Q_col;
379      mpq_t *bbar = ssx->bbar;
380      int i, k, ret;
381      ssx->tm_beg = xtime();
382      /* factorize the initial basis matrix */
383      if (ssx_factorize(ssx))
384      {  xprintf("Initial basis matrix is singular\n");
385         ret = 7;
386         goto done;
387      }
388      /* compute values of basic variables */
389      ssx_eval_bbar(ssx);
390      /* check if the initial basic solution is primal feasible */
391      for (i = 1; i <= m; i++)
392      {  int t;
393         k = Q_col[i]; /* x[k] = xB[i] */
394         t = type[k];
395         if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
396         {  /* x[k] has lower bound */
397            if (mpq_cmp(bbar[i], lb[k]) < 0)
398            {  /* which is violated */
399               break;
400            }
401         }
402         if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
403         {  /* x[k] has upper bound */
404            if (mpq_cmp(bbar[i], ub[k]) > 0)
405            {  /* which is violated */
406               break;
407            }
408         }
409      }
410      if (i > m)
411      {  /* no basic variable violates its bounds */
412         ret = 0;
413         goto skip;
414      }
415      /* phase I: find primal feasible solution */
416      ret = ssx_phase_I(ssx);
417      switch (ret)
418      {  case 0:
419            ret = 0;
420            break;
421         case 1:
422            xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
423            ret = 1;
424            break;
425         case 2:
426            xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
427            ret = 3;
428            break;
429         case 3:
430            xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
431            ret = 5;
432            break;
433         default:
434            xassert(ret != ret);
435      }
436      /* compute values of basic variables (actually only the objective
437         value needs to be computed) */
438      ssx_eval_bbar(ssx);
439skip: /* compute simplex multipliers */
440      ssx_eval_pi(ssx);
441      /* compute reduced costs of non-basic variables */
442      ssx_eval_cbar(ssx);
443      /* if phase I failed, do not start phase II */
444      if (ret != 0) goto done;
445      /* phase II: find optimal solution */
446      ret = ssx_phase_II(ssx);
447      switch (ret)
448      {  case 0:
449            xprintf("OPTIMAL SOLUTION FOUND\n");
450            ret = 0;
451            break;
452         case 1:
453            xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
454            ret = 2;
455            break;
456         case 2:
457            xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
458            ret = 4;
459            break;
460         case 3:
461            xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
462            ret = 6;
463            break;
464         default:
465            xassert(ret != ret);
466      }
467done: /* decrease the time limit by the spent amount of time */
468      if (ssx->tm_lim >= 0.0)
469#if 0
470      {  ssx->tm_lim -= utime() - ssx->tm_beg;
471#else
472      {  ssx->tm_lim -= xdifftime(xtime(), ssx->tm_beg);
473#endif
474         if (ssx->tm_lim < 0.0) ssx->tm_lim = 0.0;
475      }
476      return ret;
477}
478
479/* eof */
Note: See TracBrowser for help on using the repository browser.