glpspx02.c

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

/* glpspx02.c (dual simplex method) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[e-mail]>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glpspx.h"

struct csa
{     /* common storage area */
      /*--------------------------------------------------------------*/
      /* LP data */
      int m;
      /* number of rows (auxiliary variables), m > 0 */
      int n;
      /* number of columns (structural variables), n > 0 */
      char *type; /* char type[1+m+n]; */
      /* type[0] is not used;
         type[k], 1 <= k <= m+n, is the type of variable x[k]:
         GLP_FR - free variable
         GLP_LO - variable with lower bound
         GLP_UP - variable with upper bound
         GLP_DB - double-bounded variable
         GLP_FX - fixed variable */
      double *lb; /* double lb[1+m+n]; */
      /* lb[0] is not used;
         lb[k], 1 <= k <= m+n, is an lower bound of variable x[k];
         if x[k] has no lower bound, lb[k] is zero */
      double *ub; /* double ub[1+m+n]; */
      /* ub[0] is not used;
         ub[k], 1 <= k <= m+n, is an upper bound of variable x[k];
         if x[k] has no upper bound, ub[k] is zero;
         if x[k] is of fixed type, ub[k] is the same as lb[k] */
      double *coef; /* double coef[1+m+n]; */
      /* coef[0] is not used;
         coef[k], 1 <= k <= m+n, is an objective coefficient at
         variable x[k] */
      /*--------------------------------------------------------------*/
      /* original bounds of variables */
      char *orig_type; /* char orig_type[1+m+n]; */
      double *orig_lb; /* double orig_lb[1+m+n]; */
      double *orig_ub; /* double orig_ub[1+m+n]; */
      /*--------------------------------------------------------------*/
      /* original objective function */
      double *obj; /* double obj[1+n]; */
      /* obj[0] is a constant term of the original objective function;
         obj[j], 1 <= j <= n, is an original objective coefficient at
         structural variable x[m+j] */
      double zeta;
      /* factor used to scale original objective coefficients; its
         sign defines original optimization direction: zeta > 0 means
         minimization, zeta < 0 means maximization */
      /*--------------------------------------------------------------*/
      /* constraint matrix A; it has m rows and n columns and is stored
         by columns */
      int *A_ptr; /* int A_ptr[1+n+1]; */
      /* A_ptr[0] is not used;
         A_ptr[j], 1 <= j <= n, is starting position of j-th column in
         arrays A_ind and A_val; note that A_ptr[1] is always 1;
         A_ptr[n+1] indicates the position after the last element in
         arrays A_ind and A_val */
      int *A_ind; /* int A_ind[A_ptr[n+1]]; */
      /* row indices */
      double *A_val; /* double A_val[A_ptr[n+1]]; */
      /* non-zero element values */
      /*--------------------------------------------------------------*/
      /* basis header */
      int *head; /* int head[1+m+n]; */
      /* head[0] is not used;
         head[i], 1 <= i <= m, is the ordinal number of basic variable
         xB[i]; head[i] = k means that xB[i] = x[k] and i-th column of
         matrix B is k-th column of matrix (I|-A);
         head[m+j], 1 <= j <= n, is the ordinal number of non-basic
         variable xN[j]; head[m+j] = k means that xN[j] = x[k] and j-th
         column of matrix N is k-th column of matrix (I|-A) */
      char *stat; /* char stat[1+n]; */
      /* stat[0] is not used;
         stat[j], 1 <= j <= n, is the status of non-basic variable
         xN[j], which defines its active bound:
         GLP_NL - lower bound is active
         GLP_NU - upper bound is active
         GLP_NF - free variable
         GLP_NS - fixed variable */
      /*--------------------------------------------------------------*/
      /* matrix B is the basis matrix; it is composed from columns of
         the augmented constraint matrix (I|-A) corresponding to basic
         variables and stored in a factorized (invertable) form */
      int valid;
      /* factorization is valid only if this flag is set */
      BFD *bfd; /* BFD bfd[1:m,1:m]; */
      /* factorized (invertable) form of the basis matrix */
      /*--------------------------------------------------------------*/
      /* matrix N is a matrix composed from columns of the augmented
         constraint matrix (I|-A) corresponding to non-basic variables
         except fixed ones; it is stored by rows and changes every time
         the basis changes */
      int *N_ptr; /* int N_ptr[1+m+1]; */
      /* N_ptr[0] is not used;
         N_ptr[i], 1 <= i <= m, is starting position of i-th row in
         arrays N_ind and N_val; note that N_ptr[1] is always 1;
         N_ptr[m+1] indicates the position after the last element in
         arrays N_ind and N_val */
      int *N_len; /* int N_len[1+m]; */
      /* N_len[0] is not used;
         N_len[i], 1 <= i <= m, is length of i-th row (0 to n) */
      int *N_ind; /* int N_ind[N_ptr[m+1]]; */
      /* column indices */
      double *N_val; /* double N_val[N_ptr[m+1]]; */
      /* non-zero element values */
      /*--------------------------------------------------------------*/
      /* working parameters */
      int phase;
      /* search phase:
         0 - not determined yet
         1 - search for dual feasible solution
         2 - search for optimal solution */
      xlong_t tm_beg;
      /* time value at the beginning of the search */
      int it_beg;
      /* simplex iteration count at the beginning of the search */
      int it_cnt;
      /* simplex iteration count; it increases by one every time the
         basis changes */
      int it_dpy;
      /* simplex iteration count at the most recent display output */
      /*--------------------------------------------------------------*/
      /* basic solution components */
      double *bbar; /* double bbar[1+m]; */
      /* bbar[0] is not used on phase I; on phase II it is the current
         value of the original objective function;
         bbar[i], 1 <= i <= m, is primal value of basic variable xB[i]
         (if xB[i] is free, its primal value is not updated) */
      double *cbar; /* double cbar[1+n]; */
      /* cbar[0] is not used;
         cbar[j], 1 <= j <= n, is reduced cost of non-basic variable
         xN[j] (if xN[j] is fixed, its reduced cost is not updated) */
      /*--------------------------------------------------------------*/
      /* the following pricing technique options may be used:
         GLP_PT_STD - standard ("textbook") pricing;
         GLP_PT_PSE - projected steepest edge;
         GLP_PT_DVX - Devex pricing (not implemented yet);
         in case of GLP_PT_STD the reference space is not used, and all
         steepest edge coefficients are set to 1 */
      int refct;
      /* this count is set to an initial value when the reference space
         is defined and decreases by one every time the basis changes;
         once this count reaches zero, the reference space is redefined
         again */
      char *refsp; /* char refsp[1+m+n]; */
      /* refsp[0] is not used;
         refsp[k], 1 <= k <= m+n, is the flag which means that variable
         x[k] belongs to the current reference space */
      double *gamma; /* double gamma[1+m]; */
      /* gamma[0] is not used;
         gamma[i], 1 <= i <= n, is the steepest edge coefficient for
         basic variable xB[i]; if xB[i] is free, gamma[i] is not used
         and just set to 1 */
      /*--------------------------------------------------------------*/
      /* basic variable xB[p] chosen to leave the basis */
      int p;
      /* index of the basic variable xB[p] chosen, 1 <= p <= m;
         if the set of eligible basic variables is empty (i.e. if the
         current basic solution is primal feasible within a tolerance)
         and thus no variable has been chosen, p is set to 0 */
      double delta;
      /* change of xB[p] in the adjacent basis;
         delta > 0 means that xB[p] violates its lower bound and will
         increase to achieve it in the adjacent basis;
         delta < 0 means that xB[p] violates its upper bound and will
         decrease to achieve it in the adjacent basis */
      /*--------------------------------------------------------------*/
      /* pivot row of the simplex table corresponding to basic variable
         xB[p] chosen is the following vector:
            T' * e[p] = - N' * inv(B') * e[p] = - N' * rho,
         where B' is a matrix transposed to the current basis matrix,
         N' is a matrix, whose rows are columns of the matrix (I|-A)
         corresponding to non-basic non-fixed variables */
      int trow_nnz;
      /* number of non-zero components, 0 <= nnz <= n */
      int *trow_ind; /* int trow_ind[1+n]; */
      /* trow_ind[0] is not used;
         trow_ind[t], 1 <= t <= nnz, is an index of non-zero component,
         i.e. trow_ind[t] = j means that trow_vec[j] != 0 */
      double *trow_vec; /* int trow_vec[1+n]; */
      /* trow_vec[0] is not used;
         trow_vec[j], 1 <= j <= n, is a numeric value of j-th component
         of the row */
      double trow_max;
      /* infinity (maximum) norm of the row (max |trow_vec[j]|) */
      int trow_num;
      /* number of significant non-zero components, which means that:
         |trow_vec[j]| >= eps for j in trow_ind[1,...,num],
         |tcol_vec[j]| <  eps for j in trow_ind[num+1,...,nnz],
         where eps is a pivot tolerance */
      /*--------------------------------------------------------------*/
      /* non-basic variable xN[q] chosen to enter the basis */
      int q;
      /* index of the non-basic variable xN[q] chosen, 1 <= q <= n;
         if no variable has been chosen, q is set to 0 */
      double new_dq;
      /* reduced cost of xN[q] in the adjacent basis (it is the change
         of lambdaB[p]) */
      /*--------------------------------------------------------------*/
      /* pivot column of the simplex table corresponding to non-basic
         variable xN[q] chosen is the following vector:
            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 xN[q] */
      int tcol_nnz;
      /* number of non-zero components, 0 <= nnz <= m */
      int *tcol_ind; /* int tcol_ind[1+m]; */
      /* tcol_ind[0] is not used;
         tcol_ind[t], 1 <= t <= nnz, is an index of non-zero component,
         i.e. tcol_ind[t] = i means that tcol_vec[i] != 0 */
      double *tcol_vec; /* double tcol_vec[1+m]; */
      /* tcol_vec[0] is not used;
         tcol_vec[i], 1 <= i <= m, is a numeric value of i-th component
         of the column */
      /*--------------------------------------------------------------*/
      /* working arrays */
      double *work1; /* double work1[1+m]; */
      double *work2; /* double work2[1+m]; */
      double *work3; /* double work3[1+m]; */
      double *work4; /* double work4[1+m]; */
};

