glpspx01.c

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

            /* find element in i-th row of N */
            head = N_ptr[i];
            for (pos = head; N_ind[pos] != j; pos++) /* nop */;
            /* and remove it from the row list */
            tail = head + (--N_len[i]);
#ifdef GLP_DEBUG
            xassert(pos <= tail);
#endif
            N_ind[pos] = N_ind[tail];
            N_val[pos] = N_val[tail];
         }
      }
      return;
}

/***********************************************************************
*  build_N - build matrix N for current basis
*
*  This routine builds matrix N for the current basis from columns
*  of the augmented constraint matrix (I|-A) corresponding to non-basic
*  non-fixed variables. */

static void build_N(struct csa *csa)
{     int m = csa->m;
      int n = csa->n;
      int *head = csa->head;
      char *stat = csa->stat;
      int *N_len = csa->N_len;
      int j, k;
      /* N := empty matrix */
      memset(&N_len[1], 0, m * sizeof(int));
      /* go through non-basic columns of matrix (I|-A) */
      for (j = 1; j <= n; j++)
      {  if (stat[j] != GLP_NS)
         {  /* xN[j] is non-fixed; add j-th column to matrix N which is
               k-th column of matrix (I|-A) */
            k = head[m+j]; /* x[k] = xN[j] */
#ifdef GLP_DEBUG
            xassert(1 <= k && k <= m+n);
#endif
            add_N_col(csa, j, k);
         }
      }
      return;
}

/***********************************************************************
*  get_xN - determine current value of non-basic variable xN[j]
*
*  This routine returns the current value of non-basic variable xN[j],
*  which is a value of its active bound. */

static double get_xN(struct csa *csa, int j)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      double *lb = csa->lb;
      double *ub = csa->ub;
      int *head = csa->head;
      char *stat = csa->stat;
      int k;
      double xN;
#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);
#endif
      switch (stat[j])
      {  case GLP_NL:
            /* x[k] is on its lower bound */
            xN = lb[k]; break;
         case GLP_NU:
            /* x[k] is on its upper bound */
            xN = ub[k]; break;
         case GLP_NF:
            /* x[k] is free non-basic variable */
            xN = 0.0; break;
         case GLP_NS:
            /* x[k] is fixed non-basic variable */
            xN = lb[k]; break;
         default:
            xassert(stat != stat);
      }
      return xN;
}

/***********************************************************************
*  eval_beta - compute primal values of basic variables
*
*  This routine computes current primal values of all basic variables:
*
*     beta = - inv(B) * N * xN,
*
*  where B is the current basis matrix, N is a matrix built of columns
*  of matrix (I|-A) corresponding to non-basic variables, and xN is the
*  vector of current values of non-basic variables. */

static void eval_beta(struct csa *csa, double beta[])
{     int m = csa->m;
      int n = csa->n;
      int *A_ptr = csa->A_ptr;
      int *A_ind = csa->A_ind;
      double *A_val = csa->A_val;
      int *head = csa->head;
      double *h = csa->work2;
      int i, j, k, beg, end, ptr;
      double xN;
      /* compute the right-hand side vector:
         h := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n],
         where N[1], ..., N[n] are columns of matrix N */
      for (i = 1; i <= m; i++)
         h[i] = 0.0;
      for (j = 1; j <= n; j++)
      {  k = head[m+j]; /* x[k] = xN[j] */
#ifdef GLP_DEBUG
         xassert(1 <= k && k <= m+n);
#endif
         /* determine current value of xN[j] */
         xN = get_xN(csa, j);
         if (xN == 0.0) continue;
         if (k <= m)
         {  /* N[j] is k-th column of submatrix I */
            h[k] -= xN;
         }
         else
         {  /* N[j] is (k-m)-th column of submatrix (-A) */
            beg = A_ptr[k-m];
            end = A_ptr[k-m+1];
            for (ptr = beg; ptr < end; ptr++)
               h[A_ind[ptr]] += xN * A_val[ptr];
         }
      }
      /* solve system B * beta = h */
      memcpy(&beta[1], &h[1], m * sizeof(double));
      xassert(csa->valid);
      bfd_ftran(csa->bfd, beta);
      /* and refine the solution */
      refine_ftran(csa, h, beta);
      return;
}

/***********************************************************************
*  eval_pi - compute vector of simplex multipliers
*
*  This routine computes the vector of current simplex multipliers:
*
*     pi = inv(B') * cB,
*
*  where B' is a matrix transposed to the current basis matrix, cB is
*  a subvector of objective coefficients at basic variables. */

static void eval_pi(struct csa *csa, double pi[])
{     int m = csa->m;
      double *c = csa->coef;
      int *head = csa->head;
      double *cB = csa->work2;
      int i;
      /* construct the right-hand side vector cB */
      for (i = 1; i <= m; i++)
         cB[i] = c[head[i]];
      /* solve system B'* pi = cB */
      memcpy(&pi[1], &cB[1], m * sizeof(double));
      xassert(csa->valid);
      bfd_btran(csa->bfd, pi);
      /* and refine the solution */
      refine_btran(csa, cB, pi);
      return;
}

/***********************************************************************
*  eval_cost - compute reduced cost of non-basic variable xN[j]
*
*  This routine computes the current reduced cost of non-basic variable
*  xN[j]:
*
*     d[j] = cN[j] - N'[j] * pi,
*
*  where cN[j] is the objective coefficient at variable xN[j], N[j] is
*  a column of the augmented constraint matrix (I|-A) corresponding to
*  xN[j], pi is the vector of simplex multipliers. */

