1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpspx01.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,2954 @@
     1.4 +/* glpspx01.c (primal simplex method) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpspx.h"
    1.29 +
    1.30 +struct csa
    1.31 +{     /* common storage area */
    1.32 +      /*--------------------------------------------------------------*/
    1.33 +      /* LP data */
    1.34 +      int m;
    1.35 +      /* number of rows (auxiliary variables), m > 0 */
    1.36 +      int n;
    1.37 +      /* number of columns (structural variables), n > 0 */
    1.38 +      char *type; /* char type[1+m+n]; */
    1.39 +      /* type[0] is not used;
    1.40 +         type[k], 1 <= k <= m+n, is the type of variable x[k]:
    1.41 +         GLP_FR - free variable
    1.42 +         GLP_LO - variable with lower bound
    1.43 +         GLP_UP - variable with upper bound
    1.44 +         GLP_DB - double-bounded variable
    1.45 +         GLP_FX - fixed variable */
    1.46 +      double *lb; /* double lb[1+m+n]; */
    1.47 +      /* lb[0] is not used;
    1.48 +         lb[k], 1 <= k <= m+n, is an lower bound of variable x[k];
    1.49 +         if x[k] has no lower bound, lb[k] is zero */
    1.50 +      double *ub; /* double ub[1+m+n]; */
    1.51 +      /* ub[0] is not used;
    1.52 +         ub[k], 1 <= k <= m+n, is an upper bound of variable x[k];
    1.53 +         if x[k] has no upper bound, ub[k] is zero;
    1.54 +         if x[k] is of fixed type, ub[k] is the same as lb[k] */
    1.55 +      double *coef; /* double coef[1+m+n]; */
    1.56 +      /* coef[0] is not used;
    1.57 +         coef[k], 1 <= k <= m+n, is an objective coefficient at
    1.58 +         variable x[k] (note that on phase I auxiliary variables also
    1.59 +         may have non-zero objective coefficients) */
    1.60 +      /*--------------------------------------------------------------*/
    1.61 +      /* original objective function */
    1.62 +      double *obj; /* double obj[1+n]; */
    1.63 +      /* obj[0] is a constant term of the original objective function;
    1.64 +         obj[j], 1 <= j <= n, is an original objective coefficient at
    1.65 +         structural variable x[m+j] */
    1.66 +      double zeta;
    1.67 +      /* factor used to scale original objective coefficients; its
    1.68 +         sign defines original optimization direction: zeta > 0 means
    1.69 +         minimization, zeta < 0 means maximization */
    1.70 +      /*--------------------------------------------------------------*/
    1.71 +      /* constraint matrix A; it has m rows and n columns and is stored
    1.72 +         by columns */
    1.73 +      int *A_ptr; /* int A_ptr[1+n+1]; */
    1.74 +      /* A_ptr[0] is not used;
    1.75 +         A_ptr[j], 1 <= j <= n, is starting position of j-th column in
    1.76 +         arrays A_ind and A_val; note that A_ptr[1] is always 1;
    1.77 +         A_ptr[n+1] indicates the position after the last element in
    1.78 +         arrays A_ind and A_val */
    1.79 +      int *A_ind; /* int A_ind[A_ptr[n+1]]; */
    1.80 +      /* row indices */
    1.81 +      double *A_val; /* double A_val[A_ptr[n+1]]; */
    1.82 +      /* non-zero element values */
    1.83 +      /*--------------------------------------------------------------*/
    1.84 +      /* basis header */
    1.85 +      int *head; /* int head[1+m+n]; */
    1.86 +      /* head[0] is not used;
    1.87 +         head[i], 1 <= i <= m, is the ordinal number of basic variable
    1.88 +         xB[i]; head[i] = k means that xB[i] = x[k] and i-th column of
    1.89 +         matrix B is k-th column of matrix (I|-A);
    1.90 +         head[m+j], 1 <= j <= n, is the ordinal number of non-basic
    1.91 +         variable xN[j]; head[m+j] = k means that xN[j] = x[k] and j-th
    1.92 +         column of matrix N is k-th column of matrix (I|-A) */
    1.93 +      char *stat; /* char stat[1+n]; */
    1.94 +      /* stat[0] is not used;
    1.95 +         stat[j], 1 <= j <= n, is the status of non-basic variable
    1.96 +         xN[j], which defines its active bound:
    1.97 +         GLP_NL - lower bound is active
    1.98 +         GLP_NU - upper bound is active
    1.99 +         GLP_NF - free variable
   1.100 +         GLP_NS - fixed variable */
   1.101 +      /*--------------------------------------------------------------*/
   1.102 +      /* matrix B is the basis matrix; it is composed from columns of
   1.103 +         the augmented constraint matrix (I|-A) corresponding to basic
   1.104 +         variables and stored in a factorized (invertable) form */
   1.105 +      int valid;
   1.106 +      /* factorization is valid only if this flag is set */
   1.107 +      BFD *bfd; /* BFD bfd[1:m,1:m]; */
   1.108 +      /* factorized (invertable) form of the basis matrix */
   1.109 +      /*--------------------------------------------------------------*/
   1.110 +      /* matrix N is a matrix composed from columns of the augmented
   1.111 +         constraint matrix (I|-A) corresponding to non-basic variables
   1.112 +         except fixed ones; it is stored by rows and changes every time
   1.113 +         the basis changes */
   1.114 +      int *N_ptr; /* int N_ptr[1+m+1]; */
   1.115 +      /* N_ptr[0] is not used;
   1.116 +         N_ptr[i], 1 <= i <= m, is starting position of i-th row in
   1.117 +         arrays N_ind and N_val; note that N_ptr[1] is always 1;
   1.118 +         N_ptr[m+1] indicates the position after the last element in
   1.119 +         arrays N_ind and N_val */
   1.120 +      int *N_len; /* int N_len[1+m]; */
   1.121 +      /* N_len[0] is not used;
   1.122 +         N_len[i], 1 <= i <= m, is length of i-th row (0 to n) */
   1.123 +      int *N_ind; /* int N_ind[N_ptr[m+1]]; */
   1.124 +      /* column indices */
   1.125 +      double *N_val; /* double N_val[N_ptr[m+1]]; */
   1.126 +      /* non-zero element values */
   1.127 +      /*--------------------------------------------------------------*/
   1.128 +      /* working parameters */
   1.129 +      int phase;
   1.130 +      /* search phase:
   1.131 +         0 - not determined yet
   1.132 +         1 - search for primal feasible solution
   1.133 +         2 - search for optimal solution */
   1.134 +      glp_long tm_beg;
   1.135 +      /* time value at the beginning of the search */
   1.136 +      int it_beg;
   1.137 +      /* simplex iteration count at the beginning of the search */
   1.138 +      int it_cnt;
   1.139 +      /* simplex iteration count; it increases by one every time the
   1.140 +         basis changes (including the case when a non-basic variable
   1.141 +         jumps to its opposite bound) */
   1.142 +      int it_dpy;
   1.143 +      /* simplex iteration count at the most recent display output */
   1.144 +      /*--------------------------------------------------------------*/
   1.145 +      /* basic solution components */
   1.146 +      double *bbar; /* double bbar[1+m]; */
   1.147 +      /* bbar[0] is not used;
   1.148 +         bbar[i], 1 <= i <= m, is primal value of basic variable xB[i]
   1.149 +         (if xB[i] is free, its primal value is not updated) */
   1.150 +      double *cbar; /* double cbar[1+n]; */
   1.151 +      /* cbar[0] is not used;
   1.152 +         cbar[j], 1 <= j <= n, is reduced cost of non-basic variable
   1.153 +         xN[j] (if xN[j] is fixed, its reduced cost is not updated) */
   1.154 +      /*--------------------------------------------------------------*/
   1.155 +      /* the following pricing technique options may be used:
   1.156 +         GLP_PT_STD - standard ("textbook") pricing;
   1.157 +         GLP_PT_PSE - projected steepest edge;
   1.158 +         GLP_PT_DVX - Devex pricing (not implemented yet);
   1.159 +         in case of GLP_PT_STD the reference space is not used, and all
   1.160 +         steepest edge coefficients are set to 1 */
   1.161 +      int refct;
   1.162 +      /* this count is set to an initial value when the reference space
   1.163 +         is defined and decreases by one every time the basis changes;
   1.164 +         once this count reaches zero, the reference space is redefined
   1.165 +         again */
   1.166 +      char *refsp; /* char refsp[1+m+n]; */
   1.167 +      /* refsp[0] is not used;
   1.168 +         refsp[k], 1 <= k <= m+n, is the flag which means that variable
   1.169 +         x[k] belongs to the current reference space */
   1.170 +      double *gamma; /* double gamma[1+n]; */
   1.171 +      /* gamma[0] is not used;
   1.172 +         gamma[j], 1 <= j <= n, is the steepest edge coefficient for
   1.173 +         non-basic variable xN[j]; if xN[j] is fixed, gamma[j] is not
   1.174 +         used and just set to 1 */
   1.175 +      /*--------------------------------------------------------------*/
   1.176 +      /* non-basic variable xN[q] chosen to enter the basis */
   1.177 +      int q;
   1.178 +      /* index of the non-basic variable xN[q] chosen, 1 <= q <= n;
   1.179 +         if the set of eligible non-basic variables is empty and thus
   1.180 +         no variable has been chosen, q is set to 0 */
   1.181 +      /*--------------------------------------------------------------*/
   1.182 +      /* pivot column of the simplex table corresponding to non-basic
   1.183 +         variable xN[q] chosen is the following vector:
   1.184 +            T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q],
   1.185 +         where B is the current basis matrix, N[q] is a column of the
   1.186 +         matrix (I|-A) corresponding to xN[q] */
   1.187 +      int tcol_nnz;
   1.188 +      /* number of non-zero components, 0 <= nnz <= m */
   1.189 +      int *tcol_ind; /* int tcol_ind[1+m]; */
   1.190 +      /* tcol_ind[0] is not used;
   1.191 +         tcol_ind[t], 1 <= t <= nnz, is an index of non-zero component,
   1.192 +         i.e. tcol_ind[t] = i means that tcol_vec[i] != 0 */
   1.193 +      double *tcol_vec; /* double tcol_vec[1+m]; */
   1.194 +      /* tcol_vec[0] is not used;
   1.195 +         tcol_vec[i], 1 <= i <= m, is a numeric value of i-th component
   1.196 +         of the column */
   1.197 +      double tcol_max;
   1.198 +      /* infinity (maximum) norm of the column (max |tcol_vec[i]|) */
   1.199 +      int tcol_num;
   1.200 +      /* number of significant non-zero components, which means that:
   1.201 +         |tcol_vec[i]| >= eps for i in tcol_ind[1,...,num],
   1.202 +         |tcol_vec[i]| <  eps for i in tcol_ind[num+1,...,nnz],
   1.203 +         where eps is a pivot tolerance */
   1.204 +      /*--------------------------------------------------------------*/
   1.205 +      /* basic variable xB[p] chosen to leave the basis */
   1.206 +      int p;
   1.207 +      /* index of the basic variable xB[p] chosen, 1 <= p <= m;
   1.208 +         p = 0 means that no basic variable reaches its bound;
   1.209 +         p < 0 means that non-basic variable xN[q] reaches its opposite
   1.210 +         bound before any basic variable */
   1.211 +      int p_stat;
   1.212 +      /* new status (GLP_NL, GLP_NU, or GLP_NS) to be assigned to xB[p]
   1.213 +         once it has left the basis */
   1.214 +      double teta;
   1.215 +      /* change of non-basic variable xN[q] (see above), on which xB[p]
   1.216 +         (or, if p < 0, xN[q] itself) reaches its bound */
   1.217 +      /*--------------------------------------------------------------*/
   1.218 +      /* pivot row of the simplex table corresponding to basic variable
   1.219 +         xB[p] chosen is the following vector:
   1.220 +            T' * e[p] = - N' * inv(B') * e[p] = - N' * rho,
   1.221 +         where B' is a matrix transposed to the current basis matrix,
   1.222 +         N' is a matrix, whose rows are columns of the matrix (I|-A)
   1.223 +         corresponding to non-basic non-fixed variables */
   1.224 +      int trow_nnz;
   1.225 +      /* number of non-zero components, 0 <= nnz <= n */
   1.226 +      int *trow_ind; /* int trow_ind[1+n]; */
   1.227 +      /* trow_ind[0] is not used;
   1.228 +         trow_ind[t], 1 <= t <= nnz, is an index of non-zero component,
   1.229 +         i.e. trow_ind[t] = j means that trow_vec[j] != 0 */
   1.230 +      double *trow_vec; /* int trow_vec[1+n]; */
   1.231 +      /* trow_vec[0] is not used;
   1.232 +         trow_vec[j], 1 <= j <= n, is a numeric value of j-th component
   1.233 +         of the row */
   1.234 +      /*--------------------------------------------------------------*/
   1.235 +      /* working arrays */
   1.236 +      double *work1; /* double work1[1+m]; */
   1.237 +      double *work2; /* double work2[1+m]; */
   1.238 +      double *work3; /* double work3[1+m]; */
   1.239 +      double *work4; /* double work4[1+m]; */
   1.240 +};
   1.241 +
   1.242 +static const double kappa = 0.10;
   1.243 +
   1.244 +/***********************************************************************
   1.245 +*  alloc_csa - allocate common storage area
   1.246 +*
   1.247 +*  This routine allocates all arrays in the common storage area (CSA)
   1.248 +*  and returns a pointer to the CSA. */
   1.249 +
   1.250 +static struct csa *alloc_csa(glp_prob *lp)
   1.251 +{     struct csa *csa;
   1.252 +      int m = lp->m;
   1.253 +      int n = lp->n;
   1.254 +      int nnz = lp->nnz;
   1.255 +      csa = xmalloc(sizeof(struct csa));
   1.256 +      xassert(m > 0 && n > 0);
   1.257 +      csa->m = m;
   1.258 +      csa->n = n;
   1.259 +      csa->type = xcalloc(1+m+n, sizeof(char));
   1.260 +      csa->lb = xcalloc(1+m+n, sizeof(double));
   1.261 +      csa->ub = xcalloc(1+m+n, sizeof(double));
   1.262 +      csa->coef = xcalloc(1+m+n, sizeof(double));
   1.263 +      csa->obj = xcalloc(1+n, sizeof(double));
   1.264 +      csa->A_ptr = xcalloc(1+n+1, sizeof(int));
   1.265 +      csa->A_ind = xcalloc(1+nnz, sizeof(int));
   1.266 +      csa->A_val = xcalloc(1+nnz, sizeof(double));
   1.267 +      csa->head = xcalloc(1+m+n, sizeof(int));
   1.268 +      csa->stat = xcalloc(1+n, sizeof(char));
   1.269 +      csa->N_ptr = xcalloc(1+m+1, sizeof(int));
   1.270 +      csa->N_len = xcalloc(1+m, sizeof(int));
   1.271 +      csa->N_ind = NULL; /* will be allocated later */
   1.272 +      csa->N_val = NULL; /* will be allocated later */
   1.273 +      csa->bbar = xcalloc(1+m, sizeof(double));
   1.274 +      csa->cbar = xcalloc(1+n, sizeof(double));
   1.275 +      csa->refsp = xcalloc(1+m+n, sizeof(char));
   1.276 +      csa->gamma = xcalloc(1+n, sizeof(double));
   1.277 +      csa->tcol_ind = xcalloc(1+m, sizeof(int));
   1.278 +      csa->tcol_vec = xcalloc(1+m, sizeof(double));
   1.279 +      csa->trow_ind = xcalloc(1+n, sizeof(int));
   1.280 +      csa->trow_vec = xcalloc(1+n, sizeof(double));
   1.281 +      csa->work1 = xcalloc(1+m, sizeof(double));
   1.282 +      csa->work2 = xcalloc(1+m, sizeof(double));
   1.283 +      csa->work3 = xcalloc(1+m, sizeof(double));
   1.284 +      csa->work4 = xcalloc(1+m, sizeof(double));
   1.285 +      return csa;
   1.286 +}
   1.287 +
   1.288 +/***********************************************************************
   1.289 +*  init_csa - initialize common storage area
   1.290 +*
   1.291 +*  This routine initializes all data structures in the common storage
   1.292 +*  area (CSA). */
   1.293 +
   1.294 +static void alloc_N(struct csa *csa);
   1.295 +static void build_N(struct csa *csa);
   1.296 +
   1.297 +static void init_csa(struct csa *csa, glp_prob *lp)
   1.298 +{     int m = csa->m;
   1.299 +      int n = csa->n;
   1.300 +      char *type = csa->type;
   1.301 +      double *lb = csa->lb;
   1.302 +      double *ub = csa->ub;
   1.303 +      double *coef = csa->coef;
   1.304 +      double *obj = csa->obj;
   1.305 +      int *A_ptr = csa->A_ptr;
   1.306 +      int *A_ind = csa->A_ind;
   1.307 +      double *A_val = csa->A_val;
   1.308 +      int *head = csa->head;
   1.309 +      char *stat = csa->stat;
   1.310 +      char *refsp = csa->refsp;
   1.311 +      double *gamma = csa->gamma;
   1.312 +      int i, j, k, loc;
   1.313 +      double cmax;
   1.314 +      /* auxiliary variables */
   1.315 +      for (i = 1; i <= m; i++)
   1.316 +      {  GLPROW *row = lp->row[i];
   1.317 +         type[i] = (char)row->type;
   1.318 +         lb[i] = row->lb * row->rii;
   1.319 +         ub[i] = row->ub * row->rii;
   1.320 +         coef[i] = 0.0;
   1.321 +      }
   1.322 +      /* structural variables */
   1.323 +      for (j = 1; j <= n; j++)
   1.324 +      {  GLPCOL *col = lp->col[j];
   1.325 +         type[m+j] = (char)col->type;
   1.326 +         lb[m+j] = col->lb / col->sjj;
   1.327 +         ub[m+j] = col->ub / col->sjj;
   1.328 +         coef[m+j] = col->coef * col->sjj;
   1.329 +      }
   1.330 +      /* original objective function */
   1.331 +      obj[0] = lp->c0;
   1.332 +      memcpy(&obj[1], &coef[m+1], n * sizeof(double));
   1.333 +      /* factor used to scale original objective coefficients */
   1.334 +      cmax = 0.0;
   1.335 +      for (j = 1; j <= n; j++)
   1.336 +         if (cmax < fabs(obj[j])) cmax = fabs(obj[j]);
   1.337 +      if (cmax == 0.0) cmax = 1.0;
   1.338 +      switch (lp->dir)
   1.339 +      {  case GLP_MIN:
   1.340 +            csa->zeta = + 1.0 / cmax;
   1.341 +            break;
   1.342 +         case GLP_MAX:
   1.343 +            csa->zeta = - 1.0 / cmax;
   1.344 +            break;
   1.345 +         default:
   1.346 +            xassert(lp != lp);
   1.347 +      }
   1.348 +#if 1
   1.349 +      if (fabs(csa->zeta) < 1.0) csa->zeta *= 1000.0;
   1.350 +#endif
   1.351 +      /* matrix A (by columns) */
   1.352 +      loc = 1;
   1.353 +      for (j = 1; j <= n; j++)
   1.354 +      {  GLPAIJ *aij;
   1.355 +         A_ptr[j] = loc;
   1.356 +         for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
   1.357 +         {  A_ind[loc] = aij->row->i;
   1.358 +            A_val[loc] = aij->row->rii * aij->val * aij->col->sjj;
   1.359 +            loc++;
   1.360 +         }
   1.361 +      }
   1.362 +      A_ptr[n+1] = loc;
   1.363 +      xassert(loc == lp->nnz+1);
   1.364 +      /* basis header */
   1.365 +      xassert(lp->valid);
   1.366 +      memcpy(&head[1], &lp->head[1], m * sizeof(int));
   1.367 +      k = 0;
   1.368 +      for (i = 1; i <= m; i++)
   1.369 +      {  GLPROW *row = lp->row[i];
   1.370 +         if (row->stat != GLP_BS)
   1.371 +         {  k++;
   1.372 +            xassert(k <= n);
   1.373 +            head[m+k] = i;
   1.374 +            stat[k] = (char)row->stat;
   1.375 +         }
   1.376 +      }
   1.377 +      for (j = 1; j <= n; j++)
   1.378 +      {  GLPCOL *col = lp->col[j];
   1.379 +         if (col->stat != GLP_BS)
   1.380 +         {  k++;
   1.381 +            xassert(k <= n);
   1.382 +            head[m+k] = m + j;
   1.383 +            stat[k] = (char)col->stat;
   1.384 +         }
   1.385 +      }
   1.386 +      xassert(k == n);
   1.387 +      /* factorization of matrix B */
   1.388 +      csa->valid = 1, lp->valid = 0;
   1.389 +      csa->bfd = lp->bfd, lp->bfd = NULL;
   1.390 +      /* matrix N (by rows) */
   1.391 +      alloc_N(csa);
   1.392 +      build_N(csa);
   1.393 +      /* working parameters */
   1.394 +      csa->phase = 0;
   1.395 +      csa->tm_beg = xtime();
   1.396 +      csa->it_beg = csa->it_cnt = lp->it_cnt;
   1.397 +      csa->it_dpy = -1;
   1.398 +      /* reference space and steepest edge coefficients */
   1.399 +      csa->refct = 0;
   1.400 +      memset(&refsp[1], 0, (m+n) * sizeof(char));
   1.401 +      for (j = 1; j <= n; j++) gamma[j] = 1.0;
   1.402 +      return;
   1.403 +}
   1.404 +
   1.405 +/***********************************************************************
   1.406 +*  invert_B - compute factorization of the basis matrix
   1.407 +*
   1.408 +*  This routine computes factorization of the current basis matrix B.
   1.409 +*
   1.410 +*  If the operation is successful, the routine returns zero, otherwise
   1.411 +*  non-zero. */
   1.412 +
   1.413 +static int inv_col(void *info, int i, int ind[], double val[])
   1.414 +{     /* this auxiliary routine returns row indices and numeric values
   1.415 +         of non-zero elements of i-th column of the basis matrix */
   1.416 +      struct csa *csa = info;
   1.417 +      int m = csa->m;
   1.418 +#ifdef GLP_DEBUG
   1.419 +      int n = csa->n;
   1.420 +#endif
   1.421 +      int *A_ptr = csa->A_ptr;
   1.422 +      int *A_ind = csa->A_ind;
   1.423 +      double *A_val = csa->A_val;
   1.424 +      int *head = csa->head;
   1.425 +      int k, len, ptr, t;
   1.426 +#ifdef GLP_DEBUG
   1.427 +      xassert(1 <= i && i <= m);
   1.428 +#endif
   1.429 +      k = head[i]; /* B[i] is k-th column of (I|-A) */
   1.430 +#ifdef GLP_DEBUG
   1.431 +      xassert(1 <= k && k <= m+n);
   1.432 +#endif
   1.433 +      if (k <= m)
   1.434 +      {  /* B[i] is k-th column of submatrix I */
   1.435 +         len = 1;
   1.436 +         ind[1] = k;
   1.437 +         val[1] = 1.0;
   1.438 +      }
   1.439 +      else
   1.440 +      {  /* B[i] is (k-m)-th column of submatrix (-A) */
   1.441 +         ptr = A_ptr[k-m];
   1.442 +         len = A_ptr[k-m+1] - ptr;
   1.443 +         memcpy(&ind[1], &A_ind[ptr], len * sizeof(int));
   1.444 +         memcpy(&val[1], &A_val[ptr], len * sizeof(double));
   1.445 +         for (t = 1; t <= len; t++) val[t] = - val[t];
   1.446 +      }
   1.447 +      return len;
   1.448 +}
   1.449 +
   1.450 +static int invert_B(struct csa *csa)
   1.451 +{     int ret;
   1.452 +      ret = bfd_factorize(csa->bfd, csa->m, NULL, inv_col, csa);
   1.453 +      csa->valid = (ret == 0);
   1.454 +      return ret;
   1.455 +}
   1.456 +
   1.457 +/***********************************************************************
   1.458 +*  update_B - update factorization of the basis matrix
   1.459 +*
   1.460 +*  This routine replaces i-th column of the basis matrix B by k-th
   1.461 +*  column of the augmented constraint matrix (I|-A) and then updates
   1.462 +*  the factorization of B.
   1.463 +*
   1.464 +*  If the factorization has been successfully updated, the routine
   1.465 +*  returns zero, otherwise non-zero. */
   1.466 +
   1.467 +static int update_B(struct csa *csa, int i, int k)
   1.468 +{     int m = csa->m;
   1.469 +#ifdef GLP_DEBUG
   1.470 +      int n = csa->n;
   1.471 +#endif
   1.472 +      int ret;
   1.473 +#ifdef GLP_DEBUG
   1.474 +      xassert(1 <= i && i <= m);
   1.475 +      xassert(1 <= k && k <= m+n);
   1.476 +#endif
   1.477 +      if (k <= m)
   1.478 +      {  /* new i-th column of B is k-th column of I */
   1.479 +         int ind[1+1];
   1.480 +         double val[1+1];
   1.481 +         ind[1] = k;
   1.482 +         val[1] = 1.0;
   1.483 +         xassert(csa->valid);
   1.484 +         ret = bfd_update_it(csa->bfd, i, 0, 1, ind, val);
   1.485 +      }
   1.486 +      else
   1.487 +      {  /* new i-th column of B is (k-m)-th column of (-A) */
   1.488 +         int *A_ptr = csa->A_ptr;
   1.489 +         int *A_ind = csa->A_ind;
   1.490 +         double *A_val = csa->A_val;
   1.491 +         double *val = csa->work1;
   1.492 +         int beg, end, ptr, len;
   1.493 +         beg = A_ptr[k-m];
   1.494 +         end = A_ptr[k-m+1];
   1.495 +         len = 0;
   1.496 +         for (ptr = beg; ptr < end; ptr++)
   1.497 +            val[++len] = - A_val[ptr];
   1.498 +         xassert(csa->valid);
   1.499 +         ret = bfd_update_it(csa->bfd, i, 0, len, &A_ind[beg-1], val);
   1.500 +      }
   1.501 +      csa->valid = (ret == 0);
   1.502 +      return ret;
   1.503 +}
   1.504 +
   1.505 +/***********************************************************************
   1.506 +*  error_ftran - compute residual vector r = h - B * x
   1.507 +*
   1.508 +*  This routine computes the residual vector r = h - B * x, where B is
   1.509 +*  the current basis matrix, h is the vector of right-hand sides, x is
   1.510 +*  the solution vector. */
   1.511 +
   1.512 +static void error_ftran(struct csa *csa, double h[], double x[],
   1.513 +      double r[])
   1.514 +{     int m = csa->m;
   1.515 +#ifdef GLP_DEBUG
   1.516 +      int n = csa->n;
   1.517 +#endif
   1.518 +      int *A_ptr = csa->A_ptr;
   1.519 +      int *A_ind = csa->A_ind;
   1.520 +      double *A_val = csa->A_val;
   1.521 +      int *head = csa->head;
   1.522 +      int i, k, beg, end, ptr;
   1.523 +      double temp;
   1.524 +      /* compute the residual vector:
   1.525 +         r = h - B * x = h - B[1] * x[1] - ... - B[m] * x[m],
   1.526 +         where B[1], ..., B[m] are columns of matrix B */
   1.527 +      memcpy(&r[1], &h[1], m * sizeof(double));
   1.528 +      for (i = 1; i <= m; i++)
   1.529 +      {  temp = x[i];
   1.530 +         if (temp == 0.0) continue;
   1.531 +         k = head[i]; /* B[i] is k-th column of (I|-A) */
   1.532 +#ifdef GLP_DEBUG
   1.533 +         xassert(1 <= k && k <= m+n);
   1.534 +#endif
   1.535 +         if (k <= m)
   1.536 +         {  /* B[i] is k-th column of submatrix I */
   1.537 +            r[k] -= temp;
   1.538 +         }
   1.539 +         else
   1.540 +         {  /* B[i] is (k-m)-th column of submatrix (-A) */
   1.541 +            beg = A_ptr[k-m];
   1.542 +            end = A_ptr[k-m+1];
   1.543 +            for (ptr = beg; ptr < end; ptr++)
   1.544 +               r[A_ind[ptr]] += A_val[ptr] * temp;
   1.545 +         }
   1.546 +      }
   1.547 +      return;
   1.548 +}
   1.549 +
   1.550 +/***********************************************************************
   1.551 +*  refine_ftran - refine solution of B * x = h
   1.552 +*
   1.553 +*  This routine performs one iteration to refine the solution of
   1.554 +*  the system B * x = h, where B is the current basis matrix, h is the
   1.555 +*  vector of right-hand sides, x is the solution vector. */
   1.556 +
   1.557 +static void refine_ftran(struct csa *csa, double h[], double x[])
   1.558 +{     int m = csa->m;
   1.559 +      double *r = csa->work1;
   1.560 +      double *d = csa->work1;
   1.561 +      int i;
   1.562 +      /* compute the residual vector r = h - B * x */
   1.563 +      error_ftran(csa, h, x, r);
   1.564 +      /* compute the correction vector d = inv(B) * r */
   1.565 +      xassert(csa->valid);
   1.566 +      bfd_ftran(csa->bfd, d);
   1.567 +      /* refine the solution vector (new x) = (old x) + d */
   1.568 +      for (i = 1; i <= m; i++) x[i] += d[i];
   1.569 +      return;
   1.570 +}
   1.571 +
   1.572 +/***********************************************************************
   1.573 +*  error_btran - compute residual vector r = h - B'* x
   1.574 +*
   1.575 +*  This routine computes the residual vector r = h - B'* x, where B'
   1.576 +*  is a matrix transposed to the current basis matrix, h is the vector
   1.577 +*  of right-hand sides, x is the solution vector. */
   1.578 +
   1.579 +static void error_btran(struct csa *csa, double h[], double x[],
   1.580 +      double r[])
   1.581 +{     int m = csa->m;
   1.582 +#ifdef GLP_DEBUG
   1.583 +      int n = csa->n;
   1.584 +#endif
   1.585 +      int *A_ptr = csa->A_ptr;
   1.586 +      int *A_ind = csa->A_ind;
   1.587 +      double *A_val = csa->A_val;
   1.588 +      int *head = csa->head;
   1.589 +      int i, k, beg, end, ptr;
   1.590 +      double temp;
   1.591 +      /* compute the residual vector r = b - B'* x */
   1.592 +      for (i = 1; i <= m; i++)
   1.593 +      {  /* r[i] := b[i] - (i-th column of B)'* x */
   1.594 +         k = head[i]; /* B[i] is k-th column of (I|-A) */
   1.595 +#ifdef GLP_DEBUG
   1.596 +         xassert(1 <= k && k <= m+n);
   1.597 +#endif
   1.598 +         temp = h[i];
   1.599 +         if (k <= m)
   1.600 +         {  /* B[i] is k-th column of submatrix I */
   1.601 +            temp -= x[k];
   1.602 +         }
   1.603 +         else
   1.604 +         {  /* B[i] is (k-m)-th column of submatrix (-A) */
   1.605 +            beg = A_ptr[k-m];
   1.606 +            end = A_ptr[k-m+1];
   1.607 +            for (ptr = beg; ptr < end; ptr++)
   1.608 +               temp += A_val[ptr] * x[A_ind[ptr]];
   1.609 +         }
   1.610 +         r[i] = temp;
   1.611 +      }
   1.612 +      return;
   1.613 +}
   1.614 +
   1.615 +/***********************************************************************
   1.616 +*  refine_btran - refine solution of B'* x = h
   1.617 +*
   1.618 +*  This routine performs one iteration to refine the solution of the
   1.619 +*  system B'* x = h, where B' is a matrix transposed to the current
   1.620 +*  basis matrix, h is the vector of right-hand sides, x is the solution
   1.621 +*  vector. */
   1.622 +
   1.623 +static void refine_btran(struct csa *csa, double h[], double x[])
   1.624 +{     int m = csa->m;
   1.625 +      double *r = csa->work1;
   1.626 +      double *d = csa->work1;
   1.627 +      int i;
   1.628 +      /* compute the residual vector r = h - B'* x */
   1.629 +      error_btran(csa, h, x, r);
   1.630 +      /* compute the correction vector d = inv(B') * r */
   1.631 +      xassert(csa->valid);
   1.632 +      bfd_btran(csa->bfd, d);
   1.633 +      /* refine the solution vector (new x) = (old x) + d */
   1.634 +      for (i = 1; i <= m; i++) x[i] += d[i];
   1.635 +      return;
   1.636 +}
   1.637 +
   1.638 +/***********************************************************************
   1.639 +*  alloc_N - allocate matrix N
   1.640 +*
   1.641 +*  This routine determines maximal row lengths of matrix N, sets its
   1.642 +*  row pointers, and then allocates arrays N_ind and N_val.
   1.643 +*
   1.644 +*  Note that some fixed structural variables may temporarily become
   1.645 +*  double-bounded, so corresponding columns of matrix A should not be
   1.646 +*  ignored on calculating maximal row lengths of matrix N. */
   1.647 +
   1.648 +static void alloc_N(struct csa *csa)
   1.649 +{     int m = csa->m;
   1.650 +      int n = csa->n;
   1.651 +      int *A_ptr = csa->A_ptr;
   1.652 +      int *A_ind = csa->A_ind;
   1.653 +      int *N_ptr = csa->N_ptr;
   1.654 +      int *N_len = csa->N_len;
   1.655 +      int i, j, beg, end, ptr;
   1.656 +      /* determine number of non-zeros in each row of the augmented
   1.657 +         constraint matrix (I|-A) */
   1.658 +      for (i = 1; i <= m; i++)
   1.659 +         N_len[i] = 1;
   1.660 +      for (j = 1; j <= n; j++)
   1.661 +      {  beg = A_ptr[j];
   1.662 +         end = A_ptr[j+1];
   1.663 +         for (ptr = beg; ptr < end; ptr++)
   1.664 +            N_len[A_ind[ptr]]++;
   1.665 +      }
   1.666 +      /* determine maximal row lengths of matrix N and set its row
   1.667 +         pointers */
   1.668 +      N_ptr[1] = 1;
   1.669 +      for (i = 1; i <= m; i++)
   1.670 +      {  /* row of matrix N cannot have more than n non-zeros */
   1.671 +         if (N_len[i] > n) N_len[i] = n;
   1.672 +         N_ptr[i+1] = N_ptr[i] + N_len[i];
   1.673 +      }
   1.674 +      /* now maximal number of non-zeros in matrix N is known */
   1.675 +      csa->N_ind = xcalloc(N_ptr[m+1], sizeof(int));
   1.676 +      csa->N_val = xcalloc(N_ptr[m+1], sizeof(double));
   1.677 +      return;
   1.678 +}
   1.679 +
   1.680 +/***********************************************************************
   1.681 +*  add_N_col - add column of matrix (I|-A) to matrix N
   1.682 +*
   1.683 +*  This routine adds j-th column to matrix N which is k-th column of
   1.684 +*  the augmented constraint matrix (I|-A). (It is assumed that old j-th
   1.685 +*  column was previously removed from matrix N.) */
   1.686 +
   1.687 +static void add_N_col(struct csa *csa, int j, int k)
   1.688 +{     int m = csa->m;
   1.689 +#ifdef GLP_DEBUG
   1.690 +      int n = csa->n;
   1.691 +#endif
   1.692 +      int *N_ptr = csa->N_ptr;
   1.693 +      int *N_len = csa->N_len;
   1.694 +      int *N_ind = csa->N_ind;
   1.695 +      double *N_val = csa->N_val;
   1.696 +      int pos;
   1.697 +#ifdef GLP_DEBUG
   1.698 +      xassert(1 <= j && j <= n);
   1.699 +      xassert(1 <= k && k <= m+n);
   1.700 +#endif
   1.701 +      if (k <= m)
   1.702 +      {  /* N[j] is k-th column of submatrix I */
   1.703 +         pos = N_ptr[k] + (N_len[k]++);
   1.704 +#ifdef GLP_DEBUG
   1.705 +         xassert(pos < N_ptr[k+1]);
   1.706 +#endif
   1.707 +         N_ind[pos] = j;
   1.708 +         N_val[pos] = 1.0;
   1.709 +      }
   1.710 +      else
   1.711 +      {  /* N[j] is (k-m)-th column of submatrix (-A) */
   1.712 +         int *A_ptr = csa->A_ptr;
   1.713 +         int *A_ind = csa->A_ind;
   1.714 +         double *A_val = csa->A_val;
   1.715 +         int i, beg, end, ptr;
   1.716 +         beg = A_ptr[k-m];
   1.717 +         end = A_ptr[k-m+1];
   1.718 +         for (ptr = beg; ptr < end; ptr++)
   1.719 +         {  i = A_ind[ptr]; /* row number */
   1.720 +            pos = N_ptr[i] + (N_len[i]++);
   1.721 +#ifdef GLP_DEBUG
   1.722 +            xassert(pos < N_ptr[i+1]);
   1.723 +#endif
   1.724 +            N_ind[pos] = j;
   1.725 +            N_val[pos] = - A_val[ptr];
   1.726 +         }
   1.727 +      }
   1.728 +      return;
   1.729 +}
   1.730 +
   1.731 +/***********************************************************************
   1.732 +*  del_N_col - remove column of matrix (I|-A) from matrix N
   1.733 +*
   1.734 +*  This routine removes j-th column from matrix N which is k-th column
   1.735 +*  of the augmented constraint matrix (I|-A). */
   1.736 +
   1.737 +static void del_N_col(struct csa *csa, int j, int k)
   1.738 +{     int m = csa->m;
   1.739 +#ifdef GLP_DEBUG
   1.740 +      int n = csa->n;
   1.741 +#endif
   1.742 +      int *N_ptr = csa->N_ptr;
   1.743 +      int *N_len = csa->N_len;
   1.744 +      int *N_ind = csa->N_ind;
   1.745 +      double *N_val = csa->N_val;
   1.746 +      int pos, head, tail;
   1.747 +#ifdef GLP_DEBUG
   1.748 +      xassert(1 <= j && j <= n);
   1.749 +      xassert(1 <= k && k <= m+n);
   1.750 +#endif
   1.751 +      if (k <= m)
   1.752 +      {  /* N[j] is k-th column of submatrix I */
   1.753 +         /* find element in k-th row of N */
   1.754 +         head = N_ptr[k];
   1.755 +         for (pos = head; N_ind[pos] != j; pos++) /* nop */;
   1.756 +         /* and remove it from the row list */
   1.757 +         tail = head + (--N_len[k]);
   1.758 +#ifdef GLP_DEBUG
   1.759 +         xassert(pos <= tail);
   1.760 +#endif
   1.761 +         N_ind[pos] = N_ind[tail];
   1.762 +         N_val[pos] = N_val[tail];
   1.763 +      }
   1.764 +      else
   1.765 +      {  /* N[j] is (k-m)-th column of submatrix (-A) */
   1.766 +         int *A_ptr = csa->A_ptr;
   1.767 +         int *A_ind = csa->A_ind;
   1.768 +         int i, beg, end, ptr;
   1.769 +         beg = A_ptr[k-m];
   1.770 +         end = A_ptr[k-m+1];
   1.771 +         for (ptr = beg; ptr < end; ptr++)
   1.772 +         {  i = A_ind[ptr]; /* row number */
   1.773 +            /* find element in i-th row of N */
   1.774 +            head = N_ptr[i];
   1.775 +            for (pos = head; N_ind[pos] != j; pos++) /* nop */;
   1.776 +            /* and remove it from the row list */
   1.777 +            tail = head + (--N_len[i]);
   1.778 +#ifdef GLP_DEBUG
   1.779 +            xassert(pos <= tail);
   1.780 +#endif
   1.781 +            N_ind[pos] = N_ind[tail];
   1.782 +            N_val[pos] = N_val[tail];
   1.783 +         }
   1.784 +      }
   1.785 +      return;
   1.786 +}
   1.787 +
   1.788 +/***********************************************************************
   1.789 +*  build_N - build matrix N for current basis
   1.790 +*
   1.791 +*  This routine builds matrix N for the current basis from columns
   1.792 +*  of the augmented constraint matrix (I|-A) corresponding to non-basic
   1.793 +*  non-fixed variables. */
   1.794 +
   1.795 +static void build_N(struct csa *csa)
   1.796 +{     int m = csa->m;
   1.797 +      int n = csa->n;
   1.798 +      int *head = csa->head;
   1.799 +      char *stat = csa->stat;
   1.800 +      int *N_len = csa->N_len;
   1.801 +      int j, k;
   1.802 +      /* N := empty matrix */
   1.803 +      memset(&N_len[1], 0, m * sizeof(int));
   1.804 +      /* go through non-basic columns of matrix (I|-A) */
   1.805 +      for (j = 1; j <= n; j++)
   1.806 +      {  if (stat[j] != GLP_NS)
   1.807 +         {  /* xN[j] is non-fixed; add j-th column to matrix N which is
   1.808 +               k-th column of matrix (I|-A) */
   1.809 +            k = head[m+j]; /* x[k] = xN[j] */
   1.810 +#ifdef GLP_DEBUG
   1.811 +            xassert(1 <= k && k <= m+n);
   1.812 +#endif
   1.813 +            add_N_col(csa, j, k);
   1.814 +         }
   1.815 +      }
   1.816 +      return;
   1.817 +}
   1.818 +
   1.819 +/***********************************************************************
   1.820 +*  get_xN - determine current value of non-basic variable xN[j]
   1.821 +*
   1.822 +*  This routine returns the current value of non-basic variable xN[j],
   1.823 +*  which is a value of its active bound. */
   1.824 +
   1.825 +static double get_xN(struct csa *csa, int j)
   1.826 +{     int m = csa->m;
   1.827 +#ifdef GLP_DEBUG
   1.828 +      int n = csa->n;
   1.829 +#endif
   1.830 +      double *lb = csa->lb;
   1.831 +      double *ub = csa->ub;
   1.832 +      int *head = csa->head;
   1.833 +      char *stat = csa->stat;
   1.834 +      int k;
   1.835 +      double xN;
   1.836 +#ifdef GLP_DEBUG
   1.837 +      xassert(1 <= j && j <= n);
   1.838 +#endif
   1.839 +      k = head[m+j]; /* x[k] = xN[j] */
   1.840 +#ifdef GLP_DEBUG
   1.841 +      xassert(1 <= k && k <= m+n);
   1.842 +#endif
   1.843 +      switch (stat[j])
   1.844 +      {  case GLP_NL:
   1.845 +            /* x[k] is on its lower bound */
   1.846 +            xN = lb[k]; break;
   1.847 +         case GLP_NU:
   1.848 +            /* x[k] is on its upper bound */
   1.849 +            xN = ub[k]; break;
   1.850 +         case GLP_NF:
   1.851 +            /* x[k] is free non-basic variable */
   1.852 +            xN = 0.0; break;
   1.853 +         case GLP_NS:
   1.854 +            /* x[k] is fixed non-basic variable */
   1.855 +            xN = lb[k]; break;
   1.856 +         default:
   1.857 +            xassert(stat != stat);
   1.858 +      }
   1.859 +      return xN;
   1.860 +}
   1.861 +
   1.862 +/***********************************************************************
   1.863 +*  eval_beta - compute primal values of basic variables
   1.864 +*
   1.865 +*  This routine computes current primal values of all basic variables:
   1.866 +*
   1.867 +*     beta = - inv(B) * N * xN,
   1.868 +*
   1.869 +*  where B is the current basis matrix, N is a matrix built of columns
   1.870 +*  of matrix (I|-A) corresponding to non-basic variables, and xN is the
   1.871 +*  vector of current values of non-basic variables. */
   1.872 +
   1.873 +static void eval_beta(struct csa *csa, double beta[])
   1.874 +{     int m = csa->m;
   1.875 +      int n = csa->n;
   1.876 +      int *A_ptr = csa->A_ptr;
   1.877 +      int *A_ind = csa->A_ind;
   1.878 +      double *A_val = csa->A_val;
   1.879 +      int *head = csa->head;
   1.880 +      double *h = csa->work2;
   1.881 +      int i, j, k, beg, end, ptr;
   1.882 +      double xN;
   1.883 +      /* compute the right-hand side vector:
   1.884 +         h := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n],
   1.885 +         where N[1], ..., N[n] are columns of matrix N */
   1.886 +      for (i = 1; i <= m; i++)
   1.887 +         h[i] = 0.0;
   1.888 +      for (j = 1; j <= n; j++)
   1.889 +      {  k = head[m+j]; /* x[k] = xN[j] */
   1.890 +#ifdef GLP_DEBUG
   1.891 +         xassert(1 <= k && k <= m+n);
   1.892 +#endif
   1.893 +         /* determine current value of xN[j] */
   1.894 +         xN = get_xN(csa, j);
   1.895 +         if (xN == 0.0) continue;
   1.896 +         if (k <= m)
   1.897 +         {  /* N[j] is k-th column of submatrix I */
   1.898 +            h[k] -= xN;
   1.899 +         }
   1.900 +         else
   1.901 +         {  /* N[j] is (k-m)-th column of submatrix (-A) */
   1.902 +            beg = A_ptr[k-m];
   1.903 +            end = A_ptr[k-m+1];
   1.904 +            for (ptr = beg; ptr < end; ptr++)
   1.905 +               h[A_ind[ptr]] += xN * A_val[ptr];
   1.906 +         }
   1.907 +      }
   1.908 +      /* solve system B * beta = h */
   1.909 +      memcpy(&beta[1], &h[1], m * sizeof(double));
   1.910 +      xassert(csa->valid);
   1.911 +      bfd_ftran(csa->bfd, beta);
   1.912 +      /* and refine the solution */
   1.913 +      refine_ftran(csa, h, beta);
   1.914 +      return;
   1.915 +}
   1.916 +
   1.917 +/***********************************************************************
   1.918 +*  eval_pi - compute vector of simplex multipliers
   1.919 +*
   1.920 +*  This routine computes the vector of current simplex multipliers:
   1.921 +*
   1.922 +*     pi = inv(B') * cB,
   1.923 +*
   1.924 +*  where B' is a matrix transposed to the current basis matrix, cB is
   1.925 +*  a subvector of objective coefficients at basic variables. */
   1.926 +
   1.927 +static void eval_pi(struct csa *csa, double pi[])
   1.928 +{     int m = csa->m;
   1.929 +      double *c = csa->coef;
   1.930 +      int *head = csa->head;
   1.931 +      double *cB = csa->work2;
   1.932 +      int i;
   1.933 +      /* construct the right-hand side vector cB */
   1.934 +      for (i = 1; i <= m; i++)
   1.935 +         cB[i] = c[head[i]];
   1.936 +      /* solve system B'* pi = cB */
   1.937 +      memcpy(&pi[1], &cB[1], m * sizeof(double));
   1.938 +      xassert(csa->valid);
   1.939 +      bfd_btran(csa->bfd, pi);
   1.940 +      /* and refine the solution */
   1.941 +      refine_btran(csa, cB, pi);
   1.942 +      return;
   1.943 +}
   1.944 +
   1.945 +/***********************************************************************
   1.946 +*  eval_cost - compute reduced cost of non-basic variable xN[j]
   1.947 +*
   1.948 +*  This routine computes the current reduced cost of non-basic variable
   1.949 +*  xN[j]:
   1.950 +*
   1.951 +*     d[j] = cN[j] - N'[j] * pi,
   1.952 +*
   1.953 +*  where cN[j] is the objective coefficient at variable xN[j], N[j] is
   1.954 +*  a column of the augmented constraint matrix (I|-A) corresponding to
   1.955 +*  xN[j], pi is the vector of simplex multipliers. */
   1.956 +
   1.957 +static double eval_cost(struct csa *csa, double pi[], int j)
   1.958 +{     int m = csa->m;
   1.959 +#ifdef GLP_DEBUG
   1.960 +      int n = csa->n;
   1.961 +#endif
   1.962 +      double *coef = csa->coef;
   1.963 +      int *head = csa->head;
   1.964 +      int k;
   1.965 +      double dj;
   1.966 +#ifdef GLP_DEBUG
   1.967 +      xassert(1 <= j && j <= n);
   1.968 +#endif
   1.969 +      k = head[m+j]; /* x[k] = xN[j] */
   1.970 +#ifdef GLP_DEBUG
   1.971 +      xassert(1 <= k && k <= m+n);
   1.972 +#endif
   1.973 +      dj = coef[k];
   1.974 +      if (k <= m)
   1.975 +      {  /* N[j] is k-th column of submatrix I */
   1.976 +         dj -= pi[k];
   1.977 +      }
   1.978 +      else
   1.979 +      {  /* N[j] is (k-m)-th column of submatrix (-A) */
   1.980 +         int *A_ptr = csa->A_ptr;
   1.981 +         int *A_ind = csa->A_ind;
   1.982 +         double *A_val = csa->A_val;
   1.983 +         int beg, end, ptr;
   1.984 +         beg = A_ptr[k-m];
   1.985 +         end = A_ptr[k-m+1];
   1.986 +         for (ptr = beg; ptr < end; ptr++)
   1.987 +            dj += A_val[ptr] * pi[A_ind[ptr]];
   1.988 +      }
   1.989 +      return dj;
   1.990 +}
   1.991 +
   1.992 +/***********************************************************************
   1.993 +*  eval_bbar - compute and store primal values of basic variables
   1.994 +*
   1.995 +*  This routine computes primal values of all basic variables and then
   1.996 +*  stores them in the solution array. */
   1.997 +
   1.998 +static void eval_bbar(struct csa *csa)
   1.999 +{     eval_beta(csa, csa->bbar);
  1.1000 +      return;
  1.1001 +}
  1.1002 +
  1.1003 +/***********************************************************************
  1.1004 +*  eval_cbar - compute and store reduced costs of non-basic variables
  1.1005 +*
  1.1006 +*  This routine computes reduced costs of all non-basic variables and
  1.1007 +*  then stores them in the solution array. */
  1.1008 +
  1.1009 +static void eval_cbar(struct csa *csa)
  1.1010 +{
  1.1011 +#ifdef GLP_DEBUG
  1.1012 +      int m = csa->m;
  1.1013 +#endif
  1.1014 +      int n = csa->n;
  1.1015 +#ifdef GLP_DEBUG
  1.1016 +      int *head = csa->head;
  1.1017 +#endif
  1.1018 +      double *cbar = csa->cbar;
  1.1019 +      double *pi = csa->work3;
  1.1020 +      int j;
  1.1021 +#ifdef GLP_DEBUG
  1.1022 +      int k;
  1.1023 +#endif
  1.1024 +      /* compute simplex multipliers */
  1.1025 +      eval_pi(csa, pi);
  1.1026 +      /* compute and store reduced costs */
  1.1027 +      for (j = 1; j <= n; j++)
  1.1028 +      {
  1.1029 +#ifdef GLP_DEBUG
  1.1030 +         k = head[m+j]; /* x[k] = xN[j] */
  1.1031 +         xassert(1 <= k && k <= m+n);
  1.1032 +#endif
  1.1033 +         cbar[j] = eval_cost(csa, pi, j);
  1.1034 +      }
  1.1035 +      return;
  1.1036 +}
  1.1037 +
  1.1038 +/***********************************************************************
  1.1039 +*  reset_refsp - reset the reference space
  1.1040 +*
  1.1041 +*  This routine resets (redefines) the reference space used in the
  1.1042 +*  projected steepest edge pricing algorithm. */
  1.1043 +
  1.1044 +static void reset_refsp(struct csa *csa)
  1.1045 +{     int m = csa->m;
  1.1046 +      int n = csa->n;
  1.1047 +      int *head = csa->head;
  1.1048 +      char *refsp = csa->refsp;
  1.1049 +      double *gamma = csa->gamma;
  1.1050 +      int j, k;
  1.1051 +      xassert(csa->refct == 0);
  1.1052 +      csa->refct = 1000;
  1.1053 +      memset(&refsp[1], 0, (m+n) * sizeof(char));
  1.1054 +      for (j = 1; j <= n; j++)
  1.1055 +      {  k = head[m+j]; /* x[k] = xN[j] */
  1.1056 +         refsp[k] = 1;
  1.1057 +         gamma[j] = 1.0;
  1.1058 +      }
  1.1059 +      return;
  1.1060 +}
  1.1061 +
  1.1062 +/***********************************************************************
  1.1063 +*  eval_gamma - compute steepest edge coefficient
  1.1064 +*
  1.1065 +*  This routine computes the steepest edge coefficient for non-basic
  1.1066 +*  variable xN[j] using its direct definition:
  1.1067 +*
  1.1068 +*     gamma[j] = delta[j] +  sum   alfa[i,j]^2,
  1.1069 +*                           i in R
  1.1070 +*
  1.1071 +*  where delta[j] = 1, if xN[j] is in the current reference space,
  1.1072 +*  and 0 otherwise; R is a set of basic variables xB[i], which are in
  1.1073 +*  the current reference space; alfa[i,j] are elements of the current
  1.1074 +*  simplex table.
  1.1075 +*
  1.1076 +*  NOTE: The routine is intended only for debugginig purposes. */
  1.1077 +
  1.1078 +static double eval_gamma(struct csa *csa, int j)
  1.1079 +{     int m = csa->m;
  1.1080 +#ifdef GLP_DEBUG
  1.1081 +      int n = csa->n;
  1.1082 +#endif
  1.1083 +      int *head = csa->head;
  1.1084 +      char *refsp = csa->refsp;
  1.1085 +      double *alfa = csa->work3;
  1.1086 +      double *h = csa->work3;
  1.1087 +      int i, k;
  1.1088 +      double gamma;
  1.1089 +#ifdef GLP_DEBUG
  1.1090 +      xassert(1 <= j && j <= n);
  1.1091 +#endif
  1.1092 +      k = head[m+j]; /* x[k] = xN[j] */
  1.1093 +#ifdef GLP_DEBUG
  1.1094 +      xassert(1 <= k && k <= m+n);
  1.1095 +#endif
  1.1096 +      /* construct the right-hand side vector h = - N[j] */
  1.1097 +      for (i = 1; i <= m; i++)
  1.1098 +         h[i] = 0.0;
  1.1099 +      if (k <= m)
  1.1100 +      {  /* N[j] is k-th column of submatrix I */
  1.1101 +         h[k] = -1.0;
  1.1102 +      }
  1.1103 +      else
  1.1104 +      {  /* N[j] is (k-m)-th column of submatrix (-A) */
  1.1105 +         int *A_ptr = csa->A_ptr;
  1.1106 +         int *A_ind = csa->A_ind;
  1.1107 +         double *A_val = csa->A_val;
  1.1108 +         int beg, end, ptr;
  1.1109 +         beg = A_ptr[k-m];
  1.1110 +         end = A_ptr[k-m+1];
  1.1111 +         for (ptr = beg; ptr < end; ptr++)
  1.1112 +            h[A_ind[ptr]] = A_val[ptr];
  1.1113 +      }
  1.1114 +      /* solve system B * alfa = h */
  1.1115 +      xassert(csa->valid);
  1.1116 +      bfd_ftran(csa->bfd, alfa);
  1.1117 +      /* compute gamma */
  1.1118 +      gamma = (refsp[k] ? 1.0 : 0.0);
  1.1119 +      for (i = 1; i <= m; i++)
  1.1120 +      {  k = head[i];
  1.1121 +#ifdef GLP_DEBUG
  1.1122 +         xassert(1 <= k && k <= m+n);
  1.1123 +#endif
  1.1124 +         if (refsp[k]) gamma += alfa[i] * alfa[i];
  1.1125 +      }
  1.1126 +      return gamma;
  1.1127 +}
  1.1128 +
  1.1129 +/***********************************************************************
  1.1130 +*  chuzc - choose non-basic variable (column of the simplex table)
  1.1131 +*
  1.1132 +*  This routine chooses non-basic variable xN[q], which has largest
  1.1133 +*  weighted reduced cost:
  1.1134 +*
  1.1135 +*     |d[q]| / sqrt(gamma[q]) = max  |d[j]| / sqrt(gamma[j]),
  1.1136 +*                              j in J
  1.1137 +*
  1.1138 +*  where J is a subset of eligible non-basic variables xN[j], d[j] is
  1.1139 +*  reduced cost of xN[j], gamma[j] is the steepest edge coefficient.
  1.1140 +*
  1.1141 +*  The working objective function is always minimized, so the sign of
  1.1142 +*  d[q] determines direction, in which xN[q] has to change:
  1.1143 +*
  1.1144 +*     if d[q] < 0, xN[q] has to increase;
  1.1145 +*
  1.1146 +*     if d[q] > 0, xN[q] has to decrease.
  1.1147 +*
  1.1148 +*  If |d[j]| <= tol_dj, where tol_dj is a specified tolerance, xN[j]
  1.1149 +*  is not included in J and therefore ignored. (It is assumed that the
  1.1150 +*  working objective row is appropriately scaled, i.e. max|c[k]| = 1.)
  1.1151 +*
  1.1152 +*  If J is empty and no variable has been chosen, q is set to 0. */
  1.1153 +
  1.1154 +static void chuzc(struct csa *csa, double tol_dj)
  1.1155 +{     int n = csa->n;
  1.1156 +      char *stat = csa->stat;
  1.1157 +      double *cbar = csa->cbar;
  1.1158 +      double *gamma = csa->gamma;
  1.1159 +      int j, q;
  1.1160 +      double dj, best, temp;
  1.1161 +      /* nothing is chosen so far */
  1.1162 +      q = 0, best = 0.0;
  1.1163 +      /* look through the list of non-basic variables */
  1.1164 +      for (j = 1; j <= n; j++)
  1.1165 +      {  dj = cbar[j];
  1.1166 +         switch (stat[j])
  1.1167 +         {  case GLP_NL:
  1.1168 +               /* xN[j] can increase */
  1.1169 +               if (dj >= - tol_dj) continue;
  1.1170 +               break;
  1.1171 +            case GLP_NU:
  1.1172 +               /* xN[j] can decrease */
  1.1173 +               if (dj <= + tol_dj) continue;
  1.1174 +               break;
  1.1175 +            case GLP_NF:
  1.1176 +               /* xN[j] can change in any direction */
  1.1177 +               if (- tol_dj <= dj && dj <= + tol_dj) continue;
  1.1178 +               break;
  1.1179 +            case GLP_NS:
  1.1180 +               /* xN[j] cannot change at all */
  1.1181 +               continue;
  1.1182 +            default:
  1.1183 +               xassert(stat != stat);
  1.1184 +         }
  1.1185 +         /* xN[j] is eligible non-basic variable; choose one which has
  1.1186 +            largest weighted reduced cost */
  1.1187 +#ifdef GLP_DEBUG
  1.1188 +         xassert(gamma[j] > 0.0);
  1.1189 +#endif
  1.1190 +         temp = (dj * dj) / gamma[j];
  1.1191 +         if (best < temp)
  1.1192 +            q = j, best = temp;
  1.1193 +      }
  1.1194 +      /* store the index of non-basic variable xN[q] chosen */
  1.1195 +      csa->q = q;
  1.1196 +      return;
  1.1197 +}
  1.1198 +
  1.1199 +/***********************************************************************
  1.1200 +*  eval_tcol - compute pivot column of the simplex table
  1.1201 +*
  1.1202 +*  This routine computes the pivot column of the simplex table, which
  1.1203 +*  corresponds to non-basic variable xN[q] chosen.
  1.1204 +*
  1.1205 +*  The pivot column is the following vector:
  1.1206 +*
  1.1207 +*     tcol = T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q],
  1.1208 +*
  1.1209 +*  where B is the current basis matrix, N[q] is a column of the matrix
  1.1210 +*  (I|-A) corresponding to variable xN[q]. */
  1.1211 +
  1.1212 +static void eval_tcol(struct csa *csa)
  1.1213 +{     int m = csa->m;
  1.1214 +#ifdef GLP_DEBUG
  1.1215 +      int n = csa->n;
  1.1216 +#endif
  1.1217 +      int *head = csa->head;
  1.1218 +      int q = csa->q;
  1.1219 +      int *tcol_ind = csa->tcol_ind;
  1.1220 +      double *tcol_vec = csa->tcol_vec;
  1.1221 +      double *h = csa->tcol_vec;
  1.1222 +      int i, k, nnz;
  1.1223 +#ifdef GLP_DEBUG
  1.1224 +      xassert(1 <= q && q <= n);
  1.1225 +#endif
  1.1226 +      k = head[m+q]; /* x[k] = xN[q] */
  1.1227 +#ifdef GLP_DEBUG
  1.1228 +      xassert(1 <= k && k <= m+n);
  1.1229 +#endif
  1.1230 +      /* construct the right-hand side vector h = - N[q] */
  1.1231 +      for (i = 1; i <= m; i++)
  1.1232 +         h[i] = 0.0;
  1.1233 +      if (k <= m)
  1.1234 +      {  /* N[q] is k-th column of submatrix I */
  1.1235 +         h[k] = -1.0;
  1.1236 +      }
  1.1237 +      else
  1.1238 +      {  /* N[q] is (k-m)-th column of submatrix (-A) */
  1.1239 +         int *A_ptr = csa->A_ptr;
  1.1240 +         int *A_ind = csa->A_ind;
  1.1241 +         double *A_val = csa->A_val;
  1.1242 +         int beg, end, ptr;
  1.1243 +         beg = A_ptr[k-m];
  1.1244 +         end = A_ptr[k-m+1];
  1.1245 +         for (ptr = beg; ptr < end; ptr++)
  1.1246 +            h[A_ind[ptr]] = A_val[ptr];
  1.1247 +      }
  1.1248 +      /* solve system B * tcol = h */
  1.1249 +      xassert(csa->valid);
  1.1250 +      bfd_ftran(csa->bfd, tcol_vec);
  1.1251 +      /* construct sparse pattern of the pivot column */
  1.1252 +      nnz = 0;
  1.1253 +      for (i = 1; i <= m; i++)
  1.1254 +      {  if (tcol_vec[i] != 0.0)
  1.1255 +            tcol_ind[++nnz] = i;
  1.1256 +      }
  1.1257 +      csa->tcol_nnz = nnz;
  1.1258 +      return;
  1.1259 +}
  1.1260 +
  1.1261 +/***********************************************************************
  1.1262 +*  refine_tcol - refine pivot column of the simplex table
  1.1263 +*
  1.1264 +*  This routine refines the pivot column of the simplex table assuming
  1.1265 +*  that it was previously computed by the routine eval_tcol. */
  1.1266 +
  1.1267 +static void refine_tcol(struct csa *csa)
  1.1268 +{     int m = csa->m;
  1.1269 +#ifdef GLP_DEBUG
  1.1270 +      int n = csa->n;
  1.1271 +#endif
  1.1272 +      int *head = csa->head;
  1.1273 +      int q = csa->q;
  1.1274 +      int *tcol_ind = csa->tcol_ind;
  1.1275 +      double *tcol_vec = csa->tcol_vec;
  1.1276 +      double *h = csa->work3;
  1.1277 +      int i, k, nnz;
  1.1278 +#ifdef GLP_DEBUG
  1.1279 +      xassert(1 <= q && q <= n);
  1.1280 +#endif
  1.1281 +      k = head[m+q]; /* x[k] = xN[q] */
  1.1282 +#ifdef GLP_DEBUG
  1.1283 +      xassert(1 <= k && k <= m+n);
  1.1284 +#endif
  1.1285 +      /* construct the right-hand side vector h = - N[q] */
  1.1286 +      for (i = 1; i <= m; i++)
  1.1287 +         h[i] = 0.0;
  1.1288 +      if (k <= m)
  1.1289 +      {  /* N[q] is k-th column of submatrix I */
  1.1290 +         h[k] = -1.0;
  1.1291 +      }
  1.1292 +      else
  1.1293 +      {  /* N[q] is (k-m)-th column of submatrix (-A) */
  1.1294 +         int *A_ptr = csa->A_ptr;
  1.1295 +         int *A_ind = csa->A_ind;
  1.1296 +         double *A_val = csa->A_val;
  1.1297 +         int beg, end, ptr;
  1.1298 +         beg = A_ptr[k-m];
  1.1299 +         end = A_ptr[k-m+1];
  1.1300 +         for (ptr = beg; ptr < end; ptr++)
  1.1301 +            h[A_ind[ptr]] = A_val[ptr];
  1.1302 +      }
  1.1303 +      /* refine solution of B * tcol = h */
  1.1304 +      refine_ftran(csa, h, tcol_vec);
  1.1305 +      /* construct sparse pattern of the pivot column */
  1.1306 +      nnz = 0;
  1.1307 +      for (i = 1; i <= m; i++)
  1.1308 +      {  if (tcol_vec[i] != 0.0)
  1.1309 +            tcol_ind[++nnz] = i;
  1.1310 +      }
  1.1311 +      csa->tcol_nnz = nnz;
  1.1312 +      return;
  1.1313 +}
  1.1314 +
  1.1315 +/***********************************************************************
  1.1316 +*  sort_tcol - sort pivot column of the simplex table
  1.1317 +*
  1.1318 +*  This routine reorders the list of non-zero elements of the pivot
  1.1319 +*  column to put significant elements, whose magnitude is not less than
  1.1320 +*  a specified tolerance, in front of the list, and stores the number
  1.1321 +*  of significant elements in tcol_num. */
  1.1322 +
  1.1323 +static void sort_tcol(struct csa *csa, double tol_piv)
  1.1324 +{
  1.1325 +#ifdef GLP_DEBUG
  1.1326 +      int m = csa->m;
  1.1327 +#endif
  1.1328 +      int nnz = csa->tcol_nnz;
  1.1329 +      int *tcol_ind = csa->tcol_ind;
  1.1330 +      double *tcol_vec = csa->tcol_vec;
  1.1331 +      int i, num, pos;
  1.1332 +      double big, eps, temp;
  1.1333 +      /* compute infinity (maximum) norm of the column */
  1.1334 +      big = 0.0;
  1.1335 +      for (pos = 1; pos <= nnz; pos++)
  1.1336 +      {
  1.1337 +#ifdef GLP_DEBUG
  1.1338 +         i = tcol_ind[pos];
  1.1339 +         xassert(1 <= i && i <= m);
  1.1340 +#endif
  1.1341 +         temp = fabs(tcol_vec[tcol_ind[pos]]);
  1.1342 +         if (big < temp) big = temp;
  1.1343 +      }
  1.1344 +      csa->tcol_max = big;
  1.1345 +      /* determine absolute pivot tolerance */
  1.1346 +      eps = tol_piv * (1.0 + 0.01 * big);
  1.1347 +      /* move significant column components to front of the list */
  1.1348 +      for (num = 0; num < nnz; )
  1.1349 +      {  i = tcol_ind[nnz];
  1.1350 +         if (fabs(tcol_vec[i]) < eps)
  1.1351 +            nnz--;
  1.1352 +         else
  1.1353 +         {  num++;
  1.1354 +            tcol_ind[nnz] = tcol_ind[num];
  1.1355 +            tcol_ind[num] = i;
  1.1356 +         }
  1.1357 +      }
  1.1358 +      csa->tcol_num = num;
  1.1359 +      return;
  1.1360 +}
  1.1361 +
  1.1362 +/***********************************************************************
  1.1363 +*  chuzr - choose basic variable (row of the simplex table)
  1.1364 +*
  1.1365 +*  This routine chooses basic variable xB[p], which reaches its bound
  1.1366 +*  first on changing non-basic variable xN[q] in valid direction.
  1.1367 +*
  1.1368 +*  The parameter rtol is a relative tolerance used to relax bounds of
  1.1369 +*  basic variables. If rtol = 0, the routine implements the standard
  1.1370 +*  ratio test. Otherwise, if rtol > 0, the routine implements Harris'
  1.1371 +*  two-pass ratio test. In the latter case rtol should be about three
  1.1372 +*  times less than a tolerance used to check primal feasibility. */
  1.1373 +
  1.1374 +static void chuzr(struct csa *csa, double rtol)
  1.1375 +{     int m = csa->m;
  1.1376 +#ifdef GLP_DEBUG
  1.1377 +      int n = csa->n;
  1.1378 +#endif
  1.1379 +      char *type = csa->type;
  1.1380 +      double *lb = csa->lb;
  1.1381 +      double *ub = csa->ub;
  1.1382 +      double *coef = csa->coef;
  1.1383 +      int *head = csa->head;
  1.1384 +      int phase = csa->phase;
  1.1385 +      double *bbar = csa->bbar;
  1.1386 +      double *cbar = csa->cbar;
  1.1387 +      int q = csa->q;
  1.1388 +      int *tcol_ind = csa->tcol_ind;
  1.1389 +      double *tcol_vec = csa->tcol_vec;
  1.1390 +      int tcol_num = csa->tcol_num;
  1.1391 +      int i, i_stat, k, p, p_stat, pos;
  1.1392 +      double alfa, big, delta, s, t, teta, tmax;
  1.1393 +#ifdef GLP_DEBUG
  1.1394 +      xassert(1 <= q && q <= n);
  1.1395 +#endif
  1.1396 +      /* s := - sign(d[q]), where d[q] is reduced cost of xN[q] */
  1.1397 +#ifdef GLP_DEBUG
  1.1398 +      xassert(cbar[q] != 0.0);
  1.1399 +#endif
  1.1400 +      s = (cbar[q] > 0.0 ? -1.0 : +1.0);
  1.1401 +      /*** FIRST PASS ***/
  1.1402 +      k = head[m+q]; /* x[k] = xN[q] */
  1.1403 +#ifdef GLP_DEBUG
  1.1404 +      xassert(1 <= k && k <= m+n);
  1.1405 +#endif
  1.1406 +      if (type[k] == GLP_DB)
  1.1407 +      {  /* xN[q] has both lower and upper bounds */
  1.1408 +         p = -1, p_stat = 0, teta = ub[k] - lb[k], big = 1.0;
  1.1409 +      }
  1.1410 +      else
  1.1411 +      {  /* xN[q] has no opposite bound */
  1.1412 +         p = 0, p_stat = 0, teta = DBL_MAX, big = 0.0;
  1.1413 +      }
  1.1414 +      /* walk through significant elements of the pivot column */
  1.1415 +      for (pos = 1; pos <= tcol_num; pos++)
  1.1416 +      {  i = tcol_ind[pos];
  1.1417 +#ifdef GLP_DEBUG
  1.1418 +         xassert(1 <= i && i <= m);
  1.1419 +#endif
  1.1420 +         k = head[i]; /* x[k] = xB[i] */
  1.1421 +#ifdef GLP_DEBUG
  1.1422 +         xassert(1 <= k && k <= m+n);
  1.1423 +#endif
  1.1424 +         alfa = s * tcol_vec[i];
  1.1425 +#ifdef GLP_DEBUG
  1.1426 +         xassert(alfa != 0.0);
  1.1427 +#endif
  1.1428 +         /* xB[i] = ... + alfa * xN[q] + ..., and due to s we need to
  1.1429 +            consider the only case when xN[q] is increasing */
  1.1430 +         if (alfa > 0.0)
  1.1431 +         {  /* xB[i] is increasing */
  1.1432 +            if (phase == 1 && coef[k] < 0.0)
  1.1433 +            {  /* xB[i] violates its lower bound, which plays the role
  1.1434 +                  of an upper bound on phase I */
  1.1435 +               delta = rtol * (1.0 + kappa * fabs(lb[k]));
  1.1436 +               t = ((lb[k] + delta) - bbar[i]) / alfa;
  1.1437 +               i_stat = GLP_NL;
  1.1438 +            }
  1.1439 +            else if (phase == 1 && coef[k] > 0.0)
  1.1440 +            {  /* xB[i] violates its upper bound, which plays the role
  1.1441 +                  of an lower bound on phase I */
  1.1442 +               continue;
  1.1443 +            }
  1.1444 +            else if (type[k] == GLP_UP || type[k] == GLP_DB ||
  1.1445 +                     type[k] == GLP_FX)
  1.1446 +            {  /* xB[i] is within its bounds and has an upper bound */
  1.1447 +               delta = rtol * (1.0 + kappa * fabs(ub[k]));
  1.1448 +               t = ((ub[k] + delta) - bbar[i]) / alfa;
  1.1449 +               i_stat = GLP_NU;
  1.1450 +            }
  1.1451 +            else
  1.1452 +            {  /* xB[i] is within its bounds and has no upper bound */
  1.1453 +               continue;
  1.1454 +            }
  1.1455 +         }
  1.1456 +         else
  1.1457 +         {  /* xB[i] is decreasing */
  1.1458 +            if (phase == 1 && coef[k] > 0.0)
  1.1459 +            {  /* xB[i] violates its upper bound, which plays the role
  1.1460 +                  of an lower bound on phase I */
  1.1461 +               delta = rtol * (1.0 + kappa * fabs(ub[k]));
  1.1462 +               t = ((ub[k] - delta) - bbar[i]) / alfa;
  1.1463 +               i_stat = GLP_NU;
  1.1464 +            }
  1.1465 +            else if (phase == 1 && coef[k] < 0.0)
  1.1466 +            {  /* xB[i] violates its lower bound, which plays the role
  1.1467 +                  of an upper bound on phase I */
  1.1468 +               continue;
  1.1469 +            }
  1.1470 +            else if (type[k] == GLP_LO || type[k] == GLP_DB ||
  1.1471 +                     type[k] == GLP_FX)
  1.1472 +            {  /* xB[i] is within its bounds and has an lower bound */
  1.1473 +               delta = rtol * (1.0 + kappa * fabs(lb[k]));
  1.1474 +               t = ((lb[k] - delta) - bbar[i]) / alfa;
  1.1475 +               i_stat = GLP_NL;
  1.1476 +            }
  1.1477 +            else
  1.1478 +            {  /* xB[i] is within its bounds and has no lower bound */
  1.1479 +               continue;
  1.1480 +            }
  1.1481 +         }
  1.1482 +         /* t is a change of xN[q], on which xB[i] reaches its bound
  1.1483 +            (possibly relaxed); since the basic solution is assumed to
  1.1484 +            be primal feasible (or pseudo feasible on phase I), t has
  1.1485 +            to be non-negative by definition; however, it may happen
  1.1486 +            that xB[i] slightly (i.e. within a tolerance) violates its
  1.1487 +            bound, that leads to negative t; in the latter case, if
  1.1488 +            xB[i] is chosen, negative t means that xN[q] changes in
  1.1489 +            wrong direction; if pivot alfa[i,q] is close to zero, even
  1.1490 +            small bound violation of xB[i] may lead to a large change
  1.1491 +            of xN[q] in wrong direction; let, for example, xB[i] >= 0
  1.1492 +            and in the current basis its value be -5e-9; let also xN[q]
  1.1493 +            be on its zero bound and should increase; from the ratio
  1.1494 +            test rule it follows that the pivot alfa[i,q] < 0; however,
  1.1495 +            if alfa[i,q] is, say, -1e-9, the change of xN[q] in wrong
  1.1496 +            direction is 5e-9 / (-1e-9) = -5, and using it for updating
  1.1497 +            values of other basic variables will give absolutely wrong
  1.1498 +            results; therefore, if t is negative, we should replace it
  1.1499 +            by exact zero assuming that xB[i] is exactly on its bound,
  1.1500 +            and the violation appears due to round-off errors */
  1.1501 +         if (t < 0.0) t = 0.0;
  1.1502 +         /* apply minimal ratio test */
  1.1503 +         if (teta > t || teta == t && big < fabs(alfa))
  1.1504 +            p = i, p_stat = i_stat, teta = t, big = fabs(alfa);
  1.1505 +      }
  1.1506 +      /* the second pass is skipped in the following cases: */
  1.1507 +      /* if the standard ratio test is used */
  1.1508 +      if (rtol == 0.0) goto done;
  1.1509 +      /* if xN[q] reaches its opposite bound or if no basic variable
  1.1510 +         has been chosen on the first pass */
  1.1511 +      if (p <= 0) goto done;
  1.1512 +      /* if xB[p] is a blocking variable, i.e. if it prevents xN[q]
  1.1513 +         from any change */
  1.1514 +      if (teta == 0.0) goto done;
  1.1515 +      /*** SECOND PASS ***/
  1.1516 +      /* here tmax is a maximal change of xN[q], on which the solution
  1.1517 +         remains primal feasible (or pseudo feasible on phase I) within
  1.1518 +         a tolerance */
  1.1519 +#if 0
  1.1520 +      tmax = (1.0 + 10.0 * DBL_EPSILON) * teta;
  1.1521 +#else
  1.1522 +      tmax = teta;
  1.1523 +#endif
  1.1524 +      /* nothing is chosen so far */
  1.1525 +      p = 0, p_stat = 0, teta = DBL_MAX, big = 0.0;
  1.1526 +      /* walk through significant elements of the pivot column */
  1.1527 +      for (pos = 1; pos <= tcol_num; pos++)
  1.1528 +      {  i = tcol_ind[pos];
  1.1529 +#ifdef GLP_DEBUG
  1.1530 +         xassert(1 <= i && i <= m);
  1.1531 +#endif
  1.1532 +         k = head[i]; /* x[k] = xB[i] */
  1.1533 +#ifdef GLP_DEBUG
  1.1534 +         xassert(1 <= k && k <= m+n);
  1.1535 +#endif
  1.1536 +         alfa = s * tcol_vec[i];
  1.1537 +#ifdef GLP_DEBUG
  1.1538 +         xassert(alfa != 0.0);
  1.1539 +#endif
  1.1540 +         /* xB[i] = ... + alfa * xN[q] + ..., and due to s we need to
  1.1541 +            consider the only case when xN[q] is increasing */
  1.1542 +         if (alfa > 0.0)
  1.1543 +         {  /* xB[i] is increasing */
  1.1544 +            if (phase == 1 && coef[k] < 0.0)
  1.1545 +            {  /* xB[i] violates its lower bound, which plays the role
  1.1546 +                  of an upper bound on phase I */
  1.1547 +               t = (lb[k] - bbar[i]) / alfa;
  1.1548 +               i_stat = GLP_NL;
  1.1549 +            }
  1.1550 +            else if (phase == 1 && coef[k] > 0.0)
  1.1551 +            {  /* xB[i] violates its upper bound, which plays the role
  1.1552 +                  of an lower bound on phase I */
  1.1553 +               continue;
  1.1554 +            }
  1.1555 +            else if (type[k] == GLP_UP || type[k] == GLP_DB ||
  1.1556 +                     type[k] == GLP_FX)
  1.1557 +            {  /* xB[i] is within its bounds and has an upper bound */
  1.1558 +               t = (ub[k] - bbar[i]) / alfa;
  1.1559 +               i_stat = GLP_NU;
  1.1560 +            }
  1.1561 +            else
  1.1562 +            {  /* xB[i] is within its bounds and has no upper bound */
  1.1563 +               continue;
  1.1564 +            }
  1.1565 +         }
  1.1566 +         else
  1.1567 +         {  /* xB[i] is decreasing */
  1.1568 +            if (phase == 1 && coef[k] > 0.0)
  1.1569 +            {  /* xB[i] violates its upper bound, which plays the role
  1.1570 +                  of an lower bound on phase I */
  1.1571 +               t = (ub[k] - bbar[i]) / alfa;
  1.1572 +               i_stat = GLP_NU;
  1.1573 +            }
  1.1574 +            else if (phase == 1 && coef[k] < 0.0)
  1.1575 +            {  /* xB[i] violates its lower bound, which plays the role
  1.1576 +                  of an upper bound on phase I */
  1.1577 +               continue;
  1.1578 +            }
  1.1579 +            else if (type[k] == GLP_LO || type[k] == GLP_DB ||
  1.1580 +                     type[k] == GLP_FX)
  1.1581 +            {  /* xB[i] is within its bounds and has an lower bound */
  1.1582 +               t = (lb[k] - bbar[i]) / alfa;
  1.1583 +               i_stat = GLP_NL;
  1.1584 +            }
  1.1585 +            else
  1.1586 +            {  /* xB[i] is within its bounds and has no lower bound */
  1.1587 +               continue;
  1.1588 +            }
  1.1589 +         }
  1.1590 +         /* (see comments for the first pass) */
  1.1591 +         if (t < 0.0) t = 0.0;
  1.1592 +         /* t is a change of xN[q], on which xB[i] reaches its bound;
  1.1593 +            if t <= tmax, all basic variables can violate their bounds
  1.1594 +            only within relaxation tolerance delta; we can use this
  1.1595 +            freedom and choose basic variable having largest influence
  1.1596 +            coefficient to avoid possible numeric instability */
  1.1597 +         if (t <= tmax && big < fabs(alfa))
  1.1598 +            p = i, p_stat = i_stat, teta = t, big = fabs(alfa);
  1.1599 +      }
  1.1600 +      /* something must be chosen on the second pass */
  1.1601 +      xassert(p != 0);
  1.1602 +done: /* store the index and status of basic variable xB[p] chosen */
  1.1603 +      csa->p = p;
  1.1604 +      if (p > 0 && type[head[p]] == GLP_FX)
  1.1605 +         csa->p_stat = GLP_NS;
  1.1606 +      else
  1.1607 +         csa->p_stat = p_stat;
  1.1608 +      /* store corresponding change of non-basic variable xN[q] */
  1.1609 +#ifdef GLP_DEBUG
  1.1610 +      xassert(teta >= 0.0);
  1.1611 +#endif
  1.1612 +      csa->teta = s * teta;
  1.1613 +      return;
  1.1614 +}
  1.1615 +
  1.1616 +/***********************************************************************
  1.1617 +*  eval_rho - compute pivot row of the inverse
  1.1618 +*
  1.1619 +*  This routine computes the pivot (p-th) row of the inverse inv(B),
  1.1620 +*  which corresponds to basic variable xB[p] chosen:
  1.1621 +*
  1.1622 +*     rho = inv(B') * e[p],
  1.1623 +*
  1.1624 +*  where B' is a matrix transposed to the current basis matrix, e[p]
  1.1625 +*  is unity vector. */
  1.1626 +
  1.1627 +static void eval_rho(struct csa *csa, double rho[])
  1.1628 +{     int m = csa->m;
  1.1629 +      int p = csa->p;
  1.1630 +      double *e = rho;
  1.1631 +      int i;
  1.1632 +#ifdef GLP_DEBUG
  1.1633 +      xassert(1 <= p && p <= m);
  1.1634 +#endif
  1.1635 +      /* construct the right-hand side vector e[p] */
  1.1636 +      for (i = 1; i <= m; i++)
  1.1637 +         e[i] = 0.0;
  1.1638 +      e[p] = 1.0;
  1.1639 +      /* solve system B'* rho = e[p] */
  1.1640 +      xassert(csa->valid);
  1.1641 +      bfd_btran(csa->bfd, rho);
  1.1642 +      return;
  1.1643 +}
  1.1644 +
  1.1645 +/***********************************************************************
  1.1646 +*  refine_rho - refine pivot row of the inverse
  1.1647 +*
  1.1648 +*  This routine refines the pivot row of the inverse inv(B) assuming
  1.1649 +*  that it was previously computed by the routine eval_rho. */
  1.1650 +
  1.1651 +static void refine_rho(struct csa *csa, double rho[])
  1.1652 +{     int m = csa->m;
  1.1653 +      int p = csa->p;
  1.1654 +      double *e = csa->work3;
  1.1655 +      int i;
  1.1656 +#ifdef GLP_DEBUG
  1.1657 +      xassert(1 <= p && p <= m);
  1.1658 +#endif
  1.1659 +      /* construct the right-hand side vector e[p] */
  1.1660 +      for (i = 1; i <= m; i++)
  1.1661 +         e[i] = 0.0;
  1.1662 +      e[p] = 1.0;
  1.1663 +      /* refine solution of B'* rho = e[p] */
  1.1664 +      refine_btran(csa, e, rho);
  1.1665 +      return;
  1.1666 +}
  1.1667 +
  1.1668 +/***********************************************************************
  1.1669 +*  eval_trow - compute pivot row of the simplex table
  1.1670 +*
  1.1671 +*  This routine computes the pivot row of the simplex table, which
  1.1672 +*  corresponds to basic variable xB[p] chosen.
  1.1673 +*
  1.1674 +*  The pivot row is the following vector:
  1.1675 +*
  1.1676 +*     trow = T'* e[p] = - N'* inv(B') * e[p] = - N' * rho,
  1.1677 +*
  1.1678 +*  where rho is the pivot row of the inverse inv(B) previously computed
  1.1679 +*  by the routine eval_rho.
  1.1680 +*
  1.1681 +*  Note that elements of the pivot row corresponding to fixed non-basic
  1.1682 +*  variables are not computed. */
  1.1683 +
  1.1684 +static void eval_trow(struct csa *csa, double rho[])
  1.1685 +{     int m = csa->m;
  1.1686 +      int n = csa->n;
  1.1687 +#ifdef GLP_DEBUG
  1.1688 +      char *stat = csa->stat;
  1.1689 +#endif
  1.1690 +      int *N_ptr = csa->N_ptr;
  1.1691 +      int *N_len = csa->N_len;
  1.1692 +      int *N_ind = csa->N_ind;
  1.1693 +      double *N_val = csa->N_val;
  1.1694 +      int *trow_ind = csa->trow_ind;
  1.1695 +      double *trow_vec = csa->trow_vec;
  1.1696 +      int i, j, beg, end, ptr, nnz;
  1.1697 +      double temp;
  1.1698 +      /* clear the pivot row */
  1.1699 +      for (j = 1; j <= n; j++)
  1.1700 +         trow_vec[j] = 0.0;
  1.1701 +      /* compute the pivot row as a linear combination of rows of the
  1.1702 +         matrix N: trow = - rho[1] * N'[1] - ... - rho[m] * N'[m] */
  1.1703 +      for (i = 1; i <= m; i++)
  1.1704 +      {  temp = rho[i];
  1.1705 +         if (temp == 0.0) continue;
  1.1706 +         /* trow := trow - rho[i] * N'[i] */
  1.1707 +         beg = N_ptr[i];
  1.1708 +         end = beg + N_len[i];
  1.1709 +         for (ptr = beg; ptr < end; ptr++)
  1.1710 +         {
  1.1711 +#ifdef GLP_DEBUG
  1.1712 +            j = N_ind[ptr];
  1.1713 +            xassert(1 <= j && j <= n);
  1.1714 +            xassert(stat[j] != GLP_NS);
  1.1715 +#endif
  1.1716 +            trow_vec[N_ind[ptr]] -= temp * N_val[ptr];
  1.1717 +         }
  1.1718 +      }
  1.1719 +      /* construct sparse pattern of the pivot row */
  1.1720 +      nnz = 0;
  1.1721 +      for (j = 1; j <= n; j++)
  1.1722 +      {  if (trow_vec[j] != 0.0)
  1.1723 +            trow_ind[++nnz] = j;
  1.1724 +      }
  1.1725 +      csa->trow_nnz = nnz;
  1.1726 +      return;
  1.1727 +}
  1.1728 +
  1.1729 +/***********************************************************************
  1.1730 +*  update_bbar - update values of basic variables
  1.1731 +*
  1.1732 +*  This routine updates values of all basic variables for the adjacent
  1.1733 +*  basis. */
  1.1734 +
  1.1735 +static void update_bbar(struct csa *csa)
  1.1736 +{
  1.1737 +#ifdef GLP_DEBUG
  1.1738 +      int m = csa->m;
  1.1739 +      int n = csa->n;
  1.1740 +#endif
  1.1741 +      double *bbar = csa->bbar;
  1.1742 +      int q = csa->q;
  1.1743 +      int tcol_nnz = csa->tcol_nnz;
  1.1744 +      int *tcol_ind = csa->tcol_ind;
  1.1745 +      double *tcol_vec = csa->tcol_vec;
  1.1746 +      int p = csa->p;
  1.1747 +      double teta = csa->teta;
  1.1748 +      int i, pos;
  1.1749 +#ifdef GLP_DEBUG
  1.1750 +      xassert(1 <= q && q <= n);
  1.1751 +      xassert(p < 0 || 1 <= p && p <= m);
  1.1752 +#endif
  1.1753 +      /* if xN[q] leaves the basis, compute its value in the adjacent
  1.1754 +         basis, where it will replace xB[p] */
  1.1755 +      if (p > 0)
  1.1756 +         bbar[p] = get_xN(csa, q) + teta;
  1.1757 +      /* update values of other basic variables (except xB[p], because
  1.1758 +         it will be replaced by xN[q]) */
  1.1759 +      if (teta == 0.0) goto done;
  1.1760 +      for (pos = 1; pos <= tcol_nnz; pos++)
  1.1761 +      {  i = tcol_ind[pos];
  1.1762 +         /* skip xB[p] */
  1.1763 +         if (i == p) continue;
  1.1764 +         /* (change of xB[i]) = alfa[i,q] * (change of xN[q]) */
  1.1765 +         bbar[i] += tcol_vec[i] * teta;
  1.1766 +      }
  1.1767 +done: return;
  1.1768 +}
  1.1769 +
  1.1770 +/***********************************************************************
  1.1771 +*  reeval_cost - recompute reduced cost of non-basic variable xN[q]
  1.1772 +*
  1.1773 +*  This routine recomputes reduced cost of non-basic variable xN[q] for
  1.1774 +*  the current basis more accurately using its direct definition:
  1.1775 +*
  1.1776 +*     d[q] = cN[q] - N'[q] * pi =
  1.1777 +*
  1.1778 +*          = cN[q] - N'[q] * (inv(B') * cB) =
  1.1779 +*
  1.1780 +*          = cN[q] - (cB' * inv(B) * N[q]) =
  1.1781 +*
  1.1782 +*          = cN[q] + cB' * (pivot column).
  1.1783 +*
  1.1784 +*  It is assumed that the pivot column of the simplex table is already
  1.1785 +*  computed. */
  1.1786 +
  1.1787 +static double reeval_cost(struct csa *csa)
  1.1788 +{     int m = csa->m;
  1.1789 +#ifdef GLP_DEBUG
  1.1790 +      int n = csa->n;
  1.1791 +#endif
  1.1792 +      double *coef = csa->coef;
  1.1793 +      int *head = csa->head;
  1.1794 +      int q = csa->q;
  1.1795 +      int tcol_nnz = csa->tcol_nnz;
  1.1796 +      int *tcol_ind = csa->tcol_ind;
  1.1797 +      double *tcol_vec = csa->tcol_vec;
  1.1798 +      int i, pos;
  1.1799 +      double dq;
  1.1800 +#ifdef GLP_DEBUG
  1.1801 +      xassert(1 <= q && q <= n);
  1.1802 +#endif
  1.1803 +      dq = coef[head[m+q]];
  1.1804 +      for (pos = 1; pos <= tcol_nnz; pos++)
  1.1805 +      {  i = tcol_ind[pos];
  1.1806 +#ifdef GLP_DEBUG
  1.1807 +         xassert(1 <= i && i <= m);
  1.1808 +#endif
  1.1809 +         dq += coef[head[i]] * tcol_vec[i];
  1.1810 +      }
  1.1811 +      return dq;
  1.1812 +}
  1.1813 +
  1.1814 +/***********************************************************************
  1.1815 +*  update_cbar - update reduced costs of non-basic variables
  1.1816 +*
  1.1817 +*  This routine updates reduced costs of all (except fixed) non-basic
  1.1818 +*  variables for the adjacent basis. */
  1.1819 +
  1.1820 +static void update_cbar(struct csa *csa)
  1.1821 +{
  1.1822 +#ifdef GLP_DEBUG
  1.1823 +      int n = csa->n;
  1.1824 +#endif
  1.1825 +      double *cbar = csa->cbar;
  1.1826 +      int q = csa->q;
  1.1827 +      int trow_nnz = csa->trow_nnz;
  1.1828 +      int *trow_ind = csa->trow_ind;
  1.1829 +      double *trow_vec = csa->trow_vec;
  1.1830 +      int j, pos;
  1.1831 +      double new_dq;
  1.1832 +#ifdef GLP_DEBUG
  1.1833 +      xassert(1 <= q && q <= n);
  1.1834 +#endif
  1.1835 +      /* compute reduced cost of xB[p] in the adjacent basis, where it
  1.1836 +         will replace xN[q] */
  1.1837 +#ifdef GLP_DEBUG
  1.1838 +      xassert(trow_vec[q] != 0.0);
  1.1839 +#endif
  1.1840 +      new_dq = (cbar[q] /= trow_vec[q]);
  1.1841 +      /* update reduced costs of other non-basic variables (except
  1.1842 +         xN[q], because it will be replaced by xB[p]) */
  1.1843 +      for (pos = 1; pos <= trow_nnz; pos++)
  1.1844 +      {  j = trow_ind[pos];
  1.1845 +         /* skip xN[q] */
  1.1846 +         if (j == q) continue;
  1.1847 +         cbar[j] -= trow_vec[j] * new_dq;
  1.1848 +      }
  1.1849 +      return;
  1.1850 +}
  1.1851 +
  1.1852 +/***********************************************************************
  1.1853 +*  update_gamma - update steepest edge coefficients
  1.1854 +*
  1.1855 +*  This routine updates steepest-edge coefficients for the adjacent
  1.1856 +*  basis. */
  1.1857 +
  1.1858 +static void update_gamma(struct csa *csa)
  1.1859 +{     int m = csa->m;
  1.1860 +#ifdef GLP_DEBUG
  1.1861 +      int n = csa->n;
  1.1862 +#endif
  1.1863 +      char *type = csa->type;
  1.1864 +      int *A_ptr = csa->A_ptr;
  1.1865 +      int *A_ind = csa->A_ind;
  1.1866 +      double *A_val = csa->A_val;
  1.1867 +      int *head = csa->head;
  1.1868 +      char *refsp = csa->refsp;
  1.1869 +      double *gamma = csa->gamma;
  1.1870 +      int q = csa->q;
  1.1871 +      int tcol_nnz = csa->tcol_nnz;
  1.1872 +      int *tcol_ind = csa->tcol_ind;
  1.1873 +      double *tcol_vec = csa->tcol_vec;
  1.1874 +      int p = csa->p;
  1.1875 +      int trow_nnz = csa->trow_nnz;
  1.1876 +      int *trow_ind = csa->trow_ind;
  1.1877 +      double *trow_vec = csa->trow_vec;
  1.1878 +      double *u = csa->work3;
  1.1879 +      int i, j, k, pos, beg, end, ptr;
  1.1880 +      double gamma_q, delta_q, pivot, s, t, t1, t2;
  1.1881 +#ifdef GLP_DEBUG
  1.1882 +      xassert(1 <= p && p <= m);
  1.1883 +      xassert(1 <= q && q <= n);
  1.1884 +#endif
  1.1885 +      /* the basis changes, so decrease the count */
  1.1886 +      xassert(csa->refct > 0);
  1.1887 +      csa->refct--;
  1.1888 +      /* recompute gamma[q] for the current basis more accurately and
  1.1889 +         compute auxiliary vector u */
  1.1890 +      gamma_q = delta_q = (refsp[head[m+q]] ? 1.0 : 0.0);
  1.1891 +      for (i = 1; i <= m; i++) u[i] = 0.0;
  1.1892 +      for (pos = 1; pos <= tcol_nnz; pos++)
  1.1893 +      {  i = tcol_ind[pos];
  1.1894 +         if (refsp[head[i]])
  1.1895 +         {  u[i] = t = tcol_vec[i];
  1.1896 +            gamma_q += t * t;
  1.1897 +         }
  1.1898 +         else
  1.1899 +            u[i] = 0.0;
  1.1900 +      }
  1.1901 +      xassert(csa->valid);
  1.1902 +      bfd_btran(csa->bfd, u);
  1.1903 +      /* update gamma[k] for other non-basic variables (except fixed
  1.1904 +         variables and xN[q], because it will be replaced by xB[p]) */
  1.1905 +      pivot = trow_vec[q];
  1.1906 +#ifdef GLP_DEBUG
  1.1907 +      xassert(pivot != 0.0);
  1.1908 +#endif
  1.1909 +      for (pos = 1; pos <= trow_nnz; pos++)
  1.1910 +      {  j = trow_ind[pos];
  1.1911 +         /* skip xN[q] */
  1.1912 +         if (j == q) continue;
  1.1913 +         /* compute t */
  1.1914 +         t = trow_vec[j] / pivot;
  1.1915 +         /* compute inner product s = N'[j] * u */
  1.1916 +         k = head[m+j]; /* x[k] = xN[j] */
  1.1917 +         if (k <= m)
  1.1918 +            s = u[k];
  1.1919 +         else
  1.1920 +         {  s = 0.0;
  1.1921 +            beg = A_ptr[k-m];
  1.1922 +            end = A_ptr[k-m+1];
  1.1923 +            for (ptr = beg; ptr < end; ptr++)
  1.1924 +               s -= A_val[ptr] * u[A_ind[ptr]];
  1.1925 +         }
  1.1926 +         /* compute gamma[k] for the adjacent basis */
  1.1927 +         t1 = gamma[j] + t * t * gamma_q + 2.0 * t * s;
  1.1928 +         t2 = (refsp[k] ? 1.0 : 0.0) + delta_q * t * t;
  1.1929 +         gamma[j] = (t1 >= t2 ? t1 : t2);
  1.1930 +         if (gamma[j] < DBL_EPSILON) gamma[j] = DBL_EPSILON;
  1.1931 +      }
  1.1932 +      /* compute gamma[q] for the adjacent basis */
  1.1933 +      if (type[head[p]] == GLP_FX)
  1.1934 +         gamma[q] = 1.0;
  1.1935 +      else
  1.1936 +      {  gamma[q] = gamma_q / (pivot * pivot);
  1.1937 +         if (gamma[q] < DBL_EPSILON) gamma[q] = DBL_EPSILON;
  1.1938 +      }
  1.1939 +      return;
  1.1940 +}
  1.1941 +
  1.1942 +/***********************************************************************
  1.1943 +*  err_in_bbar - compute maximal relative error in primal solution
  1.1944 +*
  1.1945 +*  This routine returns maximal relative error:
  1.1946 +*
  1.1947 +*     max |beta[i] - bbar[i]| / (1 + |beta[i]|),
  1.1948 +*
  1.1949 +*  where beta and bbar are, respectively, directly computed and the
  1.1950 +*  current (updated) values of basic variables.
  1.1951 +*
  1.1952 +*  NOTE: The routine is intended only for debugginig purposes. */
  1.1953 +
  1.1954 +static double err_in_bbar(struct csa *csa)
  1.1955 +{     int m = csa->m;
  1.1956 +      double *bbar = csa->bbar;
  1.1957 +      int i;
  1.1958 +      double e, emax, *beta;
  1.1959 +      beta = xcalloc(1+m, sizeof(double));
  1.1960 +      eval_beta(csa, beta);
  1.1961 +      emax = 0.0;
  1.1962 +      for (i = 1; i <= m; i++)
  1.1963 +      {  e = fabs(beta[i] - bbar[i]) / (1.0 + fabs(beta[i]));
  1.1964 +         if (emax < e) emax = e;
  1.1965 +      }
  1.1966 +      xfree(beta);
  1.1967 +      return emax;
  1.1968 +}
  1.1969 +
  1.1970 +/***********************************************************************
  1.1971 +*  err_in_cbar - compute maximal relative error in dual solution
  1.1972 +*
  1.1973 +*  This routine returns maximal relative error:
  1.1974 +*
  1.1975 +*     max |cost[j] - cbar[j]| / (1 + |cost[j]|),
  1.1976 +*
  1.1977 +*  where cost and cbar are, respectively, directly computed and the
  1.1978 +*  current (updated) reduced costs of non-basic non-fixed variables.
  1.1979 +*
  1.1980 +*  NOTE: The routine is intended only for debugginig purposes. */
  1.1981 +
  1.1982 +static double err_in_cbar(struct csa *csa)
  1.1983 +{     int m = csa->m;
  1.1984 +      int n = csa->n;
  1.1985 +      char *stat = csa->stat;
  1.1986 +      double *cbar = csa->cbar;
  1.1987 +      int j;
  1.1988 +      double e, emax, cost, *pi;
  1.1989 +      pi = xcalloc(1+m, sizeof(double));
  1.1990 +      eval_pi(csa, pi);
  1.1991 +      emax = 0.0;
  1.1992 +      for (j = 1; j <= n; j++)
  1.1993 +      {  if (stat[j] == GLP_NS) continue;
  1.1994 +         cost = eval_cost(csa, pi, j);
  1.1995 +         e = fabs(cost - cbar[j]) / (1.0 + fabs(cost));
  1.1996 +         if (emax < e) emax = e;
  1.1997 +      }
  1.1998 +      xfree(pi);
  1.1999 +      return emax;
  1.2000 +}
  1.2001 +
  1.2002 +/***********************************************************************
  1.2003 +*  err_in_gamma - compute maximal relative error in steepest edge cff.
  1.2004 +*
  1.2005 +*  This routine returns maximal relative error:
  1.2006 +*
  1.2007 +*     max |gamma'[j] - gamma[j]| / (1 + |gamma'[j]),
  1.2008 +*
  1.2009 +*  where gamma'[j] and gamma[j] are, respectively, directly computed
  1.2010 +*  and the current (updated) steepest edge coefficients for non-basic
  1.2011 +*  non-fixed variable x[j].
  1.2012 +*
  1.2013 +*  NOTE: The routine is intended only for debugginig purposes. */
  1.2014 +
  1.2015 +static double err_in_gamma(struct csa *csa)
  1.2016 +{     int n = csa->n;
  1.2017 +      char *stat = csa->stat;
  1.2018 +      double *gamma = csa->gamma;
  1.2019 +      int j;
  1.2020 +      double e, emax, temp;
  1.2021 +      emax = 0.0;
  1.2022 +      for (j = 1; j <= n; j++)
  1.2023 +      {  if (stat[j] == GLP_NS)
  1.2024 +         {  xassert(gamma[j] == 1.0);
  1.2025 +            continue;
  1.2026 +         }
  1.2027 +         temp = eval_gamma(csa, j);
  1.2028 +         e = fabs(temp - gamma[j]) / (1.0 + fabs(temp));
  1.2029 +         if (emax < e) emax = e;
  1.2030 +      }
  1.2031 +      return emax;
  1.2032 +}
  1.2033 +
  1.2034 +/***********************************************************************
  1.2035 +*  change_basis - change basis header
  1.2036 +*
  1.2037 +*  This routine changes the basis header to make it corresponding to
  1.2038 +*  the adjacent basis. */
  1.2039 +
  1.2040 +static void change_basis(struct csa *csa)
  1.2041 +{     int m = csa->m;
  1.2042 +#ifdef GLP_DEBUG
  1.2043 +      int n = csa->n;
  1.2044 +      char *type = csa->type;
  1.2045 +#endif
  1.2046 +      int *head = csa->head;
  1.2047 +      char *stat = csa->stat;
  1.2048 +      int q = csa->q;
  1.2049 +      int p = csa->p;
  1.2050 +      int p_stat = csa->p_stat;
  1.2051 +      int k;
  1.2052 +#ifdef GLP_DEBUG
  1.2053 +      xassert(1 <= q && q <= n);
  1.2054 +#endif
  1.2055 +      if (p < 0)
  1.2056 +      {  /* xN[q] goes to its opposite bound */
  1.2057 +#ifdef GLP_DEBUG
  1.2058 +         k = head[m+q]; /* x[k] = xN[q] */
  1.2059 +         xassert(1 <= k && k <= m+n);
  1.2060 +         xassert(type[k] == GLP_DB);
  1.2061 +#endif
  1.2062 +         switch (stat[q])
  1.2063 +         {  case GLP_NL:
  1.2064 +               /* xN[q] increases */
  1.2065 +               stat[q] = GLP_NU;
  1.2066 +               break;
  1.2067 +            case GLP_NU:
  1.2068 +               /* xN[q] decreases */
  1.2069 +               stat[q] = GLP_NL;
  1.2070 +               break;
  1.2071 +            default:
  1.2072 +               xassert(stat != stat);
  1.2073 +         }
  1.2074 +      }
  1.2075 +      else
  1.2076 +      {  /* xB[p] leaves the basis, xN[q] enters the basis */
  1.2077 +#ifdef GLP_DEBUG
  1.2078 +         xassert(1 <= p && p <= m);
  1.2079 +         k = head[p]; /* x[k] = xB[p] */
  1.2080 +         switch (p_stat)
  1.2081 +         {  case GLP_NL:
  1.2082 +               /* xB[p] goes to its lower bound */
  1.2083 +               xassert(type[k] == GLP_LO || type[k] == GLP_DB);
  1.2084 +               break;
  1.2085 +            case GLP_NU:
  1.2086 +               /* xB[p] goes to its upper bound */
  1.2087 +               xassert(type[k] == GLP_UP || type[k] == GLP_DB);
  1.2088 +               break;
  1.2089 +            case GLP_NS:
  1.2090 +               /* xB[p] goes to its fixed value */
  1.2091 +               xassert(type[k] == GLP_NS);
  1.2092 +               break;
  1.2093 +            default:
  1.2094 +               xassert(p_stat != p_stat);
  1.2095 +         }
  1.2096 +#endif
  1.2097 +         /* xB[p] <-> xN[q] */
  1.2098 +         k = head[p], head[p] = head[m+q], head[m+q] = k;
  1.2099 +         stat[q] = (char)p_stat;
  1.2100 +      }
  1.2101 +      return;
  1.2102 +}
  1.2103 +
  1.2104 +/***********************************************************************
  1.2105 +*  set_aux_obj - construct auxiliary objective function
  1.2106 +*
  1.2107 +*  The auxiliary objective function is a separable piecewise linear
  1.2108 +*  convex function, which is the sum of primal infeasibilities:
  1.2109 +*
  1.2110 +*     z = t[1] + ... + t[m+n] -> minimize,
  1.2111 +*
  1.2112 +*  where:
  1.2113 +*
  1.2114 +*            / lb[k] - x[k], if x[k] < lb[k]
  1.2115 +*            |
  1.2116 +*     t[k] = <  0, if lb[k] <= x[k] <= ub[k]
  1.2117 +*            |
  1.2118 +*            \ x[k] - ub[k], if x[k] > ub[k]
  1.2119 +*
  1.2120 +*  This routine computes objective coefficients for the current basis
  1.2121 +*  and returns the number of non-zero terms t[k]. */
  1.2122 +
  1.2123 +static int set_aux_obj(struct csa *csa, double tol_bnd)
  1.2124 +{     int m = csa->m;
  1.2125 +      int n = csa->n;
  1.2126 +      char *type = csa->type;
  1.2127 +      double *lb = csa->lb;
  1.2128 +      double *ub = csa->ub;
  1.2129 +      double *coef = csa->coef;
  1.2130 +      int *head = csa->head;
  1.2131 +      double *bbar = csa->bbar;
  1.2132 +      int i, k, cnt = 0;
  1.2133 +      double eps;
  1.2134 +      /* use a bit more restrictive tolerance */
  1.2135 +      tol_bnd *= 0.90;
  1.2136 +      /* clear all objective coefficients */
  1.2137 +      for (k = 1; k <= m+n; k++)
  1.2138 +         coef[k] = 0.0;
  1.2139 +      /* walk through the list of basic variables */
  1.2140 +      for (i = 1; i <= m; i++)
  1.2141 +      {  k = head[i]; /* x[k] = xB[i] */
  1.2142 +         if (type[k] == GLP_LO || type[k] == GLP_DB ||
  1.2143 +             type[k] == GLP_FX)
  1.2144 +         {  /* x[k] has lower bound */
  1.2145 +            eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
  1.2146 +            if (bbar[i] < lb[k] - eps)
  1.2147 +            {  /* and violates it */
  1.2148 +               coef[k] = -1.0;
  1.2149 +               cnt++;
  1.2150 +            }
  1.2151 +         }
  1.2152 +         if (type[k] == GLP_UP || type[k] == GLP_DB ||
  1.2153 +             type[k] == GLP_FX)
  1.2154 +         {  /* x[k] has upper bound */
  1.2155 +            eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
  1.2156 +            if (bbar[i] > ub[k] + eps)
  1.2157 +            {  /* and violates it */
  1.2158 +               coef[k] = +1.0;
  1.2159 +               cnt++;
  1.2160 +            }
  1.2161 +         }
  1.2162 +      }
  1.2163 +      return cnt;
  1.2164 +}
  1.2165 +
  1.2166 +/***********************************************************************
  1.2167 +*  set_orig_obj - restore original objective function
  1.2168 +*
  1.2169 +*  This routine assigns scaled original objective coefficients to the
  1.2170 +*  working objective function. */
  1.2171 +
  1.2172 +static void set_orig_obj(struct csa *csa)
  1.2173 +{     int m = csa->m;
  1.2174 +      int n = csa->n;
  1.2175 +      double *coef = csa->coef;
  1.2176 +      double *obj = csa->obj;
  1.2177 +      double zeta = csa->zeta;
  1.2178 +      int i, j;
  1.2179 +      for (i = 1; i <= m; i++)
  1.2180 +         coef[i] = 0.0;
  1.2181 +      for (j = 1; j <= n; j++)
  1.2182 +         coef[m+j] = zeta * obj[j];
  1.2183 +      return;
  1.2184 +}
  1.2185 +
  1.2186 +/***********************************************************************
  1.2187 +*  check_stab - check numerical stability of basic solution
  1.2188 +*
  1.2189 +*  If the current basic solution is primal feasible (or pseudo feasible
  1.2190 +*  on phase I) within a tolerance, this routine returns zero, otherwise
  1.2191 +*  it returns non-zero. */
  1.2192 +
  1.2193 +static int check_stab(struct csa *csa, double tol_bnd)
  1.2194 +{     int m = csa->m;
  1.2195 +#ifdef GLP_DEBUG
  1.2196 +      int n = csa->n;
  1.2197 +#endif
  1.2198 +      char *type = csa->type;
  1.2199 +      double *lb = csa->lb;
  1.2200 +      double *ub = csa->ub;
  1.2201 +      double *coef = csa->coef;
  1.2202 +      int *head = csa->head;
  1.2203 +      int phase = csa->phase;
  1.2204 +      double *bbar = csa->bbar;
  1.2205 +      int i, k;
  1.2206 +      double eps;
  1.2207 +      /* walk through the list of basic variables */
  1.2208 +      for (i = 1; i <= m; i++)
  1.2209 +      {  k = head[i]; /* x[k] = xB[i] */
  1.2210 +#ifdef GLP_DEBUG
  1.2211 +         xassert(1 <= k && k <= m+n);
  1.2212 +#endif
  1.2213 +         if (phase == 1 && coef[k] < 0.0)
  1.2214 +         {  /* x[k] must not be greater than its lower bound */
  1.2215 +#ifdef GLP_DEBUG
  1.2216 +            xassert(type[k] == GLP_LO || type[k] == GLP_DB ||
  1.2217 +                    type[k] == GLP_FX);
  1.2218 +#endif
  1.2219 +            eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
  1.2220 +            if (bbar[i] > lb[k] + eps) return 1;
  1.2221 +         }
  1.2222 +         else if (phase == 1 && coef[k] > 0.0)
  1.2223 +         {  /* x[k] must not be less than its upper bound */
  1.2224 +#ifdef GLP_DEBUG
  1.2225 +            xassert(type[k] == GLP_UP || type[k] == GLP_DB ||
  1.2226 +                    type[k] == GLP_FX);
  1.2227 +#endif
  1.2228 +            eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
  1.2229 +            if (bbar[i] < ub[k] - eps) return 1;
  1.2230 +         }
  1.2231 +         else
  1.2232 +         {  /* either phase = 1 and coef[k] = 0, or phase = 2 */
  1.2233 +            if (type[k] == GLP_LO || type[k] == GLP_DB ||
  1.2234 +                type[k] == GLP_FX)
  1.2235 +            {  /* x[k] must not be less than its lower bound */
  1.2236 +               eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
  1.2237 +               if (bbar[i] < lb[k] - eps) return 1;
  1.2238 +            }
  1.2239 +            if (type[k] == GLP_UP || type[k] == GLP_DB ||
  1.2240 +                type[k] == GLP_FX)
  1.2241 +            {  /* x[k] must not be greater then its upper bound */
  1.2242 +               eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
  1.2243 +               if (bbar[i] > ub[k] + eps) return 1;
  1.2244 +            }
  1.2245 +         }
  1.2246 +      }
  1.2247 +      /* basic solution is primal feasible within a tolerance */
  1.2248 +      return 0;
  1.2249 +}
  1.2250 +
  1.2251 +/***********************************************************************
  1.2252 +*  check_feas - check primal feasibility of basic solution
  1.2253 +*
  1.2254 +*  If the current basic solution is primal feasible within a tolerance,
  1.2255 +*  this routine returns zero, otherwise it returns non-zero. */
  1.2256 +
  1.2257 +static int check_feas(struct csa *csa, double tol_bnd)
  1.2258 +{     int m = csa->m;
  1.2259 +#ifdef GLP_DEBUG
  1.2260 +      int n = csa->n;
  1.2261 +      char *type = csa->type;
  1.2262 +#endif
  1.2263 +      double *lb = csa->lb;
  1.2264 +      double *ub = csa->ub;
  1.2265 +      double *coef = csa->coef;
  1.2266 +      int *head = csa->head;
  1.2267 +      double *bbar = csa->bbar;
  1.2268 +      int i, k;
  1.2269 +      double eps;
  1.2270 +      xassert(csa->phase == 1);
  1.2271 +      /* walk through the list of basic variables */
  1.2272 +      for (i = 1; i <= m; i++)
  1.2273 +      {  k = head[i]; /* x[k] = xB[i] */
  1.2274 +#ifdef GLP_DEBUG
  1.2275 +         xassert(1 <= k && k <= m+n);
  1.2276 +#endif
  1.2277 +         if (coef[k] < 0.0)
  1.2278 +         {  /* check if x[k] still violates its lower bound */
  1.2279 +#ifdef GLP_DEBUG
  1.2280 +            xassert(type[k] == GLP_LO || type[k] == GLP_DB ||
  1.2281 +                    type[k] == GLP_FX);
  1.2282 +#endif
  1.2283 +            eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
  1.2284 +            if (bbar[i] < lb[k] - eps) return 1;
  1.2285 +         }
  1.2286 +         else if (coef[k] > 0.0)
  1.2287 +         {  /* check if x[k] still violates its upper bound */
  1.2288 +#ifdef GLP_DEBUG
  1.2289 +            xassert(type[k] == GLP_UP || type[k] == GLP_DB ||
  1.2290 +                    type[k] == GLP_FX);
  1.2291 +#endif
  1.2292 +            eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
  1.2293 +            if (bbar[i] > ub[k] + eps) return 1;
  1.2294 +         }
  1.2295 +      }
  1.2296 +      /* basic solution is primal feasible within a tolerance */
  1.2297 +      return 0;
  1.2298 +}
  1.2299 +
  1.2300 +/***********************************************************************
  1.2301 +*  eval_obj - compute original objective function
  1.2302 +*
  1.2303 +*  This routine computes the current value of the original objective
  1.2304 +*  function. */
  1.2305 +
  1.2306 +static double eval_obj(struct csa *csa)
  1.2307 +{     int m = csa->m;
  1.2308 +      int n = csa->n;
  1.2309 +      double *obj = csa->obj;
  1.2310 +      int *head = csa->head;
  1.2311 +      double *bbar = csa->bbar;
  1.2312 +      int i, j, k;
  1.2313 +      double sum;
  1.2314 +      sum = obj[0];
  1.2315 +      /* walk through the list of basic variables */
  1.2316 +      for (i = 1; i <= m; i++)
  1.2317 +      {  k = head[i]; /* x[k] = xB[i] */
  1.2318 +#ifdef GLP_DEBUG
  1.2319 +         xassert(1 <= k && k <= m+n);
  1.2320 +#endif
  1.2321 +         if (k > m)
  1.2322 +            sum += obj[k-m] * bbar[i];
  1.2323 +      }
  1.2324 +      /* walk through the list of non-basic variables */
  1.2325 +      for (j = 1; j <= n; j++)
  1.2326 +      {  k = head[m+j]; /* x[k] = xN[j] */
  1.2327 +#ifdef GLP_DEBUG
  1.2328 +         xassert(1 <= k && k <= m+n);
  1.2329 +#endif
  1.2330 +         if (k > m)
  1.2331 +            sum += obj[k-m] * get_xN(csa, j);
  1.2332 +      }
  1.2333 +      return sum;
  1.2334 +}
  1.2335 +
  1.2336 +/***********************************************************************
  1.2337 +*  display - display the search progress
  1.2338 +*
  1.2339 +*  This routine displays some information about the search progress
  1.2340 +*  that includes:
  1.2341 +*
  1.2342 +*  the search phase;
  1.2343 +*
  1.2344 +*  the number of simplex iterations performed by the solver;
  1.2345 +*
  1.2346 +*  the original objective value;
  1.2347 +*
  1.2348 +*  the sum of (scaled) primal infeasibilities;
  1.2349 +*
  1.2350 +*  the number of basic fixed variables. */
  1.2351 +
  1.2352 +static void display(struct csa *csa, const glp_smcp *parm, int spec)
  1.2353 +{     int m = csa->m;
  1.2354 +#ifdef GLP_DEBUG
  1.2355 +      int n = csa->n;
  1.2356 +#endif
  1.2357 +      char *type = csa->type;
  1.2358 +      double *lb = csa->lb;
  1.2359 +      double *ub = csa->ub;
  1.2360 +      int phase = csa->phase;
  1.2361 +      int *head = csa->head;
  1.2362 +      double *bbar = csa->bbar;
  1.2363 +      int i, k, cnt;
  1.2364 +      double sum;
  1.2365 +      if (parm->msg_lev < GLP_MSG_ON) goto skip;
  1.2366 +      if (parm->out_dly > 0 &&
  1.2367 +         1000.0 * xdifftime(xtime(), csa->tm_beg) < parm->out_dly)
  1.2368 +         goto skip;
  1.2369 +      if (csa->it_cnt == csa->it_dpy) goto skip;
  1.2370 +      if (!spec && csa->it_cnt % parm->out_frq != 0) goto skip;
  1.2371 +      /* compute the sum of primal infeasibilities and determine the
  1.2372 +         number of basic fixed variables */
  1.2373 +      sum = 0.0, cnt = 0;
  1.2374 +      for (i = 1; i <= m; i++)
  1.2375 +      {  k = head[i]; /* x[k] = xB[i] */
  1.2376 +#ifdef GLP_DEBUG
  1.2377 +         xassert(1 <= k && k <= m+n);
  1.2378 +#endif
  1.2379 +         if (type[k] == GLP_LO || type[k] == GLP_DB ||
  1.2380 +             type[k] == GLP_FX)
  1.2381 +         {  /* x[k] has lower bound */
  1.2382 +            if (bbar[i] < lb[k])
  1.2383 +               sum += (lb[k] - bbar[i]);
  1.2384 +         }
  1.2385 +         if (type[k] == GLP_UP || type[k] == GLP_DB ||
  1.2386 +             type[k] == GLP_FX)
  1.2387 +         {  /* x[k] has upper bound */
  1.2388 +            if (bbar[i] > ub[k])
  1.2389 +               sum += (bbar[i] - ub[k]);
  1.2390 +         }
  1.2391 +         if (type[k] == GLP_FX) cnt++;
  1.2392 +      }
  1.2393 +      xprintf("%c%6d: obj = %17.9e  infeas = %10.3e (%d)\n",
  1.2394 +         phase == 1 ? ' ' : '*', csa->it_cnt, eval_obj(csa), sum, cnt);
  1.2395 +      csa->it_dpy = csa->it_cnt;
  1.2396 +skip: return;
  1.2397 +}
  1.2398 +
  1.2399 +/***********************************************************************
  1.2400 +*  store_sol - store basic solution back to the problem object
  1.2401 +*
  1.2402 +*  This routine stores basic solution components back to the problem
  1.2403 +*  object. */
  1.2404 +
  1.2405 +static void store_sol(struct csa *csa, glp_prob *lp, int p_stat,
  1.2406 +      int d_stat, int ray)
  1.2407 +{     int m = csa->m;
  1.2408 +      int n = csa->n;
  1.2409 +      double zeta = csa->zeta;
  1.2410 +      int *head = csa->head;
  1.2411 +      char *stat = csa->stat;
  1.2412 +      double *bbar = csa->bbar;
  1.2413 +      double *cbar = csa->cbar;
  1.2414 +      int i, j, k;
  1.2415 +#ifdef GLP_DEBUG
  1.2416 +      xassert(lp->m == m);
  1.2417 +      xassert(lp->n == n);
  1.2418 +#endif
  1.2419 +      /* basis factorization */
  1.2420 +#ifdef GLP_DEBUG
  1.2421 +      xassert(!lp->valid && lp->bfd == NULL);
  1.2422 +      xassert(csa->valid && csa->bfd != NULL);
  1.2423 +#endif
  1.2424 +      lp->valid = 1, csa->valid = 0;
  1.2425 +      lp->bfd = csa->bfd, csa->bfd = NULL;
  1.2426 +      memcpy(&lp->head[1], &head[1], m * sizeof(int));
  1.2427 +      /* basic solution status */
  1.2428 +      lp->pbs_stat = p_stat;
  1.2429 +      lp->dbs_stat = d_stat;
  1.2430 +      /* objective function value */
  1.2431 +      lp->obj_val = eval_obj(csa);
  1.2432 +      /* simplex iteration count */
  1.2433 +      lp->it_cnt = csa->it_cnt;
  1.2434 +      /* unbounded ray */
  1.2435 +      lp->some = ray;
  1.2436 +      /* basic variables */
  1.2437 +      for (i = 1; i <= m; i++)
  1.2438 +      {  k = head[i]; /* x[k] = xB[i] */
  1.2439 +#ifdef GLP_DEBUG
  1.2440 +         xassert(1 <= k && k <= m+n);
  1.2441 +#endif
  1.2442 +         if (k <= m)
  1.2443 +         {  GLPROW *row = lp->row[k];
  1.2444 +            row->stat = GLP_BS;
  1.2445 +            row->bind = i;
  1.2446 +            row->prim = bbar[i] / row->rii;
  1.2447 +            row->dual = 0.0;
  1.2448 +         }
  1.2449 +         else
  1.2450 +         {  GLPCOL *col = lp->col[k-m];
  1.2451 +            col->stat = GLP_BS;
  1.2452 +            col->bind = i;
  1.2453 +            col->prim = bbar[i] * col->sjj;
  1.2454 +            col->dual = 0.0;
  1.2455 +         }
  1.2456 +      }
  1.2457 +      /* non-basic variables */
  1.2458 +      for (j = 1; j <= n; j++)
  1.2459 +      {  k = head[m+j]; /* x[k] = xN[j] */
  1.2460 +#ifdef GLP_DEBUG
  1.2461 +         xassert(1 <= k && k <= m+n);
  1.2462 +#endif
  1.2463 +         if (k <= m)
  1.2464 +         {  GLPROW *row = lp->row[k];
  1.2465 +            row->stat = stat[j];
  1.2466 +            row->bind = 0;
  1.2467 +#if 0
  1.2468 +            row->prim = get_xN(csa, j) / row->rii;
  1.2469 +#else
  1.2470 +            switch (stat[j])
  1.2471 +            {  case GLP_NL:
  1.2472 +                  row->prim = row->lb; break;
  1.2473 +               case GLP_NU:
  1.2474 +                  row->prim = row->ub; break;
  1.2475 +               case GLP_NF:
  1.2476 +                  row->prim = 0.0; break;
  1.2477 +               case GLP_NS:
  1.2478 +                  row->prim = row->lb; break;
  1.2479 +               default:
  1.2480 +                  xassert(stat != stat);
  1.2481 +            }
  1.2482 +#endif
  1.2483 +            row->dual = (cbar[j] * row->rii) / zeta;
  1.2484 +         }
  1.2485 +         else
  1.2486 +         {  GLPCOL *col = lp->col[k-m];
  1.2487 +            col->stat = stat[j];
  1.2488 +            col->bind = 0;
  1.2489 +#if 0
  1.2490 +            col->prim = get_xN(csa, j) * col->sjj;
  1.2491 +#else
  1.2492 +            switch (stat[j])
  1.2493 +            {  case GLP_NL:
  1.2494 +                  col->prim = col->lb; break;
  1.2495 +               case GLP_NU:
  1.2496 +                  col->prim = col->ub; break;
  1.2497 +               case GLP_NF:
  1.2498 +                  col->prim = 0.0; break;
  1.2499 +               case GLP_NS:
  1.2500 +                  col->prim = col->lb; break;
  1.2501 +               default:
  1.2502 +                  xassert(stat != stat);
  1.2503 +            }
  1.2504 +#endif
  1.2505 +            col->dual = (cbar[j] / col->sjj) / zeta;
  1.2506 +         }
  1.2507 +      }
  1.2508 +      return;
  1.2509 +}
  1.2510 +
  1.2511 +/***********************************************************************
  1.2512 +*  free_csa - deallocate common storage area
  1.2513 +*
  1.2514 +*  This routine frees all the memory allocated to arrays in the common
  1.2515 +*  storage area (CSA). */
  1.2516 +
  1.2517 +static void free_csa(struct csa *csa)
  1.2518 +{     xfree(csa->type);
  1.2519 +      xfree(csa->lb);
  1.2520 +      xfree(csa->ub);
  1.2521 +      xfree(csa->coef);
  1.2522 +      xfree(csa->obj);
  1.2523 +      xfree(csa->A_ptr);
  1.2524 +      xfree(csa->A_ind);
  1.2525 +      xfree(csa->A_val);
  1.2526 +      xfree(csa->head);
  1.2527 +      xfree(csa->stat);
  1.2528 +      xfree(csa->N_ptr);
  1.2529 +      xfree(csa->N_len);
  1.2530 +      xfree(csa->N_ind);
  1.2531 +      xfree(csa->N_val);
  1.2532 +      xfree(csa->bbar);
  1.2533 +      xfree(csa->cbar);
  1.2534 +      xfree(csa->refsp);
  1.2535 +      xfree(csa->gamma);
  1.2536 +      xfree(csa->tcol_ind);
  1.2537 +      xfree(csa->tcol_vec);
  1.2538 +      xfree(csa->trow_ind);
  1.2539 +      xfree(csa->trow_vec);
  1.2540 +      xfree(csa->work1);
  1.2541 +      xfree(csa->work2);
  1.2542 +      xfree(csa->work3);
  1.2543 +      xfree(csa->work4);
  1.2544 +      xfree(csa);
  1.2545 +      return;
  1.2546 +}
  1.2547 +
  1.2548 +/***********************************************************************
  1.2549 +*  spx_primal - core LP solver based on the primal simplex method
  1.2550 +*
  1.2551 +*  SYNOPSIS
  1.2552 +*
  1.2553 +*  #include "glpspx.h"
  1.2554 +*  int spx_primal(glp_prob *lp, const glp_smcp *parm);
  1.2555 +*
  1.2556 +*  DESCRIPTION
  1.2557 +*
  1.2558 +*  The routine spx_primal is a core LP solver based on the two-phase
  1.2559 +*  primal simplex method.
  1.2560 +*
  1.2561 +*  RETURNS
  1.2562 +*
  1.2563 +*  0  LP instance has been successfully solved.
  1.2564 +*
  1.2565 +*  GLP_EITLIM
  1.2566 +*     Iteration limit has been exhausted.
  1.2567 +*
  1.2568 +*  GLP_ETMLIM
  1.2569 +*     Time limit has been exhausted.
  1.2570 +*
  1.2571 +*  GLP_EFAIL
  1.2572 +*     The solver failed to solve LP instance. */
  1.2573 +
  1.2574 +int spx_primal(glp_prob *lp, const glp_smcp *parm)
  1.2575 +{     struct csa *csa;
  1.2576 +      int binv_st = 2;
  1.2577 +      /* status of basis matrix factorization:
  1.2578 +         0 - invalid; 1 - just computed; 2 - updated */
  1.2579 +      int bbar_st = 0;
  1.2580 +      /* status of primal values of basic variables:
  1.2581 +         0 - invalid; 1 - just computed; 2 - updated */
  1.2582 +      int cbar_st = 0;
  1.2583 +      /* status of reduced costs of non-basic variables:
  1.2584 +         0 - invalid; 1 - just computed; 2 - updated */
  1.2585 +      int rigorous = 0;
  1.2586 +      /* rigorous mode flag; this flag is used to enable iterative
  1.2587 +         refinement on computing pivot rows and columns of the simplex
  1.2588 +         table */
  1.2589 +      int check = 0;
  1.2590 +      int p_stat, d_stat, ret;
  1.2591 +      /* allocate and initialize the common storage area */
  1.2592 +      csa = alloc_csa(lp);
  1.2593 +      init_csa(csa, lp);
  1.2594 +      if (parm->msg_lev >= GLP_MSG_DBG)
  1.2595 +         xprintf("Objective scale factor = %g\n", csa->zeta);
  1.2596 +loop: /* main loop starts here */
  1.2597 +      /* compute factorization of the basis matrix */
  1.2598 +      if (binv_st == 0)
  1.2599 +      {  ret = invert_B(csa);
  1.2600 +         if (ret != 0)
  1.2601 +         {  if (parm->msg_lev >= GLP_MSG_ERR)
  1.2602 +            {  xprintf("Error: unable to factorize the basis matrix (%d"
  1.2603 +                  ")\n", ret);
  1.2604 +               xprintf("Sorry, basis recovery procedure not implemented"
  1.2605 +                  " yet\n");
  1.2606 +            }
  1.2607 +            xassert(!lp->valid && lp->bfd == NULL);
  1.2608 +            lp->bfd = csa->bfd, csa->bfd = NULL;
  1.2609 +            lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
  1.2610 +            lp->obj_val = 0.0;
  1.2611 +            lp->it_cnt = csa->it_cnt;
  1.2612 +            lp->some = 0;
  1.2613 +            ret = GLP_EFAIL;
  1.2614 +            goto done;
  1.2615 +         }
  1.2616 +         csa->valid = 1;
  1.2617 +         binv_st = 1; /* just computed */
  1.2618 +         /* invalidate basic solution components */
  1.2619 +         bbar_st = cbar_st = 0;
  1.2620 +      }
  1.2621 +      /* compute primal values of basic variables */
  1.2622 +      if (bbar_st == 0)
  1.2623 +      {  eval_bbar(csa);
  1.2624 +         bbar_st = 1; /* just computed */
  1.2625 +         /* determine the search phase, if not determined yet */
  1.2626 +         if (csa->phase == 0)
  1.2627 +         {  if (set_aux_obj(csa, parm->tol_bnd) > 0)
  1.2628 +            {  /* current basic solution is primal infeasible */
  1.2629 +               /* start to minimize the sum of infeasibilities */
  1.2630 +               csa->phase = 1;
  1.2631 +            }
  1.2632 +            else
  1.2633 +            {  /* current basic solution is primal feasible */
  1.2634 +               /* start to minimize the original objective function */
  1.2635 +               set_orig_obj(csa);
  1.2636 +               csa->phase = 2;
  1.2637 +            }
  1.2638 +            xassert(check_stab(csa, parm->tol_bnd) == 0);
  1.2639 +            /* working objective coefficients have been changed, so
  1.2640 +               invalidate reduced costs */
  1.2641 +            cbar_st = 0;
  1.2642 +            display(csa, parm, 1);
  1.2643 +         }
  1.2644 +         /* make sure that the current basic solution remains primal
  1.2645 +            feasible (or pseudo feasible on phase I) */
  1.2646 +         if (check_stab(csa, parm->tol_bnd))
  1.2647 +         {  /* there are excessive bound violations due to round-off
  1.2648 +               errors */
  1.2649 +            if (parm->msg_lev >= GLP_MSG_ERR)
  1.2650 +               xprintf("Warning: numerical instability (primal simplex,"
  1.2651 +                  " phase %s)\n", csa->phase == 1 ? "I" : "II");
  1.2652 +            /* restart the search */
  1.2653 +            csa->phase = 0;
  1.2654 +            binv_st = 0;
  1.2655 +            rigorous = 5;
  1.2656 +            goto loop;
  1.2657 +         }
  1.2658 +      }
  1.2659 +      xassert(csa->phase == 1 || csa->phase == 2);
  1.2660 +      /* on phase I we do not need to wait until the current basic
  1.2661 +         solution becomes dual feasible; it is sufficient to make sure
  1.2662 +         that no basic variable violates its bounds */
  1.2663 +      if (csa->phase == 1 && !check_feas(csa, parm->tol_bnd))
  1.2664 +      {  /* the current basis is primal feasible; switch to phase II */
  1.2665 +         csa->phase = 2;
  1.2666 +         set_orig_obj(csa);
  1.2667 +         cbar_st = 0;
  1.2668 +         display(csa, parm, 1);
  1.2669 +      }
  1.2670 +      /* compute reduced costs of non-basic variables */
  1.2671 +      if (cbar_st == 0)
  1.2672 +      {  eval_cbar(csa);
  1.2673 +         cbar_st = 1; /* just computed */
  1.2674 +      }
  1.2675 +      /* redefine the reference space, if required */
  1.2676 +      switch (parm->pricing)
  1.2677 +      {  case GLP_PT_STD:
  1.2678 +            break;
  1.2679 +         case GLP_PT_PSE:
  1.2680 +            if (csa->refct == 0) reset_refsp(csa);
  1.2681 +            break;
  1.2682 +         default:
  1.2683 +            xassert(parm != parm);
  1.2684 +      }
  1.2685 +      /* at this point the basis factorization and all basic solution
  1.2686 +         components are valid */
  1.2687 +      xassert(binv_st && bbar_st && cbar_st);
  1.2688 +      /* check accuracy of current basic solution components (only for
  1.2689 +         debugging) */
  1.2690 +      if (check)
  1.2691 +      {  double e_bbar = err_in_bbar(csa);
  1.2692 +         double e_cbar = err_in_cbar(csa);
  1.2693 +         double e_gamma =
  1.2694 +            (parm->pricing == GLP_PT_PSE ? err_in_gamma(csa) : 0.0);
  1.2695 +         xprintf("e_bbar = %10.3e; e_cbar = %10.3e; e_gamma = %10.3e\n",
  1.2696 +            e_bbar, e_cbar, e_gamma);
  1.2697 +         xassert(e_bbar <= 1e-5 && e_cbar <= 1e-5 && e_gamma <= 1e-3);
  1.2698 +      }
  1.2699 +      /* check if the iteration limit has been exhausted */
  1.2700 +      if (parm->it_lim < INT_MAX &&
  1.2701 +          csa->it_cnt - csa->it_beg >= parm->it_lim)
  1.2702 +      {  if (bbar_st != 1 || csa->phase == 2 && cbar_st != 1)
  1.2703 +         {  if (bbar_st != 1) bbar_st = 0;
  1.2704 +            if (csa->phase == 2 && cbar_st != 1) cbar_st = 0;
  1.2705 +            goto loop;
  1.2706 +         }
  1.2707 +         display(csa, parm, 1);
  1.2708 +         if (parm->msg_lev >= GLP_MSG_ALL)
  1.2709 +            xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n");
  1.2710 +         switch (csa->phase)
  1.2711 +         {  case 1:
  1.2712 +               p_stat = GLP_INFEAS;
  1.2713 +               set_orig_obj(csa);
  1.2714 +               eval_cbar(csa);
  1.2715 +               break;
  1.2716 +            case 2:
  1.2717 +               p_stat = GLP_FEAS;
  1.2718 +               break;
  1.2719 +            default:
  1.2720 +               xassert(csa != csa);
  1.2721 +         }
  1.2722 +         chuzc(csa, parm->tol_dj);
  1.2723 +         d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS);
  1.2724 +         store_sol(csa, lp, p_stat, d_stat, 0);
  1.2725 +         ret = GLP_EITLIM;
  1.2726 +         goto done;
  1.2727 +      }
  1.2728 +      /* check if the time limit has been exhausted */
  1.2729 +      if (parm->tm_lim < INT_MAX &&
  1.2730 +          1000.0 * xdifftime(xtime(), csa->tm_beg) >= parm->tm_lim)
  1.2731 +      {  if (bbar_st != 1 || csa->phase == 2 && cbar_st != 1)
  1.2732 +         {  if (bbar_st != 1) bbar_st = 0;
  1.2733 +            if (csa->phase == 2 && cbar_st != 1) cbar_st = 0;
  1.2734 +            goto loop;
  1.2735 +         }
  1.2736 +         display(csa, parm, 1);
  1.2737 +         if (parm->msg_lev >= GLP_MSG_ALL)
  1.2738 +            xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
  1.2739 +         switch (csa->phase)
  1.2740 +         {  case 1:
  1.2741 +               p_stat = GLP_INFEAS;
  1.2742 +               set_orig_obj(csa);
  1.2743 +               eval_cbar(csa);
  1.2744 +               break;
  1.2745 +            case 2:
  1.2746 +               p_stat = GLP_FEAS;
  1.2747 +               break;
  1.2748 +            default:
  1.2749 +               xassert(csa != csa);
  1.2750 +         }
  1.2751 +         chuzc(csa, parm->tol_dj);
  1.2752 +         d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS);
  1.2753 +         store_sol(csa, lp, p_stat, d_stat, 0);
  1.2754 +         ret = GLP_ETMLIM;
  1.2755 +         goto done;
  1.2756 +      }
  1.2757 +      /* display the search progress */
  1.2758 +      display(csa, parm, 0);
  1.2759 +      /* choose non-basic variable xN[q] */
  1.2760 +      chuzc(csa, parm->tol_dj);
  1.2761 +      if (csa->q == 0)
  1.2762 +      {  if (bbar_st != 1 || cbar_st != 1)
  1.2763 +         {  if (bbar_st != 1) bbar_st = 0;
  1.2764 +            if (cbar_st != 1) cbar_st = 0;
  1.2765 +            goto loop;
  1.2766 +         }
  1.2767 +         display(csa, parm, 1);
  1.2768 +         switch (csa->phase)
  1.2769 +         {  case 1:
  1.2770 +               if (parm->msg_lev >= GLP_MSG_ALL)
  1.2771 +                  xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
  1.2772 +               p_stat = GLP_NOFEAS;
  1.2773 +               set_orig_obj(csa);
  1.2774 +               eval_cbar(csa);
  1.2775 +               chuzc(csa, parm->tol_dj);
  1.2776 +               d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS);
  1.2777 +               break;
  1.2778 +            case 2:
  1.2779 +               if (parm->msg_lev >= GLP_MSG_ALL)
  1.2780 +                  xprintf("OPTIMAL SOLUTION FOUND\n");
  1.2781 +               p_stat = d_stat = GLP_FEAS;
  1.2782 +               break;
  1.2783 +            default:
  1.2784 +               xassert(csa != csa);
  1.2785 +         }
  1.2786 +         store_sol(csa, lp, p_stat, d_stat, 0);
  1.2787 +         ret = 0;
  1.2788 +         goto done;
  1.2789 +      }
  1.2790 +      /* compute pivot column of the simplex table */
  1.2791 +      eval_tcol(csa);
  1.2792 +      if (rigorous) refine_tcol(csa);
  1.2793 +      sort_tcol(csa, parm->tol_piv);
  1.2794 +      /* check accuracy of the reduced cost of xN[q] */
  1.2795 +      {  double d1 = csa->cbar[csa->q]; /* less accurate */
  1.2796 +         double d2 = reeval_cost(csa);  /* more accurate */
  1.2797 +         xassert(d1 != 0.0);
  1.2798 +         if (fabs(d1 - d2) > 1e-5 * (1.0 + fabs(d2)) ||
  1.2799 +             !(d1 < 0.0 && d2 < 0.0 || d1 > 0.0 && d2 > 0.0))
  1.2800 +         {  if (parm->msg_lev >= GLP_MSG_DBG)
  1.2801 +               xprintf("d1 = %.12g; d2 = %.12g\n", d1, d2);
  1.2802 +            if (cbar_st != 1 || !rigorous)
  1.2803 +            {  if (cbar_st != 1) cbar_st = 0;
  1.2804 +               rigorous = 5;
  1.2805 +               goto loop;
  1.2806 +            }
  1.2807 +         }
  1.2808 +         /* replace cbar[q] by more accurate value keeping its sign */
  1.2809 +         if (d1 > 0.0)
  1.2810 +            csa->cbar[csa->q] = (d2 > 0.0 ? d2 : +DBL_EPSILON);
  1.2811 +         else
  1.2812 +            csa->cbar[csa->q] = (d2 < 0.0 ? d2 : -DBL_EPSILON);
  1.2813 +      }
  1.2814 +      /* choose basic variable xB[p] */
  1.2815 +      switch (parm->r_test)
  1.2816 +      {  case GLP_RT_STD:
  1.2817 +            chuzr(csa, 0.0);
  1.2818 +            break;
  1.2819 +         case GLP_RT_HAR:
  1.2820 +            chuzr(csa, 0.30 * parm->tol_bnd);
  1.2821 +            break;
  1.2822 +         default:
  1.2823 +            xassert(parm != parm);
  1.2824 +      }
  1.2825 +      if (csa->p == 0)
  1.2826 +      {  if (bbar_st != 1 || cbar_st != 1 || !rigorous)
  1.2827 +         {  if (bbar_st != 1) bbar_st = 0;
  1.2828 +            if (cbar_st != 1) cbar_st = 0;
  1.2829 +            rigorous = 1;
  1.2830 +            goto loop;
  1.2831 +         }
  1.2832 +         display(csa, parm, 1);
  1.2833 +         switch (csa->phase)
  1.2834 +         {  case 1:
  1.2835 +               if (parm->msg_lev >= GLP_MSG_ERR)
  1.2836 +                  xprintf("Error: unable to choose basic variable on ph"
  1.2837 +                     "ase I\n");
  1.2838 +               xassert(!lp->valid && lp->bfd == NULL);
  1.2839 +               lp->bfd = csa->bfd, csa->bfd = NULL;
  1.2840 +               lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
  1.2841 +               lp->obj_val = 0.0;
  1.2842 +               lp->it_cnt = csa->it_cnt;
  1.2843 +               lp->some = 0;
  1.2844 +               ret = GLP_EFAIL;
  1.2845 +               break;
  1.2846 +            case 2:
  1.2847 +               if (parm->msg_lev >= GLP_MSG_ALL)
  1.2848 +                  xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
  1.2849 +               store_sol(csa, lp, GLP_FEAS, GLP_NOFEAS,
  1.2850 +                  csa->head[csa->m+csa->q]);
  1.2851 +               ret = 0;
  1.2852 +               break;
  1.2853 +            default:
  1.2854 +               xassert(csa != csa);
  1.2855 +         }
  1.2856 +         goto done;
  1.2857 +      }
  1.2858 +      /* check if the pivot element is acceptable */
  1.2859 +      if (csa->p > 0)
  1.2860 +      {  double piv = csa->tcol_vec[csa->p];
  1.2861 +         double eps = 1e-5 * (1.0 + 0.01 * csa->tcol_max);
  1.2862 +         if (fabs(piv) < eps)
  1.2863 +         {  if (parm->msg_lev >= GLP_MSG_DBG)
  1.2864 +               xprintf("piv = %.12g; eps = %g\n", piv, eps);
  1.2865 +            if (!rigorous)
  1.2866 +            {  rigorous = 5;
  1.2867 +               goto loop;
  1.2868 +            }
  1.2869 +         }
  1.2870 +      }
  1.2871 +      /* now xN[q] and xB[p] have been chosen anyhow */
  1.2872 +      /* compute pivot row of the simplex table */
  1.2873 +      if (csa->p > 0)
  1.2874 +      {  double *rho = csa->work4;
  1.2875 +         eval_rho(csa, rho);
  1.2876 +         if (rigorous) refine_rho(csa, rho);
  1.2877 +         eval_trow(csa, rho);
  1.2878 +      }
  1.2879 +      /* accuracy check based on the pivot element */
  1.2880 +      if (csa->p > 0)
  1.2881 +      {  double piv1 = csa->tcol_vec[csa->p]; /* more accurate */
  1.2882 +         double piv2 = csa->trow_vec[csa->q]; /* less accurate */
  1.2883 +         xassert(piv1 != 0.0);
  1.2884 +         if (fabs(piv1 - piv2) > 1e-8 * (1.0 + fabs(piv1)) ||
  1.2885 +             !(piv1 > 0.0 && piv2 > 0.0 || piv1 < 0.0 && piv2 < 0.0))
  1.2886 +         {  if (parm->msg_lev >= GLP_MSG_DBG)
  1.2887 +               xprintf("piv1 = %.12g; piv2 = %.12g\n", piv1, piv2);
  1.2888 +            if (binv_st != 1 || !rigorous)
  1.2889 +            {  if (binv_st != 1) binv_st = 0;
  1.2890 +               rigorous = 5;
  1.2891 +               goto loop;
  1.2892 +            }
  1.2893 +            /* use more accurate version in the pivot row */
  1.2894 +            if (csa->trow_vec[csa->q] == 0.0)
  1.2895 +            {  csa->trow_nnz++;
  1.2896 +               xassert(csa->trow_nnz <= csa->n);
  1.2897 +               csa->trow_ind[csa->trow_nnz] = csa->q;
  1.2898 +            }
  1.2899 +            csa->trow_vec[csa->q] = piv1;
  1.2900 +         }
  1.2901 +      }
  1.2902 +      /* update primal values of basic variables */
  1.2903 +      update_bbar(csa);
  1.2904 +      bbar_st = 2; /* updated */
  1.2905 +      /* update reduced costs of non-basic variables */
  1.2906 +      if (csa->p > 0)
  1.2907 +      {  update_cbar(csa);
  1.2908 +         cbar_st = 2; /* updated */
  1.2909 +         /* on phase I objective coefficient of xB[p] in the adjacent
  1.2910 +            basis becomes zero */
  1.2911 +         if (csa->phase == 1)
  1.2912 +         {  int k = csa->head[csa->p]; /* x[k] = xB[p] -> xN[q] */
  1.2913 +            csa->cbar[csa->q] -= csa->coef[k];
  1.2914 +            csa->coef[k] = 0.0;
  1.2915 +         }
  1.2916 +      }
  1.2917 +      /* update steepest edge coefficients */
  1.2918 +      if (csa->p > 0)
  1.2919 +      {  switch (parm->pricing)
  1.2920 +         {  case GLP_PT_STD:
  1.2921 +               break;
  1.2922 +            case GLP_PT_PSE:
  1.2923 +               if (csa->refct > 0) update_gamma(csa);
  1.2924 +               break;
  1.2925 +            default:
  1.2926 +               xassert(parm != parm);
  1.2927 +         }
  1.2928 +      }
  1.2929 +      /* update factorization of the basis matrix */
  1.2930 +      if (csa->p > 0)
  1.2931 +      {  ret = update_B(csa, csa->p, csa->head[csa->m+csa->q]);
  1.2932 +         if (ret == 0)
  1.2933 +            binv_st = 2; /* updated */
  1.2934 +         else
  1.2935 +         {  csa->valid = 0;
  1.2936 +            binv_st = 0; /* invalid */
  1.2937 +         }
  1.2938 +      }
  1.2939 +      /* update matrix N */
  1.2940 +      if (csa->p > 0)
  1.2941 +      {  del_N_col(csa, csa->q, csa->head[csa->m+csa->q]);
  1.2942 +         if (csa->type[csa->head[csa->p]] != GLP_FX)
  1.2943 +            add_N_col(csa, csa->q, csa->head[csa->p]);
  1.2944 +      }
  1.2945 +      /* change the basis header */
  1.2946 +      change_basis(csa);
  1.2947 +      /* iteration complete */
  1.2948 +      csa->it_cnt++;
  1.2949 +      if (rigorous > 0) rigorous--;
  1.2950 +      goto loop;
  1.2951 +done: /* deallocate the common storage area */
  1.2952 +      free_csa(csa);
  1.2953 +      /* return to the calling program */
  1.2954 +      return ret;
  1.2955 +}
  1.2956 +
  1.2957 +/* eof */