static const double kappa = 0.10;

/***********************************************************************
*  alloc_csa - allocate common storage area
*
*  This routine allocates all arrays in the common storage area (CSA)
*  and returns a pointer to the CSA. */

static struct csa *alloc_csa(glp_prob *lp)
{     struct csa *csa;
      int m = lp->m;
      int n = lp->n;
      int nnz = lp->nnz;
      csa = xmalloc(sizeof(struct csa));
      xassert(m > 0 && n > 0);
      csa->m = m;
      csa->n = n;
      csa->type = xcalloc(1+m+n, sizeof(char));
      csa->lb = xcalloc(1+m+n, sizeof(double));
      csa->ub = xcalloc(1+m+n, sizeof(double));
      csa->coef = xcalloc(1+m+n, sizeof(double));
      csa->orig_type = xcalloc(1+m+n, sizeof(char));
      csa->orig_lb = xcalloc(1+m+n, sizeof(double));
      csa->orig_ub = xcalloc(1+m+n, sizeof(double));
      csa->obj = xcalloc(1+n, sizeof(double));
      csa->A_ptr = xcalloc(1+n+1, sizeof(int));
      csa->A_ind = xcalloc(1+nnz, sizeof(int));
      csa->A_val = xcalloc(1+nnz, sizeof(double));
      csa->head = xcalloc(1+m+n, sizeof(int));
      csa->stat = xcalloc(1+n, sizeof(char));
      csa->N_ptr = xcalloc(1+m+1, sizeof(int));
      csa->N_len = xcalloc(1+m, sizeof(int));
      csa->N_ind = NULL; /* will be allocated later */
      csa->N_val = NULL; /* will be allocated later */
      csa->bbar = xcalloc(1+m, sizeof(double));
      csa->cbar = xcalloc(1+n, sizeof(double));
      csa->refsp = xcalloc(1+m+n, sizeof(char));
      csa->gamma = xcalloc(1+m, sizeof(double));
      csa->trow_ind = xcalloc(1+n, sizeof(int));
      csa->trow_vec = xcalloc(1+n, sizeof(double));
      csa->tcol_ind = xcalloc(1+m, sizeof(int));
      csa->tcol_vec = xcalloc(1+m, sizeof(double));
      csa->work1 = xcalloc(1+m, sizeof(double));
      csa->work2 = xcalloc(1+m, sizeof(double));
      csa->work3 = xcalloc(1+m, sizeof(double));
      csa->work4 = xcalloc(1+m, sizeof(double));
      return csa;
}

