alpar@1: /* glpspx01.c (primal simplex method) */
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1: *
alpar@1: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1: *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1: *  E-mail: <mao@gnu.org>.
alpar@1: *
alpar@1: *  GLPK is free software: you can redistribute it and/or modify it
alpar@1: *  under the terms of the GNU General Public License as published by
alpar@1: *  the Free Software Foundation, either version 3 of the License, or
alpar@1: *  (at your option) any later version.
alpar@1: *
alpar@1: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1: *  License for more details.
alpar@1: *
alpar@1: *  You should have received a copy of the GNU General Public License
alpar@1: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1: ***********************************************************************/
alpar@1: 
alpar@1: #include "glpspx.h"
alpar@1: 
alpar@1: struct csa
alpar@1: {     /* common storage area */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* LP data */
alpar@1:       int m;
alpar@1:       /* number of rows (auxiliary variables), m > 0 */
alpar@1:       int n;
alpar@1:       /* number of columns (structural variables), n > 0 */
alpar@1:       char *type; /* char type[1+m+n]; */
alpar@1:       /* type[0] is not used;
alpar@1:          type[k], 1 <= k <= m+n, is the type of variable x[k]:
alpar@1:          GLP_FR - free variable
alpar@1:          GLP_LO - variable with lower bound
alpar@1:          GLP_UP - variable with upper bound
alpar@1:          GLP_DB - double-bounded variable
alpar@1:          GLP_FX - fixed variable */
alpar@1:       double *lb; /* double lb[1+m+n]; */
alpar@1:       /* lb[0] is not used;
alpar@1:          lb[k], 1 <= k <= m+n, is an lower bound of variable x[k];
alpar@1:          if x[k] has no lower bound, lb[k] is zero */
alpar@1:       double *ub; /* double ub[1+m+n]; */
alpar@1:       /* ub[0] is not used;
alpar@1:          ub[k], 1 <= k <= m+n, is an upper bound of variable x[k];
alpar@1:          if x[k] has no upper bound, ub[k] is zero;
alpar@1:          if x[k] is of fixed type, ub[k] is the same as lb[k] */
alpar@1:       double *coef; /* double coef[1+m+n]; */
alpar@1:       /* coef[0] is not used;
alpar@1:          coef[k], 1 <= k <= m+n, is an objective coefficient at
alpar@1:          variable x[k] (note that on phase I auxiliary variables also
alpar@1:          may have non-zero objective coefficients) */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* original objective function */
alpar@1:       double *obj; /* double obj[1+n]; */
alpar@1:       /* obj[0] is a constant term of the original objective function;
alpar@1:          obj[j], 1 <= j <= n, is an original objective coefficient at
alpar@1:          structural variable x[m+j] */
alpar@1:       double zeta;
alpar@1:       /* factor used to scale original objective coefficients; its
alpar@1:          sign defines original optimization direction: zeta > 0 means
alpar@1:          minimization, zeta < 0 means maximization */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* constraint matrix A; it has m rows and n columns and is stored
alpar@1:          by columns */
alpar@1:       int *A_ptr; /* int A_ptr[1+n+1]; */
alpar@1:       /* A_ptr[0] is not used;
alpar@1:          A_ptr[j], 1 <= j <= n, is starting position of j-th column in
alpar@1:          arrays A_ind and A_val; note that A_ptr[1] is always 1;
alpar@1:          A_ptr[n+1] indicates the position after the last element in
alpar@1:          arrays A_ind and A_val */
alpar@1:       int *A_ind; /* int A_ind[A_ptr[n+1]]; */
alpar@1:       /* row indices */
alpar@1:       double *A_val; /* double A_val[A_ptr[n+1]]; */
alpar@1:       /* non-zero element values */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* basis header */
alpar@1:       int *head; /* int head[1+m+n]; */
alpar@1:       /* head[0] is not used;
alpar@1:          head[i], 1 <= i <= m, is the ordinal number of basic variable
alpar@1:          xB[i]; head[i] = k means that xB[i] = x[k] and i-th column of
alpar@1:          matrix B is k-th column of matrix (I|-A);
alpar@1:          head[m+j], 1 <= j <= n, is the ordinal number of non-basic
alpar@1:          variable xN[j]; head[m+j] = k means that xN[j] = x[k] and j-th
alpar@1:          column of matrix N is k-th column of matrix (I|-A) */
alpar@1:       char *stat; /* char stat[1+n]; */
alpar@1:       /* stat[0] is not used;
alpar@1:          stat[j], 1 <= j <= n, is the status of non-basic variable
alpar@1:          xN[j], which defines its active bound:
alpar@1:          GLP_NL - lower bound is active
alpar@1:          GLP_NU - upper bound is active
alpar@1:          GLP_NF - free variable
alpar@1:          GLP_NS - fixed variable */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* matrix B is the basis matrix; it is composed from columns of
alpar@1:          the augmented constraint matrix (I|-A) corresponding to basic
alpar@1:          variables and stored in a factorized (invertable) form */
alpar@1:       int valid;
alpar@1:       /* factorization is valid only if this flag is set */
alpar@1:       BFD *bfd; /* BFD bfd[1:m,1:m]; */
alpar@1:       /* factorized (invertable) form of the basis matrix */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* matrix N is a matrix composed from columns of the augmented
alpar@1:          constraint matrix (I|-A) corresponding to non-basic variables
alpar@1:          except fixed ones; it is stored by rows and changes every time
alpar@1:          the basis changes */
alpar@1:       int *N_ptr; /* int N_ptr[1+m+1]; */
alpar@1:       /* N_ptr[0] is not used;
alpar@1:          N_ptr[i], 1 <= i <= m, is starting position of i-th row in
alpar@1:          arrays N_ind and N_val; note that N_ptr[1] is always 1;
alpar@1:          N_ptr[m+1] indicates the position after the last element in
alpar@1:          arrays N_ind and N_val */
alpar@1:       int *N_len; /* int N_len[1+m]; */
alpar@1:       /* N_len[0] is not used;
alpar@1:          N_len[i], 1 <= i <= m, is length of i-th row (0 to n) */
alpar@1:       int *N_ind; /* int N_ind[N_ptr[m+1]]; */
alpar@1:       /* column indices */
alpar@1:       double *N_val; /* double N_val[N_ptr[m+1]]; */
alpar@1:       /* non-zero element values */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* working parameters */
alpar@1:       int phase;
alpar@1:       /* search phase:
alpar@1:          0 - not determined yet
alpar@1:          1 - search for primal feasible solution
alpar@1:          2 - search for optimal solution */
alpar@1:       glp_long tm_beg;
alpar@1:       /* time value at the beginning of the search */
alpar@1:       int it_beg;
alpar@1:       /* simplex iteration count at the beginning of the search */
alpar@1:       int it_cnt;
alpar@1:       /* simplex iteration count; it increases by one every time the
alpar@1:          basis changes (including the case when a non-basic variable
alpar@1:          jumps to its opposite bound) */
alpar@1:       int it_dpy;
alpar@1:       /* simplex iteration count at the most recent display output */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* basic solution components */
alpar@1:       double *bbar; /* double bbar[1+m]; */
alpar@1:       /* bbar[0] is not used;
alpar@1:          bbar[i], 1 <= i <= m, is primal value of basic variable xB[i]
alpar@1:          (if xB[i] is free, its primal value is not updated) */
alpar@1:       double *cbar; /* double cbar[1+n]; */
alpar@1:       /* cbar[0] is not used;
alpar@1:          cbar[j], 1 <= j <= n, is reduced cost of non-basic variable
alpar@1:          xN[j] (if xN[j] is fixed, its reduced cost is not updated) */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* the following pricing technique options may be used:
alpar@1:          GLP_PT_STD - standard ("textbook") pricing;
alpar@1:          GLP_PT_PSE - projected steepest edge;
alpar@1:          GLP_PT_DVX - Devex pricing (not implemented yet);
alpar@1:          in case of GLP_PT_STD the reference space is not used, and all
alpar@1:          steepest edge coefficients are set to 1 */
alpar@1:       int refct;
alpar@1:       /* this count is set to an initial value when the reference space
alpar@1:          is defined and decreases by one every time the basis changes;
alpar@1:          once this count reaches zero, the reference space is redefined
alpar@1:          again */
alpar@1:       char *refsp; /* char refsp[1+m+n]; */
alpar@1:       /* refsp[0] is not used;
alpar@1:          refsp[k], 1 <= k <= m+n, is the flag which means that variable
alpar@1:          x[k] belongs to the current reference space */
alpar@1:       double *gamma; /* double gamma[1+n]; */
alpar@1:       /* gamma[0] is not used;
alpar@1:          gamma[j], 1 <= j <= n, is the steepest edge coefficient for
alpar@1:          non-basic variable xN[j]; if xN[j] is fixed, gamma[j] is not
alpar@1:          used and just set to 1 */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* non-basic variable xN[q] chosen to enter the basis */
alpar@1:       int q;
alpar@1:       /* index of the non-basic variable xN[q] chosen, 1 <= q <= n;
alpar@1:          if the set of eligible non-basic variables is empty and thus
alpar@1:          no variable has been chosen, q is set to 0 */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* pivot column of the simplex table corresponding to non-basic
alpar@1:          variable xN[q] chosen is the following vector:
alpar@1:             T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q],
alpar@1:          where B is the current basis matrix, N[q] is a column of the
alpar@1:          matrix (I|-A) corresponding to xN[q] */
alpar@1:       int tcol_nnz;
alpar@1:       /* number of non-zero components, 0 <= nnz <= m */
alpar@1:       int *tcol_ind; /* int tcol_ind[1+m]; */
alpar@1:       /* tcol_ind[0] is not used;
alpar@1:          tcol_ind[t], 1 <= t <= nnz, is an index of non-zero component,
alpar@1:          i.e. tcol_ind[t] = i means that tcol_vec[i] != 0 */
alpar@1:       double *tcol_vec; /* double tcol_vec[1+m]; */
alpar@1:       /* tcol_vec[0] is not used;
alpar@1:          tcol_vec[i], 1 <= i <= m, is a numeric value of i-th component
alpar@1:          of the column */
alpar@1:       double tcol_max;
alpar@1:       /* infinity (maximum) norm of the column (max |tcol_vec[i]|) */
alpar@1:       int tcol_num;
alpar@1:       /* number of significant non-zero components, which means that:
alpar@1:          |tcol_vec[i]| >= eps for i in tcol_ind[1,...,num],
alpar@1:          |tcol_vec[i]| <  eps for i in tcol_ind[num+1,...,nnz],
alpar@1:          where eps is a pivot tolerance */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* basic variable xB[p] chosen to leave the basis */
alpar@1:       int p;
alpar@1:       /* index of the basic variable xB[p] chosen, 1 <= p <= m;
alpar@1:          p = 0 means that no basic variable reaches its bound;
alpar@1:          p < 0 means that non-basic variable xN[q] reaches its opposite
alpar@1:          bound before any basic variable */
alpar@1:       int p_stat;
alpar@1:       /* new status (GLP_NL, GLP_NU, or GLP_NS) to be assigned to xB[p]
alpar@1:          once it has left the basis */
alpar@1:       double teta;
alpar@1:       /* change of non-basic variable xN[q] (see above), on which xB[p]
alpar@1:          (or, if p < 0, xN[q] itself) reaches its bound */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* pivot row of the simplex table corresponding to basic variable
alpar@1:          xB[p] chosen is the following vector:
alpar@1:             T' * e[p] = - N' * inv(B') * e[p] = - N' * rho,
alpar@1:          where B' is a matrix transposed to the current basis matrix,
alpar@1:          N' is a matrix, whose rows are columns of the matrix (I|-A)
alpar@1:          corresponding to non-basic non-fixed variables */
alpar@1:       int trow_nnz;
alpar@1:       /* number of non-zero components, 0 <= nnz <= n */
alpar@1:       int *trow_ind; /* int trow_ind[1+n]; */
alpar@1:       /* trow_ind[0] is not used;
alpar@1:          trow_ind[t], 1 <= t <= nnz, is an index of non-zero component,
alpar@1:          i.e. trow_ind[t] = j means that trow_vec[j] != 0 */
alpar@1:       double *trow_vec; /* int trow_vec[1+n]; */
alpar@1:       /* trow_vec[0] is not used;
alpar@1:          trow_vec[j], 1 <= j <= n, is a numeric value of j-th component
alpar@1:          of the row */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* working arrays */
alpar@1:       double *work1; /* double work1[1+m]; */
alpar@1:       double *work2; /* double work2[1+m]; */
alpar@1:       double *work3; /* double work3[1+m]; */
alpar@1:       double *work4; /* double work4[1+m]; */
alpar@1: };
alpar@1: 
alpar@1: static const double kappa = 0.10;
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  alloc_csa - allocate common storage area
alpar@1: *
alpar@1: *  This routine allocates all arrays in the common storage area (CSA)
alpar@1: *  and returns a pointer to the CSA. */
alpar@1: 
alpar@1: static struct csa *alloc_csa(glp_prob *lp)
alpar@1: {     struct csa *csa;
alpar@1:       int m = lp->m;
alpar@1:       int n = lp->n;
alpar@1:       int nnz = lp->nnz;
alpar@1:       csa = xmalloc(sizeof(struct csa));
alpar@1:       xassert(m > 0 && n > 0);
alpar@1:       csa->m = m;
alpar@1:       csa->n = n;
alpar@1:       csa->type = xcalloc(1+m+n, sizeof(char));
alpar@1:       csa->lb = xcalloc(1+m+n, sizeof(double));
alpar@1:       csa->ub = xcalloc(1+m+n, sizeof(double));
alpar@1:       csa->coef = xcalloc(1+m+n, sizeof(double));
alpar@1:       csa->obj = xcalloc(1+n, sizeof(double));
alpar@1:       csa->A_ptr = xcalloc(1+n+1, sizeof(int));
alpar@1:       csa->A_ind = xcalloc(1+nnz, sizeof(int));
alpar@1:       csa->A_val = xcalloc(1+nnz, sizeof(double));
alpar@1:       csa->head = xcalloc(1+m+n, sizeof(int));
alpar@1:       csa->stat = xcalloc(1+n, sizeof(char));
alpar@1:       csa->N_ptr = xcalloc(1+m+1, sizeof(int));
alpar@1:       csa->N_len = xcalloc(1+m, sizeof(int));
alpar@1:       csa->N_ind = NULL; /* will be allocated later */
alpar@1:       csa->N_val = NULL; /* will be allocated later */
alpar@1:       csa->bbar = xcalloc(1+m, sizeof(double));
alpar@1:       csa->cbar = xcalloc(1+n, sizeof(double));
alpar@1:       csa->refsp = xcalloc(1+m+n, sizeof(char));
alpar@1:       csa->gamma = xcalloc(1+n, sizeof(double));
alpar@1:       csa->tcol_ind = xcalloc(1+m, sizeof(int));
alpar@1:       csa->tcol_vec = xcalloc(1+m, sizeof(double));
alpar@1:       csa->trow_ind = xcalloc(1+n, sizeof(int));
alpar@1:       csa->trow_vec = xcalloc(1+n, sizeof(double));
alpar@1:       csa->work1 = xcalloc(1+m, sizeof(double));
alpar@1:       csa->work2 = xcalloc(1+m, sizeof(double));
alpar@1:       csa->work3 = xcalloc(1+m, sizeof(double));
alpar@1:       csa->work4 = xcalloc(1+m, sizeof(double));
alpar@1:       return csa;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  init_csa - initialize common storage area
alpar@1: *
alpar@1: *  This routine initializes all data structures in the common storage
alpar@1: *  area (CSA). */
alpar@1: 
alpar@1: static void alloc_N(struct csa *csa);
alpar@1: static void build_N(struct csa *csa);
alpar@1: 
alpar@1: static void init_csa(struct csa *csa, glp_prob *lp)
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1:       char *type = csa->type;
alpar@1:       double *lb = csa->lb;
alpar@1:       double *ub = csa->ub;
alpar@1:       double *coef = csa->coef;
alpar@1:       double *obj = csa->obj;
alpar@1:       int *A_ptr = csa->A_ptr;
alpar@1:       int *A_ind = csa->A_ind;
alpar@1:       double *A_val = csa->A_val;
alpar@1:       int *head = csa->head;
alpar@1:       char *stat = csa->stat;
alpar@1:       char *refsp = csa->refsp;
alpar@1:       double *gamma = csa->gamma;
alpar@1:       int i, j, k, loc;
alpar@1:       double cmax;
alpar@1:       /* auxiliary variables */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  GLPROW *row = lp->row[i];
alpar@1:          type[i] = (char)row->type;
alpar@1:          lb[i] = row->lb * row->rii;
alpar@1:          ub[i] = row->ub * row->rii;
alpar@1:          coef[i] = 0.0;
alpar@1:       }
alpar@1:       /* structural variables */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  GLPCOL *col = lp->col[j];
alpar@1:          type[m+j] = (char)col->type;
alpar@1:          lb[m+j] = col->lb / col->sjj;
alpar@1:          ub[m+j] = col->ub / col->sjj;
alpar@1:          coef[m+j] = col->coef * col->sjj;
alpar@1:       }
alpar@1:       /* original objective function */
alpar@1:       obj[0] = lp->c0;
alpar@1:       memcpy(&obj[1], &coef[m+1], n * sizeof(double));
alpar@1:       /* factor used to scale original objective coefficients */
alpar@1:       cmax = 0.0;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:          if (cmax < fabs(obj[j])) cmax = fabs(obj[j]);
alpar@1:       if (cmax == 0.0) cmax = 1.0;
alpar@1:       switch (lp->dir)
alpar@1:       {  case GLP_MIN:
alpar@1:             csa->zeta = + 1.0 / cmax;
alpar@1:             break;
alpar@1:          case GLP_MAX:
alpar@1:             csa->zeta = - 1.0 / cmax;
alpar@1:             break;
alpar@1:          default:
alpar@1:             xassert(lp != lp);
alpar@1:       }
alpar@1: #if 1
alpar@1:       if (fabs(csa->zeta) < 1.0) csa->zeta *= 1000.0;
alpar@1: #endif
alpar@1:       /* matrix A (by columns) */
alpar@1:       loc = 1;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  GLPAIJ *aij;
alpar@1:          A_ptr[j] = loc;
alpar@1:          for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
alpar@1:          {  A_ind[loc] = aij->row->i;
alpar@1:             A_val[loc] = aij->row->rii * aij->val * aij->col->sjj;
alpar@1:             loc++;
alpar@1:          }
alpar@1:       }
alpar@1:       A_ptr[n+1] = loc;
alpar@1:       xassert(loc == lp->nnz+1);
alpar@1:       /* basis header */
alpar@1:       xassert(lp->valid);
alpar@1:       memcpy(&head[1], &lp->head[1], m * sizeof(int));
alpar@1:       k = 0;
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  GLPROW *row = lp->row[i];
alpar@1:          if (row->stat != GLP_BS)
alpar@1:          {  k++;
alpar@1:             xassert(k <= n);
alpar@1:             head[m+k] = i;
alpar@1:             stat[k] = (char)row->stat;
alpar@1:          }
alpar@1:       }
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  GLPCOL *col = lp->col[j];
alpar@1:          if (col->stat != GLP_BS)
alpar@1:          {  k++;
alpar@1:             xassert(k <= n);
alpar@1:             head[m+k] = m + j;
alpar@1:             stat[k] = (char)col->stat;
alpar@1:          }
alpar@1:       }
alpar@1:       xassert(k == n);
alpar@1:       /* factorization of matrix B */
alpar@1:       csa->valid = 1, lp->valid = 0;
alpar@1:       csa->bfd = lp->bfd, lp->bfd = NULL;
alpar@1:       /* matrix N (by rows) */
alpar@1:       alloc_N(csa);
alpar@1:       build_N(csa);
alpar@1:       /* working parameters */
alpar@1:       csa->phase = 0;
alpar@1:       csa->tm_beg = xtime();
alpar@1:       csa->it_beg = csa->it_cnt = lp->it_cnt;
alpar@1:       csa->it_dpy = -1;
alpar@1:       /* reference space and steepest edge coefficients */
alpar@1:       csa->refct = 0;
alpar@1:       memset(&refsp[1], 0, (m+n) * sizeof(char));
alpar@1:       for (j = 1; j <= n; j++) gamma[j] = 1.0;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  invert_B - compute factorization of the basis matrix
alpar@1: *
alpar@1: *  This routine computes factorization of the current basis matrix B.
alpar@1: *
alpar@1: *  If the operation is successful, the routine returns zero, otherwise
alpar@1: *  non-zero. */
alpar@1: 
alpar@1: static int inv_col(void *info, int i, int ind[], double val[])
alpar@1: {     /* this auxiliary routine returns row indices and numeric values
alpar@1:          of non-zero elements of i-th column of the basis matrix */
alpar@1:       struct csa *csa = info;
alpar@1:       int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       int *A_ptr = csa->A_ptr;
alpar@1:       int *A_ind = csa->A_ind;
alpar@1:       double *A_val = csa->A_val;
alpar@1:       int *head = csa->head;
alpar@1:       int k, len, ptr, t;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= i && i <= m);
alpar@1: #endif
alpar@1:       k = head[i]; /* B[i] is k-th column of (I|-A) */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:       if (k <= m)
alpar@1:       {  /* B[i] is k-th column of submatrix I */
alpar@1:          len = 1;
alpar@1:          ind[1] = k;
alpar@1:          val[1] = 1.0;
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* B[i] is (k-m)-th column of submatrix (-A) */
alpar@1:          ptr = A_ptr[k-m];
alpar@1:          len = A_ptr[k-m+1] - ptr;
alpar@1:          memcpy(&ind[1], &A_ind[ptr], len * sizeof(int));
alpar@1:          memcpy(&val[1], &A_val[ptr], len * sizeof(double));
alpar@1:          for (t = 1; t <= len; t++) val[t] = - val[t];
alpar@1:       }
alpar@1:       return len;
alpar@1: }
alpar@1: 
alpar@1: static int invert_B(struct csa *csa)
alpar@1: {     int ret;
alpar@1:       ret = bfd_factorize(csa->bfd, csa->m, NULL, inv_col, csa);
alpar@1:       csa->valid = (ret == 0);
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  update_B - update factorization of the basis matrix
alpar@1: *
alpar@1: *  This routine replaces i-th column of the basis matrix B by k-th
alpar@1: *  column of the augmented constraint matrix (I|-A) and then updates
alpar@1: *  the factorization of B.
alpar@1: *
alpar@1: *  If the factorization has been successfully updated, the routine
alpar@1: *  returns zero, otherwise non-zero. */
alpar@1: 
alpar@1: static int update_B(struct csa *csa, int i, int k)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       int ret;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= i && i <= m);
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:       if (k <= m)
alpar@1:       {  /* new i-th column of B is k-th column of I */
alpar@1:          int ind[1+1];
alpar@1:          double val[1+1];
alpar@1:          ind[1] = k;
alpar@1:          val[1] = 1.0;
alpar@1:          xassert(csa->valid);
alpar@1:          ret = bfd_update_it(csa->bfd, i, 0, 1, ind, val);
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* new i-th column of B is (k-m)-th column of (-A) */
alpar@1:          int *A_ptr = csa->A_ptr;
alpar@1:          int *A_ind = csa->A_ind;
alpar@1:          double *A_val = csa->A_val;
alpar@1:          double *val = csa->work1;
alpar@1:          int beg, end, ptr, len;
alpar@1:          beg = A_ptr[k-m];
alpar@1:          end = A_ptr[k-m+1];
alpar@1:          len = 0;
alpar@1:          for (ptr = beg; ptr < end; ptr++)
alpar@1:             val[++len] = - A_val[ptr];
alpar@1:          xassert(csa->valid);
alpar@1:          ret = bfd_update_it(csa->bfd, i, 0, len, &A_ind[beg-1], val);
alpar@1:       }
alpar@1:       csa->valid = (ret == 0);
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  error_ftran - compute residual vector r = h - B * x
alpar@1: *
alpar@1: *  This routine computes the residual vector r = h - B * x, where B is
alpar@1: *  the current basis matrix, h is the vector of right-hand sides, x is
alpar@1: *  the solution vector. */
alpar@1: 
alpar@1: static void error_ftran(struct csa *csa, double h[], double x[],
alpar@1:       double r[])
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       int *A_ptr = csa->A_ptr;
alpar@1:       int *A_ind = csa->A_ind;
alpar@1:       double *A_val = csa->A_val;
alpar@1:       int *head = csa->head;
alpar@1:       int i, k, beg, end, ptr;
alpar@1:       double temp;
alpar@1:       /* compute the residual vector:
alpar@1:          r = h - B * x = h - B[1] * x[1] - ... - B[m] * x[m],
alpar@1:          where B[1], ..., B[m] are columns of matrix B */
alpar@1:       memcpy(&r[1], &h[1], m * sizeof(double));
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  temp = x[i];
alpar@1:          if (temp == 0.0) continue;
alpar@1:          k = head[i]; /* B[i] is k-th column of (I|-A) */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          if (k <= m)
alpar@1:          {  /* B[i] is k-th column of submatrix I */
alpar@1:             r[k] -= temp;
alpar@1:          }
alpar@1:          else
alpar@1:          {  /* B[i] is (k-m)-th column of submatrix (-A) */
alpar@1:             beg = A_ptr[k-m];
alpar@1:             end = A_ptr[k-m+1];
alpar@1:             for (ptr = beg; ptr < end; ptr++)
alpar@1:                r[A_ind[ptr]] += A_val[ptr] * temp;
alpar@1:          }
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  refine_ftran - refine solution of B * x = h
alpar@1: *
alpar@1: *  This routine performs one iteration to refine the solution of
alpar@1: *  the system B * x = h, where B is the current basis matrix, h is the
alpar@1: *  vector of right-hand sides, x is the solution vector. */
alpar@1: 
alpar@1: static void refine_ftran(struct csa *csa, double h[], double x[])
alpar@1: {     int m = csa->m;
alpar@1:       double *r = csa->work1;
alpar@1:       double *d = csa->work1;
alpar@1:       int i;
alpar@1:       /* compute the residual vector r = h - B * x */
alpar@1:       error_ftran(csa, h, x, r);
alpar@1:       /* compute the correction vector d = inv(B) * r */
alpar@1:       xassert(csa->valid);
alpar@1:       bfd_ftran(csa->bfd, d);
alpar@1:       /* refine the solution vector (new x) = (old x) + d */
alpar@1:       for (i = 1; i <= m; i++) x[i] += d[i];
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  error_btran - compute residual vector r = h - B'* x
alpar@1: *
alpar@1: *  This routine computes the residual vector r = h - B'* x, where B'
alpar@1: *  is a matrix transposed to the current basis matrix, h is the vector
alpar@1: *  of right-hand sides, x is the solution vector. */
alpar@1: 
alpar@1: static void error_btran(struct csa *csa, double h[], double x[],
alpar@1:       double r[])
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       int *A_ptr = csa->A_ptr;
alpar@1:       int *A_ind = csa->A_ind;
alpar@1:       double *A_val = csa->A_val;
alpar@1:       int *head = csa->head;
alpar@1:       int i, k, beg, end, ptr;
alpar@1:       double temp;
alpar@1:       /* compute the residual vector r = b - B'* x */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  /* r[i] := b[i] - (i-th column of B)'* x */
alpar@1:          k = head[i]; /* B[i] is k-th column of (I|-A) */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          temp = h[i];
alpar@1:          if (k <= m)
alpar@1:          {  /* B[i] is k-th column of submatrix I */
alpar@1:             temp -= x[k];
alpar@1:          }
alpar@1:          else
alpar@1:          {  /* B[i] is (k-m)-th column of submatrix (-A) */
alpar@1:             beg = A_ptr[k-m];
alpar@1:             end = A_ptr[k-m+1];
alpar@1:             for (ptr = beg; ptr < end; ptr++)
alpar@1:                temp += A_val[ptr] * x[A_ind[ptr]];
alpar@1:          }
alpar@1:          r[i] = temp;
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  refine_btran - refine solution of B'* x = h
alpar@1: *
alpar@1: *  This routine performs one iteration to refine the solution of the
alpar@1: *  system B'* x = h, where B' is a matrix transposed to the current
alpar@1: *  basis matrix, h is the vector of right-hand sides, x is the solution
alpar@1: *  vector. */
alpar@1: 
alpar@1: static void refine_btran(struct csa *csa, double h[], double x[])
alpar@1: {     int m = csa->m;
alpar@1:       double *r = csa->work1;
alpar@1:       double *d = csa->work1;
alpar@1:       int i;
alpar@1:       /* compute the residual vector r = h - B'* x */
alpar@1:       error_btran(csa, h, x, r);
alpar@1:       /* compute the correction vector d = inv(B') * r */
alpar@1:       xassert(csa->valid);
alpar@1:       bfd_btran(csa->bfd, d);
alpar@1:       /* refine the solution vector (new x) = (old x) + d */
alpar@1:       for (i = 1; i <= m; i++) x[i] += d[i];
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  alloc_N - allocate matrix N
alpar@1: *
alpar@1: *  This routine determines maximal row lengths of matrix N, sets its
alpar@1: *  row pointers, and then allocates arrays N_ind and N_val.
alpar@1: *
alpar@1: *  Note that some fixed structural variables may temporarily become
alpar@1: *  double-bounded, so corresponding columns of matrix A should not be
alpar@1: *  ignored on calculating maximal row lengths of matrix N. */
alpar@1: 
alpar@1: static void alloc_N(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1:       int *A_ptr = csa->A_ptr;
alpar@1:       int *A_ind = csa->A_ind;
alpar@1:       int *N_ptr = csa->N_ptr;
alpar@1:       int *N_len = csa->N_len;
alpar@1:       int i, j, beg, end, ptr;
alpar@1:       /* determine number of non-zeros in each row of the augmented
alpar@1:          constraint matrix (I|-A) */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:          N_len[i] = 1;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  beg = A_ptr[j];
alpar@1:          end = A_ptr[j+1];
alpar@1:          for (ptr = beg; ptr < end; ptr++)
alpar@1:             N_len[A_ind[ptr]]++;
alpar@1:       }
alpar@1:       /* determine maximal row lengths of matrix N and set its row
alpar@1:          pointers */
alpar@1:       N_ptr[1] = 1;
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  /* row of matrix N cannot have more than n non-zeros */
alpar@1:          if (N_len[i] > n) N_len[i] = n;
alpar@1:          N_ptr[i+1] = N_ptr[i] + N_len[i];
alpar@1:       }
alpar@1:       /* now maximal number of non-zeros in matrix N is known */
alpar@1:       csa->N_ind = xcalloc(N_ptr[m+1], sizeof(int));
alpar@1:       csa->N_val = xcalloc(N_ptr[m+1], sizeof(double));
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  add_N_col - add column of matrix (I|-A) to matrix N
alpar@1: *
alpar@1: *  This routine adds j-th column to matrix N which is k-th column of
alpar@1: *  the augmented constraint matrix (I|-A). (It is assumed that old j-th
alpar@1: *  column was previously removed from matrix N.) */
alpar@1: 
alpar@1: static void add_N_col(struct csa *csa, int j, int k)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       int *N_ptr = csa->N_ptr;
alpar@1:       int *N_len = csa->N_len;
alpar@1:       int *N_ind = csa->N_ind;
alpar@1:       double *N_val = csa->N_val;
alpar@1:       int pos;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= j && j <= n);
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:       if (k <= m)
alpar@1:       {  /* N[j] is k-th column of submatrix I */
alpar@1:          pos = N_ptr[k] + (N_len[k]++);
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(pos < N_ptr[k+1]);
alpar@1: #endif
alpar@1:          N_ind[pos] = j;
alpar@1:          N_val[pos] = 1.0;
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* N[j] is (k-m)-th column of submatrix (-A) */
alpar@1:          int *A_ptr = csa->A_ptr;
alpar@1:          int *A_ind = csa->A_ind;
alpar@1:          double *A_val = csa->A_val;
alpar@1:          int i, beg, end, ptr;
alpar@1:          beg = A_ptr[k-m];
alpar@1:          end = A_ptr[k-m+1];
alpar@1:          for (ptr = beg; ptr < end; ptr++)
alpar@1:          {  i = A_ind[ptr]; /* row number */
alpar@1:             pos = N_ptr[i] + (N_len[i]++);
alpar@1: #ifdef GLP_DEBUG
alpar@1:             xassert(pos < N_ptr[i+1]);
alpar@1: #endif
alpar@1:             N_ind[pos] = j;
alpar@1:             N_val[pos] = - A_val[ptr];
alpar@1:          }
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  del_N_col - remove column of matrix (I|-A) from matrix N
alpar@1: *
alpar@1: *  This routine removes j-th column from matrix N which is k-th column
alpar@1: *  of the augmented constraint matrix (I|-A). */
alpar@1: 
alpar@1: static void del_N_col(struct csa *csa, int j, int k)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       int *N_ptr = csa->N_ptr;
alpar@1:       int *N_len = csa->N_len;
alpar@1:       int *N_ind = csa->N_ind;
alpar@1:       double *N_val = csa->N_val;
alpar@1:       int pos, head, tail;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= j && j <= n);
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:       if (k <= m)
alpar@1:       {  /* N[j] is k-th column of submatrix I */
alpar@1:          /* find element in k-th row of N */
alpar@1:          head = N_ptr[k];
alpar@1:          for (pos = head; N_ind[pos] != j; pos++) /* nop */;
alpar@1:          /* and remove it from the row list */
alpar@1:          tail = head + (--N_len[k]);
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(pos <= tail);
alpar@1: #endif
alpar@1:          N_ind[pos] = N_ind[tail];
alpar@1:          N_val[pos] = N_val[tail];
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* N[j] is (k-m)-th column of submatrix (-A) */
alpar@1:          int *A_ptr = csa->A_ptr;
alpar@1:          int *A_ind = csa->A_ind;
alpar@1:          int i, beg, end, ptr;
alpar@1:          beg = A_ptr[k-m];
alpar@1:          end = A_ptr[k-m+1];
alpar@1:          for (ptr = beg; ptr < end; ptr++)
alpar@1:          {  i = A_ind[ptr]; /* row number */
alpar@1:             /* find element in i-th row of N */
alpar@1:             head = N_ptr[i];
alpar@1:             for (pos = head; N_ind[pos] != j; pos++) /* nop */;
alpar@1:             /* and remove it from the row list */
alpar@1:             tail = head + (--N_len[i]);
alpar@1: #ifdef GLP_DEBUG
alpar@1:             xassert(pos <= tail);
alpar@1: #endif
alpar@1:             N_ind[pos] = N_ind[tail];
alpar@1:             N_val[pos] = N_val[tail];
alpar@1:          }
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  build_N - build matrix N for current basis
alpar@1: *
alpar@1: *  This routine builds matrix N for the current basis from columns
alpar@1: *  of the augmented constraint matrix (I|-A) corresponding to non-basic
alpar@1: *  non-fixed variables. */
alpar@1: 
alpar@1: static void build_N(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1:       int *head = csa->head;
alpar@1:       char *stat = csa->stat;
alpar@1:       int *N_len = csa->N_len;
alpar@1:       int j, k;
alpar@1:       /* N := empty matrix */
alpar@1:       memset(&N_len[1], 0, m * sizeof(int));
alpar@1:       /* go through non-basic columns of matrix (I|-A) */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  if (stat[j] != GLP_NS)
alpar@1:          {  /* xN[j] is non-fixed; add j-th column to matrix N which is
alpar@1:                k-th column of matrix (I|-A) */
alpar@1:             k = head[m+j]; /* x[k] = xN[j] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:             xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:             add_N_col(csa, j, k);
alpar@1:          }
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  get_xN - determine current value of non-basic variable xN[j]
alpar@1: *
alpar@1: *  This routine returns the current value of non-basic variable xN[j],
alpar@1: *  which is a value of its active bound. */
alpar@1: 
alpar@1: static double get_xN(struct csa *csa, int j)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       double *lb = csa->lb;
alpar@1:       double *ub = csa->ub;
alpar@1:       int *head = csa->head;
alpar@1:       char *stat = csa->stat;
alpar@1:       int k;
alpar@1:       double xN;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= j && j <= n);
alpar@1: #endif
alpar@1:       k = head[m+j]; /* x[k] = xN[j] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:       switch (stat[j])
alpar@1:       {  case GLP_NL:
alpar@1:             /* x[k] is on its lower bound */
alpar@1:             xN = lb[k]; break;
alpar@1:          case GLP_NU:
alpar@1:             /* x[k] is on its upper bound */
alpar@1:             xN = ub[k]; break;
alpar@1:          case GLP_NF:
alpar@1:             /* x[k] is free non-basic variable */
alpar@1:             xN = 0.0; break;
alpar@1:          case GLP_NS:
alpar@1:             /* x[k] is fixed non-basic variable */
alpar@1:             xN = lb[k]; break;
alpar@1:          default:
alpar@1:             xassert(stat != stat);
alpar@1:       }
alpar@1:       return xN;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  eval_beta - compute primal values of basic variables
alpar@1: *
alpar@1: *  This routine computes current primal values of all basic variables:
alpar@1: *
alpar@1: *     beta = - inv(B) * N * xN,
alpar@1: *
alpar@1: *  where B is the current basis matrix, N is a matrix built of columns
alpar@1: *  of matrix (I|-A) corresponding to non-basic variables, and xN is the
alpar@1: *  vector of current values of non-basic variables. */
alpar@1: 
alpar@1: static void eval_beta(struct csa *csa, double beta[])
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1:       int *A_ptr = csa->A_ptr;
alpar@1:       int *A_ind = csa->A_ind;
alpar@1:       double *A_val = csa->A_val;
alpar@1:       int *head = csa->head;
alpar@1:       double *h = csa->work2;
alpar@1:       int i, j, k, beg, end, ptr;
alpar@1:       double xN;
alpar@1:       /* compute the right-hand side vector:
alpar@1:          h := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n],
alpar@1:          where N[1], ..., N[n] are columns of matrix N */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:          h[i] = 0.0;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  k = head[m+j]; /* x[k] = xN[j] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          /* determine current value of xN[j] */
alpar@1:          xN = get_xN(csa, j);
alpar@1:          if (xN == 0.0) continue;
alpar@1:          if (k <= m)
alpar@1:          {  /* N[j] is k-th column of submatrix I */
alpar@1:             h[k] -= xN;
alpar@1:          }
alpar@1:          else
alpar@1:          {  /* N[j] is (k-m)-th column of submatrix (-A) */
alpar@1:             beg = A_ptr[k-m];
alpar@1:             end = A_ptr[k-m+1];
alpar@1:             for (ptr = beg; ptr < end; ptr++)
alpar@1:                h[A_ind[ptr]] += xN * A_val[ptr];
alpar@1:          }
alpar@1:       }
alpar@1:       /* solve system B * beta = h */
alpar@1:       memcpy(&beta[1], &h[1], m * sizeof(double));
alpar@1:       xassert(csa->valid);
alpar@1:       bfd_ftran(csa->bfd, beta);
alpar@1:       /* and refine the solution */
alpar@1:       refine_ftran(csa, h, beta);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  eval_pi - compute vector of simplex multipliers
alpar@1: *
alpar@1: *  This routine computes the vector of current simplex multipliers:
alpar@1: *
alpar@1: *     pi = inv(B') * cB,
alpar@1: *
alpar@1: *  where B' is a matrix transposed to the current basis matrix, cB is
alpar@1: *  a subvector of objective coefficients at basic variables. */
alpar@1: 
alpar@1: static void eval_pi(struct csa *csa, double pi[])
alpar@1: {     int m = csa->m;
alpar@1:       double *c = csa->coef;
alpar@1:       int *head = csa->head;
alpar@1:       double *cB = csa->work2;
alpar@1:       int i;
alpar@1:       /* construct the right-hand side vector cB */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:          cB[i] = c[head[i]];
alpar@1:       /* solve system B'* pi = cB */
alpar@1:       memcpy(&pi[1], &cB[1], m * sizeof(double));
alpar@1:       xassert(csa->valid);
alpar@1:       bfd_btran(csa->bfd, pi);
alpar@1:       /* and refine the solution */
alpar@1:       refine_btran(csa, cB, pi);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  eval_cost - compute reduced cost of non-basic variable xN[j]
alpar@1: *
alpar@1: *  This routine computes the current reduced cost of non-basic variable
alpar@1: *  xN[j]:
alpar@1: *
alpar@1: *     d[j] = cN[j] - N'[j] * pi,
alpar@1: *
alpar@1: *  where cN[j] is the objective coefficient at variable xN[j], N[j] is
alpar@1: *  a column of the augmented constraint matrix (I|-A) corresponding to
alpar@1: *  xN[j], pi is the vector of simplex multipliers. */
alpar@1: 
alpar@1: static double eval_cost(struct csa *csa, double pi[], int j)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       double *coef = csa->coef;
alpar@1:       int *head = csa->head;
alpar@1:       int k;
alpar@1:       double dj;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= j && j <= n);
alpar@1: #endif
alpar@1:       k = head[m+j]; /* x[k] = xN[j] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:       dj = coef[k];
alpar@1:       if (k <= m)
alpar@1:       {  /* N[j] is k-th column of submatrix I */
alpar@1:          dj -= pi[k];
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* N[j] is (k-m)-th column of submatrix (-A) */
alpar@1:          int *A_ptr = csa->A_ptr;
alpar@1:          int *A_ind = csa->A_ind;
alpar@1:          double *A_val = csa->A_val;
alpar@1:          int beg, end, ptr;
alpar@1:          beg = A_ptr[k-m];
alpar@1:          end = A_ptr[k-m+1];
alpar@1:          for (ptr = beg; ptr < end; ptr++)
alpar@1:             dj += A_val[ptr] * pi[A_ind[ptr]];
alpar@1:       }
alpar@1:       return dj;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  eval_bbar - compute and store primal values of basic variables
alpar@1: *
alpar@1: *  This routine computes primal values of all basic variables and then
alpar@1: *  stores them in the solution array. */
alpar@1: 
alpar@1: static void eval_bbar(struct csa *csa)
alpar@1: {     eval_beta(csa, csa->bbar);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  eval_cbar - compute and store reduced costs of non-basic variables
alpar@1: *
alpar@1: *  This routine computes reduced costs of all non-basic variables and
alpar@1: *  then stores them in the solution array. */
alpar@1: 
alpar@1: static void eval_cbar(struct csa *csa)
alpar@1: {
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int m = csa->m;
alpar@1: #endif
alpar@1:       int n = csa->n;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int *head = csa->head;
alpar@1: #endif
alpar@1:       double *cbar = csa->cbar;
alpar@1:       double *pi = csa->work3;
alpar@1:       int j;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int k;
alpar@1: #endif
alpar@1:       /* compute simplex multipliers */
alpar@1:       eval_pi(csa, pi);
alpar@1:       /* compute and store reduced costs */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {
alpar@1: #ifdef GLP_DEBUG
alpar@1:          k = head[m+j]; /* x[k] = xN[j] */
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          cbar[j] = eval_cost(csa, pi, j);
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  reset_refsp - reset the reference space
alpar@1: *
alpar@1: *  This routine resets (redefines) the reference space used in the
alpar@1: *  projected steepest edge pricing algorithm. */
alpar@1: 
alpar@1: static void reset_refsp(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1:       int *head = csa->head;
alpar@1:       char *refsp = csa->refsp;
alpar@1:       double *gamma = csa->gamma;
alpar@1:       int j, k;
alpar@1:       xassert(csa->refct == 0);
alpar@1:       csa->refct = 1000;
alpar@1:       memset(&refsp[1], 0, (m+n) * sizeof(char));
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  k = head[m+j]; /* x[k] = xN[j] */
alpar@1:          refsp[k] = 1;
alpar@1:          gamma[j] = 1.0;
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  eval_gamma - compute steepest edge coefficient
alpar@1: *
alpar@1: *  This routine computes the steepest edge coefficient for non-basic
alpar@1: *  variable xN[j] using its direct definition:
alpar@1: *
alpar@1: *     gamma[j] = delta[j] +  sum   alfa[i,j]^2,
alpar@1: *                           i in R
alpar@1: *
alpar@1: *  where delta[j] = 1, if xN[j] is in the current reference space,
alpar@1: *  and 0 otherwise; R is a set of basic variables xB[i], which are in
alpar@1: *  the current reference space; alfa[i,j] are elements of the current
alpar@1: *  simplex table.
alpar@1: *
alpar@1: *  NOTE: The routine is intended only for debugginig purposes. */
alpar@1: 
alpar@1: static double eval_gamma(struct csa *csa, int j)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       int *head = csa->head;
alpar@1:       char *refsp = csa->refsp;
alpar@1:       double *alfa = csa->work3;
alpar@1:       double *h = csa->work3;
alpar@1:       int i, k;
alpar@1:       double gamma;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= j && j <= n);
alpar@1: #endif
alpar@1:       k = head[m+j]; /* x[k] = xN[j] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:       /* construct the right-hand side vector h = - N[j] */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:          h[i] = 0.0;
alpar@1:       if (k <= m)
alpar@1:       {  /* N[j] is k-th column of submatrix I */
alpar@1:          h[k] = -1.0;
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* N[j] is (k-m)-th column of submatrix (-A) */
alpar@1:          int *A_ptr = csa->A_ptr;
alpar@1:          int *A_ind = csa->A_ind;
alpar@1:          double *A_val = csa->A_val;
alpar@1:          int beg, end, ptr;
alpar@1:          beg = A_ptr[k-m];
alpar@1:          end = A_ptr[k-m+1];
alpar@1:          for (ptr = beg; ptr < end; ptr++)
alpar@1:             h[A_ind[ptr]] = A_val[ptr];
alpar@1:       }
alpar@1:       /* solve system B * alfa = h */
alpar@1:       xassert(csa->valid);
alpar@1:       bfd_ftran(csa->bfd, alfa);
alpar@1:       /* compute gamma */
alpar@1:       gamma = (refsp[k] ? 1.0 : 0.0);
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  k = head[i];
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          if (refsp[k]) gamma += alfa[i] * alfa[i];
alpar@1:       }
alpar@1:       return gamma;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  chuzc - choose non-basic variable (column of the simplex table)
alpar@1: *
alpar@1: *  This routine chooses non-basic variable xN[q], which has largest
alpar@1: *  weighted reduced cost:
alpar@1: *
alpar@1: *     |d[q]| / sqrt(gamma[q]) = max  |d[j]| / sqrt(gamma[j]),
alpar@1: *                              j in J
alpar@1: *
alpar@1: *  where J is a subset of eligible non-basic variables xN[j], d[j] is
alpar@1: *  reduced cost of xN[j], gamma[j] is the steepest edge coefficient.
alpar@1: *
alpar@1: *  The working objective function is always minimized, so the sign of
alpar@1: *  d[q] determines direction, in which xN[q] has to change:
alpar@1: *
alpar@1: *     if d[q] < 0, xN[q] has to increase;
alpar@1: *
alpar@1: *     if d[q] > 0, xN[q] has to decrease.
alpar@1: *
alpar@1: *  If |d[j]| <= tol_dj, where tol_dj is a specified tolerance, xN[j]
alpar@1: *  is not included in J and therefore ignored. (It is assumed that the
alpar@1: *  working objective row is appropriately scaled, i.e. max|c[k]| = 1.)
alpar@1: *
alpar@1: *  If J is empty and no variable has been chosen, q is set to 0. */
alpar@1: 
alpar@1: static void chuzc(struct csa *csa, double tol_dj)
alpar@1: {     int n = csa->n;
alpar@1:       char *stat = csa->stat;
alpar@1:       double *cbar = csa->cbar;
alpar@1:       double *gamma = csa->gamma;
alpar@1:       int j, q;
alpar@1:       double dj, best, temp;
alpar@1:       /* nothing is chosen so far */
alpar@1:       q = 0, best = 0.0;
alpar@1:       /* look through the list of non-basic variables */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  dj = cbar[j];
alpar@1:          switch (stat[j])
alpar@1:          {  case GLP_NL:
alpar@1:                /* xN[j] can increase */
alpar@1:                if (dj >= - tol_dj) continue;
alpar@1:                break;
alpar@1:             case GLP_NU:
alpar@1:                /* xN[j] can decrease */
alpar@1:                if (dj <= + tol_dj) continue;
alpar@1:                break;
alpar@1:             case GLP_NF:
alpar@1:                /* xN[j] can change in any direction */
alpar@1:                if (- tol_dj <= dj && dj <= + tol_dj) continue;
alpar@1:                break;
alpar@1:             case GLP_NS:
alpar@1:                /* xN[j] cannot change at all */
alpar@1:                continue;
alpar@1:             default:
alpar@1:                xassert(stat != stat);
alpar@1:          }
alpar@1:          /* xN[j] is eligible non-basic variable; choose one which has
alpar@1:             largest weighted reduced cost */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(gamma[j] > 0.0);
alpar@1: #endif
alpar@1:          temp = (dj * dj) / gamma[j];
alpar@1:          if (best < temp)
alpar@1:             q = j, best = temp;
alpar@1:       }
alpar@1:       /* store the index of non-basic variable xN[q] chosen */
alpar@1:       csa->q = q;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  eval_tcol - compute pivot column of the simplex table
alpar@1: *
alpar@1: *  This routine computes the pivot column of the simplex table, which
alpar@1: *  corresponds to non-basic variable xN[q] chosen.
alpar@1: *
alpar@1: *  The pivot column is the following vector:
alpar@1: *
alpar@1: *     tcol = T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q],
alpar@1: *
alpar@1: *  where B is the current basis matrix, N[q] is a column of the matrix
alpar@1: *  (I|-A) corresponding to variable xN[q]. */
alpar@1: 
alpar@1: static void eval_tcol(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       int *head = csa->head;
alpar@1:       int q = csa->q;
alpar@1:       int *tcol_ind = csa->tcol_ind;
alpar@1:       double *tcol_vec = csa->tcol_vec;
alpar@1:       double *h = csa->tcol_vec;
alpar@1:       int i, k, nnz;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= q && q <= n);
alpar@1: #endif
alpar@1:       k = head[m+q]; /* x[k] = xN[q] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:       /* construct the right-hand side vector h = - N[q] */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:          h[i] = 0.0;
alpar@1:       if (k <= m)
alpar@1:       {  /* N[q] is k-th column of submatrix I */
alpar@1:          h[k] = -1.0;
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* N[q] is (k-m)-th column of submatrix (-A) */
alpar@1:          int *A_ptr = csa->A_ptr;
alpar@1:          int *A_ind = csa->A_ind;
alpar@1:          double *A_val = csa->A_val;
alpar@1:          int beg, end, ptr;
alpar@1:          beg = A_ptr[k-m];
alpar@1:          end = A_ptr[k-m+1];
alpar@1:          for (ptr = beg; ptr < end; ptr++)
alpar@1:             h[A_ind[ptr]] = A_val[ptr];
alpar@1:       }
alpar@1:       /* solve system B * tcol = h */
alpar@1:       xassert(csa->valid);
alpar@1:       bfd_ftran(csa->bfd, tcol_vec);
alpar@1:       /* construct sparse pattern of the pivot column */
alpar@1:       nnz = 0;
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  if (tcol_vec[i] != 0.0)
alpar@1:             tcol_ind[++nnz] = i;
alpar@1:       }
alpar@1:       csa->tcol_nnz = nnz;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  refine_tcol - refine pivot column of the simplex table
alpar@1: *
alpar@1: *  This routine refines the pivot column of the simplex table assuming
alpar@1: *  that it was previously computed by the routine eval_tcol. */
alpar@1: 
alpar@1: static void refine_tcol(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       int *head = csa->head;
alpar@1:       int q = csa->q;
alpar@1:       int *tcol_ind = csa->tcol_ind;
alpar@1:       double *tcol_vec = csa->tcol_vec;
alpar@1:       double *h = csa->work3;
alpar@1:       int i, k, nnz;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= q && q <= n);
alpar@1: #endif
alpar@1:       k = head[m+q]; /* x[k] = xN[q] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:       /* construct the right-hand side vector h = - N[q] */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:          h[i] = 0.0;
alpar@1:       if (k <= m)
alpar@1:       {  /* N[q] is k-th column of submatrix I */
alpar@1:          h[k] = -1.0;
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* N[q] is (k-m)-th column of submatrix (-A) */
alpar@1:          int *A_ptr = csa->A_ptr;
alpar@1:          int *A_ind = csa->A_ind;
alpar@1:          double *A_val = csa->A_val;
alpar@1:          int beg, end, ptr;
alpar@1:          beg = A_ptr[k-m];
alpar@1:          end = A_ptr[k-m+1];
alpar@1:          for (ptr = beg; ptr < end; ptr++)
alpar@1:             h[A_ind[ptr]] = A_val[ptr];
alpar@1:       }
alpar@1:       /* refine solution of B * tcol = h */
alpar@1:       refine_ftran(csa, h, tcol_vec);
alpar@1:       /* construct sparse pattern of the pivot column */
alpar@1:       nnz = 0;
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  if (tcol_vec[i] != 0.0)
alpar@1:             tcol_ind[++nnz] = i;
alpar@1:       }
alpar@1:       csa->tcol_nnz = nnz;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  sort_tcol - sort pivot column of the simplex table
alpar@1: *
alpar@1: *  This routine reorders the list of non-zero elements of the pivot
alpar@1: *  column to put significant elements, whose magnitude is not less than
alpar@1: *  a specified tolerance, in front of the list, and stores the number
alpar@1: *  of significant elements in tcol_num. */
alpar@1: 
alpar@1: static void sort_tcol(struct csa *csa, double tol_piv)
alpar@1: {
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int m = csa->m;
alpar@1: #endif
alpar@1:       int nnz = csa->tcol_nnz;
alpar@1:       int *tcol_ind = csa->tcol_ind;
alpar@1:       double *tcol_vec = csa->tcol_vec;
alpar@1:       int i, num, pos;
alpar@1:       double big, eps, temp;
alpar@1:       /* compute infinity (maximum) norm of the column */
alpar@1:       big = 0.0;
alpar@1:       for (pos = 1; pos <= nnz; pos++)
alpar@1:       {
alpar@1: #ifdef GLP_DEBUG
alpar@1:          i = tcol_ind[pos];
alpar@1:          xassert(1 <= i && i <= m);
alpar@1: #endif
alpar@1:          temp = fabs(tcol_vec[tcol_ind[pos]]);
alpar@1:          if (big < temp) big = temp;
alpar@1:       }
alpar@1:       csa->tcol_max = big;
alpar@1:       /* determine absolute pivot tolerance */
alpar@1:       eps = tol_piv * (1.0 + 0.01 * big);
alpar@1:       /* move significant column components to front of the list */
alpar@1:       for (num = 0; num < nnz; )
alpar@1:       {  i = tcol_ind[nnz];
alpar@1:          if (fabs(tcol_vec[i]) < eps)
alpar@1:             nnz--;
alpar@1:          else
alpar@1:          {  num++;
alpar@1:             tcol_ind[nnz] = tcol_ind[num];
alpar@1:             tcol_ind[num] = i;
alpar@1:          }
alpar@1:       }
alpar@1:       csa->tcol_num = num;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  chuzr - choose basic variable (row of the simplex table)
alpar@1: *
alpar@1: *  This routine chooses basic variable xB[p], which reaches its bound
alpar@1: *  first on changing non-basic variable xN[q] in valid direction.
alpar@1: *
alpar@1: *  The parameter rtol is a relative tolerance used to relax bounds of
alpar@1: *  basic variables. If rtol = 0, the routine implements the standard
alpar@1: *  ratio test. Otherwise, if rtol > 0, the routine implements Harris'
alpar@1: *  two-pass ratio test. In the latter case rtol should be about three
alpar@1: *  times less than a tolerance used to check primal feasibility. */
alpar@1: 
alpar@1: static void chuzr(struct csa *csa, double rtol)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       char *type = csa->type;
alpar@1:       double *lb = csa->lb;
alpar@1:       double *ub = csa->ub;
alpar@1:       double *coef = csa->coef;
alpar@1:       int *head = csa->head;
alpar@1:       int phase = csa->phase;
alpar@1:       double *bbar = csa->bbar;
alpar@1:       double *cbar = csa->cbar;
alpar@1:       int q = csa->q;
alpar@1:       int *tcol_ind = csa->tcol_ind;
alpar@1:       double *tcol_vec = csa->tcol_vec;
alpar@1:       int tcol_num = csa->tcol_num;
alpar@1:       int i, i_stat, k, p, p_stat, pos;
alpar@1:       double alfa, big, delta, s, t, teta, tmax;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= q && q <= n);
alpar@1: #endif
alpar@1:       /* s := - sign(d[q]), where d[q] is reduced cost of xN[q] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(cbar[q] != 0.0);
alpar@1: #endif
alpar@1:       s = (cbar[q] > 0.0 ? -1.0 : +1.0);
alpar@1:       /*** FIRST PASS ***/
alpar@1:       k = head[m+q]; /* x[k] = xN[q] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:       if (type[k] == GLP_DB)
alpar@1:       {  /* xN[q] has both lower and upper bounds */
alpar@1:          p = -1, p_stat = 0, teta = ub[k] - lb[k], big = 1.0;
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* xN[q] has no opposite bound */
alpar@1:          p = 0, p_stat = 0, teta = DBL_MAX, big = 0.0;
alpar@1:       }
alpar@1:       /* walk through significant elements of the pivot column */
alpar@1:       for (pos = 1; pos <= tcol_num; pos++)
alpar@1:       {  i = tcol_ind[pos];
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= i && i <= m);
alpar@1: #endif
alpar@1:          k = head[i]; /* x[k] = xB[i] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          alfa = s * tcol_vec[i];
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(alfa != 0.0);
alpar@1: #endif
alpar@1:          /* xB[i] = ... + alfa * xN[q] + ..., and due to s we need to
alpar@1:             consider the only case when xN[q] is increasing */
alpar@1:          if (alfa > 0.0)
alpar@1:          {  /* xB[i] is increasing */
alpar@1:             if (phase == 1 && coef[k] < 0.0)
alpar@1:             {  /* xB[i] violates its lower bound, which plays the role
alpar@1:                   of an upper bound on phase I */
alpar@1:                delta = rtol * (1.0 + kappa * fabs(lb[k]));
alpar@1:                t = ((lb[k] + delta) - bbar[i]) / alfa;
alpar@1:                i_stat = GLP_NL;
alpar@1:             }
alpar@1:             else if (phase == 1 && coef[k] > 0.0)
alpar@1:             {  /* xB[i] violates its upper bound, which plays the role
alpar@1:                   of an lower bound on phase I */
alpar@1:                continue;
alpar@1:             }
alpar@1:             else if (type[k] == GLP_UP || type[k] == GLP_DB ||
alpar@1:                      type[k] == GLP_FX)
alpar@1:             {  /* xB[i] is within its bounds and has an upper bound */
alpar@1:                delta = rtol * (1.0 + kappa * fabs(ub[k]));
alpar@1:                t = ((ub[k] + delta) - bbar[i]) / alfa;
alpar@1:                i_stat = GLP_NU;
alpar@1:             }
alpar@1:             else
alpar@1:             {  /* xB[i] is within its bounds and has no upper bound */
alpar@1:                continue;
alpar@1:             }
alpar@1:          }
alpar@1:          else
alpar@1:          {  /* xB[i] is decreasing */
alpar@1:             if (phase == 1 && coef[k] > 0.0)
alpar@1:             {  /* xB[i] violates its upper bound, which plays the role
alpar@1:                   of an lower bound on phase I */
alpar@1:                delta = rtol * (1.0 + kappa * fabs(ub[k]));
alpar@1:                t = ((ub[k] - delta) - bbar[i]) / alfa;
alpar@1:                i_stat = GLP_NU;
alpar@1:             }
alpar@1:             else if (phase == 1 && coef[k] < 0.0)
alpar@1:             {  /* xB[i] violates its lower bound, which plays the role
alpar@1:                   of an upper bound on phase I */
alpar@1:                continue;
alpar@1:             }
alpar@1:             else if (type[k] == GLP_LO || type[k] == GLP_DB ||
alpar@1:                      type[k] == GLP_FX)
alpar@1:             {  /* xB[i] is within its bounds and has an lower bound */
alpar@1:                delta = rtol * (1.0 + kappa * fabs(lb[k]));
alpar@1:                t = ((lb[k] - delta) - bbar[i]) / alfa;
alpar@1:                i_stat = GLP_NL;
alpar@1:             }
alpar@1:             else
alpar@1:             {  /* xB[i] is within its bounds and has no lower bound */
alpar@1:                continue;
alpar@1:             }
alpar@1:          }
alpar@1:          /* t is a change of xN[q], on which xB[i] reaches its bound
alpar@1:             (possibly relaxed); since the basic solution is assumed to
alpar@1:             be primal feasible (or pseudo feasible on phase I), t has
alpar@1:             to be non-negative by definition; however, it may happen
alpar@1:             that xB[i] slightly (i.e. within a tolerance) violates its
alpar@1:             bound, that leads to negative t; in the latter case, if
alpar@1:             xB[i] is chosen, negative t means that xN[q] changes in
alpar@1:             wrong direction; if pivot alfa[i,q] is close to zero, even
alpar@1:             small bound violation of xB[i] may lead to a large change
alpar@1:             of xN[q] in wrong direction; let, for example, xB[i] >= 0
alpar@1:             and in the current basis its value be -5e-9; let also xN[q]
alpar@1:             be on its zero bound and should increase; from the ratio
alpar@1:             test rule it follows that the pivot alfa[i,q] < 0; however,
alpar@1:             if alfa[i,q] is, say, -1e-9, the change of xN[q] in wrong
alpar@1:             direction is 5e-9 / (-1e-9) = -5, and using it for updating
alpar@1:             values of other basic variables will give absolutely wrong
alpar@1:             results; therefore, if t is negative, we should replace it
alpar@1:             by exact zero assuming that xB[i] is exactly on its bound,
alpar@1:             and the violation appears due to round-off errors */
alpar@1:          if (t < 0.0) t = 0.0;
alpar@1:          /* apply minimal ratio test */
alpar@1:          if (teta > t || teta == t && big < fabs(alfa))
alpar@1:             p = i, p_stat = i_stat, teta = t, big = fabs(alfa);
alpar@1:       }
alpar@1:       /* the second pass is skipped in the following cases: */
alpar@1:       /* if the standard ratio test is used */
alpar@1:       if (rtol == 0.0) goto done;
alpar@1:       /* if xN[q] reaches its opposite bound or if no basic variable
alpar@1:          has been chosen on the first pass */
alpar@1:       if (p <= 0) goto done;
alpar@1:       /* if xB[p] is a blocking variable, i.e. if it prevents xN[q]
alpar@1:          from any change */
alpar@1:       if (teta == 0.0) goto done;
alpar@1:       /*** SECOND PASS ***/
alpar@1:       /* here tmax is a maximal change of xN[q], on which the solution
alpar@1:          remains primal feasible (or pseudo feasible on phase I) within
alpar@1:          a tolerance */
alpar@1: #if 0
alpar@1:       tmax = (1.0 + 10.0 * DBL_EPSILON) * teta;
alpar@1: #else
alpar@1:       tmax = teta;
alpar@1: #endif
alpar@1:       /* nothing is chosen so far */
alpar@1:       p = 0, p_stat = 0, teta = DBL_MAX, big = 0.0;
alpar@1:       /* walk through significant elements of the pivot column */
alpar@1:       for (pos = 1; pos <= tcol_num; pos++)
alpar@1:       {  i = tcol_ind[pos];
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= i && i <= m);
alpar@1: #endif
alpar@1:          k = head[i]; /* x[k] = xB[i] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          alfa = s * tcol_vec[i];
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(alfa != 0.0);
alpar@1: #endif
alpar@1:          /* xB[i] = ... + alfa * xN[q] + ..., and due to s we need to
alpar@1:             consider the only case when xN[q] is increasing */
alpar@1:          if (alfa > 0.0)
alpar@1:          {  /* xB[i] is increasing */
alpar@1:             if (phase == 1 && coef[k] < 0.0)
alpar@1:             {  /* xB[i] violates its lower bound, which plays the role
alpar@1:                   of an upper bound on phase I */
alpar@1:                t = (lb[k] - bbar[i]) / alfa;
alpar@1:                i_stat = GLP_NL;
alpar@1:             }
alpar@1:             else if (phase == 1 && coef[k] > 0.0)
alpar@1:             {  /* xB[i] violates its upper bound, which plays the role
alpar@1:                   of an lower bound on phase I */
alpar@1:                continue;
alpar@1:             }
alpar@1:             else if (type[k] == GLP_UP || type[k] == GLP_DB ||
alpar@1:                      type[k] == GLP_FX)
alpar@1:             {  /* xB[i] is within its bounds and has an upper bound */
alpar@1:                t = (ub[k] - bbar[i]) / alfa;
alpar@1:                i_stat = GLP_NU;
alpar@1:             }
alpar@1:             else
alpar@1:             {  /* xB[i] is within its bounds and has no upper bound */
alpar@1:                continue;
alpar@1:             }
alpar@1:          }
alpar@1:          else
alpar@1:          {  /* xB[i] is decreasing */
alpar@1:             if (phase == 1 && coef[k] > 0.0)
alpar@1:             {  /* xB[i] violates its upper bound, which plays the role
alpar@1:                   of an lower bound on phase I */
alpar@1:                t = (ub[k] - bbar[i]) / alfa;
alpar@1:                i_stat = GLP_NU;
alpar@1:             }
alpar@1:             else if (phase == 1 && coef[k] < 0.0)
alpar@1:             {  /* xB[i] violates its lower bound, which plays the role
alpar@1:                   of an upper bound on phase I */
alpar@1:                continue;
alpar@1:             }
alpar@1:             else if (type[k] == GLP_LO || type[k] == GLP_DB ||
alpar@1:                      type[k] == GLP_FX)
alpar@1:             {  /* xB[i] is within its bounds and has an lower bound */
alpar@1:                t = (lb[k] - bbar[i]) / alfa;
alpar@1:                i_stat = GLP_NL;
alpar@1:             }
alpar@1:             else
alpar@1:             {  /* xB[i] is within its bounds and has no lower bound */
alpar@1:                continue;
alpar@1:             }
alpar@1:          }
alpar@1:          /* (see comments for the first pass) */
alpar@1:          if (t < 0.0) t = 0.0;
alpar@1:          /* t is a change of xN[q], on which xB[i] reaches its bound;
alpar@1:             if t <= tmax, all basic variables can violate their bounds
alpar@1:             only within relaxation tolerance delta; we can use this
alpar@1:             freedom and choose basic variable having largest influence
alpar@1:             coefficient to avoid possible numeric instability */
alpar@1:          if (t <= tmax && big < fabs(alfa))
alpar@1:             p = i, p_stat = i_stat, teta = t, big = fabs(alfa);
alpar@1:       }
alpar@1:       /* something must be chosen on the second pass */
alpar@1:       xassert(p != 0);
alpar@1: done: /* store the index and status of basic variable xB[p] chosen */
alpar@1:       csa->p = p;
alpar@1:       if (p > 0 && type[head[p]] == GLP_FX)
alpar@1:          csa->p_stat = GLP_NS;
alpar@1:       else
alpar@1:          csa->p_stat = p_stat;
alpar@1:       /* store corresponding change of non-basic variable xN[q] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(teta >= 0.0);
alpar@1: #endif
alpar@1:       csa->teta = s * teta;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  eval_rho - compute pivot row of the inverse
alpar@1: *
alpar@1: *  This routine computes the pivot (p-th) row of the inverse inv(B),
alpar@1: *  which corresponds to basic variable xB[p] chosen:
alpar@1: *
alpar@1: *     rho = inv(B') * e[p],
alpar@1: *
alpar@1: *  where B' is a matrix transposed to the current basis matrix, e[p]
alpar@1: *  is unity vector. */
alpar@1: 
alpar@1: static void eval_rho(struct csa *csa, double rho[])
alpar@1: {     int m = csa->m;
alpar@1:       int p = csa->p;
alpar@1:       double *e = rho;
alpar@1:       int i;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= p && p <= m);
alpar@1: #endif
alpar@1:       /* construct the right-hand side vector e[p] */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:          e[i] = 0.0;
alpar@1:       e[p] = 1.0;
alpar@1:       /* solve system B'* rho = e[p] */
alpar@1:       xassert(csa->valid);
alpar@1:       bfd_btran(csa->bfd, rho);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  refine_rho - refine pivot row of the inverse
alpar@1: *
alpar@1: *  This routine refines the pivot row of the inverse inv(B) assuming
alpar@1: *  that it was previously computed by the routine eval_rho. */
alpar@1: 
alpar@1: static void refine_rho(struct csa *csa, double rho[])
alpar@1: {     int m = csa->m;
alpar@1:       int p = csa->p;
alpar@1:       double *e = csa->work3;
alpar@1:       int i;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= p && p <= m);
alpar@1: #endif
alpar@1:       /* construct the right-hand side vector e[p] */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:          e[i] = 0.0;
alpar@1:       e[p] = 1.0;
alpar@1:       /* refine solution of B'* rho = e[p] */
alpar@1:       refine_btran(csa, e, rho);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  eval_trow - compute pivot row of the simplex table
alpar@1: *
alpar@1: *  This routine computes the pivot row of the simplex table, which
alpar@1: *  corresponds to basic variable xB[p] chosen.
alpar@1: *
alpar@1: *  The pivot row is the following vector:
alpar@1: *
alpar@1: *     trow = T'* e[p] = - N'* inv(B') * e[p] = - N' * rho,
alpar@1: *
alpar@1: *  where rho is the pivot row of the inverse inv(B) previously computed
alpar@1: *  by the routine eval_rho.
alpar@1: *
alpar@1: *  Note that elements of the pivot row corresponding to fixed non-basic
alpar@1: *  variables are not computed. */
alpar@1: 
alpar@1: static void eval_trow(struct csa *csa, double rho[])
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       char *stat = csa->stat;
alpar@1: #endif
alpar@1:       int *N_ptr = csa->N_ptr;
alpar@1:       int *N_len = csa->N_len;
alpar@1:       int *N_ind = csa->N_ind;
alpar@1:       double *N_val = csa->N_val;
alpar@1:       int *trow_ind = csa->trow_ind;
alpar@1:       double *trow_vec = csa->trow_vec;
alpar@1:       int i, j, beg, end, ptr, nnz;
alpar@1:       double temp;
alpar@1:       /* clear the pivot row */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:          trow_vec[j] = 0.0;
alpar@1:       /* compute the pivot row as a linear combination of rows of the
alpar@1:          matrix N: trow = - rho[1] * N'[1] - ... - rho[m] * N'[m] */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  temp = rho[i];
alpar@1:          if (temp == 0.0) continue;
alpar@1:          /* trow := trow - rho[i] * N'[i] */
alpar@1:          beg = N_ptr[i];
alpar@1:          end = beg + N_len[i];
alpar@1:          for (ptr = beg; ptr < end; ptr++)
alpar@1:          {
alpar@1: #ifdef GLP_DEBUG
alpar@1:             j = N_ind[ptr];
alpar@1:             xassert(1 <= j && j <= n);
alpar@1:             xassert(stat[j] != GLP_NS);
alpar@1: #endif
alpar@1:             trow_vec[N_ind[ptr]] -= temp * N_val[ptr];
alpar@1:          }
alpar@1:       }
alpar@1:       /* construct sparse pattern of the pivot row */
alpar@1:       nnz = 0;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  if (trow_vec[j] != 0.0)
alpar@1:             trow_ind[++nnz] = j;
alpar@1:       }
alpar@1:       csa->trow_nnz = nnz;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  update_bbar - update values of basic variables
alpar@1: *
alpar@1: *  This routine updates values of all basic variables for the adjacent
alpar@1: *  basis. */
alpar@1: 
alpar@1: static void update_bbar(struct csa *csa)
alpar@1: {
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       double *bbar = csa->bbar;
alpar@1:       int q = csa->q;
alpar@1:       int tcol_nnz = csa->tcol_nnz;
alpar@1:       int *tcol_ind = csa->tcol_ind;
alpar@1:       double *tcol_vec = csa->tcol_vec;
alpar@1:       int p = csa->p;
alpar@1:       double teta = csa->teta;
alpar@1:       int i, pos;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= q && q <= n);
alpar@1:       xassert(p < 0 || 1 <= p && p <= m);
alpar@1: #endif
alpar@1:       /* if xN[q] leaves the basis, compute its value in the adjacent
alpar@1:          basis, where it will replace xB[p] */
alpar@1:       if (p > 0)
alpar@1:          bbar[p] = get_xN(csa, q) + teta;
alpar@1:       /* update values of other basic variables (except xB[p], because
alpar@1:          it will be replaced by xN[q]) */
alpar@1:       if (teta == 0.0) goto done;
alpar@1:       for (pos = 1; pos <= tcol_nnz; pos++)
alpar@1:       {  i = tcol_ind[pos];
alpar@1:          /* skip xB[p] */
alpar@1:          if (i == p) continue;
alpar@1:          /* (change of xB[i]) = alfa[i,q] * (change of xN[q]) */
alpar@1:          bbar[i] += tcol_vec[i] * teta;
alpar@1:       }
alpar@1: done: return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  reeval_cost - recompute reduced cost of non-basic variable xN[q]
alpar@1: *
alpar@1: *  This routine recomputes reduced cost of non-basic variable xN[q] for
alpar@1: *  the current basis more accurately using its direct definition:
alpar@1: *
alpar@1: *     d[q] = cN[q] - N'[q] * pi =
alpar@1: *
alpar@1: *          = cN[q] - N'[q] * (inv(B') * cB) =
alpar@1: *
alpar@1: *          = cN[q] - (cB' * inv(B) * N[q]) =
alpar@1: *
alpar@1: *          = cN[q] + cB' * (pivot column).
alpar@1: *
alpar@1: *  It is assumed that the pivot column of the simplex table is already
alpar@1: *  computed. */
alpar@1: 
alpar@1: static double reeval_cost(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       double *coef = csa->coef;
alpar@1:       int *head = csa->head;
alpar@1:       int q = csa->q;
alpar@1:       int tcol_nnz = csa->tcol_nnz;
alpar@1:       int *tcol_ind = csa->tcol_ind;
alpar@1:       double *tcol_vec = csa->tcol_vec;
alpar@1:       int i, pos;
alpar@1:       double dq;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= q && q <= n);
alpar@1: #endif
alpar@1:       dq = coef[head[m+q]];
alpar@1:       for (pos = 1; pos <= tcol_nnz; pos++)
alpar@1:       {  i = tcol_ind[pos];
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= i && i <= m);
alpar@1: #endif
alpar@1:          dq += coef[head[i]] * tcol_vec[i];
alpar@1:       }
alpar@1:       return dq;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  update_cbar - update reduced costs of non-basic variables
alpar@1: *
alpar@1: *  This routine updates reduced costs of all (except fixed) non-basic
alpar@1: *  variables for the adjacent basis. */
alpar@1: 
alpar@1: static void update_cbar(struct csa *csa)
alpar@1: {
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       double *cbar = csa->cbar;
alpar@1:       int q = csa->q;
alpar@1:       int trow_nnz = csa->trow_nnz;
alpar@1:       int *trow_ind = csa->trow_ind;
alpar@1:       double *trow_vec = csa->trow_vec;
alpar@1:       int j, pos;
alpar@1:       double new_dq;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= q && q <= n);
alpar@1: #endif
alpar@1:       /* compute reduced cost of xB[p] in the adjacent basis, where it
alpar@1:          will replace xN[q] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(trow_vec[q] != 0.0);
alpar@1: #endif
alpar@1:       new_dq = (cbar[q] /= trow_vec[q]);
alpar@1:       /* update reduced costs of other non-basic variables (except
alpar@1:          xN[q], because it will be replaced by xB[p]) */
alpar@1:       for (pos = 1; pos <= trow_nnz; pos++)
alpar@1:       {  j = trow_ind[pos];
alpar@1:          /* skip xN[q] */
alpar@1:          if (j == q) continue;
alpar@1:          cbar[j] -= trow_vec[j] * new_dq;
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  update_gamma - update steepest edge coefficients
alpar@1: *
alpar@1: *  This routine updates steepest-edge coefficients for the adjacent
alpar@1: *  basis. */
alpar@1: 
alpar@1: static void update_gamma(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       char *type = csa->type;
alpar@1:       int *A_ptr = csa->A_ptr;
alpar@1:       int *A_ind = csa->A_ind;
alpar@1:       double *A_val = csa->A_val;
alpar@1:       int *head = csa->head;
alpar@1:       char *refsp = csa->refsp;
alpar@1:       double *gamma = csa->gamma;
alpar@1:       int q = csa->q;
alpar@1:       int tcol_nnz = csa->tcol_nnz;
alpar@1:       int *tcol_ind = csa->tcol_ind;
alpar@1:       double *tcol_vec = csa->tcol_vec;
alpar@1:       int p = csa->p;
alpar@1:       int trow_nnz = csa->trow_nnz;
alpar@1:       int *trow_ind = csa->trow_ind;
alpar@1:       double *trow_vec = csa->trow_vec;
alpar@1:       double *u = csa->work3;
alpar@1:       int i, j, k, pos, beg, end, ptr;
alpar@1:       double gamma_q, delta_q, pivot, s, t, t1, t2;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= p && p <= m);
alpar@1:       xassert(1 <= q && q <= n);
alpar@1: #endif
alpar@1:       /* the basis changes, so decrease the count */
alpar@1:       xassert(csa->refct > 0);
alpar@1:       csa->refct--;
alpar@1:       /* recompute gamma[q] for the current basis more accurately and
alpar@1:          compute auxiliary vector u */
alpar@1:       gamma_q = delta_q = (refsp[head[m+q]] ? 1.0 : 0.0);
alpar@1:       for (i = 1; i <= m; i++) u[i] = 0.0;
alpar@1:       for (pos = 1; pos <= tcol_nnz; pos++)
alpar@1:       {  i = tcol_ind[pos];
alpar@1:          if (refsp[head[i]])
alpar@1:          {  u[i] = t = tcol_vec[i];
alpar@1:             gamma_q += t * t;
alpar@1:          }
alpar@1:          else
alpar@1:             u[i] = 0.0;
alpar@1:       }
alpar@1:       xassert(csa->valid);
alpar@1:       bfd_btran(csa->bfd, u);
alpar@1:       /* update gamma[k] for other non-basic variables (except fixed
alpar@1:          variables and xN[q], because it will be replaced by xB[p]) */
alpar@1:       pivot = trow_vec[q];
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(pivot != 0.0);
alpar@1: #endif
alpar@1:       for (pos = 1; pos <= trow_nnz; pos++)
alpar@1:       {  j = trow_ind[pos];
alpar@1:          /* skip xN[q] */
alpar@1:          if (j == q) continue;
alpar@1:          /* compute t */
alpar@1:          t = trow_vec[j] / pivot;
alpar@1:          /* compute inner product s = N'[j] * u */
alpar@1:          k = head[m+j]; /* x[k] = xN[j] */
alpar@1:          if (k <= m)
alpar@1:             s = u[k];
alpar@1:          else
alpar@1:          {  s = 0.0;
alpar@1:             beg = A_ptr[k-m];
alpar@1:             end = A_ptr[k-m+1];
alpar@1:             for (ptr = beg; ptr < end; ptr++)
alpar@1:                s -= A_val[ptr] * u[A_ind[ptr]];
alpar@1:          }
alpar@1:          /* compute gamma[k] for the adjacent basis */
alpar@1:          t1 = gamma[j] + t * t * gamma_q + 2.0 * t * s;
alpar@1:          t2 = (refsp[k] ? 1.0 : 0.0) + delta_q * t * t;
alpar@1:          gamma[j] = (t1 >= t2 ? t1 : t2);
alpar@1:          if (gamma[j] < DBL_EPSILON) gamma[j] = DBL_EPSILON;
alpar@1:       }
alpar@1:       /* compute gamma[q] for the adjacent basis */
alpar@1:       if (type[head[p]] == GLP_FX)
alpar@1:          gamma[q] = 1.0;
alpar@1:       else
alpar@1:       {  gamma[q] = gamma_q / (pivot * pivot);
alpar@1:          if (gamma[q] < DBL_EPSILON) gamma[q] = DBL_EPSILON;
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  err_in_bbar - compute maximal relative error in primal solution
alpar@1: *
alpar@1: *  This routine returns maximal relative error:
alpar@1: *
alpar@1: *     max |beta[i] - bbar[i]| / (1 + |beta[i]|),
alpar@1: *
alpar@1: *  where beta and bbar are, respectively, directly computed and the
alpar@1: *  current (updated) values of basic variables.
alpar@1: *
alpar@1: *  NOTE: The routine is intended only for debugginig purposes. */
alpar@1: 
alpar@1: static double err_in_bbar(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1:       double *bbar = csa->bbar;
alpar@1:       int i;
alpar@1:       double e, emax, *beta;
alpar@1:       beta = xcalloc(1+m, sizeof(double));
alpar@1:       eval_beta(csa, beta);
alpar@1:       emax = 0.0;
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  e = fabs(beta[i] - bbar[i]) / (1.0 + fabs(beta[i]));
alpar@1:          if (emax < e) emax = e;
alpar@1:       }
alpar@1:       xfree(beta);
alpar@1:       return emax;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  err_in_cbar - compute maximal relative error in dual solution
alpar@1: *
alpar@1: *  This routine returns maximal relative error:
alpar@1: *
alpar@1: *     max |cost[j] - cbar[j]| / (1 + |cost[j]|),
alpar@1: *
alpar@1: *  where cost and cbar are, respectively, directly computed and the
alpar@1: *  current (updated) reduced costs of non-basic non-fixed variables.
alpar@1: *
alpar@1: *  NOTE: The routine is intended only for debugginig purposes. */
alpar@1: 
alpar@1: static double err_in_cbar(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1:       char *stat = csa->stat;
alpar@1:       double *cbar = csa->cbar;
alpar@1:       int j;
alpar@1:       double e, emax, cost, *pi;
alpar@1:       pi = xcalloc(1+m, sizeof(double));
alpar@1:       eval_pi(csa, pi);
alpar@1:       emax = 0.0;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  if (stat[j] == GLP_NS) continue;
alpar@1:          cost = eval_cost(csa, pi, j);
alpar@1:          e = fabs(cost - cbar[j]) / (1.0 + fabs(cost));
alpar@1:          if (emax < e) emax = e;
alpar@1:       }
alpar@1:       xfree(pi);
alpar@1:       return emax;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  err_in_gamma - compute maximal relative error in steepest edge cff.
alpar@1: *
alpar@1: *  This routine returns maximal relative error:
alpar@1: *
alpar@1: *     max |gamma'[j] - gamma[j]| / (1 + |gamma'[j]),
alpar@1: *
alpar@1: *  where gamma'[j] and gamma[j] are, respectively, directly computed
alpar@1: *  and the current (updated) steepest edge coefficients for non-basic
alpar@1: *  non-fixed variable x[j].
alpar@1: *
alpar@1: *  NOTE: The routine is intended only for debugginig purposes. */
alpar@1: 
alpar@1: static double err_in_gamma(struct csa *csa)
alpar@1: {     int n = csa->n;
alpar@1:       char *stat = csa->stat;
alpar@1:       double *gamma = csa->gamma;
alpar@1:       int j;
alpar@1:       double e, emax, temp;
alpar@1:       emax = 0.0;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  if (stat[j] == GLP_NS)
alpar@1:          {  xassert(gamma[j] == 1.0);
alpar@1:             continue;
alpar@1:          }
alpar@1:          temp = eval_gamma(csa, j);
alpar@1:          e = fabs(temp - gamma[j]) / (1.0 + fabs(temp));
alpar@1:          if (emax < e) emax = e;
alpar@1:       }
alpar@1:       return emax;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  change_basis - change basis header
alpar@1: *
alpar@1: *  This routine changes the basis header to make it corresponding to
alpar@1: *  the adjacent basis. */
alpar@1: 
alpar@1: static void change_basis(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1:       char *type = csa->type;
alpar@1: #endif
alpar@1:       int *head = csa->head;
alpar@1:       char *stat = csa->stat;
alpar@1:       int q = csa->q;
alpar@1:       int p = csa->p;
alpar@1:       int p_stat = csa->p_stat;
alpar@1:       int k;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(1 <= q && q <= n);
alpar@1: #endif
alpar@1:       if (p < 0)
alpar@1:       {  /* xN[q] goes to its opposite bound */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          k = head[m+q]; /* x[k] = xN[q] */
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1:          xassert(type[k] == GLP_DB);
alpar@1: #endif
alpar@1:          switch (stat[q])
alpar@1:          {  case GLP_NL:
alpar@1:                /* xN[q] increases */
alpar@1:                stat[q] = GLP_NU;
alpar@1:                break;
alpar@1:             case GLP_NU:
alpar@1:                /* xN[q] decreases */
alpar@1:                stat[q] = GLP_NL;
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(stat != stat);
alpar@1:          }
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* xB[p] leaves the basis, xN[q] enters the basis */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= p && p <= m);
alpar@1:          k = head[p]; /* x[k] = xB[p] */
alpar@1:          switch (p_stat)
alpar@1:          {  case GLP_NL:
alpar@1:                /* xB[p] goes to its lower bound */
alpar@1:                xassert(type[k] == GLP_LO || type[k] == GLP_DB);
alpar@1:                break;
alpar@1:             case GLP_NU:
alpar@1:                /* xB[p] goes to its upper bound */
alpar@1:                xassert(type[k] == GLP_UP || type[k] == GLP_DB);
alpar@1:                break;
alpar@1:             case GLP_NS:
alpar@1:                /* xB[p] goes to its fixed value */
alpar@1:                xassert(type[k] == GLP_NS);
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(p_stat != p_stat);
alpar@1:          }
alpar@1: #endif
alpar@1:          /* xB[p] <-> xN[q] */
alpar@1:          k = head[p], head[p] = head[m+q], head[m+q] = k;
alpar@1:          stat[q] = (char)p_stat;
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  set_aux_obj - construct auxiliary objective function
alpar@1: *
alpar@1: *  The auxiliary objective function is a separable piecewise linear
alpar@1: *  convex function, which is the sum of primal infeasibilities:
alpar@1: *
alpar@1: *     z = t[1] + ... + t[m+n] -> minimize,
alpar@1: *
alpar@1: *  where:
alpar@1: *
alpar@1: *            / lb[k] - x[k], if x[k] < lb[k]
alpar@1: *            |
alpar@1: *     t[k] = <  0, if lb[k] <= x[k] <= ub[k]
alpar@1: *            |
alpar@1: *            \ x[k] - ub[k], if x[k] > ub[k]
alpar@1: *
alpar@1: *  This routine computes objective coefficients for the current basis
alpar@1: *  and returns the number of non-zero terms t[k]. */
alpar@1: 
alpar@1: static int set_aux_obj(struct csa *csa, double tol_bnd)
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1:       char *type = csa->type;
alpar@1:       double *lb = csa->lb;
alpar@1:       double *ub = csa->ub;
alpar@1:       double *coef = csa->coef;
alpar@1:       int *head = csa->head;
alpar@1:       double *bbar = csa->bbar;
alpar@1:       int i, k, cnt = 0;
alpar@1:       double eps;
alpar@1:       /* use a bit more restrictive tolerance */
alpar@1:       tol_bnd *= 0.90;
alpar@1:       /* clear all objective coefficients */
alpar@1:       for (k = 1; k <= m+n; k++)
alpar@1:          coef[k] = 0.0;
alpar@1:       /* walk through the list of basic variables */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  k = head[i]; /* x[k] = xB[i] */
alpar@1:          if (type[k] == GLP_LO || type[k] == GLP_DB ||
alpar@1:              type[k] == GLP_FX)
alpar@1:          {  /* x[k] has lower bound */
alpar@1:             eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
alpar@1:             if (bbar[i] < lb[k] - eps)
alpar@1:             {  /* and violates it */
alpar@1:                coef[k] = -1.0;
alpar@1:                cnt++;
alpar@1:             }
alpar@1:          }
alpar@1:          if (type[k] == GLP_UP || type[k] == GLP_DB ||
alpar@1:              type[k] == GLP_FX)
alpar@1:          {  /* x[k] has upper bound */
alpar@1:             eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
alpar@1:             if (bbar[i] > ub[k] + eps)
alpar@1:             {  /* and violates it */
alpar@1:                coef[k] = +1.0;
alpar@1:                cnt++;
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       return cnt;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  set_orig_obj - restore original objective function
alpar@1: *
alpar@1: *  This routine assigns scaled original objective coefficients to the
alpar@1: *  working objective function. */
alpar@1: 
alpar@1: static void set_orig_obj(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1:       double *coef = csa->coef;
alpar@1:       double *obj = csa->obj;
alpar@1:       double zeta = csa->zeta;
alpar@1:       int i, j;
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:          coef[i] = 0.0;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:          coef[m+j] = zeta * obj[j];
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  check_stab - check numerical stability of basic solution
alpar@1: *
alpar@1: *  If the current basic solution is primal feasible (or pseudo feasible
alpar@1: *  on phase I) within a tolerance, this routine returns zero, otherwise
alpar@1: *  it returns non-zero. */
alpar@1: 
alpar@1: static int check_stab(struct csa *csa, double tol_bnd)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       char *type = csa->type;
alpar@1:       double *lb = csa->lb;
alpar@1:       double *ub = csa->ub;
alpar@1:       double *coef = csa->coef;
alpar@1:       int *head = csa->head;
alpar@1:       int phase = csa->phase;
alpar@1:       double *bbar = csa->bbar;
alpar@1:       int i, k;
alpar@1:       double eps;
alpar@1:       /* walk through the list of basic variables */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  k = head[i]; /* x[k] = xB[i] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          if (phase == 1 && coef[k] < 0.0)
alpar@1:          {  /* x[k] must not be greater than its lower bound */
alpar@1: #ifdef GLP_DEBUG
alpar@1:             xassert(type[k] == GLP_LO || type[k] == GLP_DB ||
alpar@1:                     type[k] == GLP_FX);
alpar@1: #endif
alpar@1:             eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
alpar@1:             if (bbar[i] > lb[k] + eps) return 1;
alpar@1:          }
alpar@1:          else if (phase == 1 && coef[k] > 0.0)
alpar@1:          {  /* x[k] must not be less than its upper bound */
alpar@1: #ifdef GLP_DEBUG
alpar@1:             xassert(type[k] == GLP_UP || type[k] == GLP_DB ||
alpar@1:                     type[k] == GLP_FX);
alpar@1: #endif
alpar@1:             eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
alpar@1:             if (bbar[i] < ub[k] - eps) return 1;
alpar@1:          }
alpar@1:          else
alpar@1:          {  /* either phase = 1 and coef[k] = 0, or phase = 2 */
alpar@1:             if (type[k] == GLP_LO || type[k] == GLP_DB ||
alpar@1:                 type[k] == GLP_FX)
alpar@1:             {  /* x[k] must not be less than its lower bound */
alpar@1:                eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
alpar@1:                if (bbar[i] < lb[k] - eps) return 1;
alpar@1:             }
alpar@1:             if (type[k] == GLP_UP || type[k] == GLP_DB ||
alpar@1:                 type[k] == GLP_FX)
alpar@1:             {  /* x[k] must not be greater then its upper bound */
alpar@1:                eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
alpar@1:                if (bbar[i] > ub[k] + eps) return 1;
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       /* basic solution is primal feasible within a tolerance */
alpar@1:       return 0;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  check_feas - check primal feasibility of basic solution
alpar@1: *
alpar@1: *  If the current basic solution is primal feasible within a tolerance,
alpar@1: *  this routine returns zero, otherwise it returns non-zero. */
alpar@1: 
alpar@1: static int check_feas(struct csa *csa, double tol_bnd)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1:       char *type = csa->type;
alpar@1: #endif
alpar@1:       double *lb = csa->lb;
alpar@1:       double *ub = csa->ub;
alpar@1:       double *coef = csa->coef;
alpar@1:       int *head = csa->head;
alpar@1:       double *bbar = csa->bbar;
alpar@1:       int i, k;
alpar@1:       double eps;
alpar@1:       xassert(csa->phase == 1);
alpar@1:       /* walk through the list of basic variables */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  k = head[i]; /* x[k] = xB[i] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          if (coef[k] < 0.0)
alpar@1:          {  /* check if x[k] still violates its lower bound */
alpar@1: #ifdef GLP_DEBUG
alpar@1:             xassert(type[k] == GLP_LO || type[k] == GLP_DB ||
alpar@1:                     type[k] == GLP_FX);
alpar@1: #endif
alpar@1:             eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
alpar@1:             if (bbar[i] < lb[k] - eps) return 1;
alpar@1:          }
alpar@1:          else if (coef[k] > 0.0)
alpar@1:          {  /* check if x[k] still violates its upper bound */
alpar@1: #ifdef GLP_DEBUG
alpar@1:             xassert(type[k] == GLP_UP || type[k] == GLP_DB ||
alpar@1:                     type[k] == GLP_FX);
alpar@1: #endif
alpar@1:             eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
alpar@1:             if (bbar[i] > ub[k] + eps) return 1;
alpar@1:          }
alpar@1:       }
alpar@1:       /* basic solution is primal feasible within a tolerance */
alpar@1:       return 0;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  eval_obj - compute original objective function
alpar@1: *
alpar@1: *  This routine computes the current value of the original objective
alpar@1: *  function. */
alpar@1: 
alpar@1: static double eval_obj(struct csa *csa)
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1:       double *obj = csa->obj;
alpar@1:       int *head = csa->head;
alpar@1:       double *bbar = csa->bbar;
alpar@1:       int i, j, k;
alpar@1:       double sum;
alpar@1:       sum = obj[0];
alpar@1:       /* walk through the list of basic variables */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  k = head[i]; /* x[k] = xB[i] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          if (k > m)
alpar@1:             sum += obj[k-m] * bbar[i];
alpar@1:       }
alpar@1:       /* walk through the list of non-basic variables */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  k = head[m+j]; /* x[k] = xN[j] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          if (k > m)
alpar@1:             sum += obj[k-m] * get_xN(csa, j);
alpar@1:       }
alpar@1:       return sum;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  display - display the search progress
alpar@1: *
alpar@1: *  This routine displays some information about the search progress
alpar@1: *  that includes:
alpar@1: *
alpar@1: *  the search phase;
alpar@1: *
alpar@1: *  the number of simplex iterations performed by the solver;
alpar@1: *
alpar@1: *  the original objective value;
alpar@1: *
alpar@1: *  the sum of (scaled) primal infeasibilities;
alpar@1: *
alpar@1: *  the number of basic fixed variables. */
alpar@1: 
alpar@1: static void display(struct csa *csa, const glp_smcp *parm, int spec)
alpar@1: {     int m = csa->m;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       int n = csa->n;
alpar@1: #endif
alpar@1:       char *type = csa->type;
alpar@1:       double *lb = csa->lb;
alpar@1:       double *ub = csa->ub;
alpar@1:       int phase = csa->phase;
alpar@1:       int *head = csa->head;
alpar@1:       double *bbar = csa->bbar;
alpar@1:       int i, k, cnt;
alpar@1:       double sum;
alpar@1:       if (parm->msg_lev < GLP_MSG_ON) goto skip;
alpar@1:       if (parm->out_dly > 0 &&
alpar@1:          1000.0 * xdifftime(xtime(), csa->tm_beg) < parm->out_dly)
alpar@1:          goto skip;
alpar@1:       if (csa->it_cnt == csa->it_dpy) goto skip;
alpar@1:       if (!spec && csa->it_cnt % parm->out_frq != 0) goto skip;
alpar@1:       /* compute the sum of primal infeasibilities and determine the
alpar@1:          number of basic fixed variables */
alpar@1:       sum = 0.0, cnt = 0;
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  k = head[i]; /* x[k] = xB[i] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          if (type[k] == GLP_LO || type[k] == GLP_DB ||
alpar@1:              type[k] == GLP_FX)
alpar@1:          {  /* x[k] has lower bound */
alpar@1:             if (bbar[i] < lb[k])
alpar@1:                sum += (lb[k] - bbar[i]);
alpar@1:          }
alpar@1:          if (type[k] == GLP_UP || type[k] == GLP_DB ||
alpar@1:              type[k] == GLP_FX)
alpar@1:          {  /* x[k] has upper bound */
alpar@1:             if (bbar[i] > ub[k])
alpar@1:                sum += (bbar[i] - ub[k]);
alpar@1:          }
alpar@1:          if (type[k] == GLP_FX) cnt++;
alpar@1:       }
alpar@1:       xprintf("%c%6d: obj = %17.9e  infeas = %10.3e (%d)\n",
alpar@1:          phase == 1 ? ' ' : '*', csa->it_cnt, eval_obj(csa), sum, cnt);
alpar@1:       csa->it_dpy = csa->it_cnt;
alpar@1: skip: return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  store_sol - store basic solution back to the problem object
alpar@1: *
alpar@1: *  This routine stores basic solution components back to the problem
alpar@1: *  object. */
alpar@1: 
alpar@1: static void store_sol(struct csa *csa, glp_prob *lp, int p_stat,
alpar@1:       int d_stat, int ray)
alpar@1: {     int m = csa->m;
alpar@1:       int n = csa->n;
alpar@1:       double zeta = csa->zeta;
alpar@1:       int *head = csa->head;
alpar@1:       char *stat = csa->stat;
alpar@1:       double *bbar = csa->bbar;
alpar@1:       double *cbar = csa->cbar;
alpar@1:       int i, j, k;
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(lp->m == m);
alpar@1:       xassert(lp->n == n);
alpar@1: #endif
alpar@1:       /* basis factorization */
alpar@1: #ifdef GLP_DEBUG
alpar@1:       xassert(!lp->valid && lp->bfd == NULL);
alpar@1:       xassert(csa->valid && csa->bfd != NULL);
alpar@1: #endif
alpar@1:       lp->valid = 1, csa->valid = 0;
alpar@1:       lp->bfd = csa->bfd, csa->bfd = NULL;
alpar@1:       memcpy(&lp->head[1], &head[1], m * sizeof(int));
alpar@1:       /* basic solution status */
alpar@1:       lp->pbs_stat = p_stat;
alpar@1:       lp->dbs_stat = d_stat;
alpar@1:       /* objective function value */
alpar@1:       lp->obj_val = eval_obj(csa);
alpar@1:       /* simplex iteration count */
alpar@1:       lp->it_cnt = csa->it_cnt;
alpar@1:       /* unbounded ray */
alpar@1:       lp->some = ray;
alpar@1:       /* basic variables */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  k = head[i]; /* x[k] = xB[i] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          if (k <= m)
alpar@1:          {  GLPROW *row = lp->row[k];
alpar@1:             row->stat = GLP_BS;
alpar@1:             row->bind = i;
alpar@1:             row->prim = bbar[i] / row->rii;
alpar@1:             row->dual = 0.0;
alpar@1:          }
alpar@1:          else
alpar@1:          {  GLPCOL *col = lp->col[k-m];
alpar@1:             col->stat = GLP_BS;
alpar@1:             col->bind = i;
alpar@1:             col->prim = bbar[i] * col->sjj;
alpar@1:             col->dual = 0.0;
alpar@1:          }
alpar@1:       }
alpar@1:       /* non-basic variables */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  k = head[m+j]; /* x[k] = xN[j] */
alpar@1: #ifdef GLP_DEBUG
alpar@1:          xassert(1 <= k && k <= m+n);
alpar@1: #endif
alpar@1:          if (k <= m)
alpar@1:          {  GLPROW *row = lp->row[k];
alpar@1:             row->stat = stat[j];
alpar@1:             row->bind = 0;
alpar@1: #if 0
alpar@1:             row->prim = get_xN(csa, j) / row->rii;
alpar@1: #else
alpar@1:             switch (stat[j])
alpar@1:             {  case GLP_NL:
alpar@1:                   row->prim = row->lb; break;
alpar@1:                case GLP_NU:
alpar@1:                   row->prim = row->ub; break;
alpar@1:                case GLP_NF:
alpar@1:                   row->prim = 0.0; break;
alpar@1:                case GLP_NS:
alpar@1:                   row->prim = row->lb; break;
alpar@1:                default:
alpar@1:                   xassert(stat != stat);
alpar@1:             }
alpar@1: #endif
alpar@1:             row->dual = (cbar[j] * row->rii) / zeta;
alpar@1:          }
alpar@1:          else
alpar@1:          {  GLPCOL *col = lp->col[k-m];
alpar@1:             col->stat = stat[j];
alpar@1:             col->bind = 0;
alpar@1: #if 0
alpar@1:             col->prim = get_xN(csa, j) * col->sjj;
alpar@1: #else
alpar@1:             switch (stat[j])
alpar@1:             {  case GLP_NL:
alpar@1:                   col->prim = col->lb; break;
alpar@1:                case GLP_NU:
alpar@1:                   col->prim = col->ub; break;
alpar@1:                case GLP_NF:
alpar@1:                   col->prim = 0.0; break;
alpar@1:                case GLP_NS:
alpar@1:                   col->prim = col->lb; break;
alpar@1:                default:
alpar@1:                   xassert(stat != stat);
alpar@1:             }
alpar@1: #endif
alpar@1:             col->dual = (cbar[j] / col->sjj) / zeta;
alpar@1:          }
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  free_csa - deallocate common storage area
alpar@1: *
alpar@1: *  This routine frees all the memory allocated to arrays in the common
alpar@1: *  storage area (CSA). */
alpar@1: 
alpar@1: static void free_csa(struct csa *csa)
alpar@1: {     xfree(csa->type);
alpar@1:       xfree(csa->lb);
alpar@1:       xfree(csa->ub);
alpar@1:       xfree(csa->coef);
alpar@1:       xfree(csa->obj);
alpar@1:       xfree(csa->A_ptr);
alpar@1:       xfree(csa->A_ind);
alpar@1:       xfree(csa->A_val);
alpar@1:       xfree(csa->head);
alpar@1:       xfree(csa->stat);
alpar@1:       xfree(csa->N_ptr);
alpar@1:       xfree(csa->N_len);
alpar@1:       xfree(csa->N_ind);
alpar@1:       xfree(csa->N_val);
alpar@1:       xfree(csa->bbar);
alpar@1:       xfree(csa->cbar);
alpar@1:       xfree(csa->refsp);
alpar@1:       xfree(csa->gamma);
alpar@1:       xfree(csa->tcol_ind);
alpar@1:       xfree(csa->tcol_vec);
alpar@1:       xfree(csa->trow_ind);
alpar@1:       xfree(csa->trow_vec);
alpar@1:       xfree(csa->work1);
alpar@1:       xfree(csa->work2);
alpar@1:       xfree(csa->work3);
alpar@1:       xfree(csa->work4);
alpar@1:       xfree(csa);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  spx_primal - core LP solver based on the primal simplex method
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  #include "glpspx.h"
alpar@1: *  int spx_primal(glp_prob *lp, const glp_smcp *parm);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine spx_primal is a core LP solver based on the two-phase
alpar@1: *  primal simplex method.
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  0  LP instance has been successfully solved.
alpar@1: *
alpar@1: *  GLP_EITLIM
alpar@1: *     Iteration limit has been exhausted.
alpar@1: *
alpar@1: *  GLP_ETMLIM
alpar@1: *     Time limit has been exhausted.
alpar@1: *
alpar@1: *  GLP_EFAIL
alpar@1: *     The solver failed to solve LP instance. */
alpar@1: 
alpar@1: int spx_primal(glp_prob *lp, const glp_smcp *parm)
alpar@1: {     struct csa *csa;
alpar@1:       int binv_st = 2;
alpar@1:       /* status of basis matrix factorization:
alpar@1:          0 - invalid; 1 - just computed; 2 - updated */
alpar@1:       int bbar_st = 0;
alpar@1:       /* status of primal values of basic variables:
alpar@1:          0 - invalid; 1 - just computed; 2 - updated */
alpar@1:       int cbar_st = 0;
alpar@1:       /* status of reduced costs of non-basic variables:
alpar@1:          0 - invalid; 1 - just computed; 2 - updated */
alpar@1:       int rigorous = 0;
alpar@1:       /* rigorous mode flag; this flag is used to enable iterative
alpar@1:          refinement on computing pivot rows and columns of the simplex
alpar@1:          table */
alpar@1:       int check = 0;
alpar@1:       int p_stat, d_stat, ret;
alpar@1:       /* allocate and initialize the common storage area */
alpar@1:       csa = alloc_csa(lp);
alpar@1:       init_csa(csa, lp);
alpar@1:       if (parm->msg_lev >= GLP_MSG_DBG)
alpar@1:          xprintf("Objective scale factor = %g\n", csa->zeta);
alpar@1: loop: /* main loop starts here */
alpar@1:       /* compute factorization of the basis matrix */
alpar@1:       if (binv_st == 0)
alpar@1:       {  ret = invert_B(csa);
alpar@1:          if (ret != 0)
alpar@1:          {  if (parm->msg_lev >= GLP_MSG_ERR)
alpar@1:             {  xprintf("Error: unable to factorize the basis matrix (%d"
alpar@1:                   ")\n", ret);
alpar@1:                xprintf("Sorry, basis recovery procedure not implemented"
alpar@1:                   " yet\n");
alpar@1:             }
alpar@1:             xassert(!lp->valid && lp->bfd == NULL);
alpar@1:             lp->bfd = csa->bfd, csa->bfd = NULL;
alpar@1:             lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
alpar@1:             lp->obj_val = 0.0;
alpar@1:             lp->it_cnt = csa->it_cnt;
alpar@1:             lp->some = 0;
alpar@1:             ret = GLP_EFAIL;
alpar@1:             goto done;
alpar@1:          }
alpar@1:          csa->valid = 1;
alpar@1:          binv_st = 1; /* just computed */
alpar@1:          /* invalidate basic solution components */
alpar@1:          bbar_st = cbar_st = 0;
alpar@1:       }
alpar@1:       /* compute primal values of basic variables */
alpar@1:       if (bbar_st == 0)
alpar@1:       {  eval_bbar(csa);
alpar@1:          bbar_st = 1; /* just computed */
alpar@1:          /* determine the search phase, if not determined yet */
alpar@1:          if (csa->phase == 0)
alpar@1:          {  if (set_aux_obj(csa, parm->tol_bnd) > 0)
alpar@1:             {  /* current basic solution is primal infeasible */
alpar@1:                /* start to minimize the sum of infeasibilities */
alpar@1:                csa->phase = 1;
alpar@1:             }
alpar@1:             else
alpar@1:             {  /* current basic solution is primal feasible */
alpar@1:                /* start to minimize the original objective function */
alpar@1:                set_orig_obj(csa);
alpar@1:                csa->phase = 2;
alpar@1:             }
alpar@1:             xassert(check_stab(csa, parm->tol_bnd) == 0);
alpar@1:             /* working objective coefficients have been changed, so
alpar@1:                invalidate reduced costs */
alpar@1:             cbar_st = 0;
alpar@1:             display(csa, parm, 1);
alpar@1:          }
alpar@1:          /* make sure that the current basic solution remains primal
alpar@1:             feasible (or pseudo feasible on phase I) */
alpar@1:          if (check_stab(csa, parm->tol_bnd))
alpar@1:          {  /* there are excessive bound violations due to round-off
alpar@1:                errors */
alpar@1:             if (parm->msg_lev >= GLP_MSG_ERR)
alpar@1:                xprintf("Warning: numerical instability (primal simplex,"
alpar@1:                   " phase %s)\n", csa->phase == 1 ? "I" : "II");
alpar@1:             /* restart the search */
alpar@1:             csa->phase = 0;
alpar@1:             binv_st = 0;
alpar@1:             rigorous = 5;
alpar@1:             goto loop;
alpar@1:          }
alpar@1:       }
alpar@1:       xassert(csa->phase == 1 || csa->phase == 2);
alpar@1:       /* on phase I we do not need to wait until the current basic
alpar@1:          solution becomes dual feasible; it is sufficient to make sure
alpar@1:          that no basic variable violates its bounds */
alpar@1:       if (csa->phase == 1 && !check_feas(csa, parm->tol_bnd))
alpar@1:       {  /* the current basis is primal feasible; switch to phase II */
alpar@1:          csa->phase = 2;
alpar@1:          set_orig_obj(csa);
alpar@1:          cbar_st = 0;
alpar@1:          display(csa, parm, 1);
alpar@1:       }
alpar@1:       /* compute reduced costs of non-basic variables */
alpar@1:       if (cbar_st == 0)
alpar@1:       {  eval_cbar(csa);
alpar@1:          cbar_st = 1; /* just computed */
alpar@1:       }
alpar@1:       /* redefine the reference space, if required */
alpar@1:       switch (parm->pricing)
alpar@1:       {  case GLP_PT_STD:
alpar@1:             break;
alpar@1:          case GLP_PT_PSE:
alpar@1:             if (csa->refct == 0) reset_refsp(csa);
alpar@1:             break;
alpar@1:          default:
alpar@1:             xassert(parm != parm);
alpar@1:       }
alpar@1:       /* at this point the basis factorization and all basic solution
alpar@1:          components are valid */
alpar@1:       xassert(binv_st && bbar_st && cbar_st);
alpar@1:       /* check accuracy of current basic solution components (only for
alpar@1:          debugging) */
alpar@1:       if (check)
alpar@1:       {  double e_bbar = err_in_bbar(csa);
alpar@1:          double e_cbar = err_in_cbar(csa);
alpar@1:          double e_gamma =
alpar@1:             (parm->pricing == GLP_PT_PSE ? err_in_gamma(csa) : 0.0);
alpar@1:          xprintf("e_bbar = %10.3e; e_cbar = %10.3e; e_gamma = %10.3e\n",
alpar@1:             e_bbar, e_cbar, e_gamma);
alpar@1:          xassert(e_bbar <= 1e-5 && e_cbar <= 1e-5 && e_gamma <= 1e-3);
alpar@1:       }
alpar@1:       /* check if the iteration limit has been exhausted */
alpar@1:       if (parm->it_lim < INT_MAX &&
alpar@1:           csa->it_cnt - csa->it_beg >= parm->it_lim)
alpar@1:       {  if (bbar_st != 1 || csa->phase == 2 && cbar_st != 1)
alpar@1:          {  if (bbar_st != 1) bbar_st = 0;
alpar@1:             if (csa->phase == 2 && cbar_st != 1) cbar_st = 0;
alpar@1:             goto loop;
alpar@1:          }
alpar@1:          display(csa, parm, 1);
alpar@1:          if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:             xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n");
alpar@1:          switch (csa->phase)
alpar@1:          {  case 1:
alpar@1:                p_stat = GLP_INFEAS;
alpar@1:                set_orig_obj(csa);
alpar@1:                eval_cbar(csa);
alpar@1:                break;
alpar@1:             case 2:
alpar@1:                p_stat = GLP_FEAS;
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(csa != csa);
alpar@1:          }
alpar@1:          chuzc(csa, parm->tol_dj);
alpar@1:          d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS);
alpar@1:          store_sol(csa, lp, p_stat, d_stat, 0);
alpar@1:          ret = GLP_EITLIM;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* check if the time limit has been exhausted */
alpar@1:       if (parm->tm_lim < INT_MAX &&
alpar@1:           1000.0 * xdifftime(xtime(), csa->tm_beg) >= parm->tm_lim)
alpar@1:       {  if (bbar_st != 1 || csa->phase == 2 && cbar_st != 1)
alpar@1:          {  if (bbar_st != 1) bbar_st = 0;
alpar@1:             if (csa->phase == 2 && cbar_st != 1) cbar_st = 0;
alpar@1:             goto loop;
alpar@1:          }
alpar@1:          display(csa, parm, 1);
alpar@1:          if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:             xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
alpar@1:          switch (csa->phase)
alpar@1:          {  case 1:
alpar@1:                p_stat = GLP_INFEAS;
alpar@1:                set_orig_obj(csa);
alpar@1:                eval_cbar(csa);
alpar@1:                break;
alpar@1:             case 2:
alpar@1:                p_stat = GLP_FEAS;
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(csa != csa);
alpar@1:          }
alpar@1:          chuzc(csa, parm->tol_dj);
alpar@1:          d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS);
alpar@1:          store_sol(csa, lp, p_stat, d_stat, 0);
alpar@1:          ret = GLP_ETMLIM;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* display the search progress */
alpar@1:       display(csa, parm, 0);
alpar@1:       /* choose non-basic variable xN[q] */
alpar@1:       chuzc(csa, parm->tol_dj);
alpar@1:       if (csa->q == 0)
alpar@1:       {  if (bbar_st != 1 || cbar_st != 1)
alpar@1:          {  if (bbar_st != 1) bbar_st = 0;
alpar@1:             if (cbar_st != 1) cbar_st = 0;
alpar@1:             goto loop;
alpar@1:          }
alpar@1:          display(csa, parm, 1);
alpar@1:          switch (csa->phase)
alpar@1:          {  case 1:
alpar@1:                if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:                   xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
alpar@1:                p_stat = GLP_NOFEAS;
alpar@1:                set_orig_obj(csa);
alpar@1:                eval_cbar(csa);
alpar@1:                chuzc(csa, parm->tol_dj);
alpar@1:                d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS);
alpar@1:                break;
alpar@1:             case 2:
alpar@1:                if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:                   xprintf("OPTIMAL SOLUTION FOUND\n");
alpar@1:                p_stat = d_stat = GLP_FEAS;
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(csa != csa);
alpar@1:          }
alpar@1:          store_sol(csa, lp, p_stat, d_stat, 0);
alpar@1:          ret = 0;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* compute pivot column of the simplex table */
alpar@1:       eval_tcol(csa);
alpar@1:       if (rigorous) refine_tcol(csa);
alpar@1:       sort_tcol(csa, parm->tol_piv);
alpar@1:       /* check accuracy of the reduced cost of xN[q] */
alpar@1:       {  double d1 = csa->cbar[csa->q]; /* less accurate */
alpar@1:          double d2 = reeval_cost(csa);  /* more accurate */
alpar@1:          xassert(d1 != 0.0);
alpar@1:          if (fabs(d1 - d2) > 1e-5 * (1.0 + fabs(d2)) ||
alpar@1:              !(d1 < 0.0 && d2 < 0.0 || d1 > 0.0 && d2 > 0.0))
alpar@1:          {  if (parm->msg_lev >= GLP_MSG_DBG)
alpar@1:                xprintf("d1 = %.12g; d2 = %.12g\n", d1, d2);
alpar@1:             if (cbar_st != 1 || !rigorous)
alpar@1:             {  if (cbar_st != 1) cbar_st = 0;
alpar@1:                rigorous = 5;
alpar@1:                goto loop;
alpar@1:             }
alpar@1:          }
alpar@1:          /* replace cbar[q] by more accurate value keeping its sign */
alpar@1:          if (d1 > 0.0)
alpar@1:             csa->cbar[csa->q] = (d2 > 0.0 ? d2 : +DBL_EPSILON);
alpar@1:          else
alpar@1:             csa->cbar[csa->q] = (d2 < 0.0 ? d2 : -DBL_EPSILON);
alpar@1:       }
alpar@1:       /* choose basic variable xB[p] */
alpar@1:       switch (parm->r_test)
alpar@1:       {  case GLP_RT_STD:
alpar@1:             chuzr(csa, 0.0);
alpar@1:             break;
alpar@1:          case GLP_RT_HAR:
alpar@1:             chuzr(csa, 0.30 * parm->tol_bnd);
alpar@1:             break;
alpar@1:          default:
alpar@1:             xassert(parm != parm);
alpar@1:       }
alpar@1:       if (csa->p == 0)
alpar@1:       {  if (bbar_st != 1 || cbar_st != 1 || !rigorous)
alpar@1:          {  if (bbar_st != 1) bbar_st = 0;
alpar@1:             if (cbar_st != 1) cbar_st = 0;
alpar@1:             rigorous = 1;
alpar@1:             goto loop;
alpar@1:          }
alpar@1:          display(csa, parm, 1);
alpar@1:          switch (csa->phase)
alpar@1:          {  case 1:
alpar@1:                if (parm->msg_lev >= GLP_MSG_ERR)
alpar@1:                   xprintf("Error: unable to choose basic variable on ph"
alpar@1:                      "ase I\n");
alpar@1:                xassert(!lp->valid && lp->bfd == NULL);
alpar@1:                lp->bfd = csa->bfd, csa->bfd = NULL;
alpar@1:                lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
alpar@1:                lp->obj_val = 0.0;
alpar@1:                lp->it_cnt = csa->it_cnt;
alpar@1:                lp->some = 0;
alpar@1:                ret = GLP_EFAIL;
alpar@1:                break;
alpar@1:             case 2:
alpar@1:                if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:                   xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
alpar@1:                store_sol(csa, lp, GLP_FEAS, GLP_NOFEAS,
alpar@1:                   csa->head[csa->m+csa->q]);
alpar@1:                ret = 0;
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(csa != csa);
alpar@1:          }
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* check if the pivot element is acceptable */
alpar@1:       if (csa->p > 0)
alpar@1:       {  double piv = csa->tcol_vec[csa->p];
alpar@1:          double eps = 1e-5 * (1.0 + 0.01 * csa->tcol_max);
alpar@1:          if (fabs(piv) < eps)
alpar@1:          {  if (parm->msg_lev >= GLP_MSG_DBG)
alpar@1:                xprintf("piv = %.12g; eps = %g\n", piv, eps);
alpar@1:             if (!rigorous)
alpar@1:             {  rigorous = 5;
alpar@1:                goto loop;
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       /* now xN[q] and xB[p] have been chosen anyhow */
alpar@1:       /* compute pivot row of the simplex table */
alpar@1:       if (csa->p > 0)
alpar@1:       {  double *rho = csa->work4;
alpar@1:          eval_rho(csa, rho);
alpar@1:          if (rigorous) refine_rho(csa, rho);
alpar@1:          eval_trow(csa, rho);
alpar@1:       }
alpar@1:       /* accuracy check based on the pivot element */
alpar@1:       if (csa->p > 0)
alpar@1:       {  double piv1 = csa->tcol_vec[csa->p]; /* more accurate */
alpar@1:          double piv2 = csa->trow_vec[csa->q]; /* less accurate */
alpar@1:          xassert(piv1 != 0.0);
alpar@1:          if (fabs(piv1 - piv2) > 1e-8 * (1.0 + fabs(piv1)) ||
alpar@1:              !(piv1 > 0.0 && piv2 > 0.0 || piv1 < 0.0 && piv2 < 0.0))
alpar@1:          {  if (parm->msg_lev >= GLP_MSG_DBG)
alpar@1:                xprintf("piv1 = %.12g; piv2 = %.12g\n", piv1, piv2);
alpar@1:             if (binv_st != 1 || !rigorous)
alpar@1:             {  if (binv_st != 1) binv_st = 0;
alpar@1:                rigorous = 5;
alpar@1:                goto loop;
alpar@1:             }
alpar@1:             /* use more accurate version in the pivot row */
alpar@1:             if (csa->trow_vec[csa->q] == 0.0)
alpar@1:             {  csa->trow_nnz++;
alpar@1:                xassert(csa->trow_nnz <= csa->n);
alpar@1:                csa->trow_ind[csa->trow_nnz] = csa->q;
alpar@1:             }
alpar@1:             csa->trow_vec[csa->q] = piv1;
alpar@1:          }
alpar@1:       }
alpar@1:       /* update primal values of basic variables */
alpar@1:       update_bbar(csa);
alpar@1:       bbar_st = 2; /* updated */
alpar@1:       /* update reduced costs of non-basic variables */
alpar@1:       if (csa->p > 0)
alpar@1:       {  update_cbar(csa);
alpar@1:          cbar_st = 2; /* updated */
alpar@1:          /* on phase I objective coefficient of xB[p] in the adjacent
alpar@1:             basis becomes zero */
alpar@1:          if (csa->phase == 1)
alpar@1:          {  int k = csa->head[csa->p]; /* x[k] = xB[p] -> xN[q] */
alpar@1:             csa->cbar[csa->q] -= csa->coef[k];
alpar@1:             csa->coef[k] = 0.0;
alpar@1:          }
alpar@1:       }
alpar@1:       /* update steepest edge coefficients */
alpar@1:       if (csa->p > 0)
alpar@1:       {  switch (parm->pricing)
alpar@1:          {  case GLP_PT_STD:
alpar@1:                break;
alpar@1:             case GLP_PT_PSE:
alpar@1:                if (csa->refct > 0) update_gamma(csa);
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(parm != parm);
alpar@1:          }
alpar@1:       }
alpar@1:       /* update factorization of the basis matrix */
alpar@1:       if (csa->p > 0)
alpar@1:       {  ret = update_B(csa, csa->p, csa->head[csa->m+csa->q]);
alpar@1:          if (ret == 0)
alpar@1:             binv_st = 2; /* updated */
alpar@1:          else
alpar@1:          {  csa->valid = 0;
alpar@1:             binv_st = 0; /* invalid */
alpar@1:          }
alpar@1:       }
alpar@1:       /* update matrix N */
alpar@1:       if (csa->p > 0)
alpar@1:       {  del_N_col(csa, csa->q, csa->head[csa->m+csa->q]);
alpar@1:          if (csa->type[csa->head[csa->p]] != GLP_FX)
alpar@1:             add_N_col(csa, csa->q, csa->head[csa->p]);
alpar@1:       }
alpar@1:       /* change the basis header */
alpar@1:       change_basis(csa);
alpar@1:       /* iteration complete */
alpar@1:       csa->it_cnt++;
alpar@1:       if (rigorous > 0) rigorous--;
alpar@1:       goto loop;
alpar@1: done: /* deallocate the common storage area */
alpar@1:       free_csa(csa);
alpar@1:       /* return to the calling program */
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /* eof */