1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpspx01.c Mon Dec 06 13:09:21 2010 +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 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 */