glpssx01.c

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

      k = Q_col[m+q]; /* x[k] = xN[q] */
      if (k <= m)
      {  /* N[q] is a column of the unity matrix I */
         mpq_set_si(aq[k], 1, 1);
      }
      else
      {  /* N[q] is a column of the original constraint matrix -A */
         for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
            mpq_neg(aq[A_ind[ptr]], A_val[ptr]);
      }
      /* aq := inv(B) * aq */
      bfx_ftran(ssx->binv, aq, 1);
      /* aq := - aq */
      for (i = 1; i <= m; i++) mpq_neg(aq[i], aq[i]);
      return;
}

/*----------------------------------------------------------------------
// ssx_chuzc - choose pivot column.
//
// This routine chooses non-basic variable xN[q] whose reduced cost
// indicates possible improving of the objective function to enter it
// in the basis.
//
// Currently the standard (textbook) pricing is used, i.e. that
// non-basic variable is preferred which has greatest reduced cost (in
// magnitude).
//
// If xN[q] has been chosen, the routine stores its number q and also
// sets the flag q_dir that indicates direction in which xN[q] has to
// change (+1 means increasing, -1 means decreasing).
//
// If the choice cannot be made, because the current basic solution is
// dual feasible, the routine sets the number q to 0. */

void ssx_chuzc(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      int dir = (ssx->dir == SSX_MIN ? +1 : -1);
      int *Q_col = ssx->Q_col;
      int *stat = ssx->stat;
      mpq_t *cbar = ssx->cbar;
      int j, k, s, q, q_dir;
      double best, temp;
      /* nothing is chosen so far */
      q = 0, q_dir = 0, best = 0.0;
      /* look through the list of non-basic variables */
      for (j = 1; j <= n; j++)
      {  k = Q_col[m+j]; /* x[k] = xN[j] */
         s = dir * mpq_sgn(cbar[j]);
         if ((stat[k] == SSX_NF || stat[k] == SSX_NL) && s < 0 ||
             (stat[k] == SSX_NF || stat[k] == SSX_NU) && s > 0)
         {  /* reduced cost of xN[j] indicates possible improving of
               the objective function */
            temp = fabs(mpq_get_d(cbar[j]));
            xassert(temp != 0.0);
            if (q == 0 || best < temp)
               q = j, q_dir = - s, best = temp;
         }
      }
      ssx->q = q, ssx->q_dir = q_dir;
      return;
}

/*----------------------------------------------------------------------
// ssx_chuzr - choose pivot row.
//
// This routine looks through elements of q-th column of the simplex
// table and chooses basic variable xB[p] which should leave the basis.
//
// The choice is based on the standard (textbook) ratio test.
//
// If xB[p] has been chosen, the routine stores its number p and also
// sets its non-basic status p_stat which should be assigned to xB[p]
// when it has left the basis and become xN[q].
//
// Special case p < 0 means that xN[q] is double-bounded variable and
// it reaches its opposite bound before any basic variable does that,
// so the current basis remains unchanged.
//
// If the choice cannot be made, because xN[q] can infinitely change in
// the feasible direction, the routine sets the number p to 0. */

void ssx_chuzr(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      int *type = ssx->type;
      mpq_t *lb = ssx->lb;
      mpq_t *ub = ssx->ub;
      int *Q_col = ssx->Q_col;
      mpq_t *bbar = ssx->bbar;
      int q = ssx->q;
      mpq_t *aq = ssx->aq;
      int q_dir = ssx->q_dir;
      int i, k, s, t, p, p_stat;
      mpq_t teta, temp;
      mpq_init(teta);
      mpq_init(temp);
      xassert(1 <= q && q <= n);
      xassert(q_dir == +1 || q_dir == -1);
      /* nothing is chosen so far */
      p = 0, p_stat = 0;
      /* look through the list of basic variables */
      for (i = 1; i <= m; i++)
      {  s = q_dir * mpq_sgn(aq[i]);
         if (s < 0)
         {  /* xB[i] decreases */
            k = Q_col[i]; /* x[k] = xB[i] */
            t = type[k];
            if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
            {  /* xB[i] has finite lower bound */
               mpq_sub(temp, bbar[i], lb[k]);
               mpq_div(temp, temp, aq[i]);
               mpq_abs(temp, temp);
               if (p == 0 || mpq_cmp(teta, temp) > 0)
               {  p = i;
                  p_stat = (t == SSX_FX ? SSX_NS : SSX_NL);
                  mpq_set(teta, temp);
               }
            }
         }
         else if (s > 0)
         {  /* xB[i] increases */
            k = Q_col[i]; /* x[k] = xB[i] */
            t = type[k];
            if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
            {  /* xB[i] has finite upper bound */
               mpq_sub(temp, bbar[i], ub[k]);
               mpq_div(temp, temp, aq[i]);
               mpq_abs(temp, temp);
               if (p == 0 || mpq_cmp(teta, temp) > 0)
               {  p = i;
                  p_stat = (t == SSX_FX ? SSX_NS : SSX_NU);
                  mpq_set(teta, temp);
               }
            }
         }
         /* if something has been chosen and the ratio test indicates
            exact degeneracy, the search can be finished */
         if (p != 0 && mpq_sgn(teta) == 0) break;
      }
      /* if xN[q] is double-bounded, check if it can reach its opposite
         bound before any basic variable */
      k = Q_col[m+q]; /* x[k] = xN[q] */
      if (type[k] == SSX_DB)
      {  mpq_sub(temp, ub[k], lb[k]);
         if (p == 0 || mpq_cmp(teta, temp) > 0)
         {  p = -1;
            p_stat = -1;
            mpq_set(teta, temp);
         }
      }
      ssx->p = p;
      ssx->p_stat = p_stat;
      /* if xB[p] has been chosen, determine its actual change in the
         adjacent basis (it has the same sign as q_dir) */
      if (p != 0)
      {  xassert(mpq_sgn(teta) >= 0);
         if (q_dir > 0)
            mpq_set(ssx->delta, teta);
         else
            mpq_neg(ssx->delta, teta);
      }
      mpq_clear(teta);
      mpq_clear(temp);
      return;
}

