1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpapi12.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,2219 @@
1.4 +/* glpapi12.c (basis factorization and simplex tableau routines) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpapi.h"
1.29 +
1.30 +/***********************************************************************
1.31 +* NAME
1.32 +*
1.33 +* glp_bf_exists - check if the basis factorization exists
1.34 +*
1.35 +* SYNOPSIS
1.36 +*
1.37 +* int glp_bf_exists(glp_prob *lp);
1.38 +*
1.39 +* RETURNS
1.40 +*
1.41 +* If the basis factorization for the current basis associated with
1.42 +* the specified problem object exists and therefore is available for
1.43 +* computations, the routine glp_bf_exists returns non-zero. Otherwise
1.44 +* the routine returns zero. */
1.45 +
1.46 +int glp_bf_exists(glp_prob *lp)
1.47 +{ int ret;
1.48 + ret = (lp->m == 0 || lp->valid);
1.49 + return ret;
1.50 +}
1.51 +
1.52 +/***********************************************************************
1.53 +* NAME
1.54 +*
1.55 +* glp_factorize - compute the basis factorization
1.56 +*
1.57 +* SYNOPSIS
1.58 +*
1.59 +* int glp_factorize(glp_prob *lp);
1.60 +*
1.61 +* DESCRIPTION
1.62 +*
1.63 +* The routine glp_factorize computes the basis factorization for the
1.64 +* current basis associated with the specified problem object.
1.65 +*
1.66 +* RETURNS
1.67 +*
1.68 +* 0 The basis factorization has been successfully computed.
1.69 +*
1.70 +* GLP_EBADB
1.71 +* The basis matrix is invalid, i.e. the number of basic (auxiliary
1.72 +* and structural) variables differs from the number of rows in the
1.73 +* problem object.
1.74 +*
1.75 +* GLP_ESING
1.76 +* The basis matrix is singular within the working precision.
1.77 +*
1.78 +* GLP_ECOND
1.79 +* The basis matrix is ill-conditioned. */
1.80 +
1.81 +static int b_col(void *info, int j, int ind[], double val[])
1.82 +{ glp_prob *lp = info;
1.83 + int m = lp->m;
1.84 + GLPAIJ *aij;
1.85 + int k, len;
1.86 + xassert(1 <= j && j <= m);
1.87 + /* determine the ordinal number of basic auxiliary or structural
1.88 + variable x[k] corresponding to basic variable xB[j] */
1.89 + k = lp->head[j];
1.90 + /* build j-th column of the basic matrix, which is k-th column of
1.91 + the scaled augmented matrix (I | -R*A*S) */
1.92 + if (k <= m)
1.93 + { /* x[k] is auxiliary variable */
1.94 + len = 1;
1.95 + ind[1] = k;
1.96 + val[1] = 1.0;
1.97 + }
1.98 + else
1.99 + { /* x[k] is structural variable */
1.100 + len = 0;
1.101 + for (aij = lp->col[k-m]->ptr; aij != NULL; aij = aij->c_next)
1.102 + { len++;
1.103 + ind[len] = aij->row->i;
1.104 + val[len] = - aij->row->rii * aij->val * aij->col->sjj;
1.105 + }
1.106 + }
1.107 + return len;
1.108 +}
1.109 +
1.110 +static void copy_bfcp(glp_prob *lp);
1.111 +
1.112 +int glp_factorize(glp_prob *lp)
1.113 +{ int m = lp->m;
1.114 + int n = lp->n;
1.115 + GLPROW **row = lp->row;
1.116 + GLPCOL **col = lp->col;
1.117 + int *head = lp->head;
1.118 + int j, k, stat, ret;
1.119 + /* invalidate the basis factorization */
1.120 + lp->valid = 0;
1.121 + /* build the basis header */
1.122 + j = 0;
1.123 + for (k = 1; k <= m+n; k++)
1.124 + { if (k <= m)
1.125 + { stat = row[k]->stat;
1.126 + row[k]->bind = 0;
1.127 + }
1.128 + else
1.129 + { stat = col[k-m]->stat;
1.130 + col[k-m]->bind = 0;
1.131 + }
1.132 + if (stat == GLP_BS)
1.133 + { j++;
1.134 + if (j > m)
1.135 + { /* too many basic variables */
1.136 + ret = GLP_EBADB;
1.137 + goto fini;
1.138 + }
1.139 + head[j] = k;
1.140 + if (k <= m)
1.141 + row[k]->bind = j;
1.142 + else
1.143 + col[k-m]->bind = j;
1.144 + }
1.145 + }
1.146 + if (j < m)
1.147 + { /* too few basic variables */
1.148 + ret = GLP_EBADB;
1.149 + goto fini;
1.150 + }
1.151 + /* try to factorize the basis matrix */
1.152 + if (m > 0)
1.153 + { if (lp->bfd == NULL)
1.154 + { lp->bfd = bfd_create_it();
1.155 + copy_bfcp(lp);
1.156 + }
1.157 + switch (bfd_factorize(lp->bfd, m, lp->head, b_col, lp))
1.158 + { case 0:
1.159 + /* ok */
1.160 + break;
1.161 + case BFD_ESING:
1.162 + /* singular matrix */
1.163 + ret = GLP_ESING;
1.164 + goto fini;
1.165 + case BFD_ECOND:
1.166 + /* ill-conditioned matrix */
1.167 + ret = GLP_ECOND;
1.168 + goto fini;
1.169 + default:
1.170 + xassert(lp != lp);
1.171 + }
1.172 + lp->valid = 1;
1.173 + }
1.174 + /* factorization successful */
1.175 + ret = 0;
1.176 +fini: /* bring the return code to the calling program */
1.177 + return ret;
1.178 +}
1.179 +
1.180 +/***********************************************************************
1.181 +* NAME
1.182 +*
1.183 +* glp_bf_updated - check if the basis factorization has been updated
1.184 +*
1.185 +* SYNOPSIS
1.186 +*
1.187 +* int glp_bf_updated(glp_prob *lp);
1.188 +*
1.189 +* RETURNS
1.190 +*
1.191 +* If the basis factorization has been just computed from scratch, the
1.192 +* routine glp_bf_updated returns zero. Otherwise, if the factorization
1.193 +* has been updated one or more times, the routine returns non-zero. */
1.194 +
1.195 +int glp_bf_updated(glp_prob *lp)
1.196 +{ int cnt;
1.197 + if (!(lp->m == 0 || lp->valid))
1.198 + xerror("glp_bf_update: basis factorization does not exist\n");
1.199 +#if 0 /* 15/XI-2009 */
1.200 + cnt = (lp->m == 0 ? 0 : lp->bfd->upd_cnt);
1.201 +#else
1.202 + cnt = (lp->m == 0 ? 0 : bfd_get_count(lp->bfd));
1.203 +#endif
1.204 + return cnt;
1.205 +}
1.206 +
1.207 +/***********************************************************************
1.208 +* NAME
1.209 +*
1.210 +* glp_get_bfcp - retrieve basis factorization control parameters
1.211 +*
1.212 +* SYNOPSIS
1.213 +*
1.214 +* void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm);
1.215 +*
1.216 +* DESCRIPTION
1.217 +*
1.218 +* The routine glp_get_bfcp retrieves control parameters, which are
1.219 +* used on computing and updating the basis factorization associated
1.220 +* with the specified problem object.
1.221 +*
1.222 +* Current values of control parameters are stored by the routine in
1.223 +* a glp_bfcp structure, which the parameter parm points to. */
1.224 +
1.225 +void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm)
1.226 +{ glp_bfcp *bfcp = lp->bfcp;
1.227 + if (bfcp == NULL)
1.228 + { parm->type = GLP_BF_FT;
1.229 + parm->lu_size = 0;
1.230 + parm->piv_tol = 0.10;
1.231 + parm->piv_lim = 4;
1.232 + parm->suhl = GLP_ON;
1.233 + parm->eps_tol = 1e-15;
1.234 + parm->max_gro = 1e+10;
1.235 + parm->nfs_max = 100;
1.236 + parm->upd_tol = 1e-6;
1.237 + parm->nrs_max = 100;
1.238 + parm->rs_size = 0;
1.239 + }
1.240 + else
1.241 + memcpy(parm, bfcp, sizeof(glp_bfcp));
1.242 + return;
1.243 +}
1.244 +
1.245 +/***********************************************************************
1.246 +* NAME
1.247 +*
1.248 +* glp_set_bfcp - change basis factorization control parameters
1.249 +*
1.250 +* SYNOPSIS
1.251 +*
1.252 +* void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm);
1.253 +*
1.254 +* DESCRIPTION
1.255 +*
1.256 +* The routine glp_set_bfcp changes control parameters, which are used
1.257 +* by internal GLPK routines in computing and updating the basis
1.258 +* factorization associated with the specified problem object.
1.259 +*
1.260 +* New values of the control parameters should be passed in a structure
1.261 +* glp_bfcp, which the parameter parm points to.
1.262 +*
1.263 +* The parameter parm can be specified as NULL, in which case all
1.264 +* control parameters are reset to their default values. */
1.265 +
1.266 +#if 0 /* 15/XI-2009 */
1.267 +static void copy_bfcp(glp_prob *lp)
1.268 +{ glp_bfcp _parm, *parm = &_parm;
1.269 + BFD *bfd = lp->bfd;
1.270 + glp_get_bfcp(lp, parm);
1.271 + xassert(bfd != NULL);
1.272 + bfd->type = parm->type;
1.273 + bfd->lu_size = parm->lu_size;
1.274 + bfd->piv_tol = parm->piv_tol;
1.275 + bfd->piv_lim = parm->piv_lim;
1.276 + bfd->suhl = parm->suhl;
1.277 + bfd->eps_tol = parm->eps_tol;
1.278 + bfd->max_gro = parm->max_gro;
1.279 + bfd->nfs_max = parm->nfs_max;
1.280 + bfd->upd_tol = parm->upd_tol;
1.281 + bfd->nrs_max = parm->nrs_max;
1.282 + bfd->rs_size = parm->rs_size;
1.283 + return;
1.284 +}
1.285 +#else
1.286 +static void copy_bfcp(glp_prob *lp)
1.287 +{ glp_bfcp _parm, *parm = &_parm;
1.288 + glp_get_bfcp(lp, parm);
1.289 + bfd_set_parm(lp->bfd, parm);
1.290 + return;
1.291 +}
1.292 +#endif
1.293 +
1.294 +void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm)
1.295 +{ glp_bfcp *bfcp = lp->bfcp;
1.296 + if (parm == NULL)
1.297 + { /* reset to default values */
1.298 + if (bfcp != NULL)
1.299 + xfree(bfcp), lp->bfcp = NULL;
1.300 + }
1.301 + else
1.302 + { /* set to specified values */
1.303 + if (bfcp == NULL)
1.304 + bfcp = lp->bfcp = xmalloc(sizeof(glp_bfcp));
1.305 + memcpy(bfcp, parm, sizeof(glp_bfcp));
1.306 + if (!(bfcp->type == GLP_BF_FT || bfcp->type == GLP_BF_BG ||
1.307 + bfcp->type == GLP_BF_GR))
1.308 + xerror("glp_set_bfcp: type = %d; invalid parameter\n",
1.309 + bfcp->type);
1.310 + if (bfcp->lu_size < 0)
1.311 + xerror("glp_set_bfcp: lu_size = %d; invalid parameter\n",
1.312 + bfcp->lu_size);
1.313 + if (!(0.0 < bfcp->piv_tol && bfcp->piv_tol < 1.0))
1.314 + xerror("glp_set_bfcp: piv_tol = %g; invalid parameter\n",
1.315 + bfcp->piv_tol);
1.316 + if (bfcp->piv_lim < 1)
1.317 + xerror("glp_set_bfcp: piv_lim = %d; invalid parameter\n",
1.318 + bfcp->piv_lim);
1.319 + if (!(bfcp->suhl == GLP_ON || bfcp->suhl == GLP_OFF))
1.320 + xerror("glp_set_bfcp: suhl = %d; invalid parameter\n",
1.321 + bfcp->suhl);
1.322 + if (!(0.0 <= bfcp->eps_tol && bfcp->eps_tol <= 1e-6))
1.323 + xerror("glp_set_bfcp: eps_tol = %g; invalid parameter\n",
1.324 + bfcp->eps_tol);
1.325 + if (bfcp->max_gro < 1.0)
1.326 + xerror("glp_set_bfcp: max_gro = %g; invalid parameter\n",
1.327 + bfcp->max_gro);
1.328 + if (!(1 <= bfcp->nfs_max && bfcp->nfs_max <= 32767))
1.329 + xerror("glp_set_bfcp: nfs_max = %d; invalid parameter\n",
1.330 + bfcp->nfs_max);
1.331 + if (!(0.0 < bfcp->upd_tol && bfcp->upd_tol < 1.0))
1.332 + xerror("glp_set_bfcp: upd_tol = %g; invalid parameter\n",
1.333 + bfcp->upd_tol);
1.334 + if (!(1 <= bfcp->nrs_max && bfcp->nrs_max <= 32767))
1.335 + xerror("glp_set_bfcp: nrs_max = %d; invalid parameter\n",
1.336 + bfcp->nrs_max);
1.337 + if (bfcp->rs_size < 0)
1.338 + xerror("glp_set_bfcp: rs_size = %d; invalid parameter\n",
1.339 + bfcp->nrs_max);
1.340 + if (bfcp->rs_size == 0)
1.341 + bfcp->rs_size = 20 * bfcp->nrs_max;
1.342 + }
1.343 + if (lp->bfd != NULL) copy_bfcp(lp);
1.344 + return;
1.345 +}
1.346 +
1.347 +/***********************************************************************
1.348 +* NAME
1.349 +*
1.350 +* glp_get_bhead - retrieve the basis header information
1.351 +*
1.352 +* SYNOPSIS
1.353 +*
1.354 +* int glp_get_bhead(glp_prob *lp, int k);
1.355 +*
1.356 +* DESCRIPTION
1.357 +*
1.358 +* The routine glp_get_bhead returns the basis header information for
1.359 +* the current basis associated with the specified problem object.
1.360 +*
1.361 +* RETURNS
1.362 +*
1.363 +* If xB[k], 1 <= k <= m, is i-th auxiliary variable (1 <= i <= m), the
1.364 +* routine returns i. Otherwise, if xB[k] is j-th structural variable
1.365 +* (1 <= j <= n), the routine returns m+j. Here m is the number of rows
1.366 +* and n is the number of columns in the problem object. */
1.367 +
1.368 +int glp_get_bhead(glp_prob *lp, int k)
1.369 +{ if (!(lp->m == 0 || lp->valid))
1.370 + xerror("glp_get_bhead: basis factorization does not exist\n");
1.371 + if (!(1 <= k && k <= lp->m))
1.372 + xerror("glp_get_bhead: k = %d; index out of range\n", k);
1.373 + return lp->head[k];
1.374 +}
1.375 +
1.376 +/***********************************************************************
1.377 +* NAME
1.378 +*
1.379 +* glp_get_row_bind - retrieve row index in the basis header
1.380 +*
1.381 +* SYNOPSIS
1.382 +*
1.383 +* int glp_get_row_bind(glp_prob *lp, int i);
1.384 +*
1.385 +* RETURNS
1.386 +*
1.387 +* The routine glp_get_row_bind returns the index k of basic variable
1.388 +* xB[k], 1 <= k <= m, which is i-th auxiliary variable, 1 <= i <= m,
1.389 +* in the current basis associated with the specified problem object,
1.390 +* where m is the number of rows. However, if i-th auxiliary variable
1.391 +* is non-basic, the routine returns zero. */
1.392 +
1.393 +int glp_get_row_bind(glp_prob *lp, int i)
1.394 +{ if (!(lp->m == 0 || lp->valid))
1.395 + xerror("glp_get_row_bind: basis factorization does not exist\n"
1.396 + );
1.397 + if (!(1 <= i && i <= lp->m))
1.398 + xerror("glp_get_row_bind: i = %d; row number out of range\n",
1.399 + i);
1.400 + return lp->row[i]->bind;
1.401 +}
1.402 +
1.403 +/***********************************************************************
1.404 +* NAME
1.405 +*
1.406 +* glp_get_col_bind - retrieve column index in the basis header
1.407 +*
1.408 +* SYNOPSIS
1.409 +*
1.410 +* int glp_get_col_bind(glp_prob *lp, int j);
1.411 +*
1.412 +* RETURNS
1.413 +*
1.414 +* The routine glp_get_col_bind returns the index k of basic variable
1.415 +* xB[k], 1 <= k <= m, which is j-th structural variable, 1 <= j <= n,
1.416 +* in the current basis associated with the specified problem object,
1.417 +* where m is the number of rows, n is the number of columns. However,
1.418 +* if j-th structural variable is non-basic, the routine returns zero.*/
1.419 +
1.420 +int glp_get_col_bind(glp_prob *lp, int j)
1.421 +{ if (!(lp->m == 0 || lp->valid))
1.422 + xerror("glp_get_col_bind: basis factorization does not exist\n"
1.423 + );
1.424 + if (!(1 <= j && j <= lp->n))
1.425 + xerror("glp_get_col_bind: j = %d; column number out of range\n"
1.426 + , j);
1.427 + return lp->col[j]->bind;
1.428 +}
1.429 +
1.430 +/***********************************************************************
1.431 +* NAME
1.432 +*
1.433 +* glp_ftran - perform forward transformation (solve system B*x = b)
1.434 +*
1.435 +* SYNOPSIS
1.436 +*
1.437 +* void glp_ftran(glp_prob *lp, double x[]);
1.438 +*
1.439 +* DESCRIPTION
1.440 +*
1.441 +* The routine glp_ftran performs forward transformation, i.e. solves
1.442 +* the system B*x = b, where B is the basis matrix corresponding to the
1.443 +* current basis for the specified problem object, x is the vector of
1.444 +* unknowns to be computed, b is the vector of right-hand sides.
1.445 +*
1.446 +* On entry elements of the vector b should be stored in dense format
1.447 +* in locations x[1], ..., x[m], where m is the number of rows. On exit
1.448 +* the routine stores elements of the vector x in the same locations.
1.449 +*
1.450 +* SCALING/UNSCALING
1.451 +*
1.452 +* Let A~ = (I | -A) is the augmented constraint matrix of the original
1.453 +* (unscaled) problem. In the scaled LP problem instead the matrix A the
1.454 +* scaled matrix A" = R*A*S is actually used, so
1.455 +*
1.456 +* A~" = (I | A") = (I | R*A*S) = (R*I*inv(R) | R*A*S) =
1.457 +* (1)
1.458 +* = R*(I | A)*S~ = R*A~*S~,
1.459 +*
1.460 +* is the scaled augmented constraint matrix, where R and S are diagonal
1.461 +* scaling matrices used to scale rows and columns of the matrix A, and
1.462 +*
1.463 +* S~ = diag(inv(R) | S) (2)
1.464 +*
1.465 +* is an augmented diagonal scaling matrix.
1.466 +*
1.467 +* By definition:
1.468 +*
1.469 +* A~ = (B | N), (3)
1.470 +*
1.471 +* where B is the basic matrix, which consists of basic columns of the
1.472 +* augmented constraint matrix A~, and N is a matrix, which consists of
1.473 +* non-basic columns of A~. From (1) it follows that:
1.474 +*
1.475 +* A~" = (B" | N") = (R*B*SB | R*N*SN), (4)
1.476 +*
1.477 +* where SB and SN are parts of the augmented scaling matrix S~, which
1.478 +* correspond to basic and non-basic variables, respectively. Therefore
1.479 +*
1.480 +* B" = R*B*SB, (5)
1.481 +*
1.482 +* which is the scaled basis matrix. */
1.483 +
1.484 +void glp_ftran(glp_prob *lp, double x[])
1.485 +{ int m = lp->m;
1.486 + GLPROW **row = lp->row;
1.487 + GLPCOL **col = lp->col;
1.488 + int i, k;
1.489 + /* B*x = b ===> (R*B*SB)*(inv(SB)*x) = R*b ===>
1.490 + B"*x" = b", where b" = R*b, x = SB*x" */
1.491 + if (!(m == 0 || lp->valid))
1.492 + xerror("glp_ftran: basis factorization does not exist\n");
1.493 + /* b" := R*b */
1.494 + for (i = 1; i <= m; i++)
1.495 + x[i] *= row[i]->rii;
1.496 + /* x" := inv(B")*b" */
1.497 + if (m > 0) bfd_ftran(lp->bfd, x);
1.498 + /* x := SB*x" */
1.499 + for (i = 1; i <= m; i++)
1.500 + { k = lp->head[i];
1.501 + if (k <= m)
1.502 + x[i] /= row[k]->rii;
1.503 + else
1.504 + x[i] *= col[k-m]->sjj;
1.505 + }
1.506 + return;
1.507 +}
1.508 +
1.509 +/***********************************************************************
1.510 +* NAME
1.511 +*
1.512 +* glp_btran - perform backward transformation (solve system B'*x = b)
1.513 +*
1.514 +* SYNOPSIS
1.515 +*
1.516 +* void glp_btran(glp_prob *lp, double x[]);
1.517 +*
1.518 +* DESCRIPTION
1.519 +*
1.520 +* The routine glp_btran performs backward transformation, i.e. solves
1.521 +* the system B'*x = b, where B' is a matrix transposed to the basis
1.522 +* matrix corresponding to the current basis for the specified problem
1.523 +* problem object, x is the vector of unknowns to be computed, b is the
1.524 +* vector of right-hand sides.
1.525 +*
1.526 +* On entry elements of the vector b should be stored in dense format
1.527 +* in locations x[1], ..., x[m], where m is the number of rows. On exit
1.528 +* the routine stores elements of the vector x in the same locations.
1.529 +*
1.530 +* SCALING/UNSCALING
1.531 +*
1.532 +* See comments to the routine glp_ftran. */
1.533 +
1.534 +void glp_btran(glp_prob *lp, double x[])
1.535 +{ int m = lp->m;
1.536 + GLPROW **row = lp->row;
1.537 + GLPCOL **col = lp->col;
1.538 + int i, k;
1.539 + /* B'*x = b ===> (SB*B'*R)*(inv(R)*x) = SB*b ===>
1.540 + (B")'*x" = b", where b" = SB*b, x = R*x" */
1.541 + if (!(m == 0 || lp->valid))
1.542 + xerror("glp_btran: basis factorization does not exist\n");
1.543 + /* b" := SB*b */
1.544 + for (i = 1; i <= m; i++)
1.545 + { k = lp->head[i];
1.546 + if (k <= m)
1.547 + x[i] /= row[k]->rii;
1.548 + else
1.549 + x[i] *= col[k-m]->sjj;
1.550 + }
1.551 + /* x" := inv[(B")']*b" */
1.552 + if (m > 0) bfd_btran(lp->bfd, x);
1.553 + /* x := R*x" */
1.554 + for (i = 1; i <= m; i++)
1.555 + x[i] *= row[i]->rii;
1.556 + return;
1.557 +}
1.558 +
1.559 +/***********************************************************************
1.560 +* NAME
1.561 +*
1.562 +* glp_warm_up - "warm up" LP basis
1.563 +*
1.564 +* SYNOPSIS
1.565 +*
1.566 +* int glp_warm_up(glp_prob *P);
1.567 +*
1.568 +* DESCRIPTION
1.569 +*
1.570 +* The routine glp_warm_up "warms up" the LP basis for the specified
1.571 +* problem object using current statuses assigned to rows and columns
1.572 +* (that is, to auxiliary and structural variables).
1.573 +*
1.574 +* This operation includes computing factorization of the basis matrix
1.575 +* (if it does not exist), computing primal and dual components of basic
1.576 +* solution, and determining the solution status.
1.577 +*
1.578 +* RETURNS
1.579 +*
1.580 +* 0 The operation has been successfully performed.
1.581 +*
1.582 +* GLP_EBADB
1.583 +* The basis matrix is invalid, i.e. the number of basic (auxiliary
1.584 +* and structural) variables differs from the number of rows in the
1.585 +* problem object.
1.586 +*
1.587 +* GLP_ESING
1.588 +* The basis matrix is singular within the working precision.
1.589 +*
1.590 +* GLP_ECOND
1.591 +* The basis matrix is ill-conditioned. */
1.592 +
1.593 +int glp_warm_up(glp_prob *P)
1.594 +{ GLPROW *row;
1.595 + GLPCOL *col;
1.596 + GLPAIJ *aij;
1.597 + int i, j, type, ret;
1.598 + double eps, temp, *work;
1.599 + /* invalidate basic solution */
1.600 + P->pbs_stat = P->dbs_stat = GLP_UNDEF;
1.601 + P->obj_val = 0.0;
1.602 + P->some = 0;
1.603 + for (i = 1; i <= P->m; i++)
1.604 + { row = P->row[i];
1.605 + row->prim = row->dual = 0.0;
1.606 + }
1.607 + for (j = 1; j <= P->n; j++)
1.608 + { col = P->col[j];
1.609 + col->prim = col->dual = 0.0;
1.610 + }
1.611 + /* compute the basis factorization, if necessary */
1.612 + if (!glp_bf_exists(P))
1.613 + { ret = glp_factorize(P);
1.614 + if (ret != 0) goto done;
1.615 + }
1.616 + /* allocate working array */
1.617 + work = xcalloc(1+P->m, sizeof(double));
1.618 + /* determine and store values of non-basic variables, compute
1.619 + vector (- N * xN) */
1.620 + for (i = 1; i <= P->m; i++)
1.621 + work[i] = 0.0;
1.622 + for (i = 1; i <= P->m; i++)
1.623 + { row = P->row[i];
1.624 + if (row->stat == GLP_BS)
1.625 + continue;
1.626 + else if (row->stat == GLP_NL)
1.627 + row->prim = row->lb;
1.628 + else if (row->stat == GLP_NU)
1.629 + row->prim = row->ub;
1.630 + else if (row->stat == GLP_NF)
1.631 + row->prim = 0.0;
1.632 + else if (row->stat == GLP_NS)
1.633 + row->prim = row->lb;
1.634 + else
1.635 + xassert(row != row);
1.636 + /* N[j] is i-th column of matrix (I|-A) */
1.637 + work[i] -= row->prim;
1.638 + }
1.639 + for (j = 1; j <= P->n; j++)
1.640 + { col = P->col[j];
1.641 + if (col->stat == GLP_BS)
1.642 + continue;
1.643 + else if (col->stat == GLP_NL)
1.644 + col->prim = col->lb;
1.645 + else if (col->stat == GLP_NU)
1.646 + col->prim = col->ub;
1.647 + else if (col->stat == GLP_NF)
1.648 + col->prim = 0.0;
1.649 + else if (col->stat == GLP_NS)
1.650 + col->prim = col->lb;
1.651 + else
1.652 + xassert(col != col);
1.653 + /* N[j] is (m+j)-th column of matrix (I|-A) */
1.654 + if (col->prim != 0.0)
1.655 + { for (aij = col->ptr; aij != NULL; aij = aij->c_next)
1.656 + work[aij->row->i] += aij->val * col->prim;
1.657 + }
1.658 + }
1.659 + /* compute vector of basic variables xB = - inv(B) * N * xN */
1.660 + glp_ftran(P, work);
1.661 + /* store values of basic variables, check primal feasibility */
1.662 + P->pbs_stat = GLP_FEAS;
1.663 + for (i = 1; i <= P->m; i++)
1.664 + { row = P->row[i];
1.665 + if (row->stat != GLP_BS)
1.666 + continue;
1.667 + row->prim = work[row->bind];
1.668 + type = row->type;
1.669 + if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
1.670 + { eps = 1e-6 + 1e-9 * fabs(row->lb);
1.671 + if (row->prim < row->lb - eps)
1.672 + P->pbs_stat = GLP_INFEAS;
1.673 + }
1.674 + if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
1.675 + { eps = 1e-6 + 1e-9 * fabs(row->ub);
1.676 + if (row->prim > row->ub + eps)
1.677 + P->pbs_stat = GLP_INFEAS;
1.678 + }
1.679 + }
1.680 + for (j = 1; j <= P->n; j++)
1.681 + { col = P->col[j];
1.682 + if (col->stat != GLP_BS)
1.683 + continue;
1.684 + col->prim = work[col->bind];
1.685 + type = col->type;
1.686 + if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
1.687 + { eps = 1e-6 + 1e-9 * fabs(col->lb);
1.688 + if (col->prim < col->lb - eps)
1.689 + P->pbs_stat = GLP_INFEAS;
1.690 + }
1.691 + if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
1.692 + { eps = 1e-6 + 1e-9 * fabs(col->ub);
1.693 + if (col->prim > col->ub + eps)
1.694 + P->pbs_stat = GLP_INFEAS;
1.695 + }
1.696 + }
1.697 + /* compute value of the objective function */
1.698 + P->obj_val = P->c0;
1.699 + for (j = 1; j <= P->n; j++)
1.700 + { col = P->col[j];
1.701 + P->obj_val += col->coef * col->prim;
1.702 + }
1.703 + /* build vector cB of objective coefficients at basic variables */
1.704 + for (i = 1; i <= P->m; i++)
1.705 + work[i] = 0.0;
1.706 + for (j = 1; j <= P->n; j++)
1.707 + { col = P->col[j];
1.708 + if (col->stat == GLP_BS)
1.709 + work[col->bind] = col->coef;
1.710 + }
1.711 + /* compute vector of simplex multipliers pi = inv(B') * cB */
1.712 + glp_btran(P, work);
1.713 + /* compute and store reduced costs of non-basic variables d[j] =
1.714 + c[j] - N'[j] * pi, check dual feasibility */
1.715 + P->dbs_stat = GLP_FEAS;
1.716 + for (i = 1; i <= P->m; i++)
1.717 + { row = P->row[i];
1.718 + if (row->stat == GLP_BS)
1.719 + { row->dual = 0.0;
1.720 + continue;
1.721 + }
1.722 + /* N[j] is i-th column of matrix (I|-A) */
1.723 + row->dual = - work[i];
1.724 + type = row->type;
1.725 + temp = (P->dir == GLP_MIN ? + row->dual : - row->dual);
1.726 + if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
1.727 + (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
1.728 + P->dbs_stat = GLP_INFEAS;
1.729 + }
1.730 + for (j = 1; j <= P->n; j++)
1.731 + { col = P->col[j];
1.732 + if (col->stat == GLP_BS)
1.733 + { col->dual = 0.0;
1.734 + continue;
1.735 + }
1.736 + /* N[j] is (m+j)-th column of matrix (I|-A) */
1.737 + col->dual = col->coef;
1.738 + for (aij = col->ptr; aij != NULL; aij = aij->c_next)
1.739 + col->dual += aij->val * work[aij->row->i];
1.740 + type = col->type;
1.741 + temp = (P->dir == GLP_MIN ? + col->dual : - col->dual);
1.742 + if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
1.743 + (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
1.744 + P->dbs_stat = GLP_INFEAS;
1.745 + }
1.746 + /* free working array */
1.747 + xfree(work);
1.748 + ret = 0;
1.749 +done: return ret;
1.750 +}
1.751 +
1.752 +/***********************************************************************
1.753 +* NAME
1.754 +*
1.755 +* glp_eval_tab_row - compute row of the simplex tableau
1.756 +*
1.757 +* SYNOPSIS
1.758 +*
1.759 +* int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[]);
1.760 +*
1.761 +* DESCRIPTION
1.762 +*
1.763 +* The routine glp_eval_tab_row computes a row of the current simplex
1.764 +* tableau for the basic variable, which is specified by the number k:
1.765 +* if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
1.766 +* x[k] is (k-m)-th structural variable, where m is number of rows, and
1.767 +* n is number of columns. The current basis must be available.
1.768 +*
1.769 +* The routine stores column indices and numerical values of non-zero
1.770 +* elements of the computed row using sparse format to the locations
1.771 +* ind[1], ..., ind[len] and val[1], ..., val[len], respectively, where
1.772 +* 0 <= len <= n is number of non-zeros returned on exit.
1.773 +*
1.774 +* Element indices stored in the array ind have the same sense as the
1.775 +* index k, i.e. indices 1 to m denote auxiliary variables and indices
1.776 +* m+1 to m+n denote structural ones (all these variables are obviously
1.777 +* non-basic by definition).
1.778 +*
1.779 +* The computed row shows how the specified basic variable x[k] = xB[i]
1.780 +* depends on non-basic variables:
1.781 +*
1.782 +* xB[i] = alfa[i,1]*xN[1] + alfa[i,2]*xN[2] + ... + alfa[i,n]*xN[n],
1.783 +*
1.784 +* where alfa[i,j] are elements of the simplex table row, xN[j] are
1.785 +* non-basic (auxiliary and structural) variables.
1.786 +*
1.787 +* RETURNS
1.788 +*
1.789 +* The routine returns number of non-zero elements in the simplex table
1.790 +* row stored in the arrays ind and val.
1.791 +*
1.792 +* BACKGROUND
1.793 +*
1.794 +* The system of equality constraints of the LP problem is:
1.795 +*
1.796 +* xR = A * xS, (1)
1.797 +*
1.798 +* where xR is the vector of auxliary variables, xS is the vector of
1.799 +* structural variables, A is the matrix of constraint coefficients.
1.800 +*
1.801 +* The system (1) can be written in homogenous form as follows:
1.802 +*
1.803 +* A~ * x = 0, (2)
1.804 +*
1.805 +* where A~ = (I | -A) is the augmented constraint matrix (has m rows
1.806 +* and m+n columns), x = (xR | xS) is the vector of all (auxiliary and
1.807 +* structural) variables.
1.808 +*
1.809 +* By definition for the current basis we have:
1.810 +*
1.811 +* A~ = (B | N), (3)
1.812 +*
1.813 +* where B is the basis matrix. Thus, the system (2) can be written as:
1.814 +*
1.815 +* B * xB + N * xN = 0. (4)
1.816 +*
1.817 +* From (4) it follows that:
1.818 +*
1.819 +* xB = A^ * xN, (5)
1.820 +*
1.821 +* where the matrix
1.822 +*
1.823 +* A^ = - inv(B) * N (6)
1.824 +*
1.825 +* is called the simplex table.
1.826 +*
1.827 +* It is understood that i-th row of the simplex table is:
1.828 +*
1.829 +* e * A^ = - e * inv(B) * N, (7)
1.830 +*
1.831 +* where e is a unity vector with e[i] = 1.
1.832 +*
1.833 +* To compute i-th row of the simplex table the routine first computes
1.834 +* i-th row of the inverse:
1.835 +*
1.836 +* rho = inv(B') * e, (8)
1.837 +*
1.838 +* where B' is a matrix transposed to B, and then computes elements of
1.839 +* i-th row of the simplex table as scalar products:
1.840 +*
1.841 +* alfa[i,j] = - rho * N[j] for all j, (9)
1.842 +*
1.843 +* where N[j] is a column of the augmented constraint matrix A~, which
1.844 +* corresponds to some non-basic auxiliary or structural variable. */
1.845 +
1.846 +int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[])
1.847 +{ int m = lp->m;
1.848 + int n = lp->n;
1.849 + int i, t, len, lll, *iii;
1.850 + double alfa, *rho, *vvv;
1.851 + if (!(m == 0 || lp->valid))
1.852 + xerror("glp_eval_tab_row: basis factorization does not exist\n"
1.853 + );
1.854 + if (!(1 <= k && k <= m+n))
1.855 + xerror("glp_eval_tab_row: k = %d; variable number out of range"
1.856 + , k);
1.857 + /* determine xB[i] which corresponds to x[k] */
1.858 + if (k <= m)
1.859 + i = glp_get_row_bind(lp, k);
1.860 + else
1.861 + i = glp_get_col_bind(lp, k-m);
1.862 + if (i == 0)
1.863 + xerror("glp_eval_tab_row: k = %d; variable must be basic", k);
1.864 + xassert(1 <= i && i <= m);
1.865 + /* allocate working arrays */
1.866 + rho = xcalloc(1+m, sizeof(double));
1.867 + iii = xcalloc(1+m, sizeof(int));
1.868 + vvv = xcalloc(1+m, sizeof(double));
1.869 + /* compute i-th row of the inverse; see (8) */
1.870 + for (t = 1; t <= m; t++) rho[t] = 0.0;
1.871 + rho[i] = 1.0;
1.872 + glp_btran(lp, rho);
1.873 + /* compute i-th row of the simplex table */
1.874 + len = 0;
1.875 + for (k = 1; k <= m+n; k++)
1.876 + { if (k <= m)
1.877 + { /* x[k] is auxiliary variable, so N[k] is a unity column */
1.878 + if (glp_get_row_stat(lp, k) == GLP_BS) continue;
1.879 + /* compute alfa[i,j]; see (9) */
1.880 + alfa = - rho[k];
1.881 + }
1.882 + else
1.883 + { /* x[k] is structural variable, so N[k] is a column of the
1.884 + original constraint matrix A with negative sign */
1.885 + if (glp_get_col_stat(lp, k-m) == GLP_BS) continue;
1.886 + /* compute alfa[i,j]; see (9) */
1.887 + lll = glp_get_mat_col(lp, k-m, iii, vvv);
1.888 + alfa = 0.0;
1.889 + for (t = 1; t <= lll; t++) alfa += rho[iii[t]] * vvv[t];
1.890 + }
1.891 + /* store alfa[i,j] */
1.892 + if (alfa != 0.0) len++, ind[len] = k, val[len] = alfa;
1.893 + }
1.894 + xassert(len <= n);
1.895 + /* free working arrays */
1.896 + xfree(rho);
1.897 + xfree(iii);
1.898 + xfree(vvv);
1.899 + /* return to the calling program */
1.900 + return len;
1.901 +}
1.902 +
1.903 +/***********************************************************************
1.904 +* NAME
1.905 +*
1.906 +* glp_eval_tab_col - compute column of the simplex tableau
1.907 +*
1.908 +* SYNOPSIS
1.909 +*
1.910 +* int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[]);
1.911 +*
1.912 +* DESCRIPTION
1.913 +*
1.914 +* The routine glp_eval_tab_col computes a column of the current simplex
1.915 +* table for the non-basic variable, which is specified by the number k:
1.916 +* if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
1.917 +* x[k] is (k-m)-th structural variable, where m is number of rows, and
1.918 +* n is number of columns. The current basis must be available.
1.919 +*
1.920 +* The routine stores row indices and numerical values of non-zero
1.921 +* elements of the computed column using sparse format to the locations
1.922 +* ind[1], ..., ind[len] and val[1], ..., val[len] respectively, where
1.923 +* 0 <= len <= m is number of non-zeros returned on exit.
1.924 +*
1.925 +* Element indices stored in the array ind have the same sense as the
1.926 +* index k, i.e. indices 1 to m denote auxiliary variables and indices
1.927 +* m+1 to m+n denote structural ones (all these variables are obviously
1.928 +* basic by the definition).
1.929 +*
1.930 +* The computed column shows how basic variables depend on the specified
1.931 +* non-basic variable x[k] = xN[j]:
1.932 +*
1.933 +* xB[1] = ... + alfa[1,j]*xN[j] + ...
1.934 +* xB[2] = ... + alfa[2,j]*xN[j] + ...
1.935 +* . . . . . .
1.936 +* xB[m] = ... + alfa[m,j]*xN[j] + ...
1.937 +*
1.938 +* where alfa[i,j] are elements of the simplex table column, xB[i] are
1.939 +* basic (auxiliary and structural) variables.
1.940 +*
1.941 +* RETURNS
1.942 +*
1.943 +* The routine returns number of non-zero elements in the simplex table
1.944 +* column stored in the arrays ind and val.
1.945 +*
1.946 +* BACKGROUND
1.947 +*
1.948 +* As it was explained in comments to the routine glp_eval_tab_row (see
1.949 +* above) the simplex table is the following matrix:
1.950 +*
1.951 +* A^ = - inv(B) * N. (1)
1.952 +*
1.953 +* Therefore j-th column of the simplex table is:
1.954 +*
1.955 +* A^ * e = - inv(B) * N * e = - inv(B) * N[j], (2)
1.956 +*
1.957 +* where e is a unity vector with e[j] = 1, B is the basis matrix, N[j]
1.958 +* is a column of the augmented constraint matrix A~, which corresponds
1.959 +* to the given non-basic auxiliary or structural variable. */
1.960 +
1.961 +int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[])
1.962 +{ int m = lp->m;
1.963 + int n = lp->n;
1.964 + int t, len, stat;
1.965 + double *col;
1.966 + if (!(m == 0 || lp->valid))
1.967 + xerror("glp_eval_tab_col: basis factorization does not exist\n"
1.968 + );
1.969 + if (!(1 <= k && k <= m+n))
1.970 + xerror("glp_eval_tab_col: k = %d; variable number out of range"
1.971 + , k);
1.972 + if (k <= m)
1.973 + stat = glp_get_row_stat(lp, k);
1.974 + else
1.975 + stat = glp_get_col_stat(lp, k-m);
1.976 + if (stat == GLP_BS)
1.977 + xerror("glp_eval_tab_col: k = %d; variable must be non-basic",
1.978 + k);
1.979 + /* obtain column N[k] with negative sign */
1.980 + col = xcalloc(1+m, sizeof(double));
1.981 + for (t = 1; t <= m; t++) col[t] = 0.0;
1.982 + if (k <= m)
1.983 + { /* x[k] is auxiliary variable, so N[k] is a unity column */
1.984 + col[k] = -1.0;
1.985 + }
1.986 + else
1.987 + { /* x[k] is structural variable, so N[k] is a column of the
1.988 + original constraint matrix A with negative sign */
1.989 + len = glp_get_mat_col(lp, k-m, ind, val);
1.990 + for (t = 1; t <= len; t++) col[ind[t]] = val[t];
1.991 + }
1.992 + /* compute column of the simplex table, which corresponds to the
1.993 + specified non-basic variable x[k] */
1.994 + glp_ftran(lp, col);
1.995 + len = 0;
1.996 + for (t = 1; t <= m; t++)
1.997 + { if (col[t] != 0.0)
1.998 + { len++;
1.999 + ind[len] = glp_get_bhead(lp, t);
1.1000 + val[len] = col[t];
1.1001 + }
1.1002 + }
1.1003 + xfree(col);
1.1004 + /* return to the calling program */
1.1005 + return len;
1.1006 +}
1.1007 +
1.1008 +/***********************************************************************
1.1009 +* NAME
1.1010 +*
1.1011 +* glp_transform_row - transform explicitly specified row
1.1012 +*
1.1013 +* SYNOPSIS
1.1014 +*
1.1015 +* int glp_transform_row(glp_prob *P, int len, int ind[], double val[]);
1.1016 +*
1.1017 +* DESCRIPTION
1.1018 +*
1.1019 +* The routine glp_transform_row performs the same operation as the
1.1020 +* routine glp_eval_tab_row with exception that the row to be
1.1021 +* transformed is specified explicitly as a sparse vector.
1.1022 +*
1.1023 +* The explicitly specified row may be thought as a linear form:
1.1024 +*
1.1025 +* x = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n], (1)
1.1026 +*
1.1027 +* where x is an auxiliary variable for this row, a[j] are coefficients
1.1028 +* of the linear form, x[m+j] are structural variables.
1.1029 +*
1.1030 +* On entry column indices and numerical values of non-zero elements of
1.1031 +* the row should be stored in locations ind[1], ..., ind[len] and
1.1032 +* val[1], ..., val[len], where len is the number of non-zero elements.
1.1033 +*
1.1034 +* This routine uses the system of equality constraints and the current
1.1035 +* basis in order to express the auxiliary variable x in (1) through the
1.1036 +* current non-basic variables (as if the transformed row were added to
1.1037 +* the problem object and its auxiliary variable were basic), i.e. the
1.1038 +* resultant row has the form:
1.1039 +*
1.1040 +* x = alfa[1]*xN[1] + alfa[2]*xN[2] + ... + alfa[n]*xN[n], (2)
1.1041 +*
1.1042 +* where xN[j] are non-basic (auxiliary or structural) variables, n is
1.1043 +* the number of columns in the LP problem object.
1.1044 +*
1.1045 +* On exit the routine stores indices and numerical values of non-zero
1.1046 +* elements of the resultant row (2) in locations ind[1], ..., ind[len']
1.1047 +* and val[1], ..., val[len'], where 0 <= len' <= n is the number of
1.1048 +* non-zero elements in the resultant row returned by the routine. Note
1.1049 +* that indices (numbers) of non-basic variables stored in the array ind
1.1050 +* correspond to original ordinal numbers of variables: indices 1 to m
1.1051 +* mean auxiliary variables and indices m+1 to m+n mean structural ones.
1.1052 +*
1.1053 +* RETURNS
1.1054 +*
1.1055 +* The routine returns len', which is the number of non-zero elements in
1.1056 +* the resultant row stored in the arrays ind and val.
1.1057 +*
1.1058 +* BACKGROUND
1.1059 +*
1.1060 +* The explicitly specified row (1) is transformed in the same way as it
1.1061 +* were the objective function row.
1.1062 +*
1.1063 +* From (1) it follows that:
1.1064 +*
1.1065 +* x = aB * xB + aN * xN, (3)
1.1066 +*
1.1067 +* where xB is the vector of basic variables, xN is the vector of
1.1068 +* non-basic variables.
1.1069 +*
1.1070 +* The simplex table, which corresponds to the current basis, is:
1.1071 +*
1.1072 +* xB = [-inv(B) * N] * xN. (4)
1.1073 +*
1.1074 +* Therefore substituting xB from (4) to (3) we have:
1.1075 +*
1.1076 +* x = aB * [-inv(B) * N] * xN + aN * xN =
1.1077 +* (5)
1.1078 +* = rho * (-N) * xN + aN * xN = alfa * xN,
1.1079 +*
1.1080 +* where:
1.1081 +*
1.1082 +* rho = inv(B') * aB, (6)
1.1083 +*
1.1084 +* and
1.1085 +*
1.1086 +* alfa = aN + rho * (-N) (7)
1.1087 +*
1.1088 +* is the resultant row computed by the routine. */
1.1089 +
1.1090 +int glp_transform_row(glp_prob *P, int len, int ind[], double val[])
1.1091 +{ int i, j, k, m, n, t, lll, *iii;
1.1092 + double alfa, *a, *aB, *rho, *vvv;
1.1093 + if (!glp_bf_exists(P))
1.1094 + xerror("glp_transform_row: basis factorization does not exist "
1.1095 + "\n");
1.1096 + m = glp_get_num_rows(P);
1.1097 + n = glp_get_num_cols(P);
1.1098 + /* unpack the row to be transformed to the array a */
1.1099 + a = xcalloc(1+n, sizeof(double));
1.1100 + for (j = 1; j <= n; j++) a[j] = 0.0;
1.1101 + if (!(0 <= len && len <= n))
1.1102 + xerror("glp_transform_row: len = %d; invalid row length\n",
1.1103 + len);
1.1104 + for (t = 1; t <= len; t++)
1.1105 + { j = ind[t];
1.1106 + if (!(1 <= j && j <= n))
1.1107 + xerror("glp_transform_row: ind[%d] = %d; column index out o"
1.1108 + "f range\n", t, j);
1.1109 + if (val[t] == 0.0)
1.1110 + xerror("glp_transform_row: val[%d] = 0; zero coefficient no"
1.1111 + "t allowed\n", t);
1.1112 + if (a[j] != 0.0)
1.1113 + xerror("glp_transform_row: ind[%d] = %d; duplicate column i"
1.1114 + "ndices not allowed\n", t, j);
1.1115 + a[j] = val[t];
1.1116 + }
1.1117 + /* construct the vector aB */
1.1118 + aB = xcalloc(1+m, sizeof(double));
1.1119 + for (i = 1; i <= m; i++)
1.1120 + { k = glp_get_bhead(P, i);
1.1121 + /* xB[i] is k-th original variable */
1.1122 + xassert(1 <= k && k <= m+n);
1.1123 + aB[i] = (k <= m ? 0.0 : a[k-m]);
1.1124 + }
1.1125 + /* solve the system B'*rho = aB to compute the vector rho */
1.1126 + rho = aB, glp_btran(P, rho);
1.1127 + /* compute coefficients at non-basic auxiliary variables */
1.1128 + len = 0;
1.1129 + for (i = 1; i <= m; i++)
1.1130 + { if (glp_get_row_stat(P, i) != GLP_BS)
1.1131 + { alfa = - rho[i];
1.1132 + if (alfa != 0.0)
1.1133 + { len++;
1.1134 + ind[len] = i;
1.1135 + val[len] = alfa;
1.1136 + }
1.1137 + }
1.1138 + }
1.1139 + /* compute coefficients at non-basic structural variables */
1.1140 + iii = xcalloc(1+m, sizeof(int));
1.1141 + vvv = xcalloc(1+m, sizeof(double));
1.1142 + for (j = 1; j <= n; j++)
1.1143 + { if (glp_get_col_stat(P, j) != GLP_BS)
1.1144 + { alfa = a[j];
1.1145 + lll = glp_get_mat_col(P, j, iii, vvv);
1.1146 + for (t = 1; t <= lll; t++) alfa += vvv[t] * rho[iii[t]];
1.1147 + if (alfa != 0.0)
1.1148 + { len++;
1.1149 + ind[len] = m+j;
1.1150 + val[len] = alfa;
1.1151 + }
1.1152 + }
1.1153 + }
1.1154 + xassert(len <= n);
1.1155 + xfree(iii);
1.1156 + xfree(vvv);
1.1157 + xfree(aB);
1.1158 + xfree(a);
1.1159 + return len;
1.1160 +}
1.1161 +
1.1162 +/***********************************************************************
1.1163 +* NAME
1.1164 +*
1.1165 +* glp_transform_col - transform explicitly specified column
1.1166 +*
1.1167 +* SYNOPSIS
1.1168 +*
1.1169 +* int glp_transform_col(glp_prob *P, int len, int ind[], double val[]);
1.1170 +*
1.1171 +* DESCRIPTION
1.1172 +*
1.1173 +* The routine glp_transform_col performs the same operation as the
1.1174 +* routine glp_eval_tab_col with exception that the column to be
1.1175 +* transformed is specified explicitly as a sparse vector.
1.1176 +*
1.1177 +* The explicitly specified column may be thought as if it were added
1.1178 +* to the original system of equality constraints:
1.1179 +*
1.1180 +* x[1] = a[1,1]*x[m+1] + ... + a[1,n]*x[m+n] + a[1]*x
1.1181 +* x[2] = a[2,1]*x[m+1] + ... + a[2,n]*x[m+n] + a[2]*x (1)
1.1182 +* . . . . . . . . . . . . . . .
1.1183 +* x[m] = a[m,1]*x[m+1] + ... + a[m,n]*x[m+n] + a[m]*x
1.1184 +*
1.1185 +* where x[i] are auxiliary variables, x[m+j] are structural variables,
1.1186 +* x is a structural variable for the explicitly specified column, a[i]
1.1187 +* are constraint coefficients for x.
1.1188 +*
1.1189 +* On entry row indices and numerical values of non-zero elements of
1.1190 +* the column should be stored in locations ind[1], ..., ind[len] and
1.1191 +* val[1], ..., val[len], where len is the number of non-zero elements.
1.1192 +*
1.1193 +* This routine uses the system of equality constraints and the current
1.1194 +* basis in order to express the current basic variables through the
1.1195 +* structural variable x in (1) (as if the transformed column were added
1.1196 +* to the problem object and the variable x were non-basic), i.e. the
1.1197 +* resultant column has the form:
1.1198 +*
1.1199 +* xB[1] = ... + alfa[1]*x
1.1200 +* xB[2] = ... + alfa[2]*x (2)
1.1201 +* . . . . . .
1.1202 +* xB[m] = ... + alfa[m]*x
1.1203 +*
1.1204 +* where xB are basic (auxiliary and structural) variables, m is the
1.1205 +* number of rows in the problem object.
1.1206 +*
1.1207 +* On exit the routine stores indices and numerical values of non-zero
1.1208 +* elements of the resultant column (2) in locations ind[1], ...,
1.1209 +* ind[len'] and val[1], ..., val[len'], where 0 <= len' <= m is the
1.1210 +* number of non-zero element in the resultant column returned by the
1.1211 +* routine. Note that indices (numbers) of basic variables stored in
1.1212 +* the array ind correspond to original ordinal numbers of variables:
1.1213 +* indices 1 to m mean auxiliary variables and indices m+1 to m+n mean
1.1214 +* structural ones.
1.1215 +*
1.1216 +* RETURNS
1.1217 +*
1.1218 +* The routine returns len', which is the number of non-zero elements
1.1219 +* in the resultant column stored in the arrays ind and val.
1.1220 +*
1.1221 +* BACKGROUND
1.1222 +*
1.1223 +* The explicitly specified column (1) is transformed in the same way
1.1224 +* as any other column of the constraint matrix using the formula:
1.1225 +*
1.1226 +* alfa = inv(B) * a, (3)
1.1227 +*
1.1228 +* where alfa is the resultant column computed by the routine. */
1.1229 +
1.1230 +int glp_transform_col(glp_prob *P, int len, int ind[], double val[])
1.1231 +{ int i, m, t;
1.1232 + double *a, *alfa;
1.1233 + if (!glp_bf_exists(P))
1.1234 + xerror("glp_transform_col: basis factorization does not exist "
1.1235 + "\n");
1.1236 + m = glp_get_num_rows(P);
1.1237 + /* unpack the column to be transformed to the array a */
1.1238 + a = xcalloc(1+m, sizeof(double));
1.1239 + for (i = 1; i <= m; i++) a[i] = 0.0;
1.1240 + if (!(0 <= len && len <= m))
1.1241 + xerror("glp_transform_col: len = %d; invalid column length\n",
1.1242 + len);
1.1243 + for (t = 1; t <= len; t++)
1.1244 + { i = ind[t];
1.1245 + if (!(1 <= i && i <= m))
1.1246 + xerror("glp_transform_col: ind[%d] = %d; row index out of r"
1.1247 + "ange\n", t, i);
1.1248 + if (val[t] == 0.0)
1.1249 + xerror("glp_transform_col: val[%d] = 0; zero coefficient no"
1.1250 + "t allowed\n", t);
1.1251 + if (a[i] != 0.0)
1.1252 + xerror("glp_transform_col: ind[%d] = %d; duplicate row indi"
1.1253 + "ces not allowed\n", t, i);
1.1254 + a[i] = val[t];
1.1255 + }
1.1256 + /* solve the system B*a = alfa to compute the vector alfa */
1.1257 + alfa = a, glp_ftran(P, alfa);
1.1258 + /* store resultant coefficients */
1.1259 + len = 0;
1.1260 + for (i = 1; i <= m; i++)
1.1261 + { if (alfa[i] != 0.0)
1.1262 + { len++;
1.1263 + ind[len] = glp_get_bhead(P, i);
1.1264 + val[len] = alfa[i];
1.1265 + }
1.1266 + }
1.1267 + xfree(a);
1.1268 + return len;
1.1269 +}
1.1270 +
1.1271 +/***********************************************************************
1.1272 +* NAME
1.1273 +*
1.1274 +* glp_prim_rtest - perform primal ratio test
1.1275 +*
1.1276 +* SYNOPSIS
1.1277 +*
1.1278 +* int glp_prim_rtest(glp_prob *P, int len, const int ind[],
1.1279 +* const double val[], int dir, double eps);
1.1280 +*
1.1281 +* DESCRIPTION
1.1282 +*
1.1283 +* The routine glp_prim_rtest performs the primal ratio test using an
1.1284 +* explicitly specified column of the simplex table.
1.1285 +*
1.1286 +* The current basic solution associated with the LP problem object
1.1287 +* must be primal feasible.
1.1288 +*
1.1289 +* The explicitly specified column of the simplex table shows how the
1.1290 +* basic variables xB depend on some non-basic variable x (which is not
1.1291 +* necessarily presented in the problem object):
1.1292 +*
1.1293 +* xB[1] = ... + alfa[1] * x + ...
1.1294 +* xB[2] = ... + alfa[2] * x + ... (*)
1.1295 +* . . . . . . . .
1.1296 +* xB[m] = ... + alfa[m] * x + ...
1.1297 +*
1.1298 +* The column (*) is specifed on entry to the routine using the sparse
1.1299 +* format. Ordinal numbers of basic variables xB[i] should be placed in
1.1300 +* locations ind[1], ..., ind[len], where ordinal number 1 to m denote
1.1301 +* auxiliary variables, and ordinal numbers m+1 to m+n denote structural
1.1302 +* variables. The corresponding non-zero coefficients alfa[i] should be
1.1303 +* placed in locations val[1], ..., val[len]. The arrays ind and val are
1.1304 +* not changed on exit.
1.1305 +*
1.1306 +* The parameter dir specifies direction in which the variable x changes
1.1307 +* on entering the basis: +1 means increasing, -1 means decreasing.
1.1308 +*
1.1309 +* The parameter eps is an absolute tolerance (small positive number)
1.1310 +* used by the routine to skip small alfa[j] of the row (*).
1.1311 +*
1.1312 +* The routine determines which basic variable (among specified in
1.1313 +* ind[1], ..., ind[len]) should leave the basis in order to keep primal
1.1314 +* feasibility.
1.1315 +*
1.1316 +* RETURNS
1.1317 +*
1.1318 +* The routine glp_prim_rtest returns the index piv in the arrays ind
1.1319 +* and val corresponding to the pivot element chosen, 1 <= piv <= len.
1.1320 +* If the adjacent basic solution is primal unbounded and therefore the
1.1321 +* choice cannot be made, the routine returns zero.
1.1322 +*
1.1323 +* COMMENTS
1.1324 +*
1.1325 +* If the non-basic variable x is presented in the LP problem object,
1.1326 +* the column (*) can be computed with the routine glp_eval_tab_col;
1.1327 +* otherwise it can be computed with the routine glp_transform_col. */
1.1328 +
1.1329 +int glp_prim_rtest(glp_prob *P, int len, const int ind[],
1.1330 + const double val[], int dir, double eps)
1.1331 +{ int k, m, n, piv, t, type, stat;
1.1332 + double alfa, big, beta, lb, ub, temp, teta;
1.1333 + if (glp_get_prim_stat(P) != GLP_FEAS)
1.1334 + xerror("glp_prim_rtest: basic solution is not primal feasible "
1.1335 + "\n");
1.1336 + if (!(dir == +1 || dir == -1))
1.1337 + xerror("glp_prim_rtest: dir = %d; invalid parameter\n", dir);
1.1338 + if (!(0.0 < eps && eps < 1.0))
1.1339 + xerror("glp_prim_rtest: eps = %g; invalid parameter\n", eps);
1.1340 + m = glp_get_num_rows(P);
1.1341 + n = glp_get_num_cols(P);
1.1342 + /* initial settings */
1.1343 + piv = 0, teta = DBL_MAX, big = 0.0;
1.1344 + /* walk through the entries of the specified column */
1.1345 + for (t = 1; t <= len; t++)
1.1346 + { /* get the ordinal number of basic variable */
1.1347 + k = ind[t];
1.1348 + if (!(1 <= k && k <= m+n))
1.1349 + xerror("glp_prim_rtest: ind[%d] = %d; variable number out o"
1.1350 + "f range\n", t, k);
1.1351 + /* determine type, bounds, status and primal value of basic
1.1352 + variable xB[i] = x[k] in the current basic solution */
1.1353 + if (k <= m)
1.1354 + { type = glp_get_row_type(P, k);
1.1355 + lb = glp_get_row_lb(P, k);
1.1356 + ub = glp_get_row_ub(P, k);
1.1357 + stat = glp_get_row_stat(P, k);
1.1358 + beta = glp_get_row_prim(P, k);
1.1359 + }
1.1360 + else
1.1361 + { type = glp_get_col_type(P, k-m);
1.1362 + lb = glp_get_col_lb(P, k-m);
1.1363 + ub = glp_get_col_ub(P, k-m);
1.1364 + stat = glp_get_col_stat(P, k-m);
1.1365 + beta = glp_get_col_prim(P, k-m);
1.1366 + }
1.1367 + if (stat != GLP_BS)
1.1368 + xerror("glp_prim_rtest: ind[%d] = %d; non-basic variable no"
1.1369 + "t allowed\n", t, k);
1.1370 + /* determine influence coefficient at basic variable xB[i]
1.1371 + in the explicitly specified column and turn to the case of
1.1372 + increasing the variable x in order to simplify the program
1.1373 + logic */
1.1374 + alfa = (dir > 0 ? + val[t] : - val[t]);
1.1375 + /* analyze main cases */
1.1376 + if (type == GLP_FR)
1.1377 + { /* xB[i] is free variable */
1.1378 + continue;
1.1379 + }
1.1380 + else if (type == GLP_LO)
1.1381 +lo: { /* xB[i] has an lower bound */
1.1382 + if (alfa > - eps) continue;
1.1383 + temp = (lb - beta) / alfa;
1.1384 + }
1.1385 + else if (type == GLP_UP)
1.1386 +up: { /* xB[i] has an upper bound */
1.1387 + if (alfa < + eps) continue;
1.1388 + temp = (ub - beta) / alfa;
1.1389 + }
1.1390 + else if (type == GLP_DB)
1.1391 + { /* xB[i] has both lower and upper bounds */
1.1392 + if (alfa < 0.0) goto lo; else goto up;
1.1393 + }
1.1394 + else if (type == GLP_FX)
1.1395 + { /* xB[i] is fixed variable */
1.1396 + if (- eps < alfa && alfa < + eps) continue;
1.1397 + temp = 0.0;
1.1398 + }
1.1399 + else
1.1400 + xassert(type != type);
1.1401 + /* if the value of the variable xB[i] violates its lower or
1.1402 + upper bound (slightly, because the current basis is assumed
1.1403 + to be primal feasible), temp is negative; we can think this
1.1404 + happens due to round-off errors and the value is exactly on
1.1405 + the bound; this allows replacing temp by zero */
1.1406 + if (temp < 0.0) temp = 0.0;
1.1407 + /* apply the minimal ratio test */
1.1408 + if (teta > temp || teta == temp && big < fabs(alfa))
1.1409 + piv = t, teta = temp, big = fabs(alfa);
1.1410 + }
1.1411 + /* return index of the pivot element chosen */
1.1412 + return piv;
1.1413 +}
1.1414 +
1.1415 +/***********************************************************************
1.1416 +* NAME
1.1417 +*
1.1418 +* glp_dual_rtest - perform dual ratio test
1.1419 +*
1.1420 +* SYNOPSIS
1.1421 +*
1.1422 +* int glp_dual_rtest(glp_prob *P, int len, const int ind[],
1.1423 +* const double val[], int dir, double eps);
1.1424 +*
1.1425 +* DESCRIPTION
1.1426 +*
1.1427 +* The routine glp_dual_rtest performs the dual ratio test using an
1.1428 +* explicitly specified row of the simplex table.
1.1429 +*
1.1430 +* The current basic solution associated with the LP problem object
1.1431 +* must be dual feasible.
1.1432 +*
1.1433 +* The explicitly specified row of the simplex table is a linear form
1.1434 +* that shows how some basic variable x (which is not necessarily
1.1435 +* presented in the problem object) depends on non-basic variables xN:
1.1436 +*
1.1437 +* x = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n]. (*)
1.1438 +*
1.1439 +* The row (*) is specified on entry to the routine using the sparse
1.1440 +* format. Ordinal numbers of non-basic variables xN[j] should be placed
1.1441 +* in locations ind[1], ..., ind[len], where ordinal numbers 1 to m
1.1442 +* denote auxiliary variables, and ordinal numbers m+1 to m+n denote
1.1443 +* structural variables. The corresponding non-zero coefficients alfa[j]
1.1444 +* should be placed in locations val[1], ..., val[len]. The arrays ind
1.1445 +* and val are not changed on exit.
1.1446 +*
1.1447 +* The parameter dir specifies direction in which the variable x changes
1.1448 +* on leaving the basis: +1 means that x goes to its lower bound, and -1
1.1449 +* means that x goes to its upper bound.
1.1450 +*
1.1451 +* The parameter eps is an absolute tolerance (small positive number)
1.1452 +* used by the routine to skip small alfa[j] of the row (*).
1.1453 +*
1.1454 +* The routine determines which non-basic variable (among specified in
1.1455 +* ind[1], ..., ind[len]) should enter the basis in order to keep dual
1.1456 +* feasibility.
1.1457 +*
1.1458 +* RETURNS
1.1459 +*
1.1460 +* The routine glp_dual_rtest returns the index piv in the arrays ind
1.1461 +* and val corresponding to the pivot element chosen, 1 <= piv <= len.
1.1462 +* If the adjacent basic solution is dual unbounded and therefore the
1.1463 +* choice cannot be made, the routine returns zero.
1.1464 +*
1.1465 +* COMMENTS
1.1466 +*
1.1467 +* If the basic variable x is presented in the LP problem object, the
1.1468 +* row (*) can be computed with the routine glp_eval_tab_row; otherwise
1.1469 +* it can be computed with the routine glp_transform_row. */
1.1470 +
1.1471 +int glp_dual_rtest(glp_prob *P, int len, const int ind[],
1.1472 + const double val[], int dir, double eps)
1.1473 +{ int k, m, n, piv, t, stat;
1.1474 + double alfa, big, cost, obj, temp, teta;
1.1475 + if (glp_get_dual_stat(P) != GLP_FEAS)
1.1476 + xerror("glp_dual_rtest: basic solution is not dual feasible\n")
1.1477 + ;
1.1478 + if (!(dir == +1 || dir == -1))
1.1479 + xerror("glp_dual_rtest: dir = %d; invalid parameter\n", dir);
1.1480 + if (!(0.0 < eps && eps < 1.0))
1.1481 + xerror("glp_dual_rtest: eps = %g; invalid parameter\n", eps);
1.1482 + m = glp_get_num_rows(P);
1.1483 + n = glp_get_num_cols(P);
1.1484 + /* take into account optimization direction */
1.1485 + obj = (glp_get_obj_dir(P) == GLP_MIN ? +1.0 : -1.0);
1.1486 + /* initial settings */
1.1487 + piv = 0, teta = DBL_MAX, big = 0.0;
1.1488 + /* walk through the entries of the specified row */
1.1489 + for (t = 1; t <= len; t++)
1.1490 + { /* get ordinal number of non-basic variable */
1.1491 + k = ind[t];
1.1492 + if (!(1 <= k && k <= m+n))
1.1493 + xerror("glp_dual_rtest: ind[%d] = %d; variable number out o"
1.1494 + "f range\n", t, k);
1.1495 + /* determine status and reduced cost of non-basic variable
1.1496 + x[k] = xN[j] in the current basic solution */
1.1497 + if (k <= m)
1.1498 + { stat = glp_get_row_stat(P, k);
1.1499 + cost = glp_get_row_dual(P, k);
1.1500 + }
1.1501 + else
1.1502 + { stat = glp_get_col_stat(P, k-m);
1.1503 + cost = glp_get_col_dual(P, k-m);
1.1504 + }
1.1505 + if (stat == GLP_BS)
1.1506 + xerror("glp_dual_rtest: ind[%d] = %d; basic variable not al"
1.1507 + "lowed\n", t, k);
1.1508 + /* determine influence coefficient at non-basic variable xN[j]
1.1509 + in the explicitly specified row and turn to the case of
1.1510 + increasing the variable x in order to simplify the program
1.1511 + logic */
1.1512 + alfa = (dir > 0 ? + val[t] : - val[t]);
1.1513 + /* analyze main cases */
1.1514 + if (stat == GLP_NL)
1.1515 + { /* xN[j] is on its lower bound */
1.1516 + if (alfa < + eps) continue;
1.1517 + temp = (obj * cost) / alfa;
1.1518 + }
1.1519 + else if (stat == GLP_NU)
1.1520 + { /* xN[j] is on its upper bound */
1.1521 + if (alfa > - eps) continue;
1.1522 + temp = (obj * cost) / alfa;
1.1523 + }
1.1524 + else if (stat == GLP_NF)
1.1525 + { /* xN[j] is non-basic free variable */
1.1526 + if (- eps < alfa && alfa < + eps) continue;
1.1527 + temp = 0.0;
1.1528 + }
1.1529 + else if (stat == GLP_NS)
1.1530 + { /* xN[j] is non-basic fixed variable */
1.1531 + continue;
1.1532 + }
1.1533 + else
1.1534 + xassert(stat != stat);
1.1535 + /* if the reduced cost of the variable xN[j] violates its zero
1.1536 + bound (slightly, because the current basis is assumed to be
1.1537 + dual feasible), temp is negative; we can think this happens
1.1538 + due to round-off errors and the reduced cost is exact zero;
1.1539 + this allows replacing temp by zero */
1.1540 + if (temp < 0.0) temp = 0.0;
1.1541 + /* apply the minimal ratio test */
1.1542 + if (teta > temp || teta == temp && big < fabs(alfa))
1.1543 + piv = t, teta = temp, big = fabs(alfa);
1.1544 + }
1.1545 + /* return index of the pivot element chosen */
1.1546 + return piv;
1.1547 +}
1.1548 +
1.1549 +/***********************************************************************
1.1550 +* NAME
1.1551 +*
1.1552 +* glp_analyze_row - simulate one iteration of dual simplex method
1.1553 +*
1.1554 +* SYNOPSIS
1.1555 +*
1.1556 +* int glp_analyze_row(glp_prob *P, int len, const int ind[],
1.1557 +* const double val[], int type, double rhs, double eps, int *piv,
1.1558 +* double *x, double *dx, double *y, double *dy, double *dz);
1.1559 +*
1.1560 +* DESCRIPTION
1.1561 +*
1.1562 +* Let the current basis be optimal or dual feasible, and there be
1.1563 +* specified a row (constraint), which is violated by the current basic
1.1564 +* solution. The routine glp_analyze_row simulates one iteration of the
1.1565 +* dual simplex method to determine some information on the adjacent
1.1566 +* basis (see below), where the specified row becomes active constraint
1.1567 +* (i.e. its auxiliary variable becomes non-basic).
1.1568 +*
1.1569 +* The current basic solution associated with the problem object passed
1.1570 +* to the routine must be dual feasible, and its primal components must
1.1571 +* be defined.
1.1572 +*
1.1573 +* The row to be analyzed must be previously transformed either with
1.1574 +* the routine glp_eval_tab_row (if the row is in the problem object)
1.1575 +* or with the routine glp_transform_row (if the row is external, i.e.
1.1576 +* not in the problem object). This is needed to express the row only
1.1577 +* through (auxiliary and structural) variables, which are non-basic in
1.1578 +* the current basis:
1.1579 +*
1.1580 +* y = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n],
1.1581 +*
1.1582 +* where y is an auxiliary variable of the row, alfa[j] is an influence
1.1583 +* coefficient, xN[j] is a non-basic variable.
1.1584 +*
1.1585 +* The row is passed to the routine in sparse format. Ordinal numbers
1.1586 +* of non-basic variables are stored in locations ind[1], ..., ind[len],
1.1587 +* where numbers 1 to m denote auxiliary variables while numbers m+1 to
1.1588 +* m+n denote structural variables. Corresponding non-zero coefficients
1.1589 +* alfa[j] are stored in locations val[1], ..., val[len]. The arrays
1.1590 +* ind and val are ot changed on exit.
1.1591 +*
1.1592 +* The parameters type and rhs specify the row type and its right-hand
1.1593 +* side as follows:
1.1594 +*
1.1595 +* type = GLP_LO: y = sum alfa[j] * xN[j] >= rhs
1.1596 +*
1.1597 +* type = GLP_UP: y = sum alfa[j] * xN[j] <= rhs
1.1598 +*
1.1599 +* The parameter eps is an absolute tolerance (small positive number)
1.1600 +* used by the routine to skip small coefficients alfa[j] on performing
1.1601 +* the dual ratio test.
1.1602 +*
1.1603 +* If the operation was successful, the routine stores the following
1.1604 +* information to corresponding location (if some parameter is NULL,
1.1605 +* its value is not stored):
1.1606 +*
1.1607 +* piv index in the array ind and val, 1 <= piv <= len, determining
1.1608 +* the non-basic variable, which would enter the adjacent basis;
1.1609 +*
1.1610 +* x value of the non-basic variable in the current basis;
1.1611 +*
1.1612 +* dx difference between values of the non-basic variable in the
1.1613 +* adjacent and current bases, dx = x.new - x.old;
1.1614 +*
1.1615 +* y value of the row (i.e. of its auxiliary variable) in the
1.1616 +* current basis;
1.1617 +*
1.1618 +* dy difference between values of the row in the adjacent and
1.1619 +* current bases, dy = y.new - y.old;
1.1620 +*
1.1621 +* dz difference between values of the objective function in the
1.1622 +* adjacent and current bases, dz = z.new - z.old. Note that in
1.1623 +* case of minimization dz >= 0, and in case of maximization
1.1624 +* dz <= 0, i.e. in the adjacent basis the objective function
1.1625 +* always gets worse (degrades). */
1.1626 +
1.1627 +int _glp_analyze_row(glp_prob *P, int len, const int ind[],
1.1628 + const double val[], int type, double rhs, double eps, int *_piv,
1.1629 + double *_x, double *_dx, double *_y, double *_dy, double *_dz)
1.1630 +{ int t, k, dir, piv, ret = 0;
1.1631 + double x, dx, y, dy, dz;
1.1632 + if (P->pbs_stat == GLP_UNDEF)
1.1633 + xerror("glp_analyze_row: primal basic solution components are "
1.1634 + "undefined\n");
1.1635 + if (P->dbs_stat != GLP_FEAS)
1.1636 + xerror("glp_analyze_row: basic solution is not dual feasible\n"
1.1637 + );
1.1638 + /* compute the row value y = sum alfa[j] * xN[j] in the current
1.1639 + basis */
1.1640 + if (!(0 <= len && len <= P->n))
1.1641 + xerror("glp_analyze_row: len = %d; invalid row length\n", len);
1.1642 + y = 0.0;
1.1643 + for (t = 1; t <= len; t++)
1.1644 + { /* determine value of x[k] = xN[j] in the current basis */
1.1645 + k = ind[t];
1.1646 + if (!(1 <= k && k <= P->m+P->n))
1.1647 + xerror("glp_analyze_row: ind[%d] = %d; row/column index out"
1.1648 + " of range\n", t, k);
1.1649 + if (k <= P->m)
1.1650 + { /* x[k] is auxiliary variable */
1.1651 + if (P->row[k]->stat == GLP_BS)
1.1652 + xerror("glp_analyze_row: ind[%d] = %d; basic auxiliary v"
1.1653 + "ariable is not allowed\n", t, k);
1.1654 + x = P->row[k]->prim;
1.1655 + }
1.1656 + else
1.1657 + { /* x[k] is structural variable */
1.1658 + if (P->col[k-P->m]->stat == GLP_BS)
1.1659 + xerror("glp_analyze_row: ind[%d] = %d; basic structural "
1.1660 + "variable is not allowed\n", t, k);
1.1661 + x = P->col[k-P->m]->prim;
1.1662 + }
1.1663 + y += val[t] * x;
1.1664 + }
1.1665 + /* check if the row is primal infeasible in the current basis,
1.1666 + i.e. the constraint is violated at the current point */
1.1667 + if (type == GLP_LO)
1.1668 + { if (y >= rhs)
1.1669 + { /* the constraint is not violated */
1.1670 + ret = 1;
1.1671 + goto done;
1.1672 + }
1.1673 + /* in the adjacent basis y goes to its lower bound */
1.1674 + dir = +1;
1.1675 + }
1.1676 + else if (type == GLP_UP)
1.1677 + { if (y <= rhs)
1.1678 + { /* the constraint is not violated */
1.1679 + ret = 1;
1.1680 + goto done;
1.1681 + }
1.1682 + /* in the adjacent basis y goes to its upper bound */
1.1683 + dir = -1;
1.1684 + }
1.1685 + else
1.1686 + xerror("glp_analyze_row: type = %d; invalid parameter\n",
1.1687 + type);
1.1688 + /* compute dy = y.new - y.old */
1.1689 + dy = rhs - y;
1.1690 + /* perform dual ratio test to determine which non-basic variable
1.1691 + should enter the adjacent basis to keep it dual feasible */
1.1692 + piv = glp_dual_rtest(P, len, ind, val, dir, eps);
1.1693 + if (piv == 0)
1.1694 + { /* no dual feasible adjacent basis exists */
1.1695 + ret = 2;
1.1696 + goto done;
1.1697 + }
1.1698 + /* non-basic variable x[k] = xN[j] should enter the basis */
1.1699 + k = ind[piv];
1.1700 + xassert(1 <= k && k <= P->m+P->n);
1.1701 + /* determine its value in the current basis */
1.1702 + if (k <= P->m)
1.1703 + x = P->row[k]->prim;
1.1704 + else
1.1705 + x = P->col[k-P->m]->prim;
1.1706 + /* compute dx = x.new - x.old = dy / alfa[j] */
1.1707 + xassert(val[piv] != 0.0);
1.1708 + dx = dy / val[piv];
1.1709 + /* compute dz = z.new - z.old = d[j] * dx, where d[j] is reduced
1.1710 + cost of xN[j] in the current basis */
1.1711 + if (k <= P->m)
1.1712 + dz = P->row[k]->dual * dx;
1.1713 + else
1.1714 + dz = P->col[k-P->m]->dual * dx;
1.1715 + /* store the analysis results */
1.1716 + if (_piv != NULL) *_piv = piv;
1.1717 + if (_x != NULL) *_x = x;
1.1718 + if (_dx != NULL) *_dx = dx;
1.1719 + if (_y != NULL) *_y = y;
1.1720 + if (_dy != NULL) *_dy = dy;
1.1721 + if (_dz != NULL) *_dz = dz;
1.1722 +done: return ret;
1.1723 +}
1.1724 +
1.1725 +#if 0
1.1726 +int main(void)
1.1727 +{ /* example program for the routine glp_analyze_row */
1.1728 + glp_prob *P;
1.1729 + glp_smcp parm;
1.1730 + int i, k, len, piv, ret, ind[1+100];
1.1731 + double rhs, x, dx, y, dy, dz, val[1+100];
1.1732 + P = glp_create_prob();
1.1733 + /* read plan.mps (see glpk/examples) */
1.1734 + ret = glp_read_mps(P, GLP_MPS_DECK, NULL, "plan.mps");
1.1735 + glp_assert(ret == 0);
1.1736 + /* and solve it to optimality */
1.1737 + ret = glp_simplex(P, NULL);
1.1738 + glp_assert(ret == 0);
1.1739 + glp_assert(glp_get_status(P) == GLP_OPT);
1.1740 + /* the optimal objective value is 296.217 */
1.1741 + /* we would like to know what happens if we would add a new row
1.1742 + (constraint) to plan.mps:
1.1743 + .01 * bin1 + .01 * bin2 + .02 * bin4 + .02 * bin5 <= 12 */
1.1744 + /* first, we specify this new row */
1.1745 + glp_create_index(P);
1.1746 + len = 0;
1.1747 + ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
1.1748 + ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
1.1749 + ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
1.1750 + ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
1.1751 + rhs = 12;
1.1752 + /* then we can compute value of the row (i.e. of its auxiliary
1.1753 + variable) in the current basis to see if the constraint is
1.1754 + violated */
1.1755 + y = 0.0;
1.1756 + for (k = 1; k <= len; k++)
1.1757 + y += val[k] * glp_get_col_prim(P, ind[k]);
1.1758 + glp_printf("y = %g\n", y);
1.1759 + /* this prints y = 15.1372, so the constraint is violated, since
1.1760 + we require that y <= rhs = 12 */
1.1761 + /* now we transform the row to express it only through non-basic
1.1762 + (auxiliary and artificial) variables */
1.1763 + len = glp_transform_row(P, len, ind, val);
1.1764 + /* finally, we simulate one step of the dual simplex method to
1.1765 + obtain necessary information for the adjacent basis */
1.1766 + ret = _glp_analyze_row(P, len, ind, val, GLP_UP, rhs, 1e-9, &piv,
1.1767 + &x, &dx, &y, &dy, &dz);
1.1768 + glp_assert(ret == 0);
1.1769 + glp_printf("k = %d, x = %g; dx = %g; y = %g; dy = %g; dz = %g\n",
1.1770 + ind[piv], x, dx, y, dy, dz);
1.1771 + /* this prints dz = 5.64418 and means that in the adjacent basis
1.1772 + the objective function would be 296.217 + 5.64418 = 301.861 */
1.1773 + /* now we actually include the row into the problem object; note
1.1774 + that the arrays ind and val are clobbered, so we need to build
1.1775 + them once again */
1.1776 + len = 0;
1.1777 + ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
1.1778 + ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
1.1779 + ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
1.1780 + ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
1.1781 + rhs = 12;
1.1782 + i = glp_add_rows(P, 1);
1.1783 + glp_set_row_bnds(P, i, GLP_UP, 0, rhs);
1.1784 + glp_set_mat_row(P, i, len, ind, val);
1.1785 + /* and perform one dual simplex iteration */
1.1786 + glp_init_smcp(&parm);
1.1787 + parm.meth = GLP_DUAL;
1.1788 + parm.it_lim = 1;
1.1789 + glp_simplex(P, &parm);
1.1790 + /* the current objective value is 301.861 */
1.1791 + return 0;
1.1792 +}
1.1793 +#endif
1.1794 +
1.1795 +/***********************************************************************
1.1796 +* NAME
1.1797 +*
1.1798 +* glp_analyze_bound - analyze active bound of non-basic variable
1.1799 +*
1.1800 +* SYNOPSIS
1.1801 +*
1.1802 +* void glp_analyze_bound(glp_prob *P, int k, double *limit1, int *var1,
1.1803 +* double *limit2, int *var2);
1.1804 +*
1.1805 +* DESCRIPTION
1.1806 +*
1.1807 +* The routine glp_analyze_bound analyzes the effect of varying the
1.1808 +* active bound of specified non-basic variable.
1.1809 +*
1.1810 +* The non-basic variable is specified by the parameter k, where
1.1811 +* 1 <= k <= m means auxiliary variable of corresponding row while
1.1812 +* m+1 <= k <= m+n means structural variable (column).
1.1813 +*
1.1814 +* Note that the current basic solution must be optimal, and the basis
1.1815 +* factorization must exist.
1.1816 +*
1.1817 +* Results of the analysis have the following meaning.
1.1818 +*
1.1819 +* value1 is the minimal value of the active bound, at which the basis
1.1820 +* still remains primal feasible and thus optimal. -DBL_MAX means that
1.1821 +* the active bound has no lower limit.
1.1822 +*
1.1823 +* var1 is the ordinal number of an auxiliary (1 to m) or structural
1.1824 +* (m+1 to n) basic variable, which reaches its bound first and thereby
1.1825 +* limits further decreasing the active bound being analyzed.
1.1826 +* if value1 = -DBL_MAX, var1 is set to 0.
1.1827 +*
1.1828 +* value2 is the maximal value of the active bound, at which the basis
1.1829 +* still remains primal feasible and thus optimal. +DBL_MAX means that
1.1830 +* the active bound has no upper limit.
1.1831 +*
1.1832 +* var2 is the ordinal number of an auxiliary (1 to m) or structural
1.1833 +* (m+1 to n) basic variable, which reaches its bound first and thereby
1.1834 +* limits further increasing the active bound being analyzed.
1.1835 +* if value2 = +DBL_MAX, var2 is set to 0. */
1.1836 +
1.1837 +void glp_analyze_bound(glp_prob *P, int k, double *value1, int *var1,
1.1838 + double *value2, int *var2)
1.1839 +{ GLPROW *row;
1.1840 + GLPCOL *col;
1.1841 + int m, n, stat, kase, p, len, piv, *ind;
1.1842 + double x, new_x, ll, uu, xx, delta, *val;
1.1843 + /* sanity checks */
1.1844 + if (P == NULL || P->magic != GLP_PROB_MAGIC)
1.1845 + xerror("glp_analyze_bound: P = %p; invalid problem object\n",
1.1846 + P);
1.1847 + m = P->m, n = P->n;
1.1848 + if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
1.1849 + xerror("glp_analyze_bound: optimal basic solution required\n");
1.1850 + if (!(m == 0 || P->valid))
1.1851 + xerror("glp_analyze_bound: basis factorization required\n");
1.1852 + if (!(1 <= k && k <= m+n))
1.1853 + xerror("glp_analyze_bound: k = %d; variable number out of rang"
1.1854 + "e\n", k);
1.1855 + /* retrieve information about the specified non-basic variable
1.1856 + x[k] whose active bound is to be analyzed */
1.1857 + if (k <= m)
1.1858 + { row = P->row[k];
1.1859 + stat = row->stat;
1.1860 + x = row->prim;
1.1861 + }
1.1862 + else
1.1863 + { col = P->col[k-m];
1.1864 + stat = col->stat;
1.1865 + x = col->prim;
1.1866 + }
1.1867 + if (stat == GLP_BS)
1.1868 + xerror("glp_analyze_bound: k = %d; basic variable not allowed "
1.1869 + "\n", k);
1.1870 + /* allocate working arrays */
1.1871 + ind = xcalloc(1+m, sizeof(int));
1.1872 + val = xcalloc(1+m, sizeof(double));
1.1873 + /* compute column of the simplex table corresponding to the
1.1874 + non-basic variable x[k] */
1.1875 + len = glp_eval_tab_col(P, k, ind, val);
1.1876 + xassert(0 <= len && len <= m);
1.1877 + /* perform analysis */
1.1878 + for (kase = -1; kase <= +1; kase += 2)
1.1879 + { /* kase < 0 means active bound of x[k] is decreasing;
1.1880 + kase > 0 means active bound of x[k] is increasing */
1.1881 + /* use the primal ratio test to determine some basic variable
1.1882 + x[p] which reaches its bound first */
1.1883 + piv = glp_prim_rtest(P, len, ind, val, kase, 1e-9);
1.1884 + if (piv == 0)
1.1885 + { /* nothing limits changing the active bound of x[k] */
1.1886 + p = 0;
1.1887 + new_x = (kase < 0 ? -DBL_MAX : +DBL_MAX);
1.1888 + goto store;
1.1889 + }
1.1890 + /* basic variable x[p] limits changing the active bound of
1.1891 + x[k]; determine its value in the current basis */
1.1892 + xassert(1 <= piv && piv <= len);
1.1893 + p = ind[piv];
1.1894 + if (p <= m)
1.1895 + { row = P->row[p];
1.1896 + ll = glp_get_row_lb(P, row->i);
1.1897 + uu = glp_get_row_ub(P, row->i);
1.1898 + stat = row->stat;
1.1899 + xx = row->prim;
1.1900 + }
1.1901 + else
1.1902 + { col = P->col[p-m];
1.1903 + ll = glp_get_col_lb(P, col->j);
1.1904 + uu = glp_get_col_ub(P, col->j);
1.1905 + stat = col->stat;
1.1906 + xx = col->prim;
1.1907 + }
1.1908 + xassert(stat == GLP_BS);
1.1909 + /* determine delta x[p] = bound of x[p] - value of x[p] */
1.1910 + if (kase < 0 && val[piv] > 0.0 ||
1.1911 + kase > 0 && val[piv] < 0.0)
1.1912 + { /* delta x[p] < 0, so x[p] goes toward its lower bound */
1.1913 + xassert(ll != -DBL_MAX);
1.1914 + delta = ll - xx;
1.1915 + }
1.1916 + else
1.1917 + { /* delta x[p] > 0, so x[p] goes toward its upper bound */
1.1918 + xassert(uu != +DBL_MAX);
1.1919 + delta = uu - xx;
1.1920 + }
1.1921 + /* delta x[p] = alfa[p,k] * delta x[k], so new x[k] = x[k] +
1.1922 + delta x[k] = x[k] + delta x[p] / alfa[p,k] is the value of
1.1923 + x[k] in the adjacent basis */
1.1924 + xassert(val[piv] != 0.0);
1.1925 + new_x = x + delta / val[piv];
1.1926 +store: /* store analysis results */
1.1927 + if (kase < 0)
1.1928 + { if (value1 != NULL) *value1 = new_x;
1.1929 + if (var1 != NULL) *var1 = p;
1.1930 + }
1.1931 + else
1.1932 + { if (value2 != NULL) *value2 = new_x;
1.1933 + if (var2 != NULL) *var2 = p;
1.1934 + }
1.1935 + }
1.1936 + /* free working arrays */
1.1937 + xfree(ind);
1.1938 + xfree(val);
1.1939 + return;
1.1940 +}
1.1941 +
1.1942 +/***********************************************************************
1.1943 +* NAME
1.1944 +*
1.1945 +* glp_analyze_coef - analyze objective coefficient at basic variable
1.1946 +*
1.1947 +* SYNOPSIS
1.1948 +*
1.1949 +* void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
1.1950 +* double *value1, double *coef2, int *var2, double *value2);
1.1951 +*
1.1952 +* DESCRIPTION
1.1953 +*
1.1954 +* The routine glp_analyze_coef analyzes the effect of varying the
1.1955 +* objective coefficient at specified basic variable.
1.1956 +*
1.1957 +* The basic variable is specified by the parameter k, where
1.1958 +* 1 <= k <= m means auxiliary variable of corresponding row while
1.1959 +* m+1 <= k <= m+n means structural variable (column).
1.1960 +*
1.1961 +* Note that the current basic solution must be optimal, and the basis
1.1962 +* factorization must exist.
1.1963 +*
1.1964 +* Results of the analysis have the following meaning.
1.1965 +*
1.1966 +* coef1 is the minimal value of the objective coefficient, at which
1.1967 +* the basis still remains dual feasible and thus optimal. -DBL_MAX
1.1968 +* means that the objective coefficient has no lower limit.
1.1969 +*
1.1970 +* var1 is the ordinal number of an auxiliary (1 to m) or structural
1.1971 +* (m+1 to n) non-basic variable, whose reduced cost reaches its zero
1.1972 +* bound first and thereby limits further decreasing the objective
1.1973 +* coefficient being analyzed. If coef1 = -DBL_MAX, var1 is set to 0.
1.1974 +*
1.1975 +* value1 is value of the basic variable being analyzed in an adjacent
1.1976 +* basis, which is defined as follows. Let the objective coefficient
1.1977 +* reaches its minimal value (coef1) and continues decreasing. Then the
1.1978 +* reduced cost of the limiting non-basic variable (var1) becomes dual
1.1979 +* infeasible and the current basis becomes non-optimal that forces the
1.1980 +* limiting non-basic variable to enter the basis replacing there some
1.1981 +* basic variable that leaves the basis to keep primal feasibility.
1.1982 +* Should note that on determining the adjacent basis current bounds
1.1983 +* of the basic variable being analyzed are ignored as if it were free
1.1984 +* (unbounded) variable, so it cannot leave the basis. It may happen
1.1985 +* that no dual feasible adjacent basis exists, in which case value1 is
1.1986 +* set to -DBL_MAX or +DBL_MAX.
1.1987 +*
1.1988 +* coef2 is the maximal value of the objective coefficient, at which
1.1989 +* the basis still remains dual feasible and thus optimal. +DBL_MAX
1.1990 +* means that the objective coefficient has no upper limit.
1.1991 +*
1.1992 +* var2 is the ordinal number of an auxiliary (1 to m) or structural
1.1993 +* (m+1 to n) non-basic variable, whose reduced cost reaches its zero
1.1994 +* bound first and thereby limits further increasing the objective
1.1995 +* coefficient being analyzed. If coef2 = +DBL_MAX, var2 is set to 0.
1.1996 +*
1.1997 +* value2 is value of the basic variable being analyzed in an adjacent
1.1998 +* basis, which is defined exactly in the same way as value1 above with
1.1999 +* exception that now the objective coefficient is increasing. */
1.2000 +
1.2001 +void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
1.2002 + double *value1, double *coef2, int *var2, double *value2)
1.2003 +{ GLPROW *row; GLPCOL *col;
1.2004 + int m, n, type, stat, kase, p, q, dir, clen, cpiv, rlen, rpiv,
1.2005 + *cind, *rind;
1.2006 + double lb, ub, coef, x, lim_coef, new_x, d, delta, ll, uu, xx,
1.2007 + *rval, *cval;
1.2008 + /* sanity checks */
1.2009 + if (P == NULL || P->magic != GLP_PROB_MAGIC)
1.2010 + xerror("glp_analyze_coef: P = %p; invalid problem object\n",
1.2011 + P);
1.2012 + m = P->m, n = P->n;
1.2013 + if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
1.2014 + xerror("glp_analyze_coef: optimal basic solution required\n");
1.2015 + if (!(m == 0 || P->valid))
1.2016 + xerror("glp_analyze_coef: basis factorization required\n");
1.2017 + if (!(1 <= k && k <= m+n))
1.2018 + xerror("glp_analyze_coef: k = %d; variable number out of range"
1.2019 + "\n", k);
1.2020 + /* retrieve information about the specified basic variable x[k]
1.2021 + whose objective coefficient c[k] is to be analyzed */
1.2022 + if (k <= m)
1.2023 + { row = P->row[k];
1.2024 + type = row->type;
1.2025 + lb = row->lb;
1.2026 + ub = row->ub;
1.2027 + coef = 0.0;
1.2028 + stat = row->stat;
1.2029 + x = row->prim;
1.2030 + }
1.2031 + else
1.2032 + { col = P->col[k-m];
1.2033 + type = col->type;
1.2034 + lb = col->lb;
1.2035 + ub = col->ub;
1.2036 + coef = col->coef;
1.2037 + stat = col->stat;
1.2038 + x = col->prim;
1.2039 + }
1.2040 + if (stat != GLP_BS)
1.2041 + xerror("glp_analyze_coef: k = %d; non-basic variable not allow"
1.2042 + "ed\n", k);
1.2043 + /* allocate working arrays */
1.2044 + cind = xcalloc(1+m, sizeof(int));
1.2045 + cval = xcalloc(1+m, sizeof(double));
1.2046 + rind = xcalloc(1+n, sizeof(int));
1.2047 + rval = xcalloc(1+n, sizeof(double));
1.2048 + /* compute row of the simplex table corresponding to the basic
1.2049 + variable x[k] */
1.2050 + rlen = glp_eval_tab_row(P, k, rind, rval);
1.2051 + xassert(0 <= rlen && rlen <= n);
1.2052 + /* perform analysis */
1.2053 + for (kase = -1; kase <= +1; kase += 2)
1.2054 + { /* kase < 0 means objective coefficient c[k] is decreasing;
1.2055 + kase > 0 means objective coefficient c[k] is increasing */
1.2056 + /* note that decreasing c[k] is equivalent to increasing dual
1.2057 + variable lambda[k] and vice versa; we need to correctly set
1.2058 + the dir flag as required by the routine glp_dual_rtest */
1.2059 + if (P->dir == GLP_MIN)
1.2060 + dir = - kase;
1.2061 + else if (P->dir == GLP_MAX)
1.2062 + dir = + kase;
1.2063 + else
1.2064 + xassert(P != P);
1.2065 + /* use the dual ratio test to determine non-basic variable
1.2066 + x[q] whose reduced cost d[q] reaches zero bound first */
1.2067 + rpiv = glp_dual_rtest(P, rlen, rind, rval, dir, 1e-9);
1.2068 + if (rpiv == 0)
1.2069 + { /* nothing limits changing c[k] */
1.2070 + lim_coef = (kase < 0 ? -DBL_MAX : +DBL_MAX);
1.2071 + q = 0;
1.2072 + /* x[k] keeps its current value */
1.2073 + new_x = x;
1.2074 + goto store;
1.2075 + }
1.2076 + /* non-basic variable x[q] limits changing coefficient c[k];
1.2077 + determine its status and reduced cost d[k] in the current
1.2078 + basis */
1.2079 + xassert(1 <= rpiv && rpiv <= rlen);
1.2080 + q = rind[rpiv];
1.2081 + xassert(1 <= q && q <= m+n);
1.2082 + if (q <= m)
1.2083 + { row = P->row[q];
1.2084 + stat = row->stat;
1.2085 + d = row->dual;
1.2086 + }
1.2087 + else
1.2088 + { col = P->col[q-m];
1.2089 + stat = col->stat;
1.2090 + d = col->dual;
1.2091 + }
1.2092 + /* note that delta d[q] = new d[q] - d[q] = - d[q], because
1.2093 + new d[q] = 0; delta d[q] = alfa[k,q] * delta c[k], so
1.2094 + delta c[k] = delta d[q] / alfa[k,q] = - d[q] / alfa[k,q] */
1.2095 + xassert(rval[rpiv] != 0.0);
1.2096 + delta = - d / rval[rpiv];
1.2097 + /* compute new c[k] = c[k] + delta c[k], which is the limiting
1.2098 + value of the objective coefficient c[k] */
1.2099 + lim_coef = coef + delta;
1.2100 + /* let c[k] continue decreasing/increasing that makes d[q]
1.2101 + dual infeasible and forces x[q] to enter the basis;
1.2102 + to perform the primal ratio test we need to know in which
1.2103 + direction x[q] changes on entering the basis; we determine
1.2104 + that analyzing the sign of delta d[q] (see above), since
1.2105 + d[q] may be close to zero having wrong sign */
1.2106 + /* let, for simplicity, the problem is minimization */
1.2107 + if (kase < 0 && rval[rpiv] > 0.0 ||
1.2108 + kase > 0 && rval[rpiv] < 0.0)
1.2109 + { /* delta d[q] < 0, so d[q] being non-negative will become
1.2110 + negative, so x[q] will increase */
1.2111 + dir = +1;
1.2112 + }
1.2113 + else
1.2114 + { /* delta d[q] > 0, so d[q] being non-positive will become
1.2115 + positive, so x[q] will decrease */
1.2116 + dir = -1;
1.2117 + }
1.2118 + /* if the problem is maximization, correct the direction */
1.2119 + if (P->dir == GLP_MAX) dir = - dir;
1.2120 + /* check that we didn't make a silly mistake */
1.2121 + if (dir > 0)
1.2122 + xassert(stat == GLP_NL || stat == GLP_NF);
1.2123 + else
1.2124 + xassert(stat == GLP_NU || stat == GLP_NF);
1.2125 + /* compute column of the simplex table corresponding to the
1.2126 + non-basic variable x[q] */
1.2127 + clen = glp_eval_tab_col(P, q, cind, cval);
1.2128 + /* make x[k] temporarily free (unbounded) */
1.2129 + if (k <= m)
1.2130 + { row = P->row[k];
1.2131 + row->type = GLP_FR;
1.2132 + row->lb = row->ub = 0.0;
1.2133 + }
1.2134 + else
1.2135 + { col = P->col[k-m];
1.2136 + col->type = GLP_FR;
1.2137 + col->lb = col->ub = 0.0;
1.2138 + }
1.2139 + /* use the primal ratio test to determine some basic variable
1.2140 + which leaves the basis */
1.2141 + cpiv = glp_prim_rtest(P, clen, cind, cval, dir, 1e-9);
1.2142 + /* restore original bounds of the basic variable x[k] */
1.2143 + if (k <= m)
1.2144 + { row = P->row[k];
1.2145 + row->type = type;
1.2146 + row->lb = lb, row->ub = ub;
1.2147 + }
1.2148 + else
1.2149 + { col = P->col[k-m];
1.2150 + col->type = type;
1.2151 + col->lb = lb, col->ub = ub;
1.2152 + }
1.2153 + if (cpiv == 0)
1.2154 + { /* non-basic variable x[q] can change unlimitedly */
1.2155 + if (dir < 0 && rval[rpiv] > 0.0 ||
1.2156 + dir > 0 && rval[rpiv] < 0.0)
1.2157 + { /* delta x[k] = alfa[k,q] * delta x[q] < 0 */
1.2158 + new_x = -DBL_MAX;
1.2159 + }
1.2160 + else
1.2161 + { /* delta x[k] = alfa[k,q] * delta x[q] > 0 */
1.2162 + new_x = +DBL_MAX;
1.2163 + }
1.2164 + goto store;
1.2165 + }
1.2166 + /* some basic variable x[p] limits changing non-basic variable
1.2167 + x[q] in the adjacent basis */
1.2168 + xassert(1 <= cpiv && cpiv <= clen);
1.2169 + p = cind[cpiv];
1.2170 + xassert(1 <= p && p <= m+n);
1.2171 + xassert(p != k);
1.2172 + if (p <= m)
1.2173 + { row = P->row[p];
1.2174 + xassert(row->stat == GLP_BS);
1.2175 + ll = glp_get_row_lb(P, row->i);
1.2176 + uu = glp_get_row_ub(P, row->i);
1.2177 + xx = row->prim;
1.2178 + }
1.2179 + else
1.2180 + { col = P->col[p-m];
1.2181 + xassert(col->stat == GLP_BS);
1.2182 + ll = glp_get_col_lb(P, col->j);
1.2183 + uu = glp_get_col_ub(P, col->j);
1.2184 + xx = col->prim;
1.2185 + }
1.2186 + /* determine delta x[p] = new x[p] - x[p] */
1.2187 + if (dir < 0 && cval[cpiv] > 0.0 ||
1.2188 + dir > 0 && cval[cpiv] < 0.0)
1.2189 + { /* delta x[p] < 0, so x[p] goes toward its lower bound */
1.2190 + xassert(ll != -DBL_MAX);
1.2191 + delta = ll - xx;
1.2192 + }
1.2193 + else
1.2194 + { /* delta x[p] > 0, so x[p] goes toward its upper bound */
1.2195 + xassert(uu != +DBL_MAX);
1.2196 + delta = uu - xx;
1.2197 + }
1.2198 + /* compute new x[k] = x[k] + alfa[k,q] * delta x[q], where
1.2199 + delta x[q] = delta x[p] / alfa[p,q] */
1.2200 + xassert(cval[cpiv] != 0.0);
1.2201 + new_x = x + (rval[rpiv] / cval[cpiv]) * delta;
1.2202 +store: /* store analysis results */
1.2203 + if (kase < 0)
1.2204 + { if (coef1 != NULL) *coef1 = lim_coef;
1.2205 + if (var1 != NULL) *var1 = q;
1.2206 + if (value1 != NULL) *value1 = new_x;
1.2207 + }
1.2208 + else
1.2209 + { if (coef2 != NULL) *coef2 = lim_coef;
1.2210 + if (var2 != NULL) *var2 = q;
1.2211 + if (value2 != NULL) *value2 = new_x;
1.2212 + }
1.2213 + }
1.2214 + /* free working arrays */
1.2215 + xfree(cind);
1.2216 + xfree(cval);
1.2217 + xfree(rind);
1.2218 + xfree(rval);
1.2219 + return;
1.2220 +}
1.2221 +
1.2222 +/* eof */