/***********************************************************************
*  init_csa - initialize common storage area
*
*  This routine initializes all data structures in the common storage
*  area (CSA). */

static void alloc_N(struct csa *csa);
static void build_N(struct csa *csa);

static void init_csa(struct csa *csa, glp_prob *lp)
{     int m = csa->m;
      int n = csa->n;
      char *type = csa->type;
      double *lb = csa->lb;
      double *ub = csa->ub;
      double *coef = csa->coef;
      char *orig_type = csa->orig_type;
      double *orig_lb = csa->orig_lb;
      double *orig_ub = csa->orig_ub;
      double *obj = csa->obj;
      int *A_ptr = csa->A_ptr;
      int *A_ind = csa->A_ind;
      double *A_val = csa->A_val;
      int *head = csa->head;
      char *stat = csa->stat;
      char *refsp = csa->refsp;
      double *gamma = csa->gamma;
      int i, j, k, loc;
      double cmax;
      /* auxiliary variables */
      for (i = 1; i <= m; i++)
      {  GLPROW *row = lp->row[i];
         type[i] = (char)row->type;
         lb[i] = row->lb * row->rii;
         ub[i] = row->ub * row->rii;
         coef[i] = 0.0;
      }
      /* structural variables */
      for (j = 1; j <= n; j++)
      {  GLPCOL *col = lp->col[j];
         type[m+j] = (char)col->type;
         lb[m+j] = col->lb / col->sjj;
         ub[m+j] = col->ub / col->sjj;
         coef[m+j] = col->coef * col->sjj;
      }
      /* original bounds of variables */
      memcpy(&orig_type[1], &type[1], (m+n) * sizeof(char));
      memcpy(&orig_lb[1], &lb[1], (m+n) * sizeof(double));
      memcpy(&orig_ub[1], &ub[1], (m+n) * sizeof(double));
      /* original objective function */
      obj[0] = lp->c0;
      memcpy(&obj[1], &coef[m+1], n * sizeof(double));
      /* factor used to scale original objective coefficients */
      cmax = 0.0;
      for (j = 1; j <= n; j++)
         if (cmax < fabs(obj[j])) cmax = fabs(obj[j]);
      if (cmax == 0.0) cmax = 1.0;
      switch (lp->dir)
      {  case GLP_MIN:
            csa->zeta = + 1.0 / cmax;
            break;
         case GLP_MAX:
            csa->zeta = - 1.0 / cmax;
            break;
         default:
            xassert(lp != lp);
      }
#if 1
      if (fabs(csa->zeta) < 1.0) csa->zeta *= 1000.0;
#endif
      /* scale working objective coefficients */
      for (j = 1; j <= n; j++) coef[m+j] *= csa->zeta;
      /* matrix A (by columns) */
      loc = 1;
      for (j = 1; j <= n; j++)
      {  GLPAIJ *aij;
         A_ptr[j] = loc;
         for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
         {  A_ind[loc] = aij->row->i;
            A_val[loc] = aij->row->rii * aij->val * aij->col->sjj;
            loc++;
         }
      }
      A_ptr[n+1] = loc;
      xassert(loc == lp->nnz+1);
      /* basis header */
      xassert(lp->valid);
      memcpy(&head[1], &lp->head[1], m * sizeof(int));
      k = 0;
      for (i = 1; i <= m; i++)
      {  GLPROW *row = lp->row[i];
         if (row->stat != GLP_BS)
         {  k++;
            xassert(k <= n);
            head[m+k] = i;
            stat[k] = (char)row->stat;
         }
      }