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