glpssx02.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

/* glpssx02.c */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[e-mail]>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glplib.h"
#include "glpssx.h"

static void show_progress(SSX *ssx, int phase)
{     /* this auxiliary routine displays information about progress of
         the search */
      int i, def = 0;
      for (i = 1; i <= ssx->m; i++)
         if (ssx->type[ssx->Q_col[i]] == SSX_FX) def++;
      xprintf("%s%6d:   %s = %22.15g   (%d)\n", phase == 1 ? " " : "*",
         ssx->it_cnt, phase == 1 ? "infsum" : "objval",
         mpq_get_d(ssx->bbar[0]), def);
#if 0
      ssx->tm_lag = utime();
#else
      ssx->tm_lag = xtime();
#endif
      return;
}

/*----------------------------------------------------------------------
// ssx_phase_I - find primal feasible solution.
//
// This routine implements phase I of the primal simplex method.
//
// On exit the routine returns one of the following codes:
//
// 0 - feasible solution found;
// 1 - problem has no feasible solution;
// 2 - iterations limit exceeded;
// 3 - time limit exceeded.
----------------------------------------------------------------------*/

int ssx_phase_I(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      int *type = ssx->type;
      mpq_t *lb = ssx->lb;
      mpq_t *ub = ssx->ub;
      mpq_t *coef = ssx->coef;
      int *A_ptr = ssx->A_ptr;
      int *A_ind = ssx->A_ind;
      mpq_t *A_val = ssx->A_val;
      int *Q_col = ssx->Q_col;
      mpq_t *bbar = ssx->bbar;
      mpq_t *pi = ssx->pi;
      mpq_t *cbar = ssx->cbar;
      int *orig_type, orig_dir;
      mpq_t *orig_lb, *orig_ub, *orig_coef;
      int i, k, ret;
      /* save components of the original LP problem, which are changed
         by the routine */
      orig_type = xcalloc(1+m+n, sizeof(int));
      orig_lb = xcalloc(1+m+n, sizeof(mpq_t));
      orig_ub = xcalloc(1+m+n, sizeof(mpq_t));
      orig_coef = xcalloc(1+m+n, sizeof(mpq_t));
      for (k = 1; k <= m+n; k++)
      {  orig_type[k] = type[k];
         mpq_init(orig_lb[k]);
         mpq_set(orig_lb[k], lb[k]);
         mpq_init(orig_ub[k]);
         mpq_set(orig_ub[k], ub[k]);
      }
      orig_dir = ssx->dir;
      for (k = 0; k <= m+n; k++)
      {  mpq_init(orig_coef[k]);
         mpq_set(orig_coef[k], coef[k]);
      }
      /* build an artificial basic solution, which is primal feasible,
         and also build an auxiliary objective function to minimize the
         sum of infeasibilities for the original problem */
      ssx->dir = SSX_MIN;
      for (k = 0; k <= m+n; k++) mpq_set_si(coef[k], 0, 1);
      mpq_set_si(bbar[0], 0, 1);
      for (i = 1; i <= m; i++)
      {  int t;
         k = Q_col[i]; /* x[k] = xB[i] */
         t = type[k];
         if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
         {  /* in the original problem x[k] has lower bound */
            if (mpq_cmp(bbar[i], lb[k]) < 0)
            {  /* which is violated */
               type[k] = SSX_UP;
               mpq_set(ub[k], lb[k]);
               mpq_set_si(lb[k], 0, 1);
               mpq_set_si(coef[k], -1, 1);
               mpq_add(bbar[0], bbar[0], ub[k]);
               mpq_sub(bbar[0], bbar[0], bbar[i]);
            }
         }
         if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
         {  /* in the original problem x[k] has upper bound */
            if (mpq_cmp(bbar[i], ub[k]) > 0)
            {  /* which is violated */
               type[k] = SSX_LO;
               mpq_set(lb[k], ub[k]);
               mpq_set_si(ub[k], 0, 1);
               mpq_set_si(coef[k], +1, 1);
               mpq_add(bbar[0], bbar[0], bbar[i]);
               mpq_sub(bbar[0], bbar[0], lb[k]);
            }
         }
      }
      /* now the initial basic solution should be primal feasible due
         to changes of bounds of some basic variables, which turned to
         implicit artifical variables */
      /* compute simplex multipliers and reduced costs */
      ssx_eval_pi(ssx);
      ssx_eval_cbar(ssx);
      /* display initial progress of the search */
      show_progress(ssx, 1);
      /* main loop starts here */
      for (;;)
      {  /* display current progress of the search */
#if 0
         if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
#else
         if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
#endif
            show_progress(ssx, 1);
         /* we do not need to wait until all artificial variables have
            left the basis */
         if (mpq_sgn(bbar[0]) == 0)
         {  /* the sum of infeasibilities is zero, therefore the current
               solution is primal feasible for the original problem */
            ret = 0;
            break;
         }
         /* check if the iterations limit has been exhausted */
         if (ssx->it_lim == 0)
         {  ret = 2;
            break;
         }
         /* check if the time limit has been exhausted */
#if 0
         if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
#else
         if (ssx->tm_lim >= 0.0 &&
             ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
#endif
         {  ret = 3;
            break;
         }
         /* choose non-basic variable xN[q] */
         ssx_chuzc(ssx);
         /* if xN[q] cannot be chosen, the sum of infeasibilities is
            minimal but non-zero; therefore the original problem has no
            primal feasible solution */
         if (ssx->q == 0)
         {  ret = 1;
            break;
         }
         /* compute q-th column of the simplex table */
         ssx_eval_col(ssx);
         /* choose basic variable xB[p] */
         ssx_chuzr(ssx);
         /* the sum of infeasibilities cannot be negative, therefore
            the auxiliary lp problem cannot have unbounded solution */
         xassert(ssx->p != 0);
         /* update values of basic variables */
         ssx_update_bbar(ssx);
         if (ssx->p > 0)
         {  /* compute p-th row of the inverse inv(B) */
            ssx_eval_rho(ssx);
            /* compute p-th row of the simplex table */
            ssx_eval_row(ssx);
            xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
            /* update simplex multipliers */
            ssx_update_pi(ssx);
            /* update reduced costs of non-basic variables */
            ssx_update_cbar(ssx);
         }
         /* xB[p] is leaving the basis; if it is implicit artificial
            variable, the corresponding residual vanishes; therefore
            bounds of this variable should be restored to the original
            values */
         if (ssx->p > 0)
         {  k = Q_col[ssx->p]; /* x[k] = xB[p] */
            if (type[k] != orig_type[k])
            {  /* x[k] is implicit artificial variable */
               type[k] = orig_type[k];
               mpq_set(lb[k], orig_lb[k]);
               mpq_set(ub[k], orig_ub[k]);
               xassert(ssx->p_stat == SSX_NL || ssx->p_stat == SSX_NU);
               ssx->p_stat = (ssx->p_stat == SSX_NL ? SSX_NU : SSX_NL);
               if (type[k] == SSX_FX) ssx->p_stat = SSX_NS;
               /* nullify the objective coefficient at x[k] */
               mpq_set_si(coef[k], 0, 1);
               /* since coef[k] has been changed, we need to compute
                  new reduced cost of x[k], which it will have in the
                  adjacent basis */
               /* the formula d[j] = cN[j] - pi' * N[j] is used (note
                  that the vector pi is not changed, because it depends
                  on objective coefficients at basic variables, but in
                  the adjacent basis, for which the vector pi has been
                  just recomputed, x[k] is non-basic) */
               if (k <= m)
               {  /* x[k] is auxiliary variable */
                  mpq_neg(cbar[ssx->q], pi[k]);
               }
               else
               {  /* x[k] is structural variable */
                  int ptr;
                  mpq_t temp;
                  mpq_init(temp);
                  mpq_set_si(cbar[ssx->q], 0, 1);
                  for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
                  {  mpq_mul(temp, pi[A_ind[ptr]], A_val[ptr]);
                     mpq_add(cbar[ssx->q], cbar[ssx->q], temp);
                  }
                  mpq_clear(temp);
               }
            }
         }
         /* jump to the adjacent vertex of the polyhedron */
         ssx_change_basis(ssx);
         /* one simplex iteration has been performed */