static double eval_cost(struct csa *csa, double pi[], int j)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      double *coef = csa->coef;
      int *head = csa->head;
      int k;
      double dj;
#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);
#endif
      dj = coef[k];
      if (k <= m)
      {  /* N[j] is k-th column of submatrix I */
         dj -= pi[k];
      }
      else
      {  /* N[j] 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++)
            dj += A_val[ptr] * pi[A_ind[ptr]];
      }
      return dj;
}

/***********************************************************************
*  eval_bbar - compute and store primal values of basic variables
*
*  This routine computes primal values of all basic variables and then
*  stores them in the solution array. */

static void eval_bbar(struct csa *csa)
{     eval_beta(csa, csa->bbar);
      return;
}

/***********************************************************************
*  eval_cbar - compute and store reduced costs of non-basic variables
*
*  This routine computes reduced costs of all non-basic variables and
*  then stores them in the solution array. */

static void eval_cbar(struct csa *csa)
{
#ifdef GLP_DEBUG
      int m = csa->m;
#endif
      int n = csa->n;
#ifdef GLP_DEBUG
      int *head = csa->head;
#endif
      double *cbar = csa->cbar;
      double *pi = csa->work3;
      int j;
#ifdef GLP_DEBUG
      int k;
#endif
      /* compute simplex multipliers */
      eval_pi(csa, pi);
      /* compute and store reduced costs */
      for (j = 1; j <= n; j++)
      {
#ifdef GLP_DEBUG
         k = head[m+j]; /* x[k] = xN[j] */
         xassert(1 <= k && k <= m+n);
#endif
         cbar[j] = eval_cost(csa, pi, j);
      }
      return;
}

/***********************************************************************
*  reset_refsp - reset the reference space
*
*  This routine resets (redefines) the reference space used in the
*  projected steepest edge pricing algorithm. */

static void reset_refsp(struct csa *csa)
{     int m = csa->m;
      int n = csa->n;
      int *head = csa->head;
      char *refsp = csa->refsp;
      double *gamma = csa->gamma;
      int j, k;
      xassert(csa->refct == 0);
      csa->refct = 1000;
      memset(&refsp[1], 0, (m+n) * sizeof(char));
      for (j = 1; j <= n; j++)
      {  k = head[m+j]; /* x[k] = xN[j] */
         refsp[k] = 1;
         gamma[j] = 1.0;
      }
      return;
}

/***********************************************************************
*  eval_gamma - compute steepest edge coefficient
*
*  This routine computes the steepest edge coefficient for non-basic
*  variable xN[j] using its direct definition:
*
*     gamma[j] = delta[j] +  sum   alfa[i,j]^2,
*                           i in R
*
*  where delta[j] = 1, if xN[j] is in the current reference space,
*  and 0 otherwise; R is a set of basic variables xB[i], which are in
*  the current reference space; alfa[i,j] are elements of the current
*  simplex table.
*
*  NOTE: The routine is intended only for debugginig purposes. */

static double eval_gamma(struct csa *csa, int j)
{     int m = csa->m;
#ifdef GLP_DEBUG
      int n = csa->n;
#endif
      int *head = csa->head;
      char *refsp = csa->refsp;
      double *alfa = csa->work3;
      double *h = csa->work3;
      int i, k;
      double gamma;
#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);
#endif
      /* construct the right-hand side vector h = - N[j] */
      for (i = 1; i <= m; i++)
         h[i] = 0.0;
      if (k <= m)
      {  /* N[j] is k-th column of submatrix I */
         h[k] = -1.0;
      }
      else
      {  /* N[j] 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 * alfa = h */
      xassert(csa->valid);
      bfd_ftran(csa->bfd, alfa);
      /* compute gamma */
      gamma = (refsp[k] ? 1.0 : 0.0);
      for (i = 1; i <= m; i++)
      {  k = head[i];
#ifdef GLP_DEBUG
         xassert(1 <= k && k <= m+n);
#endif
         if (refsp[k]) gamma += alfa[i] * alfa[i];
      }
      return gamma;
}

/***********************************************************************
*  chuzc - choose non-basic variable (column of the simplex table)
*
*  This routine chooses non-basic variable xN[q], which has largest
*  weighted reduced cost:
*
*     |d[q]| / sqrt(gamma[q]) = max  |d[j]| / sqrt(gamma[j]),
*                              j in J
*
*  where J is a subset of eligible non-basic variables xN[j], d[j] is
*  reduced cost of xN[j], gamma[j] is the steepest edge coefficient.
*
*  The working objective function is always minimized, so the sign of
*  d[q] determines direction, in which xN[q] has to change:
*
*     if d[q] < 0, xN[q] has to increase;
*
*     if d[q] > 0, xN[q] has to decrease.
*
*  If |d[j]| <= tol_dj, where tol_dj is a specified tolerance, xN[j]
*  is not included in J and therefore ignored. (It is assumed that the
*  working objective row is appropriately scaled, i.e. max|c[k]| = 1.)
*
*  If J is empty and no variable has been chosen, q is set to 0. */

static void chuzc(struct csa *csa, double tol_dj)
{     int n = csa->n;
      char *stat = csa->stat;