glpspx02.c

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

            for (ptr = beg; ptr < end; ptr++)
               h[A_ind[ptr]] = A_val[ptr];
         }
         /* solve system B * alfa = h */
         xassert(csa->valid);
         bfd_ftran(csa->bfd, alfa);
         /* gamma[i] := gamma[i] + alfa[i,j]^2 */
         for (i = 1; i <= m; i++)
         {  k = head[i]; /* x[k] = xB[i] */
            if (type[k] != GLP_FR)
               gamma[i] += alfa[i] * alfa[i];
         }
      }
      return;
}

/***********************************************************************
*  chuzr - choose basic variable (row of the simplex table)
*
*  This routine chooses basic variable xB[p] having largest weighted
*  bound violation:
*
*     |r[p]| / sqrt(gamma[p]) = max  |r[i]| / sqrt(gamma[i]),
*                              i in I
*
*            / lB[i] - beta[i], if beta[i] < lB[i]
*            |
*     r[i] = < 0,               if lB[i] <= beta[i] <= uB[i]
*            |
*            \ uB[i] - beta[i], if beta[i] > uB[i]
*
*  where beta[i] is primal value of xB[i] in the current basis, lB[i]
*  and uB[i] are lower and upper bounds of xB[i], I is a subset of
*  eligible basic variables, which significantly violates their bounds,
*  gamma[i] is the steepest edge coefficient.
*
*  If |r[i]| is less than a specified tolerance, xB[i] is not included
*  in I and therefore ignored.
*
*  If I is empty and no variable has been chosen, p is set to 0. */

static void chuzr(struct csa *csa, double tol_bnd)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      char *type = csa->type;
      double *lb = csa->lb;
      double *ub = csa->ub;
      int *head = csa->head;
      double *bbar = csa->bbar;
      double *gamma = csa->gamma;
      int i, k, p;
      double delta, best, eps, ri, temp;
      /* nothing is chosen so far */
      p = 0, delta = 0.0, best = 0.0;
      /* look through the list of basic variables */
      for (i = 1; i <= m; i++)
      {  k = head[i]; /* x[k] = xB[i] */
#ifdef GLP_DEBUG
         xassert(1 <= k && k <= m+n);
#endif
         /* determine bound violation ri[i] */
         ri = 0.0;
         if (type[k] == GLP_LO || type[k] == GLP_DB ||
             type[k] == GLP_FX)
         {  /* xB[i] has lower bound */
            eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
            if (bbar[i] < lb[k] - eps)
            {  /* and significantly violates it */
               ri = lb[k] - bbar[i];
            }
         }
         if (type[k] == GLP_UP || type[k] == GLP_DB ||
             type[k] == GLP_FX)
         {  /* xB[i] has upper bound */
            eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
            if (bbar[i] > ub[k] + eps)
            {  /* and significantly violates it */
               ri = ub[k] - bbar[i];
            }
         }
         /* if xB[i] is not eligible, skip it */
         if (ri == 0.0) continue;
         /* xB[i] is eligible basic variable; choose one with largest
            weighted bound violation */
#ifdef GLP_DEBUG
         xassert(gamma[i] >= 0.0);
#endif
         temp = gamma[i];
         if (temp < DBL_EPSILON) temp = DBL_EPSILON;
         temp = (ri * ri) / temp;
         if (best < temp)
            p = i, delta = ri, best = temp;
      }
      /* store the index of basic variable xB[p] chosen and its change
         in the adjacent basis */
      csa->p = p;
      csa->delta = delta;
      return;
}

#if 1 /* copied from primal */
/***********************************************************************
*  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;
}
#endif

#if 1 /* copied from primal */
/***********************************************************************
*  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;
}
#endif

#if 1 /* copied from primal */
/***********************************************************************
*  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;
}
#endif

/***********************************************************************
*  sort_trow - sort pivot row of the simplex table
*
*  This routine reorders the list of non-zero elements of the pivot
*  row to put significant elements, whose magnitude is not less than
*  a specified tolerance, in front of the list, and stores the number
*  of significant elements in trow_num. */

static void sort_trow(struct csa *csa, double tol_piv)
{
#ifdef GLP_DEBUG
      int n = csa->n;
      char *stat = csa->stat;
#endif
      int nnz = csa->trow_nnz;
      int *trow_ind = csa->trow_ind;
      double *trow_vec = csa->trow_vec;
      int j, num, pos;
      double big, eps, temp;
      /* compute infinity (maximum) norm of the row */
      big = 0.0;
      for (pos = 1; pos <= nnz; pos++)
      {
#ifdef GLP_DEBUG
         j = trow_ind[pos];
         xassert(1 <= j && j <= n);
         xassert(stat[j] != GLP_NS);
#endif
         temp = fabs(trow_vec[trow_ind[pos]]);
         if (big < temp) big = temp;
      }
      csa->trow_max = big;
      /* determine absolute pivot tolerance */
      eps = tol_piv * (1.0 + 0.01 * big);
      /* move significant row components to the front of the list */
      for (num = 0; num < nnz; )
      {  j = trow_ind[nnz];
         if (fabs(trow_vec[j]) < eps)
            nnz--;
         else
         {  num++;
            trow_ind[nnz] = trow_ind[num];
            trow_ind[num] = j;
         }
      }
      csa->trow_num = num;
      return;
}

