glpspx01.c

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

      double *cbar = csa->cbar;
      double *gamma = csa->gamma;
      int j, q;
      double dj, best, temp;
      /* nothing is chosen so far */
      q = 0, best = 0.0;
      /* look through the list of non-basic variables */
      for (j = 1; j <= n; j++)
      {  dj = cbar[j];
         switch (stat[j])
         {  case GLP_NL:
               /* xN[j] can increase */
               if (dj >= - tol_dj) continue;
               break;
            case GLP_NU:
               /* xN[j] can decrease */
               if (dj <= + tol_dj) continue;
               break;
            case GLP_NF:
               /* xN[j] can change in any direction */
               if (- tol_dj <= dj && dj <= + tol_dj) continue;
               break;
            case GLP_NS:
               /* xN[j] cannot change at all */
               continue;
            default:
               xassert(stat != stat);
         }
         /* xN[j] is eligible non-basic variable; choose one which has
            largest weighted reduced cost */
#ifdef GLP_DEBUG
         xassert(gamma[j] > 0.0);
#endif
         temp = (dj * dj) / gamma[j];
         if (best < temp)
            q = j, best = temp;
      }
      /* store the index of non-basic variable xN[q] chosen */
      csa->q = q;
      return;
}

/***********************************************************************
*  eval_tcol - compute pivot column of the simplex table
*
*  This routine computes the pivot column of the simplex table, which
*  corresponds to non-basic variable xN[q] chosen.
*
*  The pivot column is the following vector:
*
*     tcol = T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q],
*
*  where B is the current basis matrix, N[q] is a column of the matrix
*  (I|-A) corresponding to variable xN[q]. */

static void eval_tcol(struct csa *csa)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      int *head = csa->head;
      int q = csa->q;
      int *tcol_ind = csa->tcol_ind;
      double *tcol_vec = csa->tcol_vec;
      double *h = csa->tcol_vec;
      int i, k, nnz;
#ifdef GLP_DEBUG
      xassert(1 <= q && q <= n);
#endif
      k = head[m+q]; /* x[k] = xN[q] */
#ifdef GLP_DEBUG
      xassert(1 <= k && k <= m+n);
#endif
      /* construct the right-hand side vector h = - N[q] */
      for (i = 1; i <= m; i++)
         h[i] = 0.0;
      if (k <= m)
      {  /* N[q] is k-th column of submatrix I */
         h[k] = -1.0;
      }
      else
      {  /* N[q] is (k-m)-th column of submatrix (-A) */
         int *A_ptr = csa->A_ptr;
         int *A_ind = csa->A_ind;
         double *A_val = csa->A_val;
         int beg, end, ptr;
         beg = A_ptr[k-m];
         end = A_ptr[k-m+1];
         for (ptr = beg; ptr < end; ptr++)
            h[A_ind[ptr]] = A_val[ptr];
      }
      /* solve system B * tcol = h */
      xassert(csa->valid);
      bfd_ftran(csa->bfd, tcol_vec);
      /* construct sparse pattern of the pivot column */
      nnz = 0;
      for (i = 1; i <= m; i++)
      {  if (tcol_vec[i] != 0.0)
            tcol_ind[++nnz] = i;
      }
      csa->tcol_nnz = nnz;
      return;
}

/***********************************************************************
*  refine_tcol - refine pivot column of the simplex table
*
*  This routine refines the pivot column of the simplex table assuming
*  that it was previously computed by the routine eval_tcol. */

static void refine_tcol(struct csa *csa)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      int *head = csa->head;
      int q = csa->q;
      int *tcol_ind = csa->tcol_ind;
      double *tcol_vec = csa->tcol_vec;
      double *h = csa->work3;
      int i, k, nnz;
#ifdef GLP_DEBUG
      xassert(1 <= q && q <= n);
#endif
      k = head[m+q]; /* x[k] = xN[q] */
#ifdef GLP_DEBUG
      xassert(1 <= k && k <= m+n);
