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

lemon-project-template-glpk

annotate deps/glpk/src/glpspx01.c @ 9:33de93886c88

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

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