/***********************************************************************
*  chuzc - choose non-basic variable (column of the simplex table)
*
*  This routine chooses non-basic variable xN[q], which being entered
*  in the basis keeps dual feasibility of the basic solution.
*
*  The parameter rtol is a relative tolerance used to relax zero bounds
*  of reduced costs of non-basic variables. If rtol = 0, the routine
*  implements the standard ratio test. Otherwise, if rtol > 0, the
*  routine implements Harris' two-pass ratio test. In the latter case
*  rtol should be about three times less than a tolerance used to check
*  dual feasibility. */

static void chuzc(struct csa *csa, double rtol)
{
#ifdef GLP_DEBUG
      int m = csa->m;
      int n = csa->n;
#endif
      char *stat = csa->stat;
      double *cbar = csa->cbar;
#ifdef GLP_DEBUG
      int p = csa->p;
#endif
      double delta = csa->delta;
      int *trow_ind = csa->trow_ind;
      double *trow_vec = csa->trow_vec;
      int trow_num = csa->trow_num;
      int j, pos, q;
      double alfa, big, s, t, teta, tmax;
#ifdef GLP_DEBUG
      xassert(1 <= p && p <= m);
#endif
      /* delta > 0 means that xB[p] violates its lower bound and goes
         to it in the adjacent basis, so lambdaB[p] is increasing from
         its lower zero bound;
         delta < 0 means that xB[p] violates its upper bound and goes
         to it in the adjacent basis, so lambdaB[p] is decreasing from
         its upper zero bound */
#ifdef GLP_DEBUG
      xassert(delta != 0.0);
#endif
      /* s := sign(delta) */
      s = (delta > 0.0 ? +1.0 : -1.0);
      /*** FIRST PASS ***/
      /* nothing is chosen so far */
      q = 0, teta = DBL_MAX, big = 0.0;
      /* walk through significant elements of the pivot row */
      for (pos = 1; pos <= trow_num; pos++)
      {  j = trow_ind[pos];
#ifdef GLP_DEBUG
         xassert(1 <= j && j <= n);
#endif
         alfa = s * trow_vec[j];
#ifdef GLP_DEBUG
         xassert(alfa != 0.0);
#endif
         /* lambdaN[j] = ... - alfa * lambdaB[p] - ..., and due to s we
            need to consider only increasing lambdaB[p] */
         if (alfa > 0.0)
         {  /* lambdaN[j] is decreasing */
            if (stat[j] == GLP_NL || stat[j] == GLP_NF)
            {  /* lambdaN[j] has zero lower bound */
               t = (cbar[j] + rtol) / alfa;
            }
            else
            {  /* lambdaN[j] has no lower bound */
               continue;
            }
         }
         else
         {  /* lambdaN[j] is increasing */
            if (stat[j] == GLP_NU || stat[j] == GLP_NF)
            {  /* lambdaN[j] has zero upper bound */
               t = (cbar[j] - rtol) / alfa;
            }
            else
            {  /* lambdaN[j] has no upper bound */
               continue;
            }
         }
         /* t is a change of lambdaB[p], on which lambdaN[j] reaches
            its zero bound (possibly relaxed); since the basic solution
            is assumed to be dual feasible, t has to be non-negative by
            definition; however, it may happen that lambdaN[j] slightly
            (i.e. within a tolerance) violates its zero bound, that
            leads to negative t; in the latter case, if xN[j] is chosen,
            negative t means that lambdaB[p] changes in wrong direction
            that may cause wrong results on updating reduced costs;
            thus, if t is negative, we should replace it by exact zero
            assuming that lambdaN[j] is exactly on its zero bound, and
            violation appears due to round-off errors */
         if (t < 0.0) t = 0.0;
         /* apply minimal ratio test */
         if (teta > t || teta == t && big < fabs(alfa))
            q = j, teta = t, big = fabs(alfa);
      }
      /* the second pass is skipped in the following cases: */
      /* if the standard ratio test is used */
      if (rtol == 0.0) goto done;
      /* if no non-basic variable has been chosen on the first pass */
      if (q == 0) goto done;
      /* if lambdaN[q] prevents lambdaB[p] from any change */
      if (teta == 0.0) goto done;
      /*** SECOND PASS ***/
      /* here tmax is a maximal change of lambdaB[p], on which the
         solution remains dual feasible within a tolerance */
#if 0
      tmax = (1.0 + 10.0 * DBL_EPSILON) * teta;
#else
      tmax = teta;
#endif
      /* nothing is chosen so far */
      q = 0, teta = DBL_MAX, big = 0.0;
      /* walk through significant elements of the pivot row */
      for (pos = 1; pos <= trow_num; pos++)
      {  j = trow_ind[pos];
#ifdef GLP_DEBUG
         xassert(1 <= j && j <= n);
#endif