/*----------------------------------------------------------------------
// ssx_update_bbar - update values of basic variables.
//
// This routine recomputes the current values of basic variables for
// the adjacent basis.
//
// The simplex table for the current basis is the following:
//
//    xB[i] = sum{j in 1..n} alfa[i,j] * xN[q],  i = 1,...,m
//
// therefore
//
//    delta xB[i] = alfa[i,q] * delta xN[q],  i = 1,...,m
//
// where delta xN[q] = xN.new[q] - xN[q] is the change of xN[q] in the
// adjacent basis, and delta xB[i] = xB.new[i] - xB[i] is the change of
// xB[i]. This gives formulae for recomputing values of xB[i]:
//
//    xB.new[p] = xN[q] + delta xN[q]
//
// (because xN[q] becomes xB[p] in the adjacent basis), and
//
//    xB.new[i] = xB[i] + alfa[i,q] * delta xN[q],  i != p
//
// for other basic variables. */

void ssx_update_bbar(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      mpq_t *bbar = ssx->bbar;
      mpq_t *cbar = ssx->cbar;
      int p = ssx->p;
      int q = ssx->q;
      mpq_t *aq = ssx->aq;
      int i;
      mpq_t temp;
      mpq_init(temp);
      xassert(1 <= q && q <= n);
      if (p < 0)
      {  /* xN[q] is double-bounded and goes to its opposite bound */
         /* nop */;
      }
      else
      {  /* xN[q] becomes xB[p] in the adjacent basis */
         /* xB.new[p] = xN[q] + delta xN[q] */
         xassert(1 <= p && p <= m);
         ssx_get_xNj(ssx, q, temp);
         mpq_add(bbar[p], temp, ssx->delta);
      }
      /* update values of other basic variables depending on xN[q] */
      for (i = 1; i <= m; i++)
      {  if (i == p) continue;
         /* xB.new[i] = xB[i] + alfa[i,q] * delta xN[q] */
         if (mpq_sgn(aq[i]) == 0) continue;
         mpq_mul(temp, aq[i], ssx->delta);
         mpq_add(bbar[i], bbar[i], temp);
      }
#if 1
      /* update value of the objective function */
      /* z.new = z + d[q] * delta xN[q] */
      mpq_mul(temp, cbar[q], ssx->delta);
      mpq_add(bbar[0], bbar[0], temp);
#endif
      mpq_clear(temp);
      return;
}

/*----------------------------------------------------------------------
-- ssx_update_pi - update simplex multipliers.
--
-- This routine recomputes the vector of simplex multipliers for the
-- adjacent basis. */

void ssx_update_pi(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      mpq_t *pi = ssx->pi;
      mpq_t *cbar = ssx->cbar;
      int p = ssx->p;
      int q = ssx->q;
      mpq_t *aq = ssx->aq;
      mpq_t *rho = ssx->rho;
      int i;
      mpq_t new_dq, temp;
      mpq_init(new_dq);
      mpq_init(temp);
      xassert(1 <= p && p <= m);
      xassert(1 <= q && q <= n);
      /* compute d[q] in the adjacent basis */
      mpq_div(new_dq, cbar[q], aq[p]);
      /* update the vector of simplex multipliers */
      for (i = 1; i <= m; i++)
      {  if (mpq_sgn(rho[i]) == 0) continue;
         mpq_mul(temp, new_dq, rho[i]);
         mpq_sub(pi[i], pi[i], temp);
      }
      mpq_clear(new_dq);
      mpq_clear(temp);
      return;
}

/*----------------------------------------------------------------------
// ssx_update_cbar - update reduced costs of non-basic variables.
//
// This routine recomputes the vector of reduced costs of non-basic
// variables for the adjacent basis. */

