lemon-project-template-glpk: deps/glpk/src/glpspx02.c annotate

lemon-project-template-glpk

annotate deps/glpk/src/glpspx02.c @ 11:4fc6ad2fb8a6

Test GLPK in src/main.cc

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 21:43:29 +0100
parents
children

rev	line source
alpar@9	1 /* glpspx02.c (dual simplex method) */
alpar@9	2
alpar@9	3 /***********************************************************************
alpar@9	4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9	5 *
alpar@9	6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9	7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9	8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9	9 * E-mail: <mao@gnu.org>.
alpar@9	10 *
alpar@9	11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9	12 * under the terms of the GNU General Public License as published by
alpar@9	13 * the Free Software Foundation, either version 3 of the License, or
alpar@9	14 * (at your option) any later version.
alpar@9	15 *
alpar@9	16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9	17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9	18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9	19 * License for more details.
alpar@9	20 *
alpar@9	21 * You should have received a copy of the GNU General Public License
alpar@9	22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9	23 ***********************************************************************/
alpar@9	24
alpar@9	25 #include "glpspx.h"
alpar@9	26
alpar@9	27 #define GLP_DEBUG 1
alpar@9	28
alpar@9	29 #if 0
alpar@9	30 #define GLP_LONG_STEP 1
alpar@9	31 #endif
alpar@9	32
alpar@9	33 struct csa
alpar@9	34 { /* common storage area */
alpar@9	35 /--------------------------------------------------------------/
alpar@9	36 /* LP data */
alpar@9	37 int m;
alpar@9	38 /* number of rows (auxiliary variables), m > 0 */
alpar@9	39 int n;
alpar@9	40 /* number of columns (structural variables), n > 0 */
alpar@9	41 char type; / char type[1+m+n]; */
alpar@9	42 /* type[0] is not used;
alpar@9	43 type[k], 1 <= k <= m+n, is the type of variable x[k]:
alpar@9	44 GLP_FR - free variable
alpar@9	45 GLP_LO - variable with lower bound
alpar@9	46 GLP_UP - variable with upper bound
alpar@9	47 GLP_DB - double-bounded variable
alpar@9	48 GLP_FX - fixed variable */
alpar@9	49 double lb; / double lb[1+m+n]; */
alpar@9	50 /* lb[0] is not used;
alpar@9	51 lb[k], 1 <= k <= m+n, is an lower bound of variable x[k];
alpar@9	52 if x[k] has no lower bound, lb[k] is zero */
alpar@9	53 double ub; / double ub[1+m+n]; */
alpar@9	54 /* ub[0] is not used;
alpar@9	55 ub[k], 1 <= k <= m+n, is an upper bound of variable x[k];
alpar@9	56 if x[k] has no upper bound, ub[k] is zero;
alpar@9	57 if x[k] is of fixed type, ub[k] is the same as lb[k] */
alpar@9	58 double coef; / double coef[1+m+n]; */
alpar@9	59 /* coef[0] is not used;
alpar@9	60 coef[k], 1 <= k <= m+n, is an objective coefficient at
alpar@9	61 variable x[k] */
alpar@9	62 /--------------------------------------------------------------/
alpar@9	63 /* original bounds of variables */
alpar@9	64 char orig_type; / char orig_type[1+m+n]; */
alpar@9	65 double orig_lb; / double orig_lb[1+m+n]; */
alpar@9	66 double orig_ub; / double orig_ub[1+m+n]; */
alpar@9	67 /--------------------------------------------------------------/
alpar@9	68 /* original objective function */
alpar@9	69 double obj; / double obj[1+n]; */
alpar@9	70 /* obj[0] is a constant term of the original objective function;
alpar@9	71 obj[j], 1 <= j <= n, is an original objective coefficient at
alpar@9	72 structural variable x[m+j] */
alpar@9	73 double zeta;
alpar@9	74 /* factor used to scale original objective coefficients; its
alpar@9	75 sign defines original optimization direction: zeta > 0 means
alpar@9	76 minimization, zeta < 0 means maximization */
alpar@9	77 /--------------------------------------------------------------/
alpar@9	78 /* constraint matrix A; it has m rows and n columns and is stored
alpar@9	79 by columns */
alpar@9	80 int A_ptr; / int A_ptr[1+n+1]; */
alpar@9	81 /* A_ptr[0] is not used;
alpar@9	82 A_ptr[j], 1 <= j <= n, is starting position of j-th column in
alpar@9	83 arrays A_ind and A_val; note that A_ptr[1] is always 1;
alpar@9	84 A_ptr[n+1] indicates the position after the last element in
alpar@9	85 arrays A_ind and A_val */
alpar@9	86 int A_ind; / int A_ind[A_ptr[n+1]]; */
alpar@9	87 /* row indices */
alpar@9	88 double A_val; / double A_val[A_ptr[n+1]]; */
alpar@9	89 /* non-zero element values */
alpar@9	90 #if 1 /* 06/IV-2009 */
alpar@9	91 /* constraint matrix A stored by rows */
alpar@9	92 int AT_ptr; / int AT_ptr[1+m+1]; */
alpar@9	93 /* AT_ptr[0] is not used;
alpar@9	94 AT_ptr[i], 1 <= i <= m, is starting position of i-th row in
alpar@9	95 arrays AT_ind and AT_val; note that AT_ptr[1] is always 1;
alpar@9	96 AT_ptr[m+1] indicates the position after the last element in
alpar@9	97 arrays AT_ind and AT_val */
alpar@9	98 int AT_ind; / int AT_ind[AT_ptr[m+1]]; */
alpar@9	99 /* column indices */
alpar@9	100 double AT_val; / double AT_val[AT_ptr[m+1]]; */
alpar@9	101 /* non-zero element values */
alpar@9	102 #endif
alpar@9	103 /--------------------------------------------------------------/
alpar@9	104 /* basis header */
alpar@9	105 int head; / int head[1+m+n]; */
alpar@9	106 /* head[0] is not used;
alpar@9	107 head[i], 1 <= i <= m, is the ordinal number of basic variable
alpar@9	108 xB[i]; head[i] = k means that xB[i] = x[k] and i-th column of
alpar@9	109 matrix B is k-th column of matrix (I\|-A);
alpar@9	110 head[m+j], 1 <= j <= n, is the ordinal number of non-basic
alpar@9	111 variable xN[j]; head[m+j] = k means that xN[j] = x[k] and j-th
alpar@9	112 column of matrix N is k-th column of matrix (I\|-A) */
alpar@9	113 #if 1 /* 06/IV-2009 */
alpar@9	114 int bind; / int bind[1+m+n]; */
alpar@9	115 /* bind[0] is not used;
alpar@9	116 bind[k], 1 <= k <= m+n, is the position of k-th column of the
alpar@9	117 matrix (I\|-A) in the matrix (B\|N); that is, bind[k] = k' means
alpar@9	118 that head[k'] = k */
alpar@9	119 #endif
alpar@9	120 char stat; / char stat[1+n]; */
alpar@9	121 /* stat[0] is not used;
alpar@9	122 stat[j], 1 <= j <= n, is the status of non-basic variable
alpar@9	123 xN[j], which defines its active bound:
alpar@9	124 GLP_NL - lower bound is active
alpar@9	125 GLP_NU - upper bound is active
alpar@9	126 GLP_NF - free variable
alpar@9	127 GLP_NS - fixed variable */
alpar@9	128 /--------------------------------------------------------------/
alpar@9	129 /* matrix B is the basis matrix; it is composed from columns of
alpar@9	130 the augmented constraint matrix (I\|-A) corresponding to basic
alpar@9	131 variables and stored in a factorized (invertable) form */
alpar@9	132 int valid;
alpar@9	133 /* factorization is valid only if this flag is set */
alpar@9	134 BFD bfd; / BFD bfd[1:m,1:m]; */
alpar@9	135 /* factorized (invertable) form of the basis matrix */
alpar@9	136 #if 0 /* 06/IV-2009 */
alpar@9	137 /--------------------------------------------------------------/
alpar@9	138 /* matrix N is a matrix composed from columns of the augmented
alpar@9	139 constraint matrix (I\|-A) corresponding to non-basic variables
alpar@9	140 except fixed ones; it is stored by rows and changes every time
alpar@9	141 the basis changes */
alpar@9	142 int N_ptr; / int N_ptr[1+m+1]; */
alpar@9	143 /* N_ptr[0] is not used;
alpar@9	144 N_ptr[i], 1 <= i <= m, is starting position of i-th row in
alpar@9	145 arrays N_ind and N_val; note that N_ptr[1] is always 1;
alpar@9	146 N_ptr[m+1] indicates the position after the last element in
alpar@9	147 arrays N_ind and N_val */
alpar@9	148 int N_len; / int N_len[1+m]; */
alpar@9	149 /* N_len[0] is not used;
alpar@9	150 N_len[i], 1 <= i <= m, is length of i-th row (0 to n) */
alpar@9	151 int N_ind; / int N_ind[N_ptr[m+1]]; */
alpar@9	152 /* column indices */
alpar@9	153 double N_val; / double N_val[N_ptr[m+1]]; */
alpar@9	154 /* non-zero element values */
alpar@9	155 #endif
alpar@9	156 /--------------------------------------------------------------/
alpar@9	157 /* working parameters */
alpar@9	158 int phase;
alpar@9	159 /* search phase:
alpar@9	160 0 - not determined yet
alpar@9	161 1 - search for dual feasible solution
alpar@9	162 2 - search for optimal solution */
alpar@9	163 glp_long tm_beg;
alpar@9	164 /* time value at the beginning of the search */
alpar@9	165 int it_beg;
alpar@9	166 /* simplex iteration count at the beginning of the search */
alpar@9	167 int it_cnt;
alpar@9	168 /* simplex iteration count; it increases by one every time the
alpar@9	169 basis changes */
alpar@9	170 int it_dpy;
alpar@9	171 /* simplex iteration count at the most recent display output */
alpar@9	172 /--------------------------------------------------------------/
alpar@9	173 /* basic solution components */
alpar@9	174 double bbar; / double bbar[1+m]; */
alpar@9	175 /* bbar[0] is not used on phase I; on phase II it is the current
alpar@9	176 value of the original objective function;
alpar@9	177 bbar[i], 1 <= i <= m, is primal value of basic variable xB[i]
alpar@9	178 (if xB[i] is free, its primal value is not updated) */
alpar@9	179 double cbar; / double cbar[1+n]; */
alpar@9	180 /* cbar[0] is not used;
alpar@9	181 cbar[j], 1 <= j <= n, is reduced cost of non-basic variable
alpar@9	182 xN[j] (if xN[j] is fixed, its reduced cost is not updated) */
alpar@9	183 /--------------------------------------------------------------/
alpar@9	184 /* the following pricing technique options may be used:
alpar@9	185 GLP_PT_STD - standard ("textbook") pricing;
alpar@9	186 GLP_PT_PSE - projected steepest edge;
alpar@9	187 GLP_PT_DVX - Devex pricing (not implemented yet);
alpar@9	188 in case of GLP_PT_STD the reference space is not used, and all
alpar@9	189 steepest edge coefficients are set to 1 */
alpar@9	190 int refct;
alpar@9	191 /* this count is set to an initial value when the reference space
alpar@9	192 is defined and decreases by one every time the basis changes;
alpar@9	193 once this count reaches zero, the reference space is redefined
alpar@9	194 again */
alpar@9	195 char refsp; / char refsp[1+m+n]; */
alpar@9	196 /* refsp[0] is not used;
alpar@9	197 refsp[k], 1 <= k <= m+n, is the flag which means that variable
alpar@9	198 x[k] belongs to the current reference space */
alpar@9	199 double gamma; / double gamma[1+m]; */
alpar@9	200 /* gamma[0] is not used;
alpar@9	201 gamma[i], 1 <= i <= n, is the steepest edge coefficient for
alpar@9	202 basic variable xB[i]; if xB[i] is free, gamma[i] is not used
alpar@9	203 and just set to 1 */
alpar@9	204 /--------------------------------------------------------------/
alpar@9	205 /* basic variable xB[p] chosen to leave the basis */
alpar@9	206 int p;
alpar@9	207 /* index of the basic variable xB[p] chosen, 1 <= p <= m;
alpar@9	208 if the set of eligible basic variables is empty (i.e. if the
alpar@9	209 current basic solution is primal feasible within a tolerance)
alpar@9	210 and thus no variable has been chosen, p is set to 0 */
alpar@9	211 double delta;
alpar@9	212 /* change of xB[p] in the adjacent basis;
alpar@9	213 delta > 0 means that xB[p] violates its lower bound and will
alpar@9	214 increase to achieve it in the adjacent basis;
alpar@9	215 delta < 0 means that xB[p] violates its upper bound and will
alpar@9	216 decrease to achieve it in the adjacent basis */
alpar@9	217 /--------------------------------------------------------------/
alpar@9	218 /* pivot row of the simplex table corresponding to basic variable
alpar@9	219 xB[p] chosen is the following vector:
alpar@9	220 T' * e[p] = - N' * inv(B') * e[p] = - N' * rho,
alpar@9	221 where B' is a matrix transposed to the current basis matrix,
alpar@9	222 N' is a matrix, whose rows are columns of the matrix (I\|-A)
alpar@9	223 corresponding to non-basic non-fixed variables */
alpar@9	224 int trow_nnz;
alpar@9	225 /* number of non-zero components, 0 <= nnz <= n */
alpar@9	226 int trow_ind; / int trow_ind[1+n]; */
alpar@9	227 /* trow_ind[0] is not used;
alpar@9	228 trow_ind[t], 1 <= t <= nnz, is an index of non-zero component,
alpar@9	229 i.e. trow_ind[t] = j means that trow_vec[j] != 0 */
alpar@9	230 double trow_vec; / int trow_vec[1+n]; */
alpar@9	231 /* trow_vec[0] is not used;
alpar@9	232 trow_vec[j], 1 <= j <= n, is a numeric value of j-th component
alpar@9	233 of the row */
alpar@9	234 double trow_max;
alpar@9	235 /* infinity (maximum) norm of the row (max \|trow_vec[j]\|) */
alpar@9	236 int trow_num;
alpar@9	237 /* number of significant non-zero components, which means that:
alpar@9	238 \|trow_vec[j]\| >= eps for j in trow_ind[1,...,num],
alpar@9	239 \|tcol_vec[j]\| < eps for j in trow_ind[num+1,...,nnz],
alpar@9	240 where eps is a pivot tolerance */
alpar@9	241 /--------------------------------------------------------------/
alpar@9	242 #ifdef GLP_LONG_STEP /* 07/IV-2009 */
alpar@9	243 int nbps;
alpar@9	244 /* number of breakpoints, 0 <= nbps <= n */
alpar@9	245 struct bkpt
alpar@9	246 { int j;
alpar@9	247 /* index of non-basic variable xN[j], 1 <= j <= n */
alpar@9	248 double t;
alpar@9	249 /* value of dual ray parameter at breakpoint, t >= 0 */
alpar@9	250 double dz;
alpar@9	251 /* dz = zeta(t = t[k]) - zeta(t = 0) */
alpar@9	252 } bkpt; / struct bkpt bkpt[1+n]; */
alpar@9	253 /* bkpt[0] is not used;
alpar@9	254 bkpt[k], 1 <= k <= nbps, is k-th breakpoint of the dual
alpar@9	255 objective */
alpar@9	256 #endif
alpar@9	257 /--------------------------------------------------------------/
alpar@9	258 /* non-basic variable xN[q] chosen to enter the basis */
alpar@9	259 int q;
alpar@9	260 /* index of the non-basic variable xN[q] chosen, 1 <= q <= n;
alpar@9	261 if no variable has been chosen, q is set to 0 */
alpar@9	262 double new_dq;
alpar@9	263 /* reduced cost of xN[q] in the adjacent basis (it is the change
alpar@9	264 of lambdaB[p]) */
alpar@9	265 /--------------------------------------------------------------/
alpar@9	266 /* pivot column of the simplex table corresponding to non-basic
alpar@9	267 variable xN[q] chosen is the following vector:
alpar@9	268 T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q],
alpar@9	269 where B is the current basis matrix, N[q] is a column of the
alpar@9	270 matrix (I\|-A) corresponding to xN[q] */
alpar@9	271 int tcol_nnz;
alpar@9	272 /* number of non-zero components, 0 <= nnz <= m */
alpar@9	273 int tcol_ind; / int tcol_ind[1+m]; */
alpar@9	274 /* tcol_ind[0] is not used;
alpar@9	275 tcol_ind[t], 1 <= t <= nnz, is an index of non-zero component,
alpar@9	276 i.e. tcol_ind[t] = i means that tcol_vec[i] != 0 */
alpar@9	277 double tcol_vec; / double tcol_vec[1+m]; */
alpar@9	278 /* tcol_vec[0] is not used;
alpar@9	279 tcol_vec[i], 1 <= i <= m, is a numeric value of i-th component
alpar@9	280 of the column */
alpar@9	281 /--------------------------------------------------------------/
alpar@9	282 /* working arrays */
alpar@9	283 double work1; / double work1[1+m]; */
alpar@9	284 double work2; / double work2[1+m]; */
alpar@9	285 double work3; / double work3[1+m]; */
alpar@9	286 double work4; / double work4[1+m]; */
alpar@9	287 };
alpar@9	288
alpar@9	289 static const double kappa = 0.10;
alpar@9	290
alpar@9	291 /***********************************************************************
alpar@9	292 * alloc_csa - allocate common storage area
alpar@9	293 *
alpar@9	294 * This routine allocates all arrays in the common storage area (CSA)
alpar@9	295 * and returns a pointer to the CSA. */
alpar@9	296
alpar@9	297 static struct csa alloc_csa(glp_prob lp)
alpar@9	298 { struct csa *csa;
alpar@9	299 int m = lp->m;
alpar@9	300 int n = lp->n;
alpar@9	301 int nnz = lp->nnz;
alpar@9	302 csa = xmalloc(sizeof(struct csa));
alpar@9	303 xassert(m > 0 && n > 0);
alpar@9	304 csa->m = m;
alpar@9	305 csa->n = n;
alpar@9	306 csa->type = xcalloc(1+m+n, sizeof(char));
alpar@9	307 csa->lb = xcalloc(1+m+n, sizeof(double));
alpar@9	308 csa->ub = xcalloc(1+m+n, sizeof(double));
alpar@9	309 csa->coef = xcalloc(1+m+n, sizeof(double));
alpar@9	310 csa->orig_type = xcalloc(1+m+n, sizeof(char));
alpar@9	311 csa->orig_lb = xcalloc(1+m+n, sizeof(double));
alpar@9	312 csa->orig_ub = xcalloc(1+m+n, sizeof(double));
alpar@9	313 csa->obj = xcalloc(1+n, sizeof(double));
alpar@9	314 csa->A_ptr = xcalloc(1+n+1, sizeof(int));
alpar@9	315 csa->A_ind = xcalloc(1+nnz, sizeof(int));
alpar@9	316 csa->A_val = xcalloc(1+nnz, sizeof(double));
alpar@9	317 #if 1 /* 06/IV-2009 */
alpar@9	318 csa->AT_ptr = xcalloc(1+m+1, sizeof(int));
alpar@9	319 csa->AT_ind = xcalloc(1+nnz, sizeof(int));
alpar@9	320 csa->AT_val = xcalloc(1+nnz, sizeof(double));
alpar@9	321 #endif
alpar@9	322 csa->head = xcalloc(1+m+n, sizeof(int));
alpar@9	323 #if 1 /* 06/IV-2009 */
alpar@9	324 csa->bind = xcalloc(1+m+n, sizeof(int));
alpar@9	325 #endif
alpar@9	326 csa->stat = xcalloc(1+n, sizeof(char));
alpar@9	327 #if 0 /* 06/IV-2009 */
alpar@9	328 csa->N_ptr = xcalloc(1+m+1, sizeof(int));
alpar@9	329 csa->N_len = xcalloc(1+m, sizeof(int));
alpar@9	330 csa->N_ind = NULL; /* will be allocated later */
alpar@9	331 csa->N_val = NULL; /* will be allocated later */
alpar@9	332 #endif
alpar@9	333 csa->bbar = xcalloc(1+m, sizeof(double));
alpar@9	334 csa->cbar = xcalloc(1+n, sizeof(double));
alpar@9	335 csa->refsp = xcalloc(1+m+n, sizeof(char));
alpar@9	336 csa->gamma = xcalloc(1+m, sizeof(double));
alpar@9	337 csa->trow_ind = xcalloc(1+n, sizeof(int));
alpar@9	338 csa->trow_vec = xcalloc(1+n, sizeof(double));
alpar@9	339 #ifdef GLP_LONG_STEP /* 07/IV-2009 */
alpar@9	340 csa->bkpt = xcalloc(1+n, sizeof(struct bkpt));
alpar@9	341 #endif
alpar@9	342 csa->tcol_ind = xcalloc(1+m, sizeof(int));
alpar@9	343 csa->tcol_vec = xcalloc(1+m, sizeof(double));
alpar@9	344 csa->work1 = xcalloc(1+m, sizeof(double));
alpar@9	345 csa->work2 = xcalloc(1+m, sizeof(double));
alpar@9	346 csa->work3 = xcalloc(1+m, sizeof(double));
alpar@9	347 csa->work4 = xcalloc(1+m, sizeof(double));
alpar@9	348 return csa;
alpar@9	349 }
alpar@9	350
alpar@9	351 /***********************************************************************
alpar@9	352 * init_csa - initialize common storage area
alpar@9	353 *
alpar@9	354 * This routine initializes all data structures in the common storage
alpar@9	355 * area (CSA). */
alpar@9	356
alpar@9	357 static void init_csa(struct csa csa, glp_prob lp)
alpar@9	358 { int m = csa->m;
alpar@9	359 int n = csa->n;
alpar@9	360 char *type = csa->type;
alpar@9	361 double *lb = csa->lb;
alpar@9	362 double *ub = csa->ub;
alpar@9	363 double *coef = csa->coef;
alpar@9	364 char *orig_type = csa->orig_type;
alpar@9	365 double *orig_lb = csa->orig_lb;
alpar@9	366 double *orig_ub = csa->orig_ub;
alpar@9	367 double *obj = csa->obj;
alpar@9	368 int *A_ptr = csa->A_ptr;
alpar@9	369 int *A_ind = csa->A_ind;
alpar@9	370 double *A_val = csa->A_val;
alpar@9	371 #if 1 /* 06/IV-2009 */
alpar@9	372 int *AT_ptr = csa->AT_ptr;
alpar@9	373 int *AT_ind = csa->AT_ind;
alpar@9	374 double *AT_val = csa->AT_val;
alpar@9	375 #endif
alpar@9	376 int *head = csa->head;
alpar@9	377 #if 1 /* 06/IV-2009 */
alpar@9	378 int *bind = csa->bind;
alpar@9	379 #endif
alpar@9	380 char *stat = csa->stat;
alpar@9	381 char *refsp = csa->refsp;
alpar@9	382 double *gamma = csa->gamma;
alpar@9	383 int i, j, k, loc;
alpar@9	384 double cmax;
alpar@9	385 /* auxiliary variables */
alpar@9	386 for (i = 1; i <= m; i++)
alpar@9	387 { GLPROW *row = lp->row[i];
alpar@9	388 type[i] = (char)row->type;
alpar@9	389 lb[i] = row->lb * row->rii;
alpar@9	390 ub[i] = row->ub * row->rii;
alpar@9	391 coef[i] = 0.0;
alpar@9	392 }
alpar@9	393 /* structural variables */
alpar@9	394 for (j = 1; j <= n; j++)
alpar@9	395 { GLPCOL *col = lp->col[j];
alpar@9	396 type[m+j] = (char)col->type;
alpar@9	397 lb[m+j] = col->lb / col->sjj;
alpar@9	398 ub[m+j] = col->ub / col->sjj;
alpar@9	399 coef[m+j] = col->coef * col->sjj;
alpar@9	400 }
alpar@9	401 /* original bounds of variables */
alpar@9	402 memcpy(&orig_type[1], &type[1], (m+n) * sizeof(char));
alpar@9	403 memcpy(&orig_lb[1], &lb[1], (m+n) * sizeof(double));
alpar@9	404 memcpy(&orig_ub[1], &ub[1], (m+n) * sizeof(double));
alpar@9	405 /* original objective function */
alpar@9	406 obj[0] = lp->c0;
alpar@9	407 memcpy(&obj[1], &coef[m+1], n * sizeof(double));
alpar@9	408 /* factor used to scale original objective coefficients */
alpar@9	409 cmax = 0.0;
alpar@9	410 for (j = 1; j <= n; j++)
alpar@9	411 if (cmax < fabs(obj[j])) cmax = fabs(obj[j]);
alpar@9	412 if (cmax == 0.0) cmax = 1.0;
alpar@9	413 switch (lp->dir)
alpar@9	414 { case GLP_MIN:
alpar@9	415 csa->zeta = + 1.0 / cmax;
alpar@9	416 break;
alpar@9	417 case GLP_MAX:
alpar@9	418 csa->zeta = - 1.0 / cmax;
alpar@9	419 break;
alpar@9	420 default:
alpar@9	421 xassert(lp != lp);
alpar@9	422 }
alpar@9	423 #if 1
alpar@9	424 if (fabs(csa->zeta) < 1.0) csa->zeta *= 1000.0;
alpar@9	425 #endif
alpar@9	426 /* scale working objective coefficients */
alpar@9	427 for (j = 1; j <= n; j++) coef[m+j] *= csa->zeta;
alpar@9	428 /* matrix A (by columns) */
alpar@9	429 loc = 1;
alpar@9	430 for (j = 1; j <= n; j++)
alpar@9	431 { GLPAIJ *aij;
alpar@9	432 A_ptr[j] = loc;
alpar@9	433 for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
alpar@9	434 { A_ind[loc] = aij->row->i;
alpar@9	435 A_val[loc] = aij->row->rii * aij->val * aij->col->sjj;
alpar@9	436 loc++;
alpar@9	437 }
alpar@9	438 }
alpar@9	439 A_ptr[n+1] = loc;
alpar@9	440 xassert(loc-1 == lp->nnz);
alpar@9	441 #if 1 /* 06/IV-2009 */
alpar@9	442 /* matrix A (by rows) */
alpar@9	443 loc = 1;
alpar@9	444 for (i = 1; i <= m; i++)
alpar@9	445 { GLPAIJ *aij;
alpar@9	446 AT_ptr[i] = loc;
alpar@9	447 for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next)
alpar@9	448 { AT_ind[loc] = aij->col->j;
alpar@9	449 AT_val[loc] = aij->row->rii * aij->val * aij->col->sjj;
alpar@9	450 loc++;
alpar@9	451 }
alpar@9	452 }
alpar@9	453 AT_ptr[m+1] = loc;
alpar@9	454 xassert(loc-1 == lp->nnz);
alpar@9	455 #endif
alpar@9	456 /* basis header */
alpar@9	457 xassert(lp->valid);
alpar@9	458 memcpy(&head[1], &lp->head[1], m * sizeof(int));
alpar@9	459 k = 0;
alpar@9	460 for (i = 1; i <= m; i++)
alpar@9	461 { GLPROW *row = lp->row[i];
alpar@9	462 if (row->stat != GLP_BS)
alpar@9	463 { k++;
alpar@9	464 xassert(k <= n);
alpar@9	465 head[m+k] = i;
alpar@9	466 stat[k] = (char)row->stat;
alpar@9	467 }
alpar@9	468 }
alpar@9	469 for (j = 1; j <= n; j++)
alpar@9	470 { GLPCOL *col = lp->col[j];
alpar@9	471 if (col->stat != GLP_BS)
alpar@9	472 { k++;
alpar@9	473 xassert(k <= n);
alpar@9	474 head[m+k] = m + j;
alpar@9	475 stat[k] = (char)col->stat;
alpar@9	476 }
alpar@9	477 }
alpar@9	478 xassert(k == n);
alpar@9	479 #if 1 /* 06/IV-2009 */
alpar@9	480 for (k = 1; k <= m+n; k++)
alpar@9	481 bind[head[k]] = k;
alpar@9	482 #endif
alpar@9	483 /* factorization of matrix B */
alpar@9	484 csa->valid = 1, lp->valid = 0;
alpar@9	485 csa->bfd = lp->bfd, lp->bfd = NULL;
alpar@9	486 #if 0 /* 06/IV-2009 */
alpar@9	487 /* matrix N (by rows) */
alpar@9	488 alloc_N(csa);
alpar@9	489 build_N(csa);
alpar@9	490 #endif
alpar@9	491 /* working parameters */
alpar@9	492 csa->phase = 0;
alpar@9	493 csa->tm_beg = xtime();
alpar@9	494 csa->it_beg = csa->it_cnt = lp->it_cnt;
alpar@9	495 csa->it_dpy = -1;
alpar@9	496 /* reference space and steepest edge coefficients */
alpar@9	497 csa->refct = 0;
alpar@9	498 memset(&refsp[1], 0, (m+n) * sizeof(char));
alpar@9	499 for (i = 1; i <= m; i++) gamma[i] = 1.0;
alpar@9	500 return;
alpar@9	501 }
alpar@9	502
alpar@9	503 #if 1 /* copied from primal */
alpar@9	504 /***********************************************************************
alpar@9	505 * invert_B - compute factorization of the basis matrix
alpar@9	506 *
alpar@9	507 * This routine computes factorization of the current basis matrix B.
alpar@9	508 *
alpar@9	509 * If the operation is successful, the routine returns zero, otherwise
alpar@9	510 * non-zero. */
alpar@9	511
alpar@9	512 static int inv_col(void *info, int i, int ind[], double val[])
alpar@9	513 { /* this auxiliary routine returns row indices and numeric values
alpar@9	514 of non-zero elements of i-th column of the basis matrix */
alpar@9	515 struct csa *csa = info;
alpar@9	516 int m = csa->m;
alpar@9	517 #ifdef GLP_DEBUG
alpar@9	518 int n = csa->n;
alpar@9	519 #endif
alpar@9	520 int *A_ptr = csa->A_ptr;
alpar@9	521 int *A_ind = csa->A_ind;
alpar@9	522 double *A_val = csa->A_val;
alpar@9	523 int *head = csa->head;
alpar@9	524 int k, len, ptr, t;
alpar@9	525 #ifdef GLP_DEBUG
alpar@9	526 xassert(1 <= i && i <= m);
alpar@9	527 #endif
alpar@9	528 k = head[i]; /* B[i] is k-th column of (I\|-A) */
alpar@9	529 #ifdef GLP_DEBUG
alpar@9	530 xassert(1 <= k && k <= m+n);
alpar@9	531 #endif
alpar@9	532 if (k <= m)
alpar@9	533 { /* B[i] is k-th column of submatrix I */
alpar@9	534 len = 1;
alpar@9	535 ind[1] = k;
alpar@9	536 val[1] = 1.0;
alpar@9	537 }
alpar@9	538 else
alpar@9	539 { /* B[i] is (k-m)-th column of submatrix (-A) */
alpar@9	540 ptr = A_ptr[k-m];
alpar@9	541 len = A_ptr[k-m+1] - ptr;
alpar@9	542 memcpy(&ind[1], &A_ind[ptr], len * sizeof(int));
alpar@9	543 memcpy(&val[1], &A_val[ptr], len * sizeof(double));
alpar@9	544 for (t = 1; t <= len; t++) val[t] = - val[t];
alpar@9	545 }
alpar@9	546 return len;
alpar@9	547 }
alpar@9	548
alpar@9	549 static int invert_B(struct csa *csa)
alpar@9	550 { int ret;
alpar@9	551 ret = bfd_factorize(csa->bfd, csa->m, NULL, inv_col, csa);
alpar@9	552 csa->valid = (ret == 0);
alpar@9	553 return ret;
alpar@9	554 }
alpar@9	555 #endif
alpar@9	556
alpar@9	557 #if 1 /* copied from primal */
alpar@9	558 /***********************************************************************
alpar@9	559 * update_B - update factorization of the basis matrix
alpar@9	560 *
alpar@9	561 * This routine replaces i-th column of the basis matrix B by k-th
alpar@9	562 * column of the augmented constraint matrix (I\|-A) and then updates
alpar@9	563 * the factorization of B.
alpar@9	564 *
alpar@9	565 * If the factorization has been successfully updated, the routine
alpar@9	566 * returns zero, otherwise non-zero. */
alpar@9	567
alpar@9	568 static int update_B(struct csa *csa, int i, int k)
alpar@9	569 { int m = csa->m;
alpar@9	570 #ifdef GLP_DEBUG
alpar@9	571 int n = csa->n;
alpar@9	572 #endif
alpar@9	573 int ret;
alpar@9	574 #ifdef GLP_DEBUG
alpar@9	575 xassert(1 <= i && i <= m);
alpar@9	576 xassert(1 <= k && k <= m+n);
alpar@9	577 #endif
alpar@9	578 if (k <= m)
alpar@9	579 { /* new i-th column of B is k-th column of I */
alpar@9	580 int ind[1+1];
alpar@9	581 double val[1+1];
alpar@9	582 ind[1] = k;
alpar@9	583 val[1] = 1.0;
alpar@9	584 xassert(csa->valid);
alpar@9	585 ret = bfd_update_it(csa->bfd, i, 0, 1, ind, val);
alpar@9	586 }
alpar@9	587 else
alpar@9	588 { /* new i-th column of B is (k-m)-th column of (-A) */
alpar@9	589 int *A_ptr = csa->A_ptr;
alpar@9	590 int *A_ind = csa->A_ind;
alpar@9	591 double *A_val = csa->A_val;
alpar@9	592 double *val = csa->work1;
alpar@9	593 int beg, end, ptr, len;
alpar@9	594 beg = A_ptr[k-m];
alpar@9	595 end = A_ptr[k-m+1];
alpar@9	596 len = 0;
alpar@9	597 for (ptr = beg; ptr < end; ptr++)
alpar@9	598 val[++len] = - A_val[ptr];
alpar@9	599 xassert(csa->valid);
alpar@9	600 ret = bfd_update_it(csa->bfd, i, 0, len, &A_ind[beg-1], val);
alpar@9	601 }
alpar@9	602 csa->valid = (ret == 0);
alpar@9	603 return ret;
alpar@9	604 }
alpar@9	605 #endif
alpar@9	606
alpar@9	607 #if 1 /* copied from primal */
alpar@9	608 /***********************************************************************
alpar@9	609 * error_ftran - compute residual vector r = h - B * x
alpar@9	610 *
alpar@9	611 * This routine computes the residual vector r = h - B * x, where B is
alpar@9	612 * the current basis matrix, h is the vector of right-hand sides, x is
alpar@9	613 * the solution vector. */
alpar@9	614
alpar@9	615 static void error_ftran(struct csa *csa, double h[], double x[],
alpar@9	616 double r[])
alpar@9	617 { int m = csa->m;
alpar@9	618 #ifdef GLP_DEBUG
alpar@9	619 int n = csa->n;
alpar@9	620 #endif
alpar@9	621 int *A_ptr = csa->A_ptr;
alpar@9	622 int *A_ind = csa->A_ind;
alpar@9	623 double *A_val = csa->A_val;
alpar@9	624 int *head = csa->head;
alpar@9	625 int i, k, beg, end, ptr;
alpar@9	626 double temp;
alpar@9	627 /* compute the residual vector:
alpar@9	628 r = h - B * x = h - B[1] * x[1] - ... - B[m] * x[m],
alpar@9	629 where B[1], ..., B[m] are columns of matrix B */
alpar@9	630 memcpy(&r[1], &h[1], m * sizeof(double));
alpar@9	631 for (i = 1; i <= m; i++)
alpar@9	632 { temp = x[i];
alpar@9	633 if (temp == 0.0) continue;
alpar@9	634 k = head[i]; /* B[i] is k-th column of (I\|-A) */
alpar@9	635 #ifdef GLP_DEBUG
alpar@9	636 xassert(1 <= k && k <= m+n);
alpar@9	637 #endif
alpar@9	638 if (k <= m)
alpar@9	639 { /* B[i] is k-th column of submatrix I */
alpar@9	640 r[k] -= temp;
alpar@9	641 }
alpar@9	642 else
alpar@9	643 { /* B[i] is (k-m)-th column of submatrix (-A) */
alpar@9	644 beg = A_ptr[k-m];
alpar@9	645 end = A_ptr[k-m+1];
alpar@9	646 for (ptr = beg; ptr < end; ptr++)
alpar@9	647 r[A_ind[ptr]] += A_val[ptr] * temp;
alpar@9	648 }
alpar@9	649 }
alpar@9	650 return;
alpar@9	651 }
alpar@9	652 #endif
alpar@9	653
alpar@9	654 #if 1 /* copied from primal */
alpar@9	655 /***********************************************************************
alpar@9	656 * refine_ftran - refine solution of B * x = h
alpar@9	657 *
alpar@9	658 * This routine performs one iteration to refine the solution of
alpar@9	659 * the system B * x = h, where B is the current basis matrix, h is the
alpar@9	660 * vector of right-hand sides, x is the solution vector. */
alpar@9	661
alpar@9	662 static void refine_ftran(struct csa *csa, double h[], double x[])
alpar@9	663 { int m = csa->m;
alpar@9	664 double *r = csa->work1;
alpar@9	665 double *d = csa->work1;
alpar@9	666 int i;
alpar@9	667 /* compute the residual vector r = h - B * x */
alpar@9	668 error_ftran(csa, h, x, r);
alpar@9	669 /* compute the correction vector d = inv(B) * r */
alpar@9	670 xassert(csa->valid);
alpar@9	671 bfd_ftran(csa->bfd, d);
alpar@9	672 /* refine the solution vector (new x) = (old x) + d */
alpar@9	673 for (i = 1; i <= m; i++) x[i] += d[i];
alpar@9	674 return;
alpar@9	675 }
alpar@9	676 #endif
alpar@9	677
alpar@9	678 #if 1 /* copied from primal */
alpar@9	679 /***********************************************************************
alpar@9	680 * error_btran - compute residual vector r = h - B'* x
alpar@9	681 *
alpar@9	682 * This routine computes the residual vector r = h - B'* x, where B'
alpar@9	683 * is a matrix transposed to the current basis matrix, h is the vector
alpar@9	684 * of right-hand sides, x is the solution vector. */
alpar@9	685
alpar@9	686 static void error_btran(struct csa *csa, double h[], double x[],
alpar@9	687 double r[])
alpar@9	688 { int m = csa->m;
alpar@9	689 #ifdef GLP_DEBUG
alpar@9	690 int n = csa->n;
alpar@9	691 #endif
alpar@9	692 int *A_ptr = csa->A_ptr;
alpar@9	693 int *A_ind = csa->A_ind;
alpar@9	694 double *A_val = csa->A_val;
alpar@9	695 int *head = csa->head;
alpar@9	696 int i, k, beg, end, ptr;
alpar@9	697 double temp;
alpar@9	698 /* compute the residual vector r = b - B'* x */
alpar@9	699 for (i = 1; i <= m; i++)
alpar@9	700 { /* r[i] := b[i] - (i-th column of B)'* x */
alpar@9	701 k = head[i]; /* B[i] is k-th column of (I\|-A) */
alpar@9	702 #ifdef GLP_DEBUG
alpar@9	703 xassert(1 <= k && k <= m+n);
alpar@9	704 #endif
alpar@9	705 temp = h[i];
alpar@9	706 if (k <= m)
alpar@9	707 { /* B[i] is k-th column of submatrix I */
alpar@9	708 temp -= x[k];
alpar@9	709 }
alpar@9	710 else
alpar@9	711 { /* B[i] is (k-m)-th column of submatrix (-A) */
alpar@9	712 beg = A_ptr[k-m];
alpar@9	713 end = A_ptr[k-m+1];
alpar@9	714 for (ptr = beg; ptr < end; ptr++)
alpar@9	715 temp += A_val[ptr] * x[A_ind[ptr]];
alpar@9	716 }
alpar@9	717 r[i] = temp;
alpar@9	718 }
alpar@9	719 return;
alpar@9	720 }
alpar@9	721 #endif
alpar@9	722
alpar@9	723 #if 1 /* copied from primal */
alpar@9	724 /***********************************************************************
alpar@9	725 * refine_btran - refine solution of B'* x = h
alpar@9	726 *
alpar@9	727 * This routine performs one iteration to refine the solution of the
alpar@9	728 * system B'* x = h, where B' is a matrix transposed to the current
alpar@9	729 * basis matrix, h is the vector of right-hand sides, x is the solution
alpar@9	730 * vector. */
alpar@9	731
alpar@9	732 static void refine_btran(struct csa *csa, double h[], double x[])
alpar@9	733 { int m = csa->m;
alpar@9	734 double *r = csa->work1;
alpar@9	735 double *d = csa->work1;
alpar@9	736 int i;
alpar@9	737 /* compute the residual vector r = h - B'* x */
alpar@9	738 error_btran(csa, h, x, r);
alpar@9	739 /* compute the correction vector d = inv(B') * r */
alpar@9	740 xassert(csa->valid);
alpar@9	741 bfd_btran(csa->bfd, d);
alpar@9	742 /* refine the solution vector (new x) = (old x) + d */
alpar@9	743 for (i = 1; i <= m; i++) x[i] += d[i];
alpar@9	744 return;
alpar@9	745 }
alpar@9	746 #endif
alpar@9	747
alpar@9	748 #if 1 /* copied from primal */
alpar@9	749 /***********************************************************************
alpar@9	750 * get_xN - determine current value of non-basic variable xN[j]
alpar@9	751 *
alpar@9	752 * This routine returns the current value of non-basic variable xN[j],
alpar@9	753 * which is a value of its active bound. */
alpar@9	754
alpar@9	755 static double get_xN(struct csa *csa, int j)
alpar@9	756 { int m = csa->m;
alpar@9	757 #ifdef GLP_DEBUG
alpar@9	758 int n = csa->n;
alpar@9	759 #endif
alpar@9	760 double *lb = csa->lb;
alpar@9	761 double *ub = csa->ub;
alpar@9	762 int *head = csa->head;
alpar@9	763 char *stat = csa->stat;
alpar@9	764 int k;
alpar@9	765 double xN;
alpar@9	766 #ifdef GLP_DEBUG
alpar@9	767 xassert(1 <= j && j <= n);
alpar@9	768 #endif
alpar@9	769 k = head[m+j]; /* x[k] = xN[j] */
alpar@9	770 #ifdef GLP_DEBUG
alpar@9	771 xassert(1 <= k && k <= m+n);
alpar@9	772 #endif
alpar@9	773 switch (stat[j])
alpar@9	774 { case GLP_NL:
alpar@9	775 /* x[k] is on its lower bound */
alpar@9	776 xN = lb[k]; break;
alpar@9	777 case GLP_NU:
alpar@9	778 /* x[k] is on its upper bound */
alpar@9	779 xN = ub[k]; break;
alpar@9	780 case GLP_NF:
alpar@9	781 /* x[k] is free non-basic variable */
alpar@9	782 xN = 0.0; break;
alpar@9	783 case GLP_NS:
alpar@9	784 /* x[k] is fixed non-basic variable */
alpar@9	785 xN = lb[k]; break;
alpar@9	786 default:
alpar@9	787 xassert(stat != stat);
alpar@9	788 }
alpar@9	789 return xN;
alpar@9	790 }
alpar@9	791 #endif
alpar@9	792
alpar@9	793 #if 1 /* copied from primal */
alpar@9	794 /***********************************************************************
alpar@9	795 * eval_beta - compute primal values of basic variables
alpar@9	796 *
alpar@9	797 * This routine computes current primal values of all basic variables:
alpar@9	798 *
alpar@9	799 * beta = - inv(B) * N * xN,
alpar@9	800 *
alpar@9	801 * where B is the current basis matrix, N is a matrix built of columns
alpar@9	802 * of matrix (I\|-A) corresponding to non-basic variables, and xN is the
alpar@9	803 * vector of current values of non-basic variables. */
alpar@9	804
alpar@9	805 static void eval_beta(struct csa *csa, double beta[])
alpar@9	806 { int m = csa->m;
alpar@9	807 int n = csa->n;
alpar@9	808 int *A_ptr = csa->A_ptr;
alpar@9	809 int *A_ind = csa->A_ind;
alpar@9	810 double *A_val = csa->A_val;
alpar@9	811 int *head = csa->head;
alpar@9	812 double *h = csa->work2;
alpar@9	813 int i, j, k, beg, end, ptr;
alpar@9	814 double xN;
alpar@9	815 /* compute the right-hand side vector:
alpar@9	816 h := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n],
alpar@9	817 where N[1], ..., N[n] are columns of matrix N */
alpar@9	818 for (i = 1; i <= m; i++)
alpar@9	819 h[i] = 0.0;
alpar@9	820 for (j = 1; j <= n; j++)
alpar@9	821 { k = head[m+j]; /* x[k] = xN[j] */
alpar@9	822 #ifdef GLP_DEBUG
alpar@9	823 xassert(1 <= k && k <= m+n);
alpar@9	824 #endif
alpar@9	825 /* determine current value of xN[j] */
alpar@9	826 xN = get_xN(csa, j);
alpar@9	827 if (xN == 0.0) continue;
alpar@9	828 if (k <= m)
alpar@9	829 { /* N[j] is k-th column of submatrix I */
alpar@9	830 h[k] -= xN;
alpar@9	831 }
alpar@9	832 else
alpar@9	833 { /* N[j] is (k-m)-th column of submatrix (-A) */
alpar@9	834 beg = A_ptr[k-m];
alpar@9	835 end = A_ptr[k-m+1];
alpar@9	836 for (ptr = beg; ptr < end; ptr++)
alpar@9	837 h[A_ind[ptr]] += xN * A_val[ptr];
alpar@9	838 }
alpar@9	839 }
alpar@9	840 /* solve system B * beta = h */
alpar@9	841 memcpy(&beta[1], &h[1], m * sizeof(double));
alpar@9	842 xassert(csa->valid);
alpar@9	843 bfd_ftran(csa->bfd, beta);
alpar@9	844 /* and refine the solution */
alpar@9	845 refine_ftran(csa, h, beta);
alpar@9	846 return;
alpar@9	847 }
alpar@9	848 #endif
alpar@9	849
alpar@9	850 #if 1 /* copied from primal */
alpar@9	851 /***********************************************************************
alpar@9	852 * eval_pi - compute vector of simplex multipliers
alpar@9	853 *
alpar@9	854 * This routine computes the vector of current simplex multipliers:
alpar@9	855 *
alpar@9	856 * pi = inv(B') * cB,
alpar@9	857 *
alpar@9	858 * where B' is a matrix transposed to the current basis matrix, cB is
alpar@9	859 * a subvector of objective coefficients at basic variables. */
alpar@9	860
alpar@9	861 static void eval_pi(struct csa *csa, double pi[])
alpar@9	862 { int m = csa->m;
alpar@9	863 double *c = csa->coef;
alpar@9	864 int *head = csa->head;
alpar@9	865 double *cB = csa->work2;
alpar@9	866 int i;
alpar@9	867 /* construct the right-hand side vector cB */
alpar@9	868 for (i = 1; i <= m; i++)
alpar@9	869 cB[i] = c[head[i]];
alpar@9	870 /* solve system B'* pi = cB */
alpar@9	871 memcpy(&pi[1], &cB[1], m * sizeof(double));
alpar@9	872 xassert(csa->valid);
alpar@9	873 bfd_btran(csa->bfd, pi);
alpar@9	874 /* and refine the solution */
alpar@9	875 refine_btran(csa, cB, pi);
alpar@9	876 return;
alpar@9	877 }
alpar@9	878 #endif
alpar@9	879
alpar@9	880 #if 1 /* copied from primal */
alpar@9	881 /***********************************************************************
alpar@9	882 * eval_cost - compute reduced cost of non-basic variable xN[j]
alpar@9	883 *
alpar@9	884 * This routine computes the current reduced cost of non-basic variable
alpar@9	885 * xN[j]:
alpar@9	886 *
alpar@9	887 * d[j] = cN[j] - N'[j] * pi,
alpar@9	888 *
alpar@9	889 * where cN[j] is the objective coefficient at variable xN[j], N[j] is
alpar@9	890 * a column of the augmented constraint matrix (I\|-A) corresponding to
alpar@9	891 * xN[j], pi is the vector of simplex multipliers. */
alpar@9	892
alpar@9	893 static double eval_cost(struct csa *csa, double pi[], int j)
alpar@9	894 { int m = csa->m;
alpar@9	895 #ifdef GLP_DEBUG
alpar@9	896 int n = csa->n;
alpar@9	897 #endif
alpar@9	898 double *coef = csa->coef;
alpar@9	899 int *head = csa->head;
alpar@9	900 int k;
alpar@9	901 double dj;
alpar@9	902 #ifdef GLP_DEBUG
alpar@9	903 xassert(1 <= j && j <= n);
alpar@9	904 #endif
alpar@9	905 k = head[m+j]; /* x[k] = xN[j] */
alpar@9	906 #ifdef GLP_DEBUG
alpar@9	907 xassert(1 <= k && k <= m+n);
alpar@9	908 #endif
alpar@9	909 dj = coef[k];
alpar@9	910 if (k <= m)
alpar@9	911 { /* N[j] is k-th column of submatrix I */
alpar@9	912 dj -= pi[k];
alpar@9	913 }
alpar@9	914 else
alpar@9	915 { /* N[j] is (k-m)-th column of submatrix (-A) */
alpar@9	916 int *A_ptr = csa->A_ptr;
alpar@9	917 int *A_ind = csa->A_ind;
alpar@9	918 double *A_val = csa->A_val;
alpar@9	919 int beg, end, ptr;
alpar@9	920 beg = A_ptr[k-m];
alpar@9	921 end = A_ptr[k-m+1];
alpar@9	922 for (ptr = beg; ptr < end; ptr++)
alpar@9	923 dj += A_val[ptr] * pi[A_ind[ptr]];
alpar@9	924 }
alpar@9	925 return dj;
alpar@9	926 }
alpar@9	927 #endif
alpar@9	928
alpar@9	929 #if 1 /* copied from primal */
alpar@9	930 /***********************************************************************
alpar@9	931 * eval_bbar - compute and store primal values of basic variables
alpar@9	932 *
alpar@9	933 * This routine computes primal values of all basic variables and then
alpar@9	934 * stores them in the solution array. */
alpar@9	935
alpar@9	936 static void eval_bbar(struct csa *csa)
alpar@9	937 { eval_beta(csa, csa->bbar);
alpar@9	938 return;
alpar@9	939 }
alpar@9	940 #endif
alpar@9	941
alpar@9	942 #if 1 /* copied from primal */
alpar@9	943 /***********************************************************************
alpar@9	944 * eval_cbar - compute and store reduced costs of non-basic variables
alpar@9	945 *
alpar@9	946 * This routine computes reduced costs of all non-basic variables and
alpar@9	947 * then stores them in the solution array. */
alpar@9	948
alpar@9	949 static void eval_cbar(struct csa *csa)
alpar@9	950 {
alpar@9	951 #ifdef GLP_DEBUG
alpar@9	952 int m = csa->m;
alpar@9	953 #endif
alpar@9	954 int n = csa->n;
alpar@9	955 #ifdef GLP_DEBUG
alpar@9	956 int *head = csa->head;
alpar@9	957 #endif
alpar@9	958 double *cbar = csa->cbar;
alpar@9	959 double *pi = csa->work3;
alpar@9	960 int j;
alpar@9	961 #ifdef GLP_DEBUG
alpar@9	962 int k;
alpar@9	963 #endif
alpar@9	964 /* compute simplex multipliers */
alpar@9	965 eval_pi(csa, pi);
alpar@9	966 /* compute and store reduced costs */
alpar@9	967 for (j = 1; j <= n; j++)
alpar@9	968 {
alpar@9	969 #ifdef GLP_DEBUG
alpar@9	970 k = head[m+j]; /* x[k] = xN[j] */
alpar@9	971 xassert(1 <= k && k <= m+n);
alpar@9	972 #endif
alpar@9	973 cbar[j] = eval_cost(csa, pi, j);
alpar@9	974 }
alpar@9	975 return;
alpar@9	976 }
alpar@9	977 #endif
alpar@9	978
alpar@9	979 /***********************************************************************
alpar@9	980 * reset_refsp - reset the reference space
alpar@9	981 *
alpar@9	982 * This routine resets (redefines) the reference space used in the
alpar@9	983 * projected steepest edge pricing algorithm. */
alpar@9	984
alpar@9	985 static void reset_refsp(struct csa *csa)
alpar@9	986 { int m = csa->m;
alpar@9	987 int n = csa->n;
alpar@9	988 int *head = csa->head;
alpar@9	989 char *refsp = csa->refsp;
alpar@9	990 double *gamma = csa->gamma;
alpar@9	991 int i, k;
alpar@9	992 xassert(csa->refct == 0);
alpar@9	993 csa->refct = 1000;
alpar@9	994 memset(&refsp[1], 0, (m+n) * sizeof(char));
alpar@9	995 for (i = 1; i <= m; i++)
alpar@9	996 { k = head[i]; /* x[k] = xB[i] */
alpar@9	997 refsp[k] = 1;
alpar@9	998 gamma[i] = 1.0;
alpar@9	999 }
alpar@9	1000 return;
alpar@9	1001 }
alpar@9	1002
alpar@9	1003 /***********************************************************************
alpar@9	1004 * eval_gamma - compute steepest edge coefficients
alpar@9	1005 *
alpar@9	1006 * This routine computes the vector of steepest edge coefficients for
alpar@9	1007 * all basic variables (except free ones) using its direct definition:
alpar@9	1008 *
alpar@9	1009 * gamma[i] = eta[i] + sum alfa[i,j]^2, i = 1,...,m,
alpar@9	1010 * j in C
alpar@9	1011 *
alpar@9	1012 * where eta[i] = 1 means that xB[i] is in the current reference space,
alpar@9	1013 * and 0 otherwise; C is a set of non-basic non-fixed variables xN[j],
alpar@9	1014 * which are in the current reference space; alfa[i,j] are elements of
alpar@9	1015 * the current simplex table.
alpar@9	1016 *
alpar@9	1017 * NOTE: The routine is intended only for debugginig purposes. */
alpar@9	1018
alpar@9	1019 static void eval_gamma(struct csa *csa, double gamma[])
alpar@9	1020 { int m = csa->m;
alpar@9	1021 int n = csa->n;
alpar@9	1022 char *type = csa->type;
alpar@9	1023 int *head = csa->head;
alpar@9	1024 char *refsp = csa->refsp;
alpar@9	1025 double *alfa = csa->work3;
alpar@9	1026 double *h = csa->work3;
alpar@9	1027 int i, j, k;
alpar@9	1028 /* gamma[i] := eta[i] (or 1, if xB[i] is free) */
alpar@9	1029 for (i = 1; i <= m; i++)
alpar@9	1030 { k = head[i]; /* x[k] = xB[i] */
alpar@9	1031 #ifdef GLP_DEBUG
alpar@9	1032 xassert(1 <= k && k <= m+n);
alpar@9	1033 #endif
alpar@9	1034 if (type[k] == GLP_FR)
alpar@9	1035 gamma[i] = 1.0;
alpar@9	1036 else
alpar@9	1037 gamma[i] = (refsp[k] ? 1.0 : 0.0);
alpar@9	1038 }
alpar@9	1039 /* compute columns of the current simplex table */
alpar@9	1040 for (j = 1; j <= n; j++)
alpar@9	1041 { k = head[m+j]; /* x[k] = xN[j] */
alpar@9	1042 #ifdef GLP_DEBUG
alpar@9	1043 xassert(1 <= k && k <= m+n);
alpar@9	1044 #endif
alpar@9	1045 /* skip column, if xN[j] is not in C */
alpar@9	1046 if (!refsp[k]) continue;
alpar@9	1047 #ifdef GLP_DEBUG
alpar@9	1048 /* set C must not contain fixed variables */
alpar@9	1049 xassert(type[k] != GLP_FX);
alpar@9	1050 #endif
alpar@9	1051 /* construct the right-hand side vector h = - N[j] */
alpar@9	1052 for (i = 1; i <= m; i++)
alpar@9	1053 h[i] = 0.0;
alpar@9	1054 if (k <= m)
alpar@9	1055 { /* N[j] is k-th column of submatrix I */
alpar@9	1056 h[k] = -1.0;
alpar@9	1057 }
alpar@9	1058 else
alpar@9	1059 { /* N[j] is (k-m)-th column of submatrix (-A) */
alpar@9	1060 int *A_ptr = csa->A_ptr;
alpar@9	1061 int *A_ind = csa->A_ind;
alpar@9	1062 double *A_val = csa->A_val;
alpar@9	1063 int beg, end, ptr;
alpar@9	1064 beg = A_ptr[k-m];
alpar@9	1065 end = A_ptr[k-m+1];
alpar@9	1066 for (ptr = beg; ptr < end; ptr++)
alpar@9	1067 h[A_ind[ptr]] = A_val[ptr];
alpar@9	1068 }
alpar@9	1069 /* solve system B * alfa = h */
alpar@9	1070 xassert(csa->valid);
alpar@9	1071 bfd_ftran(csa->bfd, alfa);
alpar@9	1072 /* gamma[i] := gamma[i] + alfa[i,j]^2 */
alpar@9	1073 for (i = 1; i <= m; i++)
alpar@9	1074 { k = head[i]; /* x[k] = xB[i] */
alpar@9	1075 if (type[k] != GLP_FR)
alpar@9	1076 gamma[i] += alfa[i] * alfa[i];
alpar@9	1077 }
alpar@9	1078 }
alpar@9	1079 return;
alpar@9	1080 }
alpar@9	1081
alpar@9	1082 /***********************************************************************
alpar@9	1083 * chuzr - choose basic variable (row of the simplex table)
alpar@9	1084 *
alpar@9	1085 * This routine chooses basic variable xB[p] having largest weighted
alpar@9	1086 * bound violation:
alpar@9	1087 *
alpar@9	1088 * \|r[p]\| / sqrt(gamma[p]) = max \|r[i]\| / sqrt(gamma[i]),
alpar@9	1089 * i in I
alpar@9	1090 *
alpar@9	1091 * / lB[i] - beta[i], if beta[i] < lB[i]
alpar@9	1092 * \|
alpar@9	1093 * r[i] = < 0, if lB[i] <= beta[i] <= uB[i]
alpar@9	1094 * \|
alpar@9	1095 * \ uB[i] - beta[i], if beta[i] > uB[i]
alpar@9	1096 *
alpar@9	1097 * where beta[i] is primal value of xB[i] in the current basis, lB[i]
alpar@9	1098 * and uB[i] are lower and upper bounds of xB[i], I is a subset of
alpar@9	1099 * eligible basic variables, which significantly violates their bounds,
alpar@9	1100 * gamma[i] is the steepest edge coefficient.
alpar@9	1101 *
alpar@9	1102 * If \|r[i]\| is less than a specified tolerance, xB[i] is not included
alpar@9	1103 * in I and therefore ignored.
alpar@9	1104 *
alpar@9	1105 * If I is empty and no variable has been chosen, p is set to 0. */
alpar@9	1106
alpar@9	1107 static void chuzr(struct csa *csa, double tol_bnd)
alpar@9	1108 { int m = csa->m;
alpar@9	1109 #ifdef GLP_DEBUG
alpar@9	1110 int n = csa->n;
alpar@9	1111 #endif
alpar@9	1112 char *type = csa->type;
alpar@9	1113 double *lb = csa->lb;
alpar@9	1114 double *ub = csa->ub;
alpar@9	1115 int *head = csa->head;
alpar@9	1116 double *bbar = csa->bbar;
alpar@9	1117 double *gamma = csa->gamma;
alpar@9	1118 int i, k, p;
alpar@9	1119 double delta, best, eps, ri, temp;
alpar@9	1120 /* nothing is chosen so far */
alpar@9	1121 p = 0, delta = 0.0, best = 0.0;
alpar@9	1122 /* look through the list of basic variables */
alpar@9	1123 for (i = 1; i <= m; i++)
alpar@9	1124 { k = head[i]; /* x[k] = xB[i] */
alpar@9	1125 #ifdef GLP_DEBUG
alpar@9	1126 xassert(1 <= k && k <= m+n);
alpar@9	1127 #endif
alpar@9	1128 /* determine bound violation ri[i] */
alpar@9	1129 ri = 0.0;
alpar@9	1130 if (type[k] == GLP_LO \|\| type[k] == GLP_DB \|\|
alpar@9	1131 type[k] == GLP_FX)
alpar@9	1132 { /* xB[i] has lower bound */
alpar@9	1133 eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
alpar@9	1134 if (bbar[i] < lb[k] - eps)
alpar@9	1135 { /* and significantly violates it */
alpar@9	1136 ri = lb[k] - bbar[i];
alpar@9	1137 }
alpar@9	1138 }
alpar@9	1139 if (type[k] == GLP_UP \|\| type[k] == GLP_DB \|\|
alpar@9	1140 type[k] == GLP_FX)
alpar@9	1141 { /* xB[i] has upper bound */
alpar@9	1142 eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
alpar@9	1143 if (bbar[i] > ub[k] + eps)
alpar@9	1144 { /* and significantly violates it */
alpar@9	1145 ri = ub[k] - bbar[i];
alpar@9	1146 }
alpar@9	1147 }
alpar@9	1148 /* if xB[i] is not eligible, skip it */
alpar@9	1149 if (ri == 0.0) continue;
alpar@9	1150 /* xB[i] is eligible basic variable; choose one with largest
alpar@9	1151 weighted bound violation */
alpar@9	1152 #ifdef GLP_DEBUG
alpar@9	1153 xassert(gamma[i] >= 0.0);
alpar@9	1154 #endif
alpar@9	1155 temp = gamma[i];
alpar@9	1156 if (temp < DBL_EPSILON) temp = DBL_EPSILON;
alpar@9	1157 temp = (ri * ri) / temp;
alpar@9	1158 if (best < temp)
alpar@9	1159 p = i, delta = ri, best = temp;
alpar@9	1160 }
alpar@9	1161 /* store the index of basic variable xB[p] chosen and its change
alpar@9	1162 in the adjacent basis */
alpar@9	1163 csa->p = p;
alpar@9	1164 csa->delta = delta;
alpar@9	1165 return;
alpar@9	1166 }
alpar@9	1167
alpar@9	1168 #if 1 /* copied from primal */
alpar@9	1169 /***********************************************************************
alpar@9	1170 * eval_rho - compute pivot row of the inverse
alpar@9	1171 *
alpar@9	1172 * This routine computes the pivot (p-th) row of the inverse inv(B),
alpar@9	1173 * which corresponds to basic variable xB[p] chosen:
alpar@9	1174 *
alpar@9	1175 * rho = inv(B') * e[p],
alpar@9	1176 *
alpar@9	1177 * where B' is a matrix transposed to the current basis matrix, e[p]
alpar@9	1178 * is unity vector. */
alpar@9	1179
alpar@9	1180 static void eval_rho(struct csa *csa, double rho[])
alpar@9	1181 { int m = csa->m;
alpar@9	1182 int p = csa->p;
alpar@9	1183 double *e = rho;
alpar@9	1184 int i;
alpar@9	1185 #ifdef GLP_DEBUG
alpar@9	1186 xassert(1 <= p && p <= m);
alpar@9	1187 #endif
alpar@9	1188 /* construct the right-hand side vector e[p] */
alpar@9	1189 for (i = 1; i <= m; i++)
alpar@9	1190 e[i] = 0.0;
alpar@9	1191 e[p] = 1.0;
alpar@9	1192 /* solve system B'* rho = e[p] */
alpar@9	1193 xassert(csa->valid);
alpar@9	1194 bfd_btran(csa->bfd, rho);
alpar@9	1195 return;
alpar@9	1196 }
alpar@9	1197 #endif
alpar@9	1198
alpar@9	1199 #if 1 /* copied from primal */
alpar@9	1200 /***********************************************************************
alpar@9	1201 * refine_rho - refine pivot row of the inverse
alpar@9	1202 *
alpar@9	1203 * This routine refines the pivot row of the inverse inv(B) assuming
alpar@9	1204 * that it was previously computed by the routine eval_rho. */
alpar@9	1205
alpar@9	1206 static void refine_rho(struct csa *csa, double rho[])
alpar@9	1207 { int m = csa->m;
alpar@9	1208 int p = csa->p;
alpar@9	1209 double *e = csa->work3;
alpar@9	1210 int i;
alpar@9	1211 #ifdef GLP_DEBUG
alpar@9	1212 xassert(1 <= p && p <= m);
alpar@9	1213 #endif
alpar@9	1214 /* construct the right-hand side vector e[p] */
alpar@9	1215 for (i = 1; i <= m; i++)
alpar@9	1216 e[i] = 0.0;
alpar@9	1217 e[p] = 1.0;
alpar@9	1218 /* refine solution of B'* rho = e[p] */
alpar@9	1219 refine_btran(csa, e, rho);
alpar@9	1220 return;
alpar@9	1221 }
alpar@9	1222 #endif
alpar@9	1223
alpar@9	1224 #if 1 /* 06/IV-2009 */
alpar@9	1225 /***********************************************************************
alpar@9	1226 * eval_trow - compute pivot row of the simplex table
alpar@9	1227 *
alpar@9	1228 * This routine computes the pivot row of the simplex table, which
alpar@9	1229 * corresponds to basic variable xB[p] chosen.
alpar@9	1230 *
alpar@9	1231 * The pivot row is the following vector:
alpar@9	1232 *
alpar@9	1233 * trow = T'* e[p] = - N'* inv(B') * e[p] = - N' * rho,
alpar@9	1234 *
alpar@9	1235 * where rho is the pivot row of the inverse inv(B) previously computed
alpar@9	1236 * by the routine eval_rho.
alpar@9	1237 *
alpar@9	1238 * Note that elements of the pivot row corresponding to fixed non-basic
alpar@9	1239 * variables are not computed.
alpar@9	1240 *
alpar@9	1241 * NOTES
alpar@9	1242 *
alpar@9	1243 * Computing pivot row of the simplex table is one of the most time
alpar@9	1244 * consuming operations, and for some instances it may take more than
alpar@9	1245 * 50% of the total solution time.
alpar@9	1246 *
alpar@9	1247 * In the current implementation there are two routines to compute the
alpar@9	1248 * pivot row. The routine eval_trow1 computes elements of the pivot row
alpar@9	1249 * as inner products of columns of the matrix N and the vector rho; it
alpar@9	1250 * is used when the vector rho is relatively dense. The routine
alpar@9	1251 * eval_trow2 computes the pivot row as a linear combination of rows of
alpar@9	1252 * the matrix N; it is used when the vector rho is relatively sparse. */
alpar@9	1253
alpar@9	1254 static void eval_trow1(struct csa *csa, double rho[])
alpar@9	1255 { int m = csa->m;
alpar@9	1256 int n = csa->n;
alpar@9	1257 int *A_ptr = csa->A_ptr;
alpar@9	1258 int *A_ind = csa->A_ind;
alpar@9	1259 double *A_val = csa->A_val;
alpar@9	1260 int *head = csa->head;
alpar@9	1261 char *stat = csa->stat;
alpar@9	1262 int *trow_ind = csa->trow_ind;
alpar@9	1263 double *trow_vec = csa->trow_vec;
alpar@9	1264 int j, k, beg, end, ptr, nnz;
alpar@9	1265 double temp;
alpar@9	1266 /* compute the pivot row as inner products of columns of the
alpar@9	1267 matrix N and vector rho: trow[j] = - rho * N[j] */
alpar@9	1268 nnz = 0;
alpar@9	1269 for (j = 1; j <= n; j++)
alpar@9	1270 { if (stat[j] == GLP_NS)
alpar@9	1271 { /* xN[j] is fixed */
alpar@9	1272 trow_vec[j] = 0.0;
alpar@9	1273 continue;
alpar@9	1274 }
alpar@9	1275 k = head[m+j]; /* x[k] = xN[j] */
alpar@9	1276 if (k <= m)
alpar@9	1277 { /* N[j] is k-th column of submatrix I */
alpar@9	1278 temp = - rho[k];
alpar@9	1279 }
alpar@9	1280 else
alpar@9	1281 { /* N[j] is (k-m)-th column of submatrix (-A) */
alpar@9	1282 beg = A_ptr[k-m], end = A_ptr[k-m+1];
alpar@9	1283 temp = 0.0;
alpar@9	1284 for (ptr = beg; ptr < end; ptr++)
alpar@9	1285 temp += rho[A_ind[ptr]] * A_val[ptr];
alpar@9	1286 }
alpar@9	1287 if (temp != 0.0)
alpar@9	1288 trow_ind[++nnz] = j;
alpar@9	1289 trow_vec[j] = temp;
alpar@9	1290 }
alpar@9	1291 csa->trow_nnz = nnz;
alpar@9	1292 return;
alpar@9	1293 }
alpar@9	1294
alpar@9	1295 static void eval_trow2(struct csa *csa, double rho[])
alpar@9	1296 { int m = csa->m;
alpar@9	1297 int n = csa->n;
alpar@9	1298 int *AT_ptr = csa->AT_ptr;
alpar@9	1299 int *AT_ind = csa->AT_ind;
alpar@9	1300 double *AT_val = csa->AT_val;
alpar@9	1301 int *bind = csa->bind;
alpar@9	1302 char *stat = csa->stat;
alpar@9	1303 int *trow_ind = csa->trow_ind;
alpar@9	1304 double *trow_vec = csa->trow_vec;
alpar@9	1305 int i, j, beg, end, ptr, nnz;
alpar@9	1306 double temp;
alpar@9	1307 /* clear the pivot row */
alpar@9	1308 for (j = 1; j <= n; j++)
alpar@9	1309 trow_vec[j] = 0.0;
alpar@9	1310 /* compute the pivot row as a linear combination of rows of the
alpar@9	1311 matrix N: trow = - rho[1] * N'[1] - ... - rho[m] * N'[m] */
alpar@9	1312 for (i = 1; i <= m; i++)
alpar@9	1313 { temp = rho[i];
alpar@9	1314 if (temp == 0.0) continue;
alpar@9	1315 /* trow := trow - rho[i] * N'[i] */
alpar@9	1316 j = bind[i] - m; /* x[i] = xN[j] */
alpar@9	1317 if (j >= 1 && stat[j] != GLP_NS)
alpar@9	1318 trow_vec[j] -= temp;
alpar@9	1319 beg = AT_ptr[i], end = AT_ptr[i+1];
alpar@9	1320 for (ptr = beg; ptr < end; ptr++)
alpar@9	1321 { j = bind[m + AT_ind[ptr]] - m; /* x[k] = xN[j] */
alpar@9	1322 if (j >= 1 && stat[j] != GLP_NS)
alpar@9	1323 trow_vec[j] += temp * AT_val[ptr];
alpar@9	1324 }
alpar@9	1325 }
alpar@9	1326 /* construct sparse pattern of the pivot row */
alpar@9	1327 nnz = 0;
alpar@9	1328 for (j = 1; j <= n; j++)
alpar@9	1329 { if (trow_vec[j] != 0.0)
alpar@9	1330 trow_ind[++nnz] = j;
alpar@9	1331 }
alpar@9	1332 csa->trow_nnz = nnz;
alpar@9	1333 return;
alpar@9	1334 }
alpar@9	1335
alpar@9	1336 static void eval_trow(struct csa *csa, double rho[])
alpar@9	1337 { int m = csa->m;
alpar@9	1338 int i, nnz;
alpar@9	1339 double dens;
alpar@9	1340 /* determine the density of the vector rho */
alpar@9	1341 nnz = 0;
alpar@9	1342 for (i = 1; i <= m; i++)
alpar@9	1343 if (rho[i] != 0.0) nnz++;
alpar@9	1344 dens = (double)nnz / (double)m;
alpar@9	1345 if (dens >= 0.20)
alpar@9	1346 { /* rho is relatively dense */
alpar@9	1347 eval_trow1(csa, rho);
alpar@9	1348 }
alpar@9	1349 else
alpar@9	1350 { /* rho is relatively sparse */
alpar@9	1351 eval_trow2(csa, rho);
alpar@9	1352 }
alpar@9	1353 return;
alpar@9	1354 }
alpar@9	1355 #endif
alpar@9	1356
alpar@9	1357 /***********************************************************************
alpar@9	1358 * sort_trow - sort pivot row of the simplex table
alpar@9	1359 *
alpar@9	1360 * This routine reorders the list of non-zero elements of the pivot
alpar@9	1361 * row to put significant elements, whose magnitude is not less than
alpar@9	1362 * a specified tolerance, in front of the list, and stores the number
alpar@9	1363 * of significant elements in trow_num. */
alpar@9	1364
alpar@9	1365 static void sort_trow(struct csa *csa, double tol_piv)
alpar@9	1366 {
alpar@9	1367 #ifdef GLP_DEBUG
alpar@9	1368 int n = csa->n;
alpar@9	1369 char *stat = csa->stat;
alpar@9	1370 #endif
alpar@9	1371 int nnz = csa->trow_nnz;
alpar@9	1372 int *trow_ind = csa->trow_ind;
alpar@9	1373 double *trow_vec = csa->trow_vec;
alpar@9	1374 int j, num, pos;
alpar@9	1375 double big, eps, temp;
alpar@9	1376 /* compute infinity (maximum) norm of the row */
alpar@9	1377 big = 0.0;
alpar@9	1378 for (pos = 1; pos <= nnz; pos++)
alpar@9	1379 {
alpar@9	1380 #ifdef GLP_DEBUG
alpar@9	1381 j = trow_ind[pos];
alpar@9	1382 xassert(1 <= j && j <= n);
alpar@9	1383 xassert(stat[j] != GLP_NS);
alpar@9	1384 #endif
alpar@9	1385 temp = fabs(trow_vec[trow_ind[pos]]);
alpar@9	1386 if (big < temp) big = temp;
alpar@9	1387 }
alpar@9	1388 csa->trow_max = big;
alpar@9	1389 /* determine absolute pivot tolerance */
alpar@9	1390 eps = tol_piv * (1.0 + 0.01 * big);
alpar@9	1391 /* move significant row components to the front of the list */
alpar@9	1392 for (num = 0; num < nnz; )
alpar@9	1393 { j = trow_ind[nnz];
alpar@9	1394 if (fabs(trow_vec[j]) < eps)
alpar@9	1395 nnz--;
alpar@9	1396 else
alpar@9	1397 { num++;
alpar@9	1398 trow_ind[nnz] = trow_ind[num];
alpar@9	1399 trow_ind[num] = j;
alpar@9	1400 }
alpar@9	1401 }
alpar@9	1402 csa->trow_num = num;
alpar@9	1403 return;
alpar@9	1404 }
alpar@9	1405
alpar@9	1406 #ifdef GLP_LONG_STEP /* 07/IV-2009 */
alpar@9	1407 static int ls_func(const void p1_, const void p2_)
alpar@9	1408 { const struct bkpt p1 = p1_, p2 = p2_;
alpar@9	1409 if (p1->t < p2->t) return -1;
alpar@9	1410 if (p1->t > p2->t) return +1;
alpar@9	1411 return 0;
alpar@9	1412 }
alpar@9	1413
alpar@9	1414 static int ls_func1(const void p1_, const void p2_)
alpar@9	1415 { const struct bkpt p1 = p1_, p2 = p2_;
alpar@9	1416 if (p1->dz < p2->dz) return -1;
alpar@9	1417 if (p1->dz > p2->dz) return +1;
alpar@9	1418 return 0;
alpar@9	1419 }
alpar@9	1420
alpar@9	1421 static void long_step(struct csa *csa)
alpar@9	1422 { int m = csa->m;
alpar@9	1423 #ifdef GLP_DEBUG
alpar@9	1424 int n = csa->n;
alpar@9	1425 #endif
alpar@9	1426 char *type = csa->type;
alpar@9	1427 double *lb = csa->lb;
alpar@9	1428 double *ub = csa->ub;
alpar@9	1429 int *head = csa->head;
alpar@9	1430 char *stat = csa->stat;
alpar@9	1431 double *cbar = csa->cbar;
alpar@9	1432 double delta = csa->delta;
alpar@9	1433 int *trow_ind = csa->trow_ind;
alpar@9	1434 double *trow_vec = csa->trow_vec;
alpar@9	1435 int trow_num = csa->trow_num;
alpar@9	1436 struct bkpt *bkpt = csa->bkpt;
alpar@9	1437 int j, k, kk, nbps, pos;
alpar@9	1438 double alfa, s, slope, dzmax;
alpar@9	1439 /* delta > 0 means that xB[p] violates its lower bound, so to
alpar@9	1440 increase the dual objective lambdaB[p] must increase;
alpar@9	1441 delta < 0 means that xB[p] violates its upper bound, so to
alpar@9	1442 increase the dual objective lambdaB[p] must decrease */
alpar@9	1443 /* s := sign(delta) */
alpar@9	1444 s = (delta > 0.0 ? +1.0 : -1.0);
alpar@9	1445 /* determine breakpoints of the dual objective */
alpar@9	1446 nbps = 0;
alpar@9	1447 for (pos = 1; pos <= trow_num; pos++)
alpar@9	1448 { j = trow_ind[pos];
alpar@9	1449 #ifdef GLP_DEBUG
alpar@9	1450 xassert(1 <= j && j <= n);
alpar@9	1451 xassert(stat[j] != GLP_NS);
alpar@9	1452 #endif
alpar@9	1453 /* if there is free non-basic variable, switch to the standard
alpar@9	1454 ratio test */
alpar@9	1455 if (stat[j] == GLP_NF)
alpar@9	1456 { nbps = 0;
alpar@9	1457 goto done;
alpar@9	1458 }
alpar@9	1459 /* lambdaN[j] = ... - alfa * t - ..., where t = s * lambdaB[i]
alpar@9	1460 is the dual ray parameter, t >= 0 */
alpar@9	1461 alfa = s * trow_vec[j];
alpar@9	1462 #ifdef GLP_DEBUG
alpar@9	1463 xassert(alfa != 0.0);
alpar@9	1464 xassert(stat[j] == GLP_NL \|\| stat[j] == GLP_NU);
alpar@9	1465 #endif
alpar@9	1466 if (alfa > 0.0 && stat[j] == GLP_NL \|\|
alpar@9	1467 alfa < 0.0 && stat[j] == GLP_NU)
alpar@9	1468 { /* either lambdaN[j] >= 0 (if stat = GLP_NL) and decreases
alpar@9	1469 or lambdaN[j] <= 0 (if stat = GLP_NU) and increases; in
alpar@9	1470 both cases we have a breakpoint */
alpar@9	1471 nbps++;
alpar@9	1472 #ifdef GLP_DEBUG
alpar@9	1473 xassert(nbps <= n);
alpar@9	1474 #endif
alpar@9	1475 bkpt[nbps].j = j;
alpar@9	1476 bkpt[nbps].t = cbar[j] / alfa;
alpar@9	1477 /*
alpar@9	1478 if (stat[j] == GLP_NL && cbar[j] < 0.0 \|\|
alpar@9	1479 stat[j] == GLP_NU && cbar[j] > 0.0)
alpar@9	1480 xprintf("%d %g\n", stat[j], cbar[j]);
alpar@9	1481 */
alpar@9	1482 /* if t is negative, replace it by exact zero (see comments
alpar@9	1483 in the routine chuzc) */
alpar@9	1484 if (bkpt[nbps].t < 0.0) bkpt[nbps].t = 0.0;
alpar@9	1485 }
alpar@9	1486 }
alpar@9	1487 /* if there are less than two breakpoints, switch to the standard
alpar@9	1488 ratio test */
alpar@9	1489 if (nbps < 2)
alpar@9	1490 { nbps = 0;
alpar@9	1491 goto done;
alpar@9	1492 }
alpar@9	1493 /* sort breakpoints by ascending the dual ray parameter, t */
alpar@9	1494 qsort(&bkpt[1], nbps, sizeof(struct bkpt), ls_func);
alpar@9	1495 /* determine last breakpoint, at which the dual objective still
alpar@9	1496 greater than at t = 0 */
alpar@9	1497 dzmax = 0.0;
alpar@9	1498 slope = fabs(delta); /* initial slope */
alpar@9	1499 for (kk = 1; kk <= nbps; kk++)
alpar@9	1500 { if (kk == 1)
alpar@9	1501 bkpt[kk].dz =
alpar@9	1502 0.0 + slope * (bkpt[kk].t - 0.0);
alpar@9	1503 else
alpar@9	1504 bkpt[kk].dz =
alpar@9	1505 bkpt[kk-1].dz + slope * (bkpt[kk].t - bkpt[kk-1].t);
alpar@9	1506 if (dzmax < bkpt[kk].dz)
alpar@9	1507 dzmax = bkpt[kk].dz;
alpar@9	1508 else if (bkpt[kk].dz < 0.05 * (1.0 + dzmax))
alpar@9	1509 { nbps = kk - 1;
alpar@9	1510 break;
alpar@9	1511 }
alpar@9	1512 j = bkpt[kk].j;
alpar@9	1513 k = head[m+j]; /* x[k] = xN[j] */
alpar@9	1514 if (type[k] == GLP_DB)
alpar@9	1515 slope -= fabs(trow_vec[j]) * (ub[k] - lb[k]);
alpar@9	1516 else
alpar@9	1517 { nbps = kk;
alpar@9	1518 break;
alpar@9	1519 }
alpar@9	1520 }
alpar@9	1521 /* if there are less than two breakpoints, switch to the standard
alpar@9	1522 ratio test */
alpar@9	1523 if (nbps < 2)
alpar@9	1524 { nbps = 0;
alpar@9	1525 goto done;
alpar@9	1526 }
alpar@9	1527 /* sort breakpoints by ascending the dual change, dz */
alpar@9	1528 qsort(&bkpt[1], nbps, sizeof(struct bkpt), ls_func1);
alpar@9	1529 /*
alpar@9	1530 for (kk = 1; kk <= nbps; kk++)
alpar@9	1531 xprintf("%d; t = %g; dz = %g\n", kk, bkpt[kk].t, bkpt[kk].dz);
alpar@9	1532 */
alpar@9	1533 done: csa->nbps = nbps;
alpar@9	1534 return;
alpar@9	1535 }
alpar@9	1536 #endif
alpar@9	1537
alpar@9	1538 /***********************************************************************
alpar@9	1539 * chuzc - choose non-basic variable (column of the simplex table)
alpar@9	1540 *
alpar@9	1541 * This routine chooses non-basic variable xN[q], which being entered
alpar@9	1542 * in the basis keeps dual feasibility of the basic solution.
alpar@9	1543 *
alpar@9	1544 * The parameter rtol is a relative tolerance used to relax zero bounds
alpar@9	1545 * of reduced costs of non-basic variables. If rtol = 0, the routine
alpar@9	1546 * implements the standard ratio test. Otherwise, if rtol > 0, the
alpar@9	1547 * routine implements Harris' two-pass ratio test. In the latter case
alpar@9	1548 * rtol should be about three times less than a tolerance used to check
alpar@9	1549 * dual feasibility. */
alpar@9	1550
alpar@9	1551 static void chuzc(struct csa *csa, double rtol)
alpar@9	1552 {
alpar@9	1553 #ifdef GLP_DEBUG
alpar@9	1554 int m = csa->m;
alpar@9	1555 int n = csa->n;
alpar@9	1556 #endif
alpar@9	1557 char *stat = csa->stat;
alpar@9	1558 double *cbar = csa->cbar;
alpar@9	1559 #ifdef GLP_DEBUG
alpar@9	1560 int p = csa->p;
alpar@9	1561 #endif
alpar@9	1562 double delta = csa->delta;
alpar@9	1563 int *trow_ind = csa->trow_ind;
alpar@9	1564 double *trow_vec = csa->trow_vec;
alpar@9	1565 int trow_num = csa->trow_num;
alpar@9	1566 int j, pos, q;
alpar@9	1567 double alfa, big, s, t, teta, tmax;
alpar@9	1568 #ifdef GLP_DEBUG
alpar@9	1569 xassert(1 <= p && p <= m);
alpar@9	1570 #endif
alpar@9	1571 /* delta > 0 means that xB[p] violates its lower bound and goes
alpar@9	1572 to it in the adjacent basis, so lambdaB[p] is increasing from
alpar@9	1573 its lower zero bound;
alpar@9	1574 delta < 0 means that xB[p] violates its upper bound and goes
alpar@9	1575 to it in the adjacent basis, so lambdaB[p] is decreasing from
alpar@9	1576 its upper zero bound */
alpar@9	1577 #ifdef GLP_DEBUG
alpar@9	1578 xassert(delta != 0.0);
alpar@9	1579 #endif
alpar@9	1580 /* s := sign(delta) */
alpar@9	1581 s = (delta > 0.0 ? +1.0 : -1.0);
alpar@9	1582 /* FIRST PASS */
alpar@9	1583 /* nothing is chosen so far */
alpar@9	1584 q = 0, teta = DBL_MAX, big = 0.0;
alpar@9	1585 /* walk through significant elements of the pivot row */
alpar@9	1586 for (pos = 1; pos <= trow_num; pos++)
alpar@9	1587 { j = trow_ind[pos];
alpar@9	1588 #ifdef GLP_DEBUG
alpar@9	1589 xassert(1 <= j && j <= n);
alpar@9	1590 #endif
alpar@9	1591 alfa = s * trow_vec[j];
alpar@9	1592 #ifdef GLP_DEBUG
alpar@9	1593 xassert(alfa != 0.0);
alpar@9	1594 #endif
alpar@9	1595 /* lambdaN[j] = ... - alfa * lambdaB[p] - ..., and due to s we
alpar@9	1596 need to consider only increasing lambdaB[p] */
alpar@9	1597 if (alfa > 0.0)
alpar@9	1598 { /* lambdaN[j] is decreasing */
alpar@9	1599 if (stat[j] == GLP_NL \|\| stat[j] == GLP_NF)
alpar@9	1600 { /* lambdaN[j] has zero lower bound */
alpar@9	1601 t = (cbar[j] + rtol) / alfa;
alpar@9	1602 }
alpar@9	1603 else
alpar@9	1604 { /* lambdaN[j] has no lower bound */
alpar@9	1605 continue;
alpar@9	1606 }
alpar@9	1607 }
alpar@9	1608 else
alpar@9	1609 { /* lambdaN[j] is increasing */
alpar@9	1610 if (stat[j] == GLP_NU \|\| stat[j] == GLP_NF)
alpar@9	1611 { /* lambdaN[j] has zero upper bound */
alpar@9	1612 t = (cbar[j] - rtol) / alfa;
alpar@9	1613 }
alpar@9	1614 else
alpar@9	1615 { /* lambdaN[j] has no upper bound */
alpar@9	1616 continue;
alpar@9	1617 }
alpar@9	1618 }
alpar@9	1619 /* t is a change of lambdaB[p], on which lambdaN[j] reaches
alpar@9	1620 its zero bound (possibly relaxed); since the basic solution
alpar@9	1621 is assumed to be dual feasible, t has to be non-negative by
alpar@9	1622 definition; however, it may happen that lambdaN[j] slightly
alpar@9	1623 (i.e. within a tolerance) violates its zero bound, that
alpar@9	1624 leads to negative t; in the latter case, if xN[j] is chosen,
alpar@9	1625 negative t means that lambdaB[p] changes in wrong direction
alpar@9	1626 that may cause wrong results on updating reduced costs;
alpar@9	1627 thus, if t is negative, we should replace it by exact zero
alpar@9	1628 assuming that lambdaN[j] is exactly on its zero bound, and
alpar@9	1629 violation appears due to round-off errors */
alpar@9	1630 if (t < 0.0) t = 0.0;
alpar@9	1631 /* apply minimal ratio test */
alpar@9	1632 if (teta > t \|\| teta == t && big < fabs(alfa))
alpar@9	1633 q = j, teta = t, big = fabs(alfa);
alpar@9	1634 }
alpar@9	1635 /* the second pass is skipped in the following cases: */
alpar@9	1636 /* if the standard ratio test is used */
alpar@9	1637 if (rtol == 0.0) goto done;
alpar@9	1638 /* if no non-basic variable has been chosen on the first pass */
alpar@9	1639 if (q == 0) goto done;
alpar@9	1640 /* if lambdaN[q] prevents lambdaB[p] from any change */
alpar@9	1641 if (teta == 0.0) goto done;
alpar@9	1642 /* SECOND PASS */
alpar@9	1643 /* here tmax is a maximal change of lambdaB[p], on which the
alpar@9	1644 solution remains dual feasible within a tolerance */
alpar@9	1645 #if 0
alpar@9	1646 tmax = (1.0 + 10.0 * DBL_EPSILON) * teta;
alpar@9	1647 #else
alpar@9	1648 tmax = teta;
alpar@9	1649 #endif
alpar@9	1650 /* nothing is chosen so far */
alpar@9	1651 q = 0, teta = DBL_MAX, big = 0.0;
alpar@9	1652 /* walk through significant elements of the pivot row */
alpar@9	1653 for (pos = 1; pos <= trow_num; pos++)
alpar@9	1654 { j = trow_ind[pos];
alpar@9	1655 #ifdef GLP_DEBUG
alpar@9	1656 xassert(1 <= j && j <= n);
alpar@9	1657 #endif
alpar@9	1658 alfa = s * trow_vec[j];
alpar@9	1659 #ifdef GLP_DEBUG
alpar@9	1660 xassert(alfa != 0.0);
alpar@9	1661 #endif
alpar@9	1662 /* lambdaN[j] = ... - alfa * lambdaB[p] - ..., and due to s we
alpar@9	1663 need to consider only increasing lambdaB[p] */
alpar@9	1664 if (alfa > 0.0)
alpar@9	1665 { /* lambdaN[j] is decreasing */
alpar@9	1666 if (stat[j] == GLP_NL \|\| stat[j] == GLP_NF)
alpar@9	1667 { /* lambdaN[j] has zero lower bound */
alpar@9	1668 t = cbar[j] / alfa;
alpar@9	1669 }
alpar@9	1670 else
alpar@9	1671 { /* lambdaN[j] has no lower bound */
alpar@9	1672 continue;
alpar@9	1673 }
alpar@9	1674 }
alpar@9	1675 else
alpar@9	1676 { /* lambdaN[j] is increasing */
alpar@9	1677 if (stat[j] == GLP_NU \|\| stat[j] == GLP_NF)
alpar@9	1678 { /* lambdaN[j] has zero upper bound */
alpar@9	1679 t = cbar[j] / alfa;
alpar@9	1680 }
alpar@9	1681 else
alpar@9	1682 { /* lambdaN[j] has no upper bound */
alpar@9	1683 continue;
alpar@9	1684 }
alpar@9	1685 }
alpar@9	1686 /* (see comments for the first pass) */
alpar@9	1687 if (t < 0.0) t = 0.0;
alpar@9	1688 /* t is a change of lambdaB[p], on which lambdaN[j] reaches
alpar@9	1689 its zero (lower or upper) bound; if t <= tmax, all reduced
alpar@9	1690 costs can violate their zero bounds only within relaxation
alpar@9	1691 tolerance rtol, so we can choose non-basic variable having
alpar@9	1692 largest influence coefficient to avoid possible numerical
alpar@9	1693 instability */
alpar@9	1694 if (t <= tmax && big < fabs(alfa))
alpar@9	1695 q = j, teta = t, big = fabs(alfa);
alpar@9	1696 }
alpar@9	1697 /* something must be chosen on the second pass */
alpar@9	1698 xassert(q != 0);
alpar@9	1699 done: /* store the index of non-basic variable xN[q] chosen */
alpar@9	1700 csa->q = q;
alpar@9	1701 /* store reduced cost of xN[q] in the adjacent basis */
alpar@9	1702 csa->new_dq = s * teta;
alpar@9	1703 return;
alpar@9	1704 }
alpar@9	1705
alpar@9	1706 #if 1 /* copied from primal */
alpar@9	1707 /***********************************************************************
alpar@9	1708 * eval_tcol - compute pivot column of the simplex table
alpar@9	1709 *
alpar@9	1710 * This routine computes the pivot column of the simplex table, which
alpar@9	1711 * corresponds to non-basic variable xN[q] chosen.
alpar@9	1712 *
alpar@9	1713 * The pivot column is the following vector:
alpar@9	1714 *
alpar@9	1715 * tcol = T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q],
alpar@9	1716 *
alpar@9	1717 * where B is the current basis matrix, N[q] is a column of the matrix
alpar@9	1718 * (I\|-A) corresponding to variable xN[q]. */
alpar@9	1719
alpar@9	1720 static void eval_tcol(struct csa *csa)
alpar@9	1721 { int m = csa->m;
alpar@9	1722 #ifdef GLP_DEBUG
alpar@9	1723 int n = csa->n;
alpar@9	1724 #endif
alpar@9	1725 int *head = csa->head;
alpar@9	1726 int q = csa->q;
alpar@9	1727 int *tcol_ind = csa->tcol_ind;
alpar@9	1728 double *tcol_vec = csa->tcol_vec;
alpar@9	1729 double *h = csa->tcol_vec;
alpar@9	1730 int i, k, nnz;
alpar@9	1731 #ifdef GLP_DEBUG
alpar@9	1732 xassert(1 <= q && q <= n);
alpar@9	1733 #endif
alpar@9	1734 k = head[m+q]; /* x[k] = xN[q] */
alpar@9	1735 #ifdef GLP_DEBUG
alpar@9	1736 xassert(1 <= k && k <= m+n);
alpar@9	1737 #endif
alpar@9	1738 /* construct the right-hand side vector h = - N[q] */
alpar@9	1739 for (i = 1; i <= m; i++)
alpar@9	1740 h[i] = 0.0;
alpar@9	1741 if (k <= m)
alpar@9	1742 { /* N[q] is k-th column of submatrix I */
alpar@9	1743 h[k] = -1.0;
alpar@9	1744 }
alpar@9	1745 else
alpar@9	1746 { /* N[q] is (k-m)-th column of submatrix (-A) */
alpar@9	1747 int *A_ptr = csa->A_ptr;
alpar@9	1748 int *A_ind = csa->A_ind;
alpar@9	1749 double *A_val = csa->A_val;
alpar@9	1750 int beg, end, ptr;
alpar@9	1751 beg = A_ptr[k-m];
alpar@9	1752 end = A_ptr[k-m+1];
alpar@9	1753 for (ptr = beg; ptr < end; ptr++)
alpar@9	1754 h[A_ind[ptr]] = A_val[ptr];
alpar@9	1755 }
alpar@9	1756 /* solve system B * tcol = h */
alpar@9	1757 xassert(csa->valid);
alpar@9	1758 bfd_ftran(csa->bfd, tcol_vec);
alpar@9	1759 /* construct sparse pattern of the pivot column */
alpar@9	1760 nnz = 0;
alpar@9	1761 for (i = 1; i <= m; i++)
alpar@9	1762 { if (tcol_vec[i] != 0.0)
alpar@9	1763 tcol_ind[++nnz] = i;
alpar@9	1764 }
alpar@9	1765 csa->tcol_nnz = nnz;
alpar@9	1766 return;
alpar@9	1767 }
alpar@9	1768 #endif
alpar@9	1769
alpar@9	1770 #if 1 /* copied from primal */
alpar@9	1771 /***********************************************************************
alpar@9	1772 * refine_tcol - refine pivot column of the simplex table
alpar@9	1773 *
alpar@9	1774 * This routine refines the pivot column of the simplex table assuming
alpar@9	1775 * that it was previously computed by the routine eval_tcol. */
alpar@9	1776
alpar@9	1777 static void refine_tcol(struct csa *csa)
alpar@9	1778 { int m = csa->m;
alpar@9	1779 #ifdef GLP_DEBUG
alpar@9	1780 int n = csa->n;
alpar@9	1781 #endif
alpar@9	1782 int *head = csa->head;
alpar@9	1783 int q = csa->q;
alpar@9	1784 int *tcol_ind = csa->tcol_ind;
alpar@9	1785 double *tcol_vec = csa->tcol_vec;
alpar@9	1786 double *h = csa->work3;
alpar@9	1787 int i, k, nnz;
alpar@9	1788 #ifdef GLP_DEBUG
alpar@9	1789 xassert(1 <= q && q <= n);
alpar@9	1790 #endif
alpar@9	1791 k = head[m+q]; /* x[k] = xN[q] */
alpar@9	1792 #ifdef GLP_DEBUG
alpar@9	1793 xassert(1 <= k && k <= m+n);
alpar@9	1794 #endif
alpar@9	1795 /* construct the right-hand side vector h = - N[q] */
alpar@9	1796 for (i = 1; i <= m; i++)
alpar@9	1797 h[i] = 0.0;
alpar@9	1798 if (k <= m)
alpar@9	1799 { /* N[q] is k-th column of submatrix I */
alpar@9	1800 h[k] = -1.0;
alpar@9	1801 }
alpar@9	1802 else
alpar@9	1803 { /* N[q] is (k-m)-th column of submatrix (-A) */
alpar@9	1804 int *A_ptr = csa->A_ptr;
alpar@9	1805 int *A_ind = csa->A_ind;
alpar@9	1806 double *A_val = csa->A_val;
alpar@9	1807 int beg, end, ptr;
alpar@9	1808 beg = A_ptr[k-m];
alpar@9	1809 end = A_ptr[k-m+1];
alpar@9	1810 for (ptr = beg; ptr < end; ptr++)
alpar@9	1811 h[A_ind[ptr]] = A_val[ptr];
alpar@9	1812 }
alpar@9	1813 /* refine solution of B * tcol = h */
alpar@9	1814 refine_ftran(csa, h, tcol_vec);
alpar@9	1815 /* construct sparse pattern of the pivot column */
alpar@9	1816 nnz = 0;
alpar@9	1817 for (i = 1; i <= m; i++)
alpar@9	1818 { if (tcol_vec[i] != 0.0)
alpar@9	1819 tcol_ind[++nnz] = i;
alpar@9	1820 }
alpar@9	1821 csa->tcol_nnz = nnz;
alpar@9	1822 return;
alpar@9	1823 }
alpar@9	1824 #endif
alpar@9	1825
alpar@9	1826 /***********************************************************************
alpar@9	1827 * update_cbar - update reduced costs of non-basic variables
alpar@9	1828 *
alpar@9	1829 * This routine updates reduced costs of all (except fixed) non-basic
alpar@9	1830 * variables for the adjacent basis. */
alpar@9	1831
alpar@9	1832 static void update_cbar(struct csa *csa)
alpar@9	1833 {
alpar@9	1834 #ifdef GLP_DEBUG
alpar@9	1835 int n = csa->n;
alpar@9	1836 #endif
alpar@9	1837 double *cbar = csa->cbar;
alpar@9	1838 int trow_nnz = csa->trow_nnz;
alpar@9	1839 int *trow_ind = csa->trow_ind;
alpar@9	1840 double *trow_vec = csa->trow_vec;
alpar@9	1841 int q = csa->q;
alpar@9	1842 double new_dq = csa->new_dq;
alpar@9	1843 int j, pos;
alpar@9	1844 #ifdef GLP_DEBUG
alpar@9	1845 xassert(1 <= q && q <= n);
alpar@9	1846 #endif
alpar@9	1847 /* set new reduced cost of xN[q] */
alpar@9	1848 cbar[q] = new_dq;
alpar@9	1849 /* update reduced costs of other non-basic variables */
alpar@9	1850 if (new_dq == 0.0) goto done;
alpar@9	1851 for (pos = 1; pos <= trow_nnz; pos++)
alpar@9	1852 { j = trow_ind[pos];
alpar@9	1853 #ifdef GLP_DEBUG
alpar@9	1854 xassert(1 <= j && j <= n);
alpar@9	1855 #endif
alpar@9	1856 if (j != q)
alpar@9	1857 cbar[j] -= trow_vec[j] * new_dq;
alpar@9	1858 }
alpar@9	1859 done: return;
alpar@9	1860 }
alpar@9	1861
alpar@9	1862 /***********************************************************************
alpar@9	1863 * update_bbar - update values of basic variables
alpar@9	1864 *
alpar@9	1865 * This routine updates values of all basic variables for the adjacent
alpar@9	1866 * basis. */
alpar@9	1867
alpar@9	1868 static void update_bbar(struct csa *csa)
alpar@9	1869 {
alpar@9	1870 #ifdef GLP_DEBUG
alpar@9	1871 int m = csa->m;
alpar@9	1872 int n = csa->n;
alpar@9	1873 #endif
alpar@9	1874 double *bbar = csa->bbar;
alpar@9	1875 int p = csa->p;
alpar@9	1876 double delta = csa->delta;
alpar@9	1877 int q = csa->q;
alpar@9	1878 int tcol_nnz = csa->tcol_nnz;
alpar@9	1879 int *tcol_ind = csa->tcol_ind;
alpar@9	1880 double *tcol_vec = csa->tcol_vec;
alpar@9	1881 int i, pos;
alpar@9	1882 double teta;
alpar@9	1883 #ifdef GLP_DEBUG
alpar@9	1884 xassert(1 <= p && p <= m);
alpar@9	1885 xassert(1 <= q && q <= n);
alpar@9	1886 #endif
alpar@9	1887 /* determine the change of xN[q] in the adjacent basis */
alpar@9	1888 #ifdef GLP_DEBUG
alpar@9	1889 xassert(tcol_vec[p] != 0.0);
alpar@9	1890 #endif
alpar@9	1891 teta = delta / tcol_vec[p];
alpar@9	1892 /* set new primal value of xN[q] */
alpar@9	1893 bbar[p] = get_xN(csa, q) + teta;
alpar@9	1894 /* update primal values of other basic variables */
alpar@9	1895 if (teta == 0.0) goto done;
alpar@9	1896 for (pos = 1; pos <= tcol_nnz; pos++)
alpar@9	1897 { i = tcol_ind[pos];
alpar@9	1898 #ifdef GLP_DEBUG
alpar@9	1899 xassert(1 <= i && i <= m);
alpar@9	1900 #endif
alpar@9	1901 if (i != p)
alpar@9	1902 bbar[i] += tcol_vec[i] * teta;
alpar@9	1903 }
alpar@9	1904 done: return;
alpar@9	1905 }
alpar@9	1906
alpar@9	1907 /***********************************************************************
alpar@9	1908 * update_gamma - update steepest edge coefficients
alpar@9	1909 *
alpar@9	1910 * This routine updates steepest-edge coefficients for the adjacent
alpar@9	1911 * basis. */
alpar@9	1912
alpar@9	1913 static void update_gamma(struct csa *csa)
alpar@9	1914 { int m = csa->m;
alpar@9	1915 #ifdef GLP_DEBUG
alpar@9	1916 int n = csa->n;
alpar@9	1917 #endif
alpar@9	1918 char *type = csa->type;
alpar@9	1919 int *head = csa->head;
alpar@9	1920 char *refsp = csa->refsp;
alpar@9	1921 double *gamma = csa->gamma;
alpar@9	1922 int p = csa->p;
alpar@9	1923 int trow_nnz = csa->trow_nnz;
alpar@9	1924 int *trow_ind = csa->trow_ind;
alpar@9	1925 double *trow_vec = csa->trow_vec;
alpar@9	1926 int q = csa->q;
alpar@9	1927 int tcol_nnz = csa->tcol_nnz;
alpar@9	1928 int *tcol_ind = csa->tcol_ind;
alpar@9	1929 double *tcol_vec = csa->tcol_vec;
alpar@9	1930 double *u = csa->work3;
alpar@9	1931 int i, j, k,pos;
alpar@9	1932 double gamma_p, eta_p, pivot, t, t1, t2;
alpar@9	1933 #ifdef GLP_DEBUG
alpar@9	1934 xassert(1 <= p && p <= m);
alpar@9	1935 xassert(1 <= q && q <= n);
alpar@9	1936 #endif
alpar@9	1937 /* the basis changes, so decrease the count */
alpar@9	1938 xassert(csa->refct > 0);
alpar@9	1939 csa->refct--;
alpar@9	1940 /* recompute gamma[p] for the current basis more accurately and
alpar@9	1941 compute auxiliary vector u */
alpar@9	1942 #ifdef GLP_DEBUG
alpar@9	1943 xassert(type[head[p]] != GLP_FR);
alpar@9	1944 #endif
alpar@9	1945 gamma_p = eta_p = (refsp[head[p]] ? 1.0 : 0.0);
alpar@9	1946 for (i = 1; i <= m; i++) u[i] = 0.0;
alpar@9	1947 for (pos = 1; pos <= trow_nnz; pos++)
alpar@9	1948 { j = trow_ind[pos];
alpar@9	1949 #ifdef GLP_DEBUG
alpar@9	1950 xassert(1 <= j && j <= n);
alpar@9	1951 #endif
alpar@9	1952 k = head[m+j]; /* x[k] = xN[j] */
alpar@9	1953 #ifdef GLP_DEBUG
alpar@9	1954 xassert(1 <= k && k <= m+n);
alpar@9	1955 xassert(type[k] != GLP_FX);
alpar@9	1956 #endif
alpar@9	1957 if (!refsp[k]) continue;
alpar@9	1958 t = trow_vec[j];
alpar@9	1959 gamma_p += t * t;
alpar@9	1960 /* u := u + N[j] * delta[j] * trow[j] */
alpar@9	1961 if (k <= m)
alpar@9	1962 { /* N[k] = k-j stolbec submatrix I */
alpar@9	1963 u[k] += t;
alpar@9	1964 }
alpar@9	1965 else
alpar@9	1966 { /* N[k] = k-m-k stolbec (-A) */
alpar@9	1967 int *A_ptr = csa->A_ptr;
alpar@9	1968 int *A_ind = csa->A_ind;
alpar@9	1969 double *A_val = csa->A_val;
alpar@9	1970 int beg, end, ptr;
alpar@9	1971 beg = A_ptr[k-m];
alpar@9	1972 end = A_ptr[k-m+1];
alpar@9	1973 for (ptr = beg; ptr < end; ptr++)
alpar@9	1974 u[A_ind[ptr]] -= t * A_val[ptr];
alpar@9	1975 }
alpar@9	1976 }
alpar@9	1977 xassert(csa->valid);
alpar@9	1978 bfd_ftran(csa->bfd, u);
alpar@9	1979 /* update gamma[i] for other basic variables (except xB[p] and
alpar@9	1980 free variables) */
alpar@9	1981 pivot = tcol_vec[p];
alpar@9	1982 #ifdef GLP_DEBUG
alpar@9	1983 xassert(pivot != 0.0);
alpar@9	1984 #endif
alpar@9	1985 for (pos = 1; pos <= tcol_nnz; pos++)
alpar@9	1986 { i = tcol_ind[pos];
alpar@9	1987 #ifdef GLP_DEBUG
alpar@9	1988 xassert(1 <= i && i <= m);
alpar@9	1989 #endif
alpar@9	1990 k = head[i];
alpar@9	1991 #ifdef GLP_DEBUG
alpar@9	1992 xassert(1 <= k && k <= m+n);
alpar@9	1993 #endif
alpar@9	1994 /* skip xB[p] */
alpar@9	1995 if (i == p) continue;
alpar@9	1996 /* skip free basic variable */
alpar@9	1997 if (type[head[i]] == GLP_FR)
alpar@9	1998 {
alpar@9	1999 #ifdef GLP_DEBUG
alpar@9	2000 xassert(gamma[i] == 1.0);
alpar@9	2001 #endif
alpar@9	2002 continue;
alpar@9	2003 }
alpar@9	2004 /* compute gamma[i] for the adjacent basis */
alpar@9	2005 t = tcol_vec[i] / pivot;
alpar@9	2006 t1 = gamma[i] + t * t * gamma_p + 2.0 * t * u[i];
alpar@9	2007 t2 = (refsp[k] ? 1.0 : 0.0) + eta_p * t * t;
alpar@9	2008 gamma[i] = (t1 >= t2 ? t1 : t2);
alpar@9	2009 /* (though gamma[i] can be exact zero, because the reference
alpar@9	2010 space does not include non-basic fixed variables) */
alpar@9	2011 if (gamma[i] < DBL_EPSILON) gamma[i] = DBL_EPSILON;
alpar@9	2012 }
alpar@9	2013 /* compute gamma[p] for the adjacent basis */
alpar@9	2014 if (type[head[m+q]] == GLP_FR)
alpar@9	2015 gamma[p] = 1.0;
alpar@9	2016 else
alpar@9	2017 { gamma[p] = gamma_p / (pivot * pivot);
alpar@9	2018 if (gamma[p] < DBL_EPSILON) gamma[p] = DBL_EPSILON;
alpar@9	2019 }
alpar@9	2020 /* if xB[p], which becomes xN[q] in the adjacent basis, is fixed
alpar@9	2021 and belongs to the reference space, remove it from there, and
alpar@9	2022 change all gamma's appropriately */
alpar@9	2023 k = head[p];
alpar@9	2024 if (type[k] == GLP_FX && refsp[k])
alpar@9	2025 { refsp[k] = 0;
alpar@9	2026 for (pos = 1; pos <= tcol_nnz; pos++)
alpar@9	2027 { i = tcol_ind[pos];
alpar@9	2028 if (i == p)
alpar@9	2029 { if (type[head[m+q]] == GLP_FR) continue;
alpar@9	2030 t = 1.0 / tcol_vec[p];
alpar@9	2031 }
alpar@9	2032 else
alpar@9	2033 { if (type[head[i]] == GLP_FR) continue;
alpar@9	2034 t = tcol_vec[i] / tcol_vec[p];
alpar@9	2035 }
alpar@9	2036 gamma[i] -= t * t;
alpar@9	2037 if (gamma[i] < DBL_EPSILON) gamma[i] = DBL_EPSILON;
alpar@9	2038 }
alpar@9	2039 }
alpar@9	2040 return;
alpar@9	2041 }
alpar@9	2042
alpar@9	2043 #if 1 /* copied from primal */
alpar@9	2044 /***********************************************************************
alpar@9	2045 * err_in_bbar - compute maximal relative error in primal solution
alpar@9	2046 *
alpar@9	2047 * This routine returns maximal relative error:
alpar@9	2048 *
alpar@9	2049 * max \|beta[i] - bbar[i]\| / (1 + \|beta[i]\|),
alpar@9	2050 *
alpar@9	2051 * where beta and bbar are, respectively, directly computed and the
alpar@9	2052 * current (updated) values of basic variables.
alpar@9	2053 *
alpar@9	2054 * NOTE: The routine is intended only for debugginig purposes. */
alpar@9	2055
alpar@9	2056 static double err_in_bbar(struct csa *csa)
alpar@9	2057 { int m = csa->m;
alpar@9	2058 double *bbar = csa->bbar;
alpar@9	2059 int i;
alpar@9	2060 double e, emax, *beta;
alpar@9	2061 beta = xcalloc(1+m, sizeof(double));
alpar@9	2062 eval_beta(csa, beta);
alpar@9	2063 emax = 0.0;
alpar@9	2064 for (i = 1; i <= m; i++)
alpar@9	2065 { e = fabs(beta[i] - bbar[i]) / (1.0 + fabs(beta[i]));
alpar@9	2066 if (emax < e) emax = e;
alpar@9	2067 }
alpar@9	2068 xfree(beta);
alpar@9	2069 return emax;
alpar@9	2070 }
alpar@9	2071 #endif
alpar@9	2072
alpar@9	2073 #if 1 /* copied from primal */
alpar@9	2074 /***********************************************************************
alpar@9	2075 * err_in_cbar - compute maximal relative error in dual solution
alpar@9	2076 *
alpar@9	2077 * This routine returns maximal relative error:
alpar@9	2078 *
alpar@9	2079 * max \|cost[j] - cbar[j]\| / (1 + \|cost[j]\|),
alpar@9	2080 *
alpar@9	2081 * where cost and cbar are, respectively, directly computed and the
alpar@9	2082 * current (updated) reduced costs of non-basic non-fixed variables.
alpar@9	2083 *
alpar@9	2084 * NOTE: The routine is intended only for debugginig purposes. */
alpar@9	2085
alpar@9	2086 static double err_in_cbar(struct csa *csa)
alpar@9	2087 { int m = csa->m;
alpar@9	2088 int n = csa->n;
alpar@9	2089 char *stat = csa->stat;
alpar@9	2090 double *cbar = csa->cbar;
alpar@9	2091 int j;
alpar@9	2092 double e, emax, cost, *pi;
alpar@9	2093 pi = xcalloc(1+m, sizeof(double));
alpar@9	2094 eval_pi(csa, pi);
alpar@9	2095 emax = 0.0;
alpar@9	2096 for (j = 1; j <= n; j++)
alpar@9	2097 { if (stat[j] == GLP_NS) continue;
alpar@9	2098 cost = eval_cost(csa, pi, j);
alpar@9	2099 e = fabs(cost - cbar[j]) / (1.0 + fabs(cost));
alpar@9	2100 if (emax < e) emax = e;
alpar@9	2101 }
alpar@9	2102 xfree(pi);
alpar@9	2103 return emax;
alpar@9	2104 }
alpar@9	2105 #endif
alpar@9	2106
alpar@9	2107 /***********************************************************************
alpar@9	2108 * err_in_gamma - compute maximal relative error in steepest edge cff.
alpar@9	2109 *
alpar@9	2110 * This routine returns maximal relative error:
alpar@9	2111 *
alpar@9	2112 * max \|gamma'[j] - gamma[j]\| / (1 + \|gamma'[j]),
alpar@9	2113 *
alpar@9	2114 * where gamma'[j] and gamma[j] are, respectively, directly computed
alpar@9	2115 * and the current (updated) steepest edge coefficients for non-basic
alpar@9	2116 * non-fixed variable x[j].
alpar@9	2117 *
alpar@9	2118 * NOTE: The routine is intended only for debugginig purposes. */
alpar@9	2119
alpar@9	2120 static double err_in_gamma(struct csa *csa)
alpar@9	2121 { int m = csa->m;
alpar@9	2122 char *type = csa->type;
alpar@9	2123 int *head = csa->head;
alpar@9	2124 double *gamma = csa->gamma;
alpar@9	2125 double *exact = csa->work4;
alpar@9	2126 int i;
alpar@9	2127 double e, emax, temp;
alpar@9	2128 eval_gamma(csa, exact);
alpar@9	2129 emax = 0.0;
alpar@9	2130 for (i = 1; i <= m; i++)
alpar@9	2131 { if (type[head[i]] == GLP_FR)
alpar@9	2132 { xassert(gamma[i] == 1.0);
alpar@9	2133 xassert(exact[i] == 1.0);
alpar@9	2134 continue;
alpar@9	2135 }
alpar@9	2136 temp = exact[i];
alpar@9	2137 e = fabs(temp - gamma[i]) / (1.0 + fabs(temp));
alpar@9	2138 if (emax < e) emax = e;
alpar@9	2139 }
alpar@9	2140 return emax;
alpar@9	2141 }
alpar@9	2142
alpar@9	2143 /***********************************************************************
alpar@9	2144 * change_basis - change basis header
alpar@9	2145 *
alpar@9	2146 * This routine changes the basis header to make it corresponding to
alpar@9	2147 * the adjacent basis. */
alpar@9	2148
alpar@9	2149 static void change_basis(struct csa *csa)
alpar@9	2150 { int m = csa->m;
alpar@9	2151 #ifdef GLP_DEBUG
alpar@9	2152 int n = csa->n;
alpar@9	2153 #endif
alpar@9	2154 char *type = csa->type;
alpar@9	2155 int *head = csa->head;
alpar@9	2156 #if 1 /* 06/IV-2009 */
alpar@9	2157 int *bind = csa->bind;
alpar@9	2158 #endif
alpar@9	2159 char *stat = csa->stat;
alpar@9	2160 int p = csa->p;
alpar@9	2161 double delta = csa->delta;
alpar@9	2162 int q = csa->q;
alpar@9	2163 int k;
alpar@9	2164 /* xB[p] leaves the basis, xN[q] enters the basis */
alpar@9	2165 #ifdef GLP_DEBUG
alpar@9	2166 xassert(1 <= p && p <= m);
alpar@9	2167 xassert(1 <= q && q <= n);
alpar@9	2168 #endif
alpar@9	2169 /* xB[p] <-> xN[q] */
alpar@9	2170 k = head[p], head[p] = head[m+q], head[m+q] = k;
alpar@9	2171 #if 1 /* 06/IV-2009 */
alpar@9	2172 bind[head[p]] = p, bind[head[m+q]] = m + q;
alpar@9	2173 #endif
alpar@9	2174 if (type[k] == GLP_FX)
alpar@9	2175 stat[q] = GLP_NS;
alpar@9	2176 else if (delta > 0.0)
alpar@9	2177 {
alpar@9	2178 #ifdef GLP_DEBUG
alpar@9	2179 xassert(type[k] == GLP_LO \|\| type[k] == GLP_DB);
alpar@9	2180 #endif
alpar@9	2181 stat[q] = GLP_NL;
alpar@9	2182 }
alpar@9	2183 else /* delta < 0.0 */
alpar@9	2184 {
alpar@9	2185 #ifdef GLP_DEBUG
alpar@9	2186 xassert(type[k] == GLP_UP \|\| type[k] == GLP_DB);
alpar@9	2187 #endif
alpar@9	2188 stat[q] = GLP_NU;
alpar@9	2189 }
alpar@9	2190 return;
alpar@9	2191 }
alpar@9	2192
alpar@9	2193 /***********************************************************************
alpar@9	2194 * check_feas - check dual feasibility of basic solution
alpar@9	2195 *
alpar@9	2196 * If the current basic solution is dual feasible within a tolerance,
alpar@9	2197 * this routine returns zero, otherwise it returns non-zero. */
alpar@9	2198
alpar@9	2199 static int check_feas(struct csa *csa, double tol_dj)
alpar@9	2200 { int m = csa->m;
alpar@9	2201 int n = csa->n;
alpar@9	2202 char *orig_type = csa->orig_type;
alpar@9	2203 int *head = csa->head;
alpar@9	2204 double *cbar = csa->cbar;
alpar@9	2205 int j, k;
alpar@9	2206 for (j = 1; j <= n; j++)
alpar@9	2207 { k = head[m+j]; /* x[k] = xN[j] */
alpar@9	2208 #ifdef GLP_DEBUG
alpar@9	2209 xassert(1 <= k && k <= m+n);
alpar@9	2210 #endif
alpar@9	2211 if (cbar[j] < - tol_dj)
alpar@9	2212 if (orig_type[k] == GLP_LO \|\| orig_type[k] == GLP_FR)
alpar@9	2213 return 1;
alpar@9	2214 if (cbar[j] > + tol_dj)
alpar@9	2215 if (orig_type[k] == GLP_UP \|\| orig_type[k] == GLP_FR)
alpar@9	2216 return 1;
alpar@9	2217 }
alpar@9	2218 return 0;
alpar@9	2219 }
alpar@9	2220
alpar@9	2221 /***********************************************************************
alpar@9	2222 * set_aux_bnds - assign auxiliary bounds to variables
alpar@9	2223 *
alpar@9	2224 * This routine assigns auxiliary bounds to variables to construct an
alpar@9	2225 * LP problem solved on phase I. */
alpar@9	2226
alpar@9	2227 static void set_aux_bnds(struct csa *csa)
alpar@9	2228 { int m = csa->m;
alpar@9	2229 int n = csa->n;
alpar@9	2230 char *type = csa->type;
alpar@9	2231 double *lb = csa->lb;
alpar@9	2232 double *ub = csa->ub;
alpar@9	2233 char *orig_type = csa->orig_type;
alpar@9	2234 int *head = csa->head;
alpar@9	2235 char *stat = csa->stat;
alpar@9	2236 double *cbar = csa->cbar;
alpar@9	2237 int j, k;
alpar@9	2238 for (k = 1; k <= m+n; k++)
alpar@9	2239 { switch (orig_type[k])
alpar@9	2240 { case GLP_FR:
alpar@9	2241 #if 0
alpar@9	2242 type[k] = GLP_DB, lb[k] = -1.0, ub[k] = +1.0;
alpar@9	2243 #else
alpar@9	2244 /* to force free variables to enter the basis */
alpar@9	2245 type[k] = GLP_DB, lb[k] = -1e3, ub[k] = +1e3;
alpar@9	2246 #endif
alpar@9	2247 break;
alpar@9	2248 case GLP_LO:
alpar@9	2249 type[k] = GLP_DB, lb[k] = 0.0, ub[k] = +1.0;
alpar@9	2250 break;
alpar@9	2251 case GLP_UP:
alpar@9	2252 type[k] = GLP_DB, lb[k] = -1.0, ub[k] = 0.0;
alpar@9	2253 break;
alpar@9	2254 case GLP_DB:
alpar@9	2255 case GLP_FX:
alpar@9	2256 type[k] = GLP_FX, lb[k] = ub[k] = 0.0;
alpar@9	2257 break;
alpar@9	2258 default:
alpar@9	2259 xassert(orig_type != orig_type);
alpar@9	2260 }
alpar@9	2261 }
alpar@9	2262 for (j = 1; j <= n; j++)
alpar@9	2263 { k = head[m+j]; /* x[k] = xN[j] */
alpar@9	2264 #ifdef GLP_DEBUG
alpar@9	2265 xassert(1 <= k && k <= m+n);
alpar@9	2266 #endif
alpar@9	2267 if (type[k] == GLP_FX)
alpar@9	2268 stat[j] = GLP_NS;
alpar@9	2269 else if (cbar[j] >= 0.0)
alpar@9	2270 stat[j] = GLP_NL;
alpar@9	2271 else
alpar@9	2272 stat[j] = GLP_NU;
alpar@9	2273 }
alpar@9	2274 return;
alpar@9	2275 }
alpar@9	2276
alpar@9	2277 /***********************************************************************
alpar@9	2278 * set_orig_bnds - restore original bounds of variables
alpar@9	2279 *
alpar@9	2280 * This routine restores original types and bounds of variables and
alpar@9	2281 * determines statuses of non-basic variables assuming that the current
alpar@9	2282 * basis is dual feasible. */
alpar@9	2283
alpar@9	2284 static void set_orig_bnds(struct csa *csa)
alpar@9	2285 { int m = csa->m;
alpar@9	2286 int n = csa->n;
alpar@9	2287 char *type = csa->type;
alpar@9	2288 double *lb = csa->lb;
alpar@9	2289 double *ub = csa->ub;
alpar@9	2290 char *orig_type = csa->orig_type;
alpar@9	2291 double *orig_lb = csa->orig_lb;
alpar@9	2292 double *orig_ub = csa->orig_ub;
alpar@9	2293 int *head = csa->head;
alpar@9	2294 char *stat = csa->stat;
alpar@9	2295 double *cbar = csa->cbar;
alpar@9	2296 int j, k;
alpar@9	2297 memcpy(&type[1], &orig_type[1], (m+n) * sizeof(char));
alpar@9	2298 memcpy(&lb[1], &orig_lb[1], (m+n) * sizeof(double));
alpar@9	2299 memcpy(&ub[1], &orig_ub[1], (m+n) * sizeof(double));
alpar@9	2300 for (j = 1; j <= n; j++)
alpar@9	2301 { k = head[m+j]; /* x[k] = xN[j] */
alpar@9	2302 #ifdef GLP_DEBUG
alpar@9	2303 xassert(1 <= k && k <= m+n);
alpar@9	2304 #endif
alpar@9	2305 switch (type[k])
alpar@9	2306 { case GLP_FR:
alpar@9	2307 stat[j] = GLP_NF;
alpar@9	2308 break;
alpar@9	2309 case GLP_LO:
alpar@9	2310 stat[j] = GLP_NL;
alpar@9	2311 break;
alpar@9	2312 case GLP_UP:
alpar@9	2313 stat[j] = GLP_NU;
alpar@9	2314 break;
alpar@9	2315 case GLP_DB:
alpar@9	2316 if (cbar[j] >= +DBL_EPSILON)
alpar@9	2317 stat[j] = GLP_NL;
alpar@9	2318 else if (cbar[j] <= -DBL_EPSILON)
alpar@9	2319 stat[j] = GLP_NU;
alpar@9	2320 else if (fabs(lb[k]) <= fabs(ub[k]))
alpar@9	2321 stat[j] = GLP_NL;
alpar@9	2322 else
alpar@9	2323 stat[j] = GLP_NU;
alpar@9	2324 break;
alpar@9	2325 case GLP_FX:
alpar@9	2326 stat[j] = GLP_NS;
alpar@9	2327 break;
alpar@9	2328 default:
alpar@9	2329 xassert(type != type);
alpar@9	2330 }
alpar@9	2331 }
alpar@9	2332 return;
alpar@9	2333 }
alpar@9	2334
alpar@9	2335 /***********************************************************************
alpar@9	2336 * check_stab - check numerical stability of basic solution
alpar@9	2337 *
alpar@9	2338 * If the current basic solution is dual feasible within a tolerance,
alpar@9	2339 * this routine returns zero, otherwise it returns non-zero. */
alpar@9	2340
alpar@9	2341 static int check_stab(struct csa *csa, double tol_dj)
alpar@9	2342 { int n = csa->n;
alpar@9	2343 char *stat = csa->stat;
alpar@9	2344 double *cbar = csa->cbar;
alpar@9	2345 int j;
alpar@9	2346 for (j = 1; j <= n; j++)
alpar@9	2347 { if (cbar[j] < - tol_dj)
alpar@9	2348 if (stat[j] == GLP_NL \|\| stat[j] == GLP_NF) return 1;
alpar@9	2349 if (cbar[j] > + tol_dj)
alpar@9	2350 if (stat[j] == GLP_NU \|\| stat[j] == GLP_NF) return 1;
alpar@9	2351 }
alpar@9	2352 return 0;
alpar@9	2353 }
alpar@9	2354
alpar@9	2355 #if 1 /* copied from primal */
alpar@9	2356 /***********************************************************************
alpar@9	2357 * eval_obj - compute original objective function
alpar@9	2358 *
alpar@9	2359 * This routine computes the current value of the original objective
alpar@9	2360 * function. */
alpar@9	2361
alpar@9	2362 static double eval_obj(struct csa *csa)
alpar@9	2363 { int m = csa->m;
alpar@9	2364 int n = csa->n;
alpar@9	2365 double *obj = csa->obj;
alpar@9	2366 int *head = csa->head;
alpar@9	2367 double *bbar = csa->bbar;
alpar@9	2368 int i, j, k;
alpar@9	2369 double sum;
alpar@9	2370 sum = obj[0];
alpar@9	2371 /* walk through the list of basic variables */
alpar@9	2372 for (i = 1; i <= m; i++)
alpar@9	2373 { k = head[i]; /* x[k] = xB[i] */
alpar@9	2374 #ifdef GLP_DEBUG
alpar@9	2375 xassert(1 <= k && k <= m+n);
alpar@9	2376 #endif
alpar@9	2377 if (k > m)
alpar@9	2378 sum += obj[k-m] * bbar[i];
alpar@9	2379 }
alpar@9	2380 /* walk through the list of non-basic variables */
alpar@9	2381 for (j = 1; j <= n; j++)
alpar@9	2382 { k = head[m+j]; /* x[k] = xN[j] */
alpar@9	2383 #ifdef GLP_DEBUG
alpar@9	2384 xassert(1 <= k && k <= m+n);
alpar@9	2385 #endif
alpar@9	2386 if (k > m)
alpar@9	2387 sum += obj[k-m] * get_xN(csa, j);
alpar@9	2388 }
alpar@9	2389 return sum;
alpar@9	2390 }
alpar@9	2391 #endif
alpar@9	2392
alpar@9	2393 /***********************************************************************
alpar@9	2394 * display - display the search progress
alpar@9	2395 *
alpar@9	2396 * This routine displays some information about the search progress. */
alpar@9	2397
alpar@9	2398 static void display(struct csa csa, const glp_smcp parm, int spec)
alpar@9	2399 { int m = csa->m;
alpar@9	2400 int n = csa->n;
alpar@9	2401 double *coef = csa->coef;
alpar@9	2402 char *orig_type = csa->orig_type;
alpar@9	2403 int *head = csa->head;
alpar@9	2404 char *stat = csa->stat;
alpar@9	2405 int phase = csa->phase;
alpar@9	2406 double *bbar = csa->bbar;
alpar@9	2407 double *cbar = csa->cbar;
alpar@9	2408 int i, j, cnt;
alpar@9	2409 double sum;
alpar@9	2410 if (parm->msg_lev < GLP_MSG_ON) goto skip;
alpar@9	2411 if (parm->out_dly > 0 &&
alpar@9	2412 1000.0 * xdifftime(xtime(), csa->tm_beg) < parm->out_dly)
alpar@9	2413 goto skip;
alpar@9	2414 if (csa->it_cnt == csa->it_dpy) goto skip;
alpar@9	2415 if (!spec && csa->it_cnt % parm->out_frq != 0) goto skip;
alpar@9	2416 /* compute the sum of dual infeasibilities */
alpar@9	2417 sum = 0.0;
alpar@9	2418 if (phase == 1)
alpar@9	2419 { for (i = 1; i <= m; i++)
alpar@9	2420 sum -= coef[head[i]] * bbar[i];
alpar@9	2421 for (j = 1; j <= n; j++)
alpar@9	2422 sum -= coef[head[m+j]] * get_xN(csa, j);
alpar@9	2423 }
alpar@9	2424 else
alpar@9	2425 { for (j = 1; j <= n; j++)
alpar@9	2426 { if (cbar[j] < 0.0)
alpar@9	2427 if (stat[j] == GLP_NL \|\| stat[j] == GLP_NF)
alpar@9	2428 sum -= cbar[j];
alpar@9	2429 if (cbar[j] > 0.0)
alpar@9	2430 if (stat[j] == GLP_NU \|\| stat[j] == GLP_NF)
alpar@9	2431 sum += cbar[j];
alpar@9	2432 }
alpar@9	2433 }
alpar@9	2434 /* determine the number of basic fixed variables */
alpar@9	2435 cnt = 0;
alpar@9	2436 for (i = 1; i <= m; i++)
alpar@9	2437 if (orig_type[head[i]] == GLP_FX) cnt++;
alpar@9	2438 if (csa->phase == 1)
alpar@9	2439 xprintf(" %6d: %24s infeas = %10.3e (%d)\n",
alpar@9	2440 csa->it_cnt, "", sum, cnt);
alpar@9	2441 else
alpar@9	2442 xprintf("\|%6d: obj = %17.9e infeas = %10.3e (%d)\n",
alpar@9	2443 csa->it_cnt, eval_obj(csa), sum, cnt);
alpar@9	2444 csa->it_dpy = csa->it_cnt;
alpar@9	2445 skip: return;
alpar@9	2446 }
alpar@9	2447
alpar@9	2448 #if 1 /* copied from primal */
alpar@9	2449 /***********************************************************************
alpar@9	2450 * store_sol - store basic solution back to the problem object
alpar@9	2451 *
alpar@9	2452 * This routine stores basic solution components back to the problem
alpar@9	2453 * object. */
alpar@9	2454
alpar@9	2455 static void store_sol(struct csa csa, glp_prob lp, int p_stat,
alpar@9	2456 int d_stat, int ray)
alpar@9	2457 { int m = csa->m;
alpar@9	2458 int n = csa->n;
alpar@9	2459 double zeta = csa->zeta;
alpar@9	2460 int *head = csa->head;
alpar@9	2461 char *stat = csa->stat;
alpar@9	2462 double *bbar = csa->bbar;
alpar@9	2463 double *cbar = csa->cbar;
alpar@9	2464 int i, j, k;
alpar@9	2465 #ifdef GLP_DEBUG
alpar@9	2466 xassert(lp->m == m);
alpar@9	2467 xassert(lp->n == n);
alpar@9	2468 #endif
alpar@9	2469 /* basis factorization */
alpar@9	2470 #ifdef GLP_DEBUG
alpar@9	2471 xassert(!lp->valid && lp->bfd == NULL);
alpar@9	2472 xassert(csa->valid && csa->bfd != NULL);
alpar@9	2473 #endif
alpar@9	2474 lp->valid = 1, csa->valid = 0;
alpar@9	2475 lp->bfd = csa->bfd, csa->bfd = NULL;
alpar@9	2476 memcpy(&lp->head[1], &head[1], m * sizeof(int));
alpar@9	2477 /* basic solution status */
alpar@9	2478 lp->pbs_stat = p_stat;
alpar@9	2479 lp->dbs_stat = d_stat;
alpar@9	2480 /* objective function value */
alpar@9	2481 lp->obj_val = eval_obj(csa);
alpar@9	2482 /* simplex iteration count */
alpar@9	2483 lp->it_cnt = csa->it_cnt;
alpar@9	2484 /* unbounded ray */
alpar@9	2485 lp->some = ray;
alpar@9	2486 /* basic variables */
alpar@9	2487 for (i = 1; i <= m; i++)
alpar@9	2488 { k = head[i]; /* x[k] = xB[i] */
alpar@9	2489 #ifdef GLP_DEBUG
alpar@9	2490 xassert(1 <= k && k <= m+n);
alpar@9	2491 #endif
alpar@9	2492 if (k <= m)
alpar@9	2493 { GLPROW *row = lp->row[k];
alpar@9	2494 row->stat = GLP_BS;
alpar@9	2495 row->bind = i;
alpar@9	2496 row->prim = bbar[i] / row->rii;
alpar@9	2497 row->dual = 0.0;
alpar@9	2498 }
alpar@9	2499 else
alpar@9	2500 { GLPCOL *col = lp->col[k-m];
alpar@9	2501 col->stat = GLP_BS;
alpar@9	2502 col->bind = i;
alpar@9	2503 col->prim = bbar[i] * col->sjj;
alpar@9	2504 col->dual = 0.0;
alpar@9	2505 }
alpar@9	2506 }
alpar@9	2507 /* non-basic variables */
alpar@9	2508 for (j = 1; j <= n; j++)
alpar@9	2509 { k = head[m+j]; /* x[k] = xN[j] */
alpar@9	2510 #ifdef GLP_DEBUG
alpar@9	2511 xassert(1 <= k && k <= m+n);
alpar@9	2512 #endif
alpar@9	2513 if (k <= m)
alpar@9	2514 { GLPROW *row = lp->row[k];
alpar@9	2515 row->stat = stat[j];
alpar@9	2516 row->bind = 0;
alpar@9	2517 #if 0
alpar@9	2518 row->prim = get_xN(csa, j) / row->rii;
alpar@9	2519 #else
alpar@9	2520 switch (stat[j])
alpar@9	2521 { case GLP_NL:
alpar@9	2522 row->prim = row->lb; break;
alpar@9	2523 case GLP_NU:
alpar@9	2524 row->prim = row->ub; break;
alpar@9	2525 case GLP_NF:
alpar@9	2526 row->prim = 0.0; break;
alpar@9	2527 case GLP_NS:
alpar@9	2528 row->prim = row->lb; break;
alpar@9	2529 default:
alpar@9	2530 xassert(stat != stat);
alpar@9	2531 }
alpar@9	2532 #endif
alpar@9	2533 row->dual = (cbar[j] * row->rii) / zeta;
alpar@9	2534 }
alpar@9	2535 else
alpar@9	2536 { GLPCOL *col = lp->col[k-m];
alpar@9	2537 col->stat = stat[j];
alpar@9	2538 col->bind = 0;
alpar@9	2539 #if 0
alpar@9	2540 col->prim = get_xN(csa, j) * col->sjj;
alpar@9	2541 #else
alpar@9	2542 switch (stat[j])
alpar@9	2543 { case GLP_NL:
alpar@9	2544 col->prim = col->lb; break;
alpar@9	2545 case GLP_NU:
alpar@9	2546 col->prim = col->ub; break;
alpar@9	2547 case GLP_NF:
alpar@9	2548 col->prim = 0.0; break;
alpar@9	2549 case GLP_NS:
alpar@9	2550 col->prim = col->lb; break;
alpar@9	2551 default:
alpar@9	2552 xassert(stat != stat);
alpar@9	2553 }
alpar@9	2554 #endif
alpar@9	2555 col->dual = (cbar[j] / col->sjj) / zeta;
alpar@9	2556 }
alpar@9	2557 }
alpar@9	2558 return;
alpar@9	2559 }
alpar@9	2560 #endif
alpar@9	2561
alpar@9	2562 /***********************************************************************
alpar@9	2563 * free_csa - deallocate common storage area
alpar@9	2564 *
alpar@9	2565 * This routine frees all the memory allocated to arrays in the common
alpar@9	2566 * storage area (CSA). */
alpar@9	2567
alpar@9	2568 static void free_csa(struct csa *csa)
alpar@9	2569 { xfree(csa->type);
alpar@9	2570 xfree(csa->lb);
alpar@9	2571 xfree(csa->ub);
alpar@9	2572 xfree(csa->coef);
alpar@9	2573 xfree(csa->orig_type);
alpar@9	2574 xfree(csa->orig_lb);
alpar@9	2575 xfree(csa->orig_ub);
alpar@9	2576 xfree(csa->obj);
alpar@9	2577 xfree(csa->A_ptr);
alpar@9	2578 xfree(csa->A_ind);
alpar@9	2579 xfree(csa->A_val);
alpar@9	2580 #if 1 /* 06/IV-2009 */
alpar@9	2581 xfree(csa->AT_ptr);
alpar@9	2582 xfree(csa->AT_ind);
alpar@9	2583 xfree(csa->AT_val);
alpar@9	2584 #endif
alpar@9	2585 xfree(csa->head);
alpar@9	2586 #if 1 /* 06/IV-2009 */
alpar@9	2587 xfree(csa->bind);
alpar@9	2588 #endif
alpar@9	2589 xfree(csa->stat);
alpar@9	2590 #if 0 /* 06/IV-2009 */
alpar@9	2591 xfree(csa->N_ptr);
alpar@9	2592 xfree(csa->N_len);
alpar@9	2593 xfree(csa->N_ind);
alpar@9	2594 xfree(csa->N_val);
alpar@9	2595 #endif
alpar@9	2596 xfree(csa->bbar);
alpar@9	2597 xfree(csa->cbar);
alpar@9	2598 xfree(csa->refsp);
alpar@9	2599 xfree(csa->gamma);
alpar@9	2600 xfree(csa->trow_ind);
alpar@9	2601 xfree(csa->trow_vec);
alpar@9	2602 #ifdef GLP_LONG_STEP /* 07/IV-2009 */
alpar@9	2603 xfree(csa->bkpt);
alpar@9	2604 #endif
alpar@9	2605 xfree(csa->tcol_ind);
alpar@9	2606 xfree(csa->tcol_vec);
alpar@9	2607 xfree(csa->work1);
alpar@9	2608 xfree(csa->work2);
alpar@9	2609 xfree(csa->work3);
alpar@9	2610 xfree(csa->work4);
alpar@9	2611 xfree(csa);
alpar@9	2612 return;
alpar@9	2613 }
alpar@9	2614
alpar@9	2615 /***********************************************************************
alpar@9	2616 * spx_dual - core LP solver based on the dual simplex method
alpar@9	2617 *
alpar@9	2618 * SYNOPSIS
alpar@9	2619 *
alpar@9	2620 * #include "glpspx.h"
alpar@9	2621 * int spx_dual(glp_prob lp, const glp_smcp parm);
alpar@9	2622 *
alpar@9	2623 * DESCRIPTION
alpar@9	2624 *
alpar@9	2625 * The routine spx_dual is a core LP solver based on the two-phase dual
alpar@9	2626 * simplex method.
alpar@9	2627 *
alpar@9	2628 * RETURNS
alpar@9	2629 *
alpar@9	2630 * 0 LP instance has been successfully solved.
alpar@9	2631 *
alpar@9	2632 * GLP_EOBJLL
alpar@9	2633 * Objective lower limit has been reached (maximization).
alpar@9	2634 *
alpar@9	2635 * GLP_EOBJUL
alpar@9	2636 * Objective upper limit has been reached (minimization).
alpar@9	2637 *
alpar@9	2638 * GLP_EITLIM
alpar@9	2639 * Iteration limit has been exhausted.
alpar@9	2640 *
alpar@9	2641 * GLP_ETMLIM
alpar@9	2642 * Time limit has been exhausted.
alpar@9	2643 *
alpar@9	2644 * GLP_EFAIL
alpar@9	2645 * The solver failed to solve LP instance. */
alpar@9	2646
alpar@9	2647 int spx_dual(glp_prob lp, const glp_smcp parm)
alpar@9	2648 { struct csa *csa;
alpar@9	2649 int binv_st = 2;
alpar@9	2650 /* status of basis matrix factorization:
alpar@9	2651 0 - invalid; 1 - just computed; 2 - updated */
alpar@9	2652 int bbar_st = 0;
alpar@9	2653 /* status of primal values of basic variables:
alpar@9	2654 0 - invalid; 1 - just computed; 2 - updated */
alpar@9	2655 int cbar_st = 0;
alpar@9	2656 /* status of reduced costs of non-basic variables:
alpar@9	2657 0 - invalid; 1 - just computed; 2 - updated */
alpar@9	2658 int rigorous = 0;
alpar@9	2659 /* rigorous mode flag; this flag is used to enable iterative
alpar@9	2660 refinement on computing pivot rows and columns of the simplex
alpar@9	2661 table */
alpar@9	2662 int check = 0;
alpar@9	2663 int p_stat, d_stat, ret;
alpar@9	2664 /* allocate and initialize the common storage area */
alpar@9	2665 csa = alloc_csa(lp);
alpar@9	2666 init_csa(csa, lp);
alpar@9	2667 if (parm->msg_lev >= GLP_MSG_DBG)
alpar@9	2668 xprintf("Objective scale factor = %g\n", csa->zeta);
alpar@9	2669 loop: /* main loop starts here */
alpar@9	2670 /* compute factorization of the basis matrix */
alpar@9	2671 if (binv_st == 0)
alpar@9	2672 { ret = invert_B(csa);
alpar@9	2673 if (ret != 0)
alpar@9	2674 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9	2675 { xprintf("Error: unable to factorize the basis matrix (%d"
alpar@9	2676 ")\n", ret);
alpar@9	2677 xprintf("Sorry, basis recovery procedure not implemented"
alpar@9	2678 " yet\n");
alpar@9	2679 }
alpar@9	2680 xassert(!lp->valid && lp->bfd == NULL);
alpar@9	2681 lp->bfd = csa->bfd, csa->bfd = NULL;
alpar@9	2682 lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
alpar@9	2683 lp->obj_val = 0.0;
alpar@9	2684 lp->it_cnt = csa->it_cnt;
alpar@9	2685 lp->some = 0;
alpar@9	2686 ret = GLP_EFAIL;
alpar@9	2687 goto done;
alpar@9	2688 }
alpar@9	2689 csa->valid = 1;
alpar@9	2690 binv_st = 1; /* just computed */
alpar@9	2691 /* invalidate basic solution components */
alpar@9	2692 bbar_st = cbar_st = 0;
alpar@9	2693 }
alpar@9	2694 /* compute reduced costs of non-basic variables */
alpar@9	2695 if (cbar_st == 0)
alpar@9	2696 { eval_cbar(csa);
alpar@9	2697 cbar_st = 1; /* just computed */
alpar@9	2698 /* determine the search phase, if not determined yet */
alpar@9	2699 if (csa->phase == 0)
alpar@9	2700 { if (check_feas(csa, 0.90 * parm->tol_dj) != 0)
alpar@9	2701 { /* current basic solution is dual infeasible */
alpar@9	2702 /* start searching for dual feasible solution */
alpar@9	2703 csa->phase = 1;
alpar@9	2704 set_aux_bnds(csa);
alpar@9	2705 }
alpar@9	2706 else
alpar@9	2707 { /* current basic solution is dual feasible */
alpar@9	2708 /* start searching for optimal solution */
alpar@9	2709 csa->phase = 2;
alpar@9	2710 set_orig_bnds(csa);
alpar@9	2711 }
alpar@9	2712 xassert(check_stab(csa, parm->tol_dj) == 0);
alpar@9	2713 /* some non-basic double-bounded variables might become
alpar@9	2714 fixed (on phase I) or vice versa (on phase II) */
alpar@9	2715 #if 0 /* 06/IV-2009 */
alpar@9	2716 build_N(csa);
alpar@9	2717 #endif
alpar@9	2718 csa->refct = 0;
alpar@9	2719 /* bounds of non-basic variables have been changed, so
alpar@9	2720 invalidate primal values */
alpar@9	2721 bbar_st = 0;
alpar@9	2722 }
alpar@9	2723 /* make sure that the current basic solution remains dual
alpar@9	2724 feasible */
alpar@9	2725 if (check_stab(csa, parm->tol_dj) != 0)
alpar@9	2726 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9	2727 xprintf("Warning: numerical instability (dual simplex, p"
alpar@9	2728 "hase %s)\n", csa->phase == 1 ? "I" : "II");
alpar@9	2729 #if 1
alpar@9	2730 if (parm->meth == GLP_DUALP)
alpar@9	2731 { store_sol(csa, lp, GLP_UNDEF, GLP_UNDEF, 0);
alpar@9	2732 ret = GLP_EFAIL;
alpar@9	2733 goto done;
alpar@9	2734 }
alpar@9	2735 #endif
alpar@9	2736 /* restart the search */
alpar@9	2737 csa->phase = 0;
alpar@9	2738 binv_st = 0;
alpar@9	2739 rigorous = 5;
alpar@9	2740 goto loop;
alpar@9	2741 }
alpar@9	2742 }
alpar@9	2743 xassert(csa->phase == 1 \|\| csa->phase == 2);
alpar@9	2744 /* on phase I we do not need to wait until the current basic
alpar@9	2745 solution becomes primal feasible; it is sufficient to make
alpar@9	2746 sure that all reduced costs have correct signs */
alpar@9	2747 if (csa->phase == 1 && check_feas(csa, parm->tol_dj) == 0)
alpar@9	2748 { /* the current basis is dual feasible; switch to phase II */
alpar@9	2749 display(csa, parm, 1);
alpar@9	2750 csa->phase = 2;
alpar@9	2751 if (cbar_st != 1)
alpar@9	2752 { eval_cbar(csa);
alpar@9	2753 cbar_st = 1;
alpar@9	2754 }
alpar@9	2755 set_orig_bnds(csa);
alpar@9	2756 #if 0 /* 06/IV-2009 */
alpar@9	2757 build_N(csa);
alpar@9	2758 #endif
alpar@9	2759 csa->refct = 0;
alpar@9	2760 bbar_st = 0;
alpar@9	2761 }
alpar@9	2762 /* compute primal values of basic variables */
alpar@9	2763 if (bbar_st == 0)
alpar@9	2764 { eval_bbar(csa);
alpar@9	2765 if (csa->phase == 2)
alpar@9	2766 csa->bbar[0] = eval_obj(csa);
alpar@9	2767 bbar_st = 1; /* just computed */
alpar@9	2768 }
alpar@9	2769 /* redefine the reference space, if required */
alpar@9	2770 switch (parm->pricing)
alpar@9	2771 { case GLP_PT_STD:
alpar@9	2772 break;
alpar@9	2773 case GLP_PT_PSE:
alpar@9	2774 if (csa->refct == 0) reset_refsp(csa);
alpar@9	2775 break;
alpar@9	2776 default:
alpar@9	2777 xassert(parm != parm);
alpar@9	2778 }
alpar@9	2779 /* at this point the basis factorization and all basic solution
alpar@9	2780 components are valid */
alpar@9	2781 xassert(binv_st && bbar_st && cbar_st);
alpar@9	2782 /* check accuracy of current basic solution components (only for
alpar@9	2783 debugging) */
alpar@9	2784 if (check)
alpar@9	2785 { double e_bbar = err_in_bbar(csa);
alpar@9	2786 double e_cbar = err_in_cbar(csa);
alpar@9	2787 double e_gamma =
alpar@9	2788 (parm->pricing == GLP_PT_PSE ? err_in_gamma(csa) : 0.0);
alpar@9	2789 xprintf("e_bbar = %10.3e; e_cbar = %10.3e; e_gamma = %10.3e\n",
alpar@9	2790 e_bbar, e_cbar, e_gamma);
alpar@9	2791 xassert(e_bbar <= 1e-5 && e_cbar <= 1e-5 && e_gamma <= 1e-3);
alpar@9	2792 }
alpar@9	2793 /* if the objective has to be maximized, check if it has reached
alpar@9	2794 its lower limit */
alpar@9	2795 if (csa->phase == 2 && csa->zeta < 0.0 &&
alpar@9	2796 parm->obj_ll > -DBL_MAX && csa->bbar[0] <= parm->obj_ll)
alpar@9	2797 { if (bbar_st != 1 \|\| cbar_st != 1)
alpar@9	2798 { if (bbar_st != 1) bbar_st = 0;
alpar@9	2799 if (cbar_st != 1) cbar_st = 0;
alpar@9	2800 goto loop;
alpar@9	2801 }
alpar@9	2802 display(csa, parm, 1);
alpar@9	2803 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9	2804 xprintf("OBJECTIVE LOWER LIMIT REACHED; SEARCH TERMINATED\n"
alpar@9	2805 );
alpar@9	2806 store_sol(csa, lp, GLP_INFEAS, GLP_FEAS, 0);
alpar@9	2807 ret = GLP_EOBJLL;
alpar@9	2808 goto done;
alpar@9	2809 }
alpar@9	2810 /* if the objective has to be minimized, check if it has reached
alpar@9	2811 its upper limit */
alpar@9	2812 if (csa->phase == 2 && csa->zeta > 0.0 &&
alpar@9	2813 parm->obj_ul < +DBL_MAX && csa->bbar[0] >= parm->obj_ul)
alpar@9	2814 { if (bbar_st != 1 \|\| cbar_st != 1)
alpar@9	2815 { if (bbar_st != 1) bbar_st = 0;
alpar@9	2816 if (cbar_st != 1) cbar_st = 0;
alpar@9	2817 goto loop;
alpar@9	2818 }
alpar@9	2819 display(csa, parm, 1);
alpar@9	2820 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9	2821 xprintf("OBJECTIVE UPPER LIMIT REACHED; SEARCH TERMINATED\n"
alpar@9	2822 );
alpar@9	2823 store_sol(csa, lp, GLP_INFEAS, GLP_FEAS, 0);
alpar@9	2824 ret = GLP_EOBJUL;
alpar@9	2825 goto done;
alpar@9	2826 }
alpar@9	2827 /* check if the iteration limit has been exhausted */
alpar@9	2828 if (parm->it_lim < INT_MAX &&
alpar@9	2829 csa->it_cnt - csa->it_beg >= parm->it_lim)
alpar@9	2830 { if (csa->phase == 2 && bbar_st != 1 \|\| cbar_st != 1)
alpar@9	2831 { if (csa->phase == 2 && bbar_st != 1) bbar_st = 0;
alpar@9	2832 if (cbar_st != 1) cbar_st = 0;
alpar@9	2833 goto loop;
alpar@9	2834 }
alpar@9	2835 display(csa, parm, 1);
alpar@9	2836 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9	2837 xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n");
alpar@9	2838 switch (csa->phase)
alpar@9	2839 { case 1:
alpar@9	2840 d_stat = GLP_INFEAS;
alpar@9	2841 set_orig_bnds(csa);
alpar@9	2842 eval_bbar(csa);
alpar@9	2843 break;
alpar@9	2844 case 2:
alpar@9	2845 d_stat = GLP_FEAS;
alpar@9	2846 break;
alpar@9	2847 default:
alpar@9	2848 xassert(csa != csa);
alpar@9	2849 }
alpar@9	2850 store_sol(csa, lp, GLP_INFEAS, d_stat, 0);
alpar@9	2851 ret = GLP_EITLIM;
alpar@9	2852 goto done;
alpar@9	2853 }
alpar@9	2854 /* check if the time limit has been exhausted */
alpar@9	2855 if (parm->tm_lim < INT_MAX &&
alpar@9	2856 1000.0 * xdifftime(xtime(), csa->tm_beg) >= parm->tm_lim)
alpar@9	2857 { if (csa->phase == 2 && bbar_st != 1 \|\| cbar_st != 1)
alpar@9	2858 { if (csa->phase == 2 && bbar_st != 1) bbar_st = 0;
alpar@9	2859 if (cbar_st != 1) cbar_st = 0;
alpar@9	2860 goto loop;
alpar@9	2861 }
alpar@9	2862 display(csa, parm, 1);
alpar@9	2863 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9	2864 xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
alpar@9	2865 switch (csa->phase)
alpar@9	2866 { case 1:
alpar@9	2867 d_stat = GLP_INFEAS;
alpar@9	2868 set_orig_bnds(csa);
alpar@9	2869 eval_bbar(csa);
alpar@9	2870 break;
alpar@9	2871 case 2:
alpar@9	2872 d_stat = GLP_FEAS;
alpar@9	2873 break;
alpar@9	2874 default:
alpar@9	2875 xassert(csa != csa);
alpar@9	2876 }
alpar@9	2877 store_sol(csa, lp, GLP_INFEAS, d_stat, 0);
alpar@9	2878 ret = GLP_ETMLIM;
alpar@9	2879 goto done;
alpar@9	2880 }
alpar@9	2881 /* display the search progress */
alpar@9	2882 display(csa, parm, 0);
alpar@9	2883 /* choose basic variable xB[p] */
alpar@9	2884 chuzr(csa, parm->tol_bnd);
alpar@9	2885 if (csa->p == 0)
alpar@9	2886 { if (bbar_st != 1 \|\| cbar_st != 1)
alpar@9	2887 { if (bbar_st != 1) bbar_st = 0;
alpar@9	2888 if (cbar_st != 1) cbar_st = 0;
alpar@9	2889 goto loop;
alpar@9	2890 }
alpar@9	2891 display(csa, parm, 1);
alpar@9	2892 switch (csa->phase)
alpar@9	2893 { case 1:
alpar@9	2894 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9	2895 xprintf("PROBLEM HAS NO DUAL FEASIBLE SOLUTION\n");
alpar@9	2896 set_orig_bnds(csa);
alpar@9	2897 eval_bbar(csa);
alpar@9	2898 p_stat = GLP_INFEAS, d_stat = GLP_NOFEAS;
alpar@9	2899 break;
alpar@9	2900 case 2:
alpar@9	2901 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9	2902 xprintf("OPTIMAL SOLUTION FOUND\n");
alpar@9	2903 p_stat = d_stat = GLP_FEAS;
alpar@9	2904 break;
alpar@9	2905 default:
alpar@9	2906 xassert(csa != csa);
alpar@9	2907 }
alpar@9	2908 store_sol(csa, lp, p_stat, d_stat, 0);
alpar@9	2909 ret = 0;
alpar@9	2910 goto done;
alpar@9	2911 }
alpar@9	2912 /* compute pivot row of the simplex table */
alpar@9	2913 { double *rho = csa->work4;
alpar@9	2914 eval_rho(csa, rho);
alpar@9	2915 if (rigorous) refine_rho(csa, rho);
alpar@9	2916 eval_trow(csa, rho);
alpar@9	2917 sort_trow(csa, parm->tol_bnd);
alpar@9	2918 }
alpar@9	2919 /* unlike primal simplex there is no need to check accuracy of
alpar@9	2920 the primal value of xB[p] (which might be computed using the
alpar@9	2921 pivot row), since bbar is a result of FTRAN */
alpar@9	2922 #ifdef GLP_LONG_STEP /* 07/IV-2009 */
alpar@9	2923 long_step(csa);
alpar@9	2924 if (csa->nbps > 0)
alpar@9	2925 { csa->q = csa->bkpt[csa->nbps].j;
alpar@9	2926 if (csa->delta > 0.0)
alpar@9	2927 csa->new_dq = + csa->bkpt[csa->nbps].t;
alpar@9	2928 else
alpar@9	2929 csa->new_dq = - csa->bkpt[csa->nbps].t;
alpar@9	2930 }
alpar@9	2931 else
alpar@9	2932 #endif
alpar@9	2933 /* choose non-basic variable xN[q] */
alpar@9	2934 switch (parm->r_test)
alpar@9	2935 { case GLP_RT_STD:
alpar@9	2936 chuzc(csa, 0.0);
alpar@9	2937 break;
alpar@9	2938 case GLP_RT_HAR:
alpar@9	2939 chuzc(csa, 0.30 * parm->tol_dj);
alpar@9	2940 break;
alpar@9	2941 default:
alpar@9	2942 xassert(parm != parm);
alpar@9	2943 }
alpar@9	2944 if (csa->q == 0)
alpar@9	2945 { if (bbar_st != 1 \|\| cbar_st != 1 \|\| !rigorous)
alpar@9	2946 { if (bbar_st != 1) bbar_st = 0;
alpar@9	2947 if (cbar_st != 1) cbar_st = 0;
alpar@9	2948 rigorous = 1;
alpar@9	2949 goto loop;
alpar@9	2950 }
alpar@9	2951 display(csa, parm, 1);
alpar@9	2952 switch (csa->phase)
alpar@9	2953 { case 1:
alpar@9	2954 if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9	2955 xprintf("Error: unable to choose basic variable on ph"
alpar@9	2956 "ase I\n");
alpar@9	2957 xassert(!lp->valid && lp->bfd == NULL);
alpar@9	2958 lp->bfd = csa->bfd, csa->bfd = NULL;
alpar@9	2959 lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
alpar@9	2960 lp->obj_val = 0.0;
alpar@9	2961 lp->it_cnt = csa->it_cnt;
alpar@9	2962 lp->some = 0;
alpar@9	2963 ret = GLP_EFAIL;
alpar@9	2964 break;
alpar@9	2965 case 2:
alpar@9	2966 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9	2967 xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
alpar@9	2968 store_sol(csa, lp, GLP_NOFEAS, GLP_FEAS,
alpar@9	2969 csa->head[csa->p]);
alpar@9	2970 ret = 0;
alpar@9	2971 break;
alpar@9	2972 default:
alpar@9	2973 xassert(csa != csa);
alpar@9	2974 }
alpar@9	2975 goto done;
alpar@9	2976 }
alpar@9	2977 /* check if the pivot element is acceptable */
alpar@9	2978 { double piv = csa->trow_vec[csa->q];
alpar@9	2979 double eps = 1e-5 * (1.0 + 0.01 * csa->trow_max);
alpar@9	2980 if (fabs(piv) < eps)
alpar@9	2981 { if (parm->msg_lev >= GLP_MSG_DBG)
alpar@9	2982 xprintf("piv = %.12g; eps = %g\n", piv, eps);
alpar@9	2983 if (!rigorous)
alpar@9	2984 { rigorous = 5;
alpar@9	2985 goto loop;
alpar@9	2986 }
alpar@9	2987 }
alpar@9	2988 }
alpar@9	2989 /* now xN[q] and xB[p] have been chosen anyhow */
alpar@9	2990 /* compute pivot column of the simplex table */
alpar@9	2991 eval_tcol(csa);
alpar@9	2992 if (rigorous) refine_tcol(csa);
alpar@9	2993 /* accuracy check based on the pivot element */
alpar@9	2994 { double piv1 = csa->tcol_vec[csa->p]; /* more accurate */
alpar@9	2995 double piv2 = csa->trow_vec[csa->q]; /* less accurate */
alpar@9	2996 xassert(piv1 != 0.0);
alpar@9	2997 if (fabs(piv1 - piv2) > 1e-8 * (1.0 + fabs(piv1)) \|\|
alpar@9	2998 !(piv1 > 0.0 && piv2 > 0.0 \|\| piv1 < 0.0 && piv2 < 0.0))
alpar@9	2999 { if (parm->msg_lev >= GLP_MSG_DBG)
alpar@9	3000 xprintf("piv1 = %.12g; piv2 = %.12g\n", piv1, piv2);
alpar@9	3001 if (binv_st != 1 \|\| !rigorous)
alpar@9	3002 { if (binv_st != 1) binv_st = 0;
alpar@9	3003 rigorous = 5;
alpar@9	3004 goto loop;
alpar@9	3005 }
alpar@9	3006 /* (not a good idea; should be revised later) */
alpar@9	3007 if (csa->tcol_vec[csa->p] == 0.0)
alpar@9	3008 { csa->tcol_nnz++;
alpar@9	3009 xassert(csa->tcol_nnz <= csa->m);
alpar@9	3010 csa->tcol_ind[csa->tcol_nnz] = csa->p;
alpar@9	3011 }
alpar@9	3012 csa->tcol_vec[csa->p] = piv2;
alpar@9	3013 }
alpar@9	3014 }
alpar@9	3015 /* update primal values of basic variables */
alpar@9	3016 #ifdef GLP_LONG_STEP /* 07/IV-2009 */
alpar@9	3017 if (csa->nbps > 0)
alpar@9	3018 { int kk, j, k;
alpar@9	3019 for (kk = 1; kk < csa->nbps; kk++)
alpar@9	3020 { if (csa->bkpt[kk].t >= csa->bkpt[csa->nbps].t) continue;
alpar@9	3021 j = csa->bkpt[kk].j;
alpar@9	3022 k = csa->head[csa->m + j];
alpar@9	3023 xassert(csa->type[k] == GLP_DB);
alpar@9	3024 if (csa->stat[j] == GLP_NL)
alpar@9	3025 csa->stat[j] = GLP_NU;
alpar@9	3026 else
alpar@9	3027 csa->stat[j] = GLP_NL;
alpar@9	3028 }
alpar@9	3029 }
alpar@9	3030 bbar_st = 0;
alpar@9	3031 #else
alpar@9	3032 update_bbar(csa);
alpar@9	3033 if (csa->phase == 2)
alpar@9	3034 csa->bbar[0] += (csa->cbar[csa->q] / csa->zeta) *
alpar@9	3035 (csa->delta / csa->tcol_vec[csa->p]);
alpar@9	3036 bbar_st = 2; /* updated */
alpar@9	3037 #endif
alpar@9	3038 /* update reduced costs of non-basic variables */
alpar@9	3039 update_cbar(csa);
alpar@9	3040 cbar_st = 2; /* updated */
alpar@9	3041 /* update steepest edge coefficients */
alpar@9	3042 switch (parm->pricing)
alpar@9	3043 { case GLP_PT_STD:
alpar@9	3044 break;
alpar@9	3045 case GLP_PT_PSE:
alpar@9	3046 if (csa->refct > 0) update_gamma(csa);
alpar@9	3047 break;
alpar@9	3048 default:
alpar@9	3049 xassert(parm != parm);
alpar@9	3050 }
alpar@9	3051 /* update factorization of the basis matrix */
alpar@9	3052 ret = update_B(csa, csa->p, csa->head[csa->m+csa->q]);
alpar@9	3053 if (ret == 0)
alpar@9	3054 binv_st = 2; /* updated */
alpar@9	3055 else
alpar@9	3056 { csa->valid = 0;
alpar@9	3057 binv_st = 0; /* invalid */
alpar@9	3058 }
alpar@9	3059 #if 0 /* 06/IV-2009 */
alpar@9	3060 /* update matrix N */
alpar@9	3061 del_N_col(csa, csa->q, csa->head[csa->m+csa->q]);
alpar@9	3062 if (csa->type[csa->head[csa->p]] != GLP_FX)
alpar@9	3063 add_N_col(csa, csa->q, csa->head[csa->p]);
alpar@9	3064 #endif
alpar@9	3065 /* change the basis header */
alpar@9	3066 change_basis(csa);
alpar@9	3067 /* iteration complete */
alpar@9	3068 csa->it_cnt++;
alpar@9	3069 if (rigorous > 0) rigorous--;
alpar@9	3070 goto loop;
alpar@9	3071 done: /* deallocate the common storage area */
alpar@9	3072 free_csa(csa);
alpar@9	3073 /* return to the calling program */
alpar@9	3074 return ret;
alpar@9	3075 }
alpar@9	3076
alpar@9	3077 /* eof */