#endif
      /* construct the right-hand side vector h = - N[q] */
      for (i = 1; i <= m; i++)
         h[i] = 0.0;
      if (k <= m)
      {  /* N[q] is k-th column of submatrix I */
         h[k] = -1.0;
      }
      else
      {  /* N[q] is (k-m)-th column of submatrix (-A) */
         int *A_ptr = csa->A_ptr;
         int *A_ind = csa->A_ind;
         double *A_val = csa->A_val;
         int beg, end, ptr;
         beg = A_ptr[k-m];
         end = A_ptr[k-m+1];
         for (ptr = beg; ptr < end; ptr++)
            h[A_ind[ptr]] = A_val[ptr];
      }
      /* refine solution of B * tcol = h */
      refine_ftran(csa, h, tcol_vec);
      /* construct sparse pattern of the pivot column */
      nnz = 0;
      for (i = 1; i <= m; i++)
      {  if (tcol_vec[i] != 0.0)
            tcol_ind[++nnz] = i;
      }
      csa->tcol_nnz = nnz;
      return;
}

/***********************************************************************
*  sort_tcol - sort pivot column of the simplex table
*
*  This routine reorders the list of non-zero elements of the pivot
*  column 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 tcol_num. */

static void sort_tcol(struct csa *csa, double tol_piv)
{
#ifdef GLP_DEBUG
      int m = csa->m;
#endif
      int nnz = csa->tcol_nnz;
      int *tcol_ind = csa->tcol_ind;
      double *tcol_vec = csa->tcol_vec;
      int i, num, pos;
      double big, eps, temp;
      /* compute infinity (maximum) norm of the column */
      big = 0.0;
      for (pos = 1; pos <= nnz; pos++)
      {
#ifdef GLP_DEBUG
         i = tcol_ind[pos];
         xassert(1 <= i && i <= m);
#endif
         temp = fabs(tcol_vec[tcol_ind[pos]]);
         if (big < temp) big = temp;
      }
      csa->tcol_max = big;
      /* determine absolute pivot tolerance */
      eps = tol_piv * (1.0 + 0.01 * big);
      /* move significant column components to front of the list */
      for (num = 0; num < nnz; )
      {  i = tcol_ind[nnz];
         if (fabs(tcol_vec[i]) < eps)
            nnz--;
         else
         {  num++;
            tcol_ind[nnz] = tcol_ind[num];
            tcol_ind[num] = i;
         }
      }
      csa->tcol_num = num;
      return;
}

/***********************************************************************
*  chuzr - choose basic variable (row of the simplex table)
*
*  This routine chooses basic variable xB[p], which reaches its bound
*  first on changing non-basic variable xN[q] in valid direction.
*
*  The parameter rtol is a relative tolerance used to relax bounds of
*  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 primal feasibility. */

static void chuzr(struct csa *csa, double rtol)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      char *type = csa->type;
      double *lb = csa->lb;
      double *ub = csa->ub;
      double *coef = csa->coef;
      int *head = csa->head;
      int phase = csa->phase;
      double *bbar = csa->bbar;
      double *cbar = csa->cbar;
      int q = csa->q;
      int *tcol_ind = csa->tcol_ind;
      double *tcol_vec = csa->tcol_vec;
      int tcol_num = csa->tcol_num;
      int i, i_stat, k, p, p_stat, pos;
      double alfa, big, delta, s, t, teta, tmax;
#ifdef GLP_DEBUG
      xassert(1 <= q && q <= n);
#endif
      /* s := - sign(d[q]), where d[q] is reduced cost of xN[q] */
#ifdef GLP_DEBUG
      xassert(cbar[q] != 0.0);
#endif
      s = (cbar[q] > 0.0 ? -1.0 : +1.0);
      /*** FIRST PASS ***/
      k = head[m+q]; /* x[k] = xN[q] */
#ifdef GLP_DEBUG
      xassert(1 <= k && k <= m+n);
