glpspx02.c

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

         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] / 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] / alfa;
            }
            else
            {  /* lambdaN[j] has no upper bound */
               continue;
            }
         }
         /* (see comments for the first pass) */
         if (t < 0.0) t = 0.0;
         /* t is a change of lambdaB[p], on which lambdaN[j] reaches
            its zero (lower or upper) bound; if t <= tmax, all reduced
            costs can violate their zero bounds only within relaxation
            tolerance rtol, so we can choose non-basic variable having
            largest influence coefficient to avoid possible numerical
            instability */
         if (t <= tmax && big < fabs(alfa))
            q = j, teta = t, big = fabs(alfa);
      }
      /* something must be chosen on the second pass */
      xassert(q != 0);
done: /* store the index of non-basic variable xN[q] chosen */
      csa->q = q;
      /* store reduced cost of xN[q] in the adjacent basis */
      csa->new_dq = s * teta;
      return;
}

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

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

/***********************************************************************
*  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 trow_nnz = csa->trow_nnz;
      int *trow_ind = csa->trow_ind;
      double *trow_vec = csa->trow_vec;
      int q = csa->q;
      double new_dq = csa->new_dq;
      int j, pos;
#ifdef GLP_DEBUG
      xassert(1 <= q && q <= n);
#endif
      /* set new reduced cost of xN[q] */
      cbar[q] = new_dq;
      /* update reduced costs of other non-basic variables */
      if (new_dq == 0.0) goto done;
      for (pos = 1; pos <= trow_nnz; pos++)
      {  j = trow_ind[pos];
#ifdef GLP_DEBUG
         xassert(1 <= j && j <= n);
#endif
         if (j != q)
            cbar[j] -= trow_vec[j] * new_dq;
      }
done: 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 p = csa->p;
      double delta = csa->delta;
      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 teta;
#ifdef GLP_DEBUG
      xassert(1 <= p && p <= m);
      xassert(1 <= q && q <= n);
#endif
      /* determine the change of xN[q] in the adjacent basis */
#ifdef GLP_DEBUG
      xassert(tcol_vec[p] != 0.0);
#endif
      teta = delta / tcol_vec[p];
      /* set new primal value of xN[q] */
      bbar[p] = get_xN(csa, q) + teta;
      /* update primal values of other basic variables */
      if (teta == 0.0) goto done;
      for (pos = 1; pos <= tcol_nnz; pos++)
      {  i = tcol_ind[pos];
#ifdef GLP_DEBUG
         xassert(1 <= i && i <= m);
#endif
         if (i != p)
            bbar[i] += tcol_vec[i] * teta;
      }
done: 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 *head = csa->head;
      char *refsp = csa->refsp;
      double *gamma = csa->gamma;
      int p = csa->p;
      int trow_nnz = csa->trow_nnz;
      int *trow_ind = csa->trow_ind;
      double *trow_vec = csa->trow_vec;
      int q = csa->q;
      int tcol_nnz = csa->tcol_nnz;
      int *tcol_ind = csa->tcol_ind;
      double *tcol_vec = csa->tcol_vec;
      double *u = csa->work3;
      int i, j, k,pos;
      double gamma_p, eta_p, pivot, 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[p] for the current basis more accurately and
         compute auxiliary vector u */
#ifdef GLP_DEBUG
      xassert(type[head[p]] != GLP_FR);
#endif
      gamma_p = eta_p = (refsp[head[p]] ? 1.0 : 0.0);
      for (i = 1; i <= m; i++) u[i] = 0.0;
      for (pos = 1; pos <= trow_nnz; pos++)
      {  j = trow_ind[pos];
#ifdef GLP_DEBUG
         xassert(1 <= j && j <= n);
#endif
         k = head[m+j]; /* x[k] = xN[j] */
#ifdef GLP_DEBUG
         xassert(1 <= k && k <= m+n);
         xassert(type[k] != GLP_FX);
#endif
         if (!refsp[k]) continue;
         t = trow_vec[j];
         gamma_p += t * t;
         /* u := u + N[j] * delta[j] * trow[j] */
         if (k <= m)
         {  /* N[k] = k-j stolbec submatrix I */
            u[k] += t;
         }
         else
         {  /* N[k] = k-m-k stolbec (-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++)
               u[A_ind[ptr]] -= t * A_val[ptr];
         }
      }
      xassert(csa->valid);
      bfd_ftran(csa->bfd, u);
      /* update gamma[i] for other basic variables (except xB[p] and
         free variables) */
      pivot = tcol_vec[p];
#ifdef GLP_DEBUG
      xassert(pivot != 0.0);
#endif
      for (pos = 1; pos <= tcol_nnz; pos++)
      {  i = tcol_ind[pos];
#ifdef GLP_DEBUG
         xassert(1 <= i && i <= m);
#endif
         k = head[i];
#ifdef GLP_DEBUG
         xassert(1 <= k && k <= m+n);
#endif
         /* skip xB[p] */
         if (i == p) continue;
         /* skip free basic variable */
         if (type[head[i]] == GLP_FR)
         {
#ifdef GLP_DEBUG
            xassert(gamma[i] == 1.0);
#endif
            continue;
         }
         /* compute gamma[i] for the adjacent basis */
         t = tcol_vec[i] / pivot;
         t1 = gamma[i] + t * t * gamma_p + 2.0 * t * u[i];
         t2 = (refsp[k] ? 1.0 : 0.0) + eta_p * t * t;
         gamma[i] = (t1 >= t2 ? t1 : t2);
         /* (though gamma[i] can be exact zero, because the reference
            space does not include non-basic fixed variables) */
         if (gamma[i] < DBL_EPSILON) gamma[i] = DBL_EPSILON;
      }
      /* compute gamma[p] for the adjacent basis */
      if (type[head[m+q]] == GLP_FR)
         gamma[p] = 1.0;
      else
      {  gamma[p] = gamma_p / (pivot * pivot);
         if (gamma[p] < DBL_EPSILON) gamma[p] = DBL_EPSILON;
      }
      /* if xB[p], which becomes xN[q] in the adjacent basis, is fixed
         and belongs to the reference space, remove it from there, and
         change all gamma's appropriately */
      k = head[p];
      if (type[k] == GLP_FX && refsp[k])
      {  refsp[k] = 0;
         for (pos = 1; pos <= tcol_nnz; pos++)
         {  i = tcol_ind[pos];
            if (i == p)
            {  if (type[head[m+q]] == GLP_FR) continue;
               t = 1.0 / tcol_vec[p];
            }
            else
            {  if (type[head[i]] == GLP_FR) continue;
               t = tcol_vec[i] / tcol_vec[p];
            }
            gamma[i] -= t * t;
            if (gamma[i] < DBL_EPSILON) gamma[i] = DBL_EPSILON;
         }
      }
      return;
}

#if 1 /* copied from primal */
/***********************************************************************
*  err_in_bbar - compute maximal relative error in primal solut