glpspx01.c

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

            consider the only case when xN[q] is increasing */
         if (alfa > 0.0)
         {  /* xB[i] is increasing */
            if (phase == 1 && coef[k] < 0.0)
            {  /* xB[i] violates its lower bound, which plays the role
                  of an upper bound on phase I */
               t = (lb[k] - bbar[i]) / alfa;
               i_stat = GLP_NL;
            }
            else if (phase == 1 && coef[k] > 0.0)
            {  /* xB[i] violates its upper bound, which plays the role
                  of an lower bound on phase I */
               continue;
            }
            else if (type[k] == GLP_UP || type[k] == GLP_DB ||
                     type[k] == GLP_FX)
            {  /* xB[i] is within its bounds and has an upper bound */
               t = (ub[k] - bbar[i]) / alfa;
               i_stat = GLP_NU;
            }
            else
            {  /* xB[i] is within its bounds and has no upper bound */
               continue;
            }
         }
         else
         {  /* xB[i] is decreasing */
            if (phase == 1 && coef[k] > 0.0)
            {  /* xB[i] violates its upper bound, which plays the role
                  of an lower bound on phase I */
               t = (ub[k] - bbar[i]) / alfa;
               i_stat = GLP_NU;
            }
            else if (phase == 1 && coef[k] < 0.0)
            {  /* xB[i] violates its lower bound, which plays the role
                  of an upper bound on phase I */
               continue;
            }
            else if (type[k] == GLP_LO || type[k] == GLP_DB ||
                     type[k] == GLP_FX)
            {  /* xB[i] is within its bounds and has an lower bound */
               t = (lb[k] - bbar[i]) / alfa;
               i_stat = GLP_NL;
            }
            else
            {  /* xB[i] is within its bounds and has no lower bound */
               continue;
            }
         }
         /* (see comments for the first pass) */
         if (t < 0.0) t = 0.0;
         /* t is a change of xN[q], on which xB[i] reaches its bound;
            if t <= tmax, all basic variables can violate their bounds
            only within relaxation tolerance delta; we can use this
            freedom and choose basic variable having largest influence
            coefficient to avoid possible numeric instability */
         if (t <= tmax && big < fabs(alfa))
            p = i, p_stat = i_stat, teta = t, big = fabs(alfa);
      }
      /* something must be chosen on the second pass */
      xassert(p != 0);
done: /* store the index and status of basic variable xB[p] chosen */
      csa->p = p;
      if (p > 0 && type[head[p]] == GLP_FX)
         csa->p_stat = GLP_NS;
      else
         csa->p_stat = p_stat;
      /* store corresponding change of non-basic variable xN[q] */
#ifdef GLP_DEBUG
      xassert(teta >= 0.0);
#endif
      csa->teta = s * teta;
      return;
}

/***********************************************************************
*  eval_rho - compute pivot row of the inverse
*
*  This routine computes the pivot (p-th) row of the inverse inv(B),
*  which corresponds to basic variable xB[p] chosen:
*
*     rho = inv(B') * e[p],
*
*  where B' is a matrix transposed to the current basis matrix, e[p]
*  is unity vector. */

static void eval_rho(struct csa *csa, double rho[])
{     int m = csa->m;
      int p = csa->p;
      double *e = rho;
      int i;
#ifdef GLP_DEBUG
      xassert(1 <= p && p <= m);
#endif
      /* construct the right-hand side vector e[p] */
      for (i = 1; i <= m; i++)
         e[i] = 0.0;
      e[p] = 1.0;
      /* solve system B'* rho = e[p] */
      xassert(csa->valid);
      bfd_btran(csa->bfd, rho);
      return;
}

/***********************************************************************
*  refine_rho - refine pivot row of the inverse
*
*  This routine refines the pivot row of the inverse inv(B) assuming
*  that it was previously computed by the routine eval_rho. */

static void refine_rho(struct csa *csa, double rho[])
{     int m = csa->m;
      int p = csa->p;
      double *e = csa->work3;
      int i;
#ifdef GLP_DEBUG
      xassert(1 <= p && p <= m);
#endif
      /* construct the right-hand side vector e[p] */
      for (i = 1; i <= m; i++)
         e[i] = 0.0;
      e[p] = 1.0;
      /* refine solution of B'* rho = e[p] */
      refine_btran(csa, e, rho);
      return;
}