void ssx_update_cbar(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      mpq_t *cbar = ssx->cbar;
      int p = ssx->p;
      int q = ssx->q;
      mpq_t *ap = ssx->ap;
      int j;
      mpq_t temp;
      mpq_init(temp);
      xassert(1 <= p && p <= m);
      xassert(1 <= q && q <= n);
      /* compute d[q] in the adjacent basis */
      /* d.new[q] = d[q] / alfa[p,q] */
      mpq_div(cbar[q], cbar[q], ap[q]);
      /* update reduced costs of other non-basic variables */
      for (j = 1; j <= n; j++)
      {  if (j == q) continue;
         /* d.new[j] = d[j] - (alfa[p,j] / alfa[p,q]) * d[q] */
         if (mpq_sgn(ap[j]) == 0) continue;
         mpq_mul(temp, ap[j], cbar[q]);
         mpq_sub(cbar[j], cbar[j], temp);
      }
      mpq_clear(temp);
      return;
}

/*----------------------------------------------------------------------
// ssx_change_basis - change current basis to adjacent one.
//
// This routine changes the current basis to the adjacent one swapping
// basic variable xB[p] and non-basic variable xN[q]. */

void ssx_change_basis(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      int *type = ssx->type;
      int *stat = ssx->stat;
      int *Q_row = ssx->Q_row;
      int *Q_col = ssx->Q_col;
      int p = ssx->p;
      int q = ssx->q;
      int p_stat = ssx->p_stat;
      int k, kp, kq;
      if (p < 0)
      {  /* special case: xN[q] goes to its opposite bound */
         xassert(1 <= q && q <= n);
         k = Q_col[m+q]; /* x[k] = xN[q] */
         xassert(type[k] == SSX_DB);
         switch (stat[k])
         {  case SSX_NL:
               stat[k] = SSX_NU;
               break;
            case SSX_NU:
               stat[k] = SSX_NL;
               break;
            default:
               xassert(stat != stat);
         }
      }
      else
      {  /* xB[p] leaves the basis, xN[q] enters the basis */
         xassert(1 <= p && p <= m);
         xassert(1 <= q && q <= n);
         kp = Q_col[p];   /* x[kp] = xB[p] */
         kq = Q_col[m+q]; /* x[kq] = xN[q] */
         /* check non-basic status of xB[p] which becomes xN[q] */
         switch (type[kp])
         {  case SSX_FR:
               xassert(p_stat == SSX_NF);
               break;
            case SSX_LO:
               xassert(p_stat == SSX_NL);
               break;
            case SSX_UP:
               xassert(p_stat == SSX_NU);
               break;
            case SSX_DB:
               xassert(p_stat == SSX_NL || p_stat == SSX_NU);
               break;
            case SSX_FX:
               xassert(p_stat == SSX_NS);
               break;
            default:
               xassert(type != type);
         }
         /* swap xB[p] and xN[q] */
         stat[kp] = (char)p_stat, stat[kq] = SSX_BS;
         Q_row[kp] = m+q, Q_row[kq] = p;
         Q_col[p] = kq, Q_col[m+q] = kp;
         /* update factorization of the basis matrix */
         if (bfx_update(ssx->binv, p))
         {  if (ssx_factorize(ssx))
               xassert(("Internal error: basis matrix is singular", 0));
         }
      }
      return;
}

/*----------------------------------------------------------------------
// ssx_delete - delete simplex solver workspace.
//
// This routine deletes the simplex solver workspace freeing all the
// memory allocated to this object. */

void ssx_delete(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      int nnz = ssx->A_ptr[n+1]-1;
      int i, j, k;
      xfree(ssx->type);
      for (k = 1; k <= m+n; k++) mpq_clear(ssx->lb[k]);
      xfree(ssx->lb);
      for (k = 1; k <= m+n; k++) mpq_clear(ssx->ub[k]);
      xfree(ssx->ub);
      for (k = 0; k <= m+n; k++) mpq_clear(ssx->coef[k]);
      xfree(ssx->coef);
      xfree(ssx->A_ptr);
      xfree(ssx->A_ind);
      for (k = 1; k <= nnz; k++) mpq_clear(ssx->A_val[k]);
      xfree(ssx->A_val);
      xfree(ssx->stat);
      xfree(ssx->Q_row);
      xfree(ssx->Q_col);
      bfx_delete_binv(ssx->binv);
      for (i = 0; i <= m; i++) mpq_clear(ssx->bbar[i]);
      xfree(ssx->bbar);
      for (i = 1; i <= m; i++) mpq_clear(ssx->pi[i]);
      xfree(ssx->pi);
      for (j = 1; j <= n; j++) mpq_clear(ssx->cbar[j]);
      xfree(ssx->cbar);
      for (i = 1; i <= m; i++) mpq_clear(ssx->rho[i]);
      xfree(ssx->rho);
      for (j = 1; j <= n; j++) mpq_clear(ssx->ap[j]);
      xfree(ssx->ap);
      for (i = 1; i <= m; i++) mpq_clear(ssx->aq[i]);
      xfree(ssx->aq);
      mpq_clear(ssx->delta);
      xfree(ssx);
      return;
}

/* eof */