PageRenderTime 37ms CodeModel.GetById 10ms app.highlight 23ms RepoModel.GetById 1ms app.codeStats 0ms

/farmR/src/glpssx02.c

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