lemon-project-template-glpk

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

Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
rev   line source
alpar@9 1 /* 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 */