/***********************************************************************
*  eval_trow - compute pivot row of the simplex table
*
*  This routine computes the pivot row of the simplex table, which
*  corresponds to basic variable xB[p] chosen.
*
*  The pivot row is the following vector:
*
*     trow = T'* e[p] = - N'* inv(B') * e[p] = - N' * rho,
*
*  where rho is the pivot row of the inverse inv(B) previously computed
*  by the routine eval_rho.
*
*  Note that elements of the pivot row corresponding to fixed non-basic
*  variables are not computed. */

static void eval_trow(struct csa *csa, double rho[])
{     int m = csa->m;
      int n = csa->n;
#ifdef GLP_DEBUG
      char *stat = csa->stat;
#endif
      int *N_ptr = csa->N_ptr;
      int *N_len = csa->N_len;
      int *N_ind = csa->N_ind;
      double *N_val = csa->N_val;
      int *trow_ind = csa->trow_ind;
      double *trow_vec = csa->trow_vec;
      int i, j, beg, end, ptr, nnz;
      double temp;
      /* clear the pivot row */
      for (j = 1; j <= n; j++)
         trow_vec[j] = 0.0;
      /* compute the pivot row as a linear combination of rows of the
         matrix N: trow = - rho[1] * N'[1] - ... - rho[m] * N'[m] */
      for (i = 1; i <= m; i++)
      {  temp = rho[i];
         if (temp == 0.0) continue;
         /* trow := trow - rho[i] * N'[i] */
         beg = N_ptr[i];
         end = beg + N_len[i];
         for (ptr = beg; ptr < end; ptr++)
         {
#ifdef GLP_DEBUG
            j = N_ind[ptr];
            xassert(1 <= j && j <= n);
            xassert(stat[j] != GLP_NS);
#endif
            trow_vec[N_ind[ptr]] -= temp * N_val[ptr];
         }
      }
      /* construct sparse pattern of the pivot row */
      nnz = 0;
      for (j = 1; j <= n; j++)
      {  if (trow_vec[j] != 0.0)
            trow_ind[++nnz] = j;
      }
      csa->trow_nnz = nnz;
      return;
}

/***********************************************************************
*  update_bbar - update values of basic variables
*
*  This routine updates values of all basic variables for the adjacent
*  basis. */

static void update_bbar(struct csa *csa)
{
#ifdef GLP_DEBUG
      int m = csa->m;
      int n = csa->n;
#endif
      double *bbar = csa->bbar;
      int q = csa->q;
      int tcol_nnz = csa->tcol_nnz;
      int *tcol_ind = csa->tcol_ind;
      double *tcol_vec = csa->tcol_vec;
      int p = csa->p;
      double teta = csa->teta;
      int i, pos;
#ifdef GLP_DEBUG
      xassert(1 <= q && q <= n);
      xassert(p < 0 || 1 <= p && p <= m);
#endif
      /* if xN[q] leaves the basis, compute its value in the adjacent
         basis, where it will replace xB[p] */
      if (p > 0)
         bbar[p] = get_xN(csa, q) + teta;
      /* update values of other basic variables (except xB[p], because
         it will be replaced by xN[q]) */
      if (teta == 0.0) goto done;
      for (pos = 1; pos <= tcol_nnz; pos++)
      {  i = tcol_ind[pos];
         /* skip xB[p] */
         if (i == p) continue;
         /* (change of xB[i]) = alfa[i,q] * (change of xN[q]) */
         bbar[i] += tcol_vec[i] * teta;
      }
done: return;
}

/***********************************************************************
*  reeval_cost - recompute reduced cost of non-basic variable xN[q]
*
*  This routine recomputes reduced cost of non-basic variable xN[q] for
*  the current basis more accurately using its direct definition:
*
*     d[q] = cN[q] - N'[q] * pi =
*
*          = cN[q] - N'[q] * (inv(B') * cB) =
*
*          = cN[q] - (cB' * inv(B) * N[q]) =
*
*          = cN[q] + cB' * (pivot column).
*
*  It is assumed that the pivot column of the simplex table is already
*  computed. */

static double reeval_cost(struct csa *csa)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      double *coef = csa->coef;
      int *head = csa->head;
      int q = csa->q;
      int tcol_nnz = csa->tcol_nnz;
      int *tcol_ind = csa->tcol_ind;
      double *tcol_vec = csa->tcol_vec;
      int i, pos;
      double dq;