#endif
      if (type[k] == GLP_DB)
      {  /* xN[q] has both lower and upper bounds */
         p = -1, p_stat = 0, teta = ub[k] - lb[k], big = 1.0;
      }
      else
      {  /* xN[q] has no opposite bound */
         p = 0, p_stat = 0, teta = DBL_MAX, big = 0.0;
      }
      /* walk through significant elements of the pivot column */
      for (pos = 1; pos <= tcol_num; pos++)
      {  i = tcol_ind[pos];
#ifdef GLP_DEBUG
         xassert(1 <= i && i <= m);
#endif
         k = head[i]; /* x[k] = xB[i] */
#ifdef GLP_DEBUG
         xassert(1 <= k && k <= m+n);
#endif
         alfa = s * tcol_vec[i];
#ifdef GLP_DEBUG
         xassert(alfa != 0.0);
#endif
         /* xB[i] = ... + alfa * xN[q] + ..., and due to s we need to
            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 */
               delta = rtol * (1.0 + kappa * fabs(lb[k]));
               t = ((lb[k] + delta) - 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 */
               delta = rtol * (1.0 + kappa * fabs(ub[k]));
               t = ((ub[k] + delta) - 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 */
               delta = rtol * (1.0 + kappa * fabs(ub[k]));
               t = ((ub[k] - delta) - 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 */
               delta = rtol * (1.0 + kappa * fabs(lb[k]));
               t = ((lb[k] - delta) - bbar[i]) / alfa;
               i_stat = GLP_NL;
            }
            else
            {  /* xB[i] is within its bounds and has no lower bound */
               continue;
            }
         }
         /* t is a change of xN[q], on which xB[i] reaches its bound
            (possibly relaxed); since the basic solution is assumed to
            be primal feasible (or pseudo feasible on phase I), t has
            to be non-negative by definition; however, it may happen
            that xB[i] slightly (i.e. within a tolerance) violates its
            bound, that leads to negative t; in the latter case, if
            xB[i] is chosen, negative t means that xN[q] changes in
            wrong direction; if pivot alfa[i,q] is close to zero, even
            small bound violation of xB[i] may lead to a large change
            of xN[q] in wrong direction; let, for example, xB[i] >= 0
            and in the current basis its value be -5e-9; let also xN[q]
            be on its zero bound and should increase; from the ratio
            test rule it follows that the pivot alfa[i,q] < 0; however,
            if alfa[i,q] is, say, -1e-9, the change of xN[q] in wrong
            direction is 5e-9 / (-1e-9) = -5, and using it for updating
            values of other basic variables will give absolutely wrong
            results; therefore, if t is negative, we should replace it
            by exact zero assuming that xB[i] is exactly on its bound,
            and the 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))
            p = i, p_stat = i_stat, 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 xN[q] reaches its opposite bound or if no basic variable
         has been chosen on the first pass */
      if (p <= 0) goto done;
      /* if xB[p] is a blocking variable, i.e. if it prevents xN[q]
         from any change */
      if (teta == 0.0) goto done;
      /*** SECOND PASS ***/
      /* here tmax is a maximal change of xN[q], on which the solution
         remains primal feasible (or pseudo feasible on phase I) within
         a tolerance */
#if 0
      tmax = (1.0 + 10.0 * DBL_EPSILON) * teta;
#else
      tmax = teta;
#endif
      /* nothing is chosen so far */
      p = 0, p_stat = 0, teta = DBL_MAX, big = 0.0;
      /* walk through significant elements of the pivot column */
      for (pos = 1; pos <= tcol_num; pos++)
      {  i = tcol_ind[pos];
#ifdef GLP_DEBUG
         xassert(1 <= i && i <= m);
#endif
         k = head[i]; /* x[k] = xB[i] */
#ifdef GLP_DEBUG
         xassert(1 <= k && k <= m+n);
#endif
         alfa = s * tcol_vec[i];
#ifdef GLP_DEBUG
         xassert(alfa != 0.0);
#endif
         /* xB[i] = ... + alfa * xN[q] + ..., and due to s we need to