#ifdef GLP_DEBUG
      xassert(1 <= q && q <= n);
#endif
      dq = coef[head[m+q]];
      for (pos = 1; pos <= tcol_nnz; pos++)
      {  i = tcol_ind[pos];
#ifdef GLP_DEBUG
         xassert(1 <= i && i <= m);
#endif
         dq += coef[head[i]] * tcol_vec[i];
      }
      return dq;
}

/***********************************************************************
*  update_cbar - update reduced costs of non-basic variables
*
*  This routine updates reduced costs of all (except fixed) non-basic
*  variables for the adjacent basis. */

static void update_cbar(struct csa *csa)
{
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      double *cbar = csa->cbar;
      int q = csa->q;
      int trow_nnz = csa->trow_nnz;
      int *trow_ind = csa->trow_ind;
      double *trow_vec = csa->trow_vec;
      int j, pos;
      double new_dq;
#ifdef GLP_DEBUG
      xassert(1 <= q && q <= n);
#endif
      /* compute reduced cost of xB[p] in the adjacent basis, where it
         will replace xN[q] */
#ifdef GLP_DEBUG
      xassert(trow_vec[q] != 0.0);
#endif
      new_dq = (cbar[q] /= trow_vec[q]);
      /* update reduced costs of other non-basic variables (except
         xN[q], because it will be replaced by xB[p]) */
      for (pos = 1; pos <= trow_nnz; pos++)
      {  j = trow_ind[pos];
         /* skip xN[q] */
         if (j == q) continue;
         cbar[j] -= trow_vec[j] * new_dq;
      }
      return;
}

/***********************************************************************
*  update_gamma - update steepest edge coefficients
*
*  This routine updates steepest-edge coefficients for the adjacent
*  basis. */

static void update_gamma(struct csa *csa)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      char *type = csa->type;
      int *A_ptr = csa->A_ptr;
      int *A_ind = csa->A_ind;
      double *A_val = csa->A_val;
      int *head = csa->head;
      char *refsp = csa->refsp;
      double *gamma = csa->gamma;
      int q = csa->q;
      int tcol_nnz = csa->tcol_nnz;
      int *tcol_ind = csa->tcol_ind;
      double *tcol_vec = csa->tcol_vec;
      int p = csa->p;
      int trow_nnz = csa->trow_nnz;
      int *trow_ind = csa->trow_ind;
      double *trow_vec = csa->trow_vec;
      double *u = csa->work3;
      int i, j, k, pos, beg, end, ptr;
      double gamma_q, delta_q, pivot, s, t, t1, t2;
#ifdef GLP_DEBUG
      xassert(1 <= p && p <= m);
      xassert(1 <= q && q <= n);
#endif
      /* the basis changes, so decrease the count */
      xassert(csa->refct > 0);
      csa->refct--;
      /* recompute gamma[q] for the current basis more accurately and
         compute auxiliary vector u */
      gamma_q = delta_q = (refsp[head[m+q]] ? 1.0 : 0.0);
      for (i = 1; i <= m; i++) u[i] = 0.0;
      for (pos = 1; pos <= tcol_nnz; pos++)
      {  i = tcol_ind[pos];
         if (refsp[head[i]])
         {  u[i] = t = tcol_vec[i];
            gamma_q += t * t;
         }
         else
            u[i] = 0.0;
      }
      xassert(csa->valid);
      bfd_btran(csa->bfd, u);
      /* update gamma[k] for other non-basic variables (except fixed
         variables and xN[q], because it will be replaced by xB[p]) */
      pivot = trow_vec[q];
#ifdef GLP_DEBUG
      xassert(pivot != 0.0);
#endif
      for (pos = 1; pos <= trow_nnz; pos++)
      {  j = trow_ind[pos];
         /* skip xN[q] */
         if (j == q) continue;
         /* compute t */
         t = trow_vec[j] / pivot;
         /* compute inner product s = N'[j] * u */
         k = head[m+j]; /* x[k] = xN[j] */
         if (k <= m)
            s = u[k];
         else
         {  s = 0.0;
            beg = A_ptr[k-m];
            end = A_ptr[k-m+1];
            for (ptr = beg