lemon-project-template-glpk

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

find changesets by author, revision, files, or words in the commit message
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 /* glpapi12.c (basis factorization and simplex tableau routines) */
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 "glpapi.h"
alpar@9 26
alpar@9 27 /***********************************************************************
alpar@9 28 * NAME
alpar@9 29 *
alpar@9 30 * glp_bf_exists - check if the basis factorization exists
alpar@9 31 *
alpar@9 32 * SYNOPSIS
alpar@9 33 *
alpar@9 34 * int glp_bf_exists(glp_prob *lp);
alpar@9 35 *
alpar@9 36 * RETURNS
alpar@9 37 *
alpar@9 38 * If the basis factorization for the current basis associated with
alpar@9 39 * the specified problem object exists and therefore is available for
alpar@9 40 * computations, the routine glp_bf_exists returns non-zero. Otherwise
alpar@9 41 * the routine returns zero. */
alpar@9 42
alpar@9 43 int glp_bf_exists(glp_prob *lp)
alpar@9 44 { int ret;
alpar@9 45 ret = (lp->m == 0 || lp->valid);
alpar@9 46 return ret;
alpar@9 47 }
alpar@9 48
alpar@9 49 /***********************************************************************
alpar@9 50 * NAME
alpar@9 51 *
alpar@9 52 * glp_factorize - compute the basis factorization
alpar@9 53 *
alpar@9 54 * SYNOPSIS
alpar@9 55 *
alpar@9 56 * int glp_factorize(glp_prob *lp);
alpar@9 57 *
alpar@9 58 * DESCRIPTION
alpar@9 59 *
alpar@9 60 * The routine glp_factorize computes the basis factorization for the
alpar@9 61 * current basis associated with the specified problem object.
alpar@9 62 *
alpar@9 63 * RETURNS
alpar@9 64 *
alpar@9 65 * 0 The basis factorization has been successfully computed.
alpar@9 66 *
alpar@9 67 * GLP_EBADB
alpar@9 68 * The basis matrix is invalid, i.e. the number of basic (auxiliary
alpar@9 69 * and structural) variables differs from the number of rows in the
alpar@9 70 * problem object.
alpar@9 71 *
alpar@9 72 * GLP_ESING
alpar@9 73 * The basis matrix is singular within the working precision.
alpar@9 74 *
alpar@9 75 * GLP_ECOND
alpar@9 76 * The basis matrix is ill-conditioned. */
alpar@9 77
alpar@9 78 static int b_col(void *info, int j, int ind[], double val[])
alpar@9 79 { glp_prob *lp = info;
alpar@9 80 int m = lp->m;
alpar@9 81 GLPAIJ *aij;
alpar@9 82 int k, len;
alpar@9 83 xassert(1 <= j && j <= m);
alpar@9 84 /* determine the ordinal number of basic auxiliary or structural
alpar@9 85 variable x[k] corresponding to basic variable xB[j] */
alpar@9 86 k = lp->head[j];
alpar@9 87 /* build j-th column of the basic matrix, which is k-th column of
alpar@9 88 the scaled augmented matrix (I | -R*A*S) */
alpar@9 89 if (k <= m)
alpar@9 90 { /* x[k] is auxiliary variable */
alpar@9 91 len = 1;
alpar@9 92 ind[1] = k;
alpar@9 93 val[1] = 1.0;
alpar@9 94 }
alpar@9 95 else
alpar@9 96 { /* x[k] is structural variable */
alpar@9 97 len = 0;
alpar@9 98 for (aij = lp->col[k-m]->ptr; aij != NULL; aij = aij->c_next)
alpar@9 99 { len++;
alpar@9 100 ind[len] = aij->row->i;
alpar@9 101 val[len] = - aij->row->rii * aij->val * aij->col->sjj;
alpar@9 102 }
alpar@9 103 }
alpar@9 104 return len;
alpar@9 105 }
alpar@9 106
alpar@9 107 static void copy_bfcp(glp_prob *lp);
alpar@9 108
alpar@9 109 int glp_factorize(glp_prob *lp)
alpar@9 110 { int m = lp->m;
alpar@9 111 int n = lp->n;
alpar@9 112 GLPROW **row = lp->row;
alpar@9 113 GLPCOL **col = lp->col;
alpar@9 114 int *head = lp->head;
alpar@9 115 int j, k, stat, ret;
alpar@9 116 /* invalidate the basis factorization */
alpar@9 117 lp->valid = 0;
alpar@9 118 /* build the basis header */
alpar@9 119 j = 0;
alpar@9 120 for (k = 1; k <= m+n; k++)
alpar@9 121 { if (k <= m)
alpar@9 122 { stat = row[k]->stat;
alpar@9 123 row[k]->bind = 0;
alpar@9 124 }
alpar@9 125 else
alpar@9 126 { stat = col[k-m]->stat;
alpar@9 127 col[k-m]->bind = 0;
alpar@9 128 }
alpar@9 129 if (stat == GLP_BS)
alpar@9 130 { j++;
alpar@9 131 if (j > m)
alpar@9 132 { /* too many basic variables */
alpar@9 133 ret = GLP_EBADB;
alpar@9 134 goto fini;
alpar@9 135 }
alpar@9 136 head[j] = k;
alpar@9 137 if (k <= m)
alpar@9 138 row[k]->bind = j;
alpar@9 139 else
alpar@9 140 col[k-m]->bind = j;
alpar@9 141 }
alpar@9 142 }
alpar@9 143 if (j < m)
alpar@9 144 { /* too few basic variables */
alpar@9 145 ret = GLP_EBADB;
alpar@9 146 goto fini;
alpar@9 147 }
alpar@9 148 /* try to factorize the basis matrix */
alpar@9 149 if (m > 0)
alpar@9 150 { if (lp->bfd == NULL)
alpar@9 151 { lp->bfd = bfd_create_it();
alpar@9 152 copy_bfcp(lp);
alpar@9 153 }
alpar@9 154 switch (bfd_factorize(lp->bfd, m, lp->head, b_col, lp))
alpar@9 155 { case 0:
alpar@9 156 /* ok */
alpar@9 157 break;
alpar@9 158 case BFD_ESING:
alpar@9 159 /* singular matrix */
alpar@9 160 ret = GLP_ESING;
alpar@9 161 goto fini;
alpar@9 162 case BFD_ECOND:
alpar@9 163 /* ill-conditioned matrix */
alpar@9 164 ret = GLP_ECOND;
alpar@9 165 goto fini;
alpar@9 166 default:
alpar@9 167 xassert(lp != lp);
alpar@9 168 }
alpar@9 169 lp->valid = 1;
alpar@9 170 }
alpar@9 171 /* factorization successful */
alpar@9 172 ret = 0;
alpar@9 173 fini: /* bring the return code to the calling program */
alpar@9 174 return ret;
alpar@9 175 }
alpar@9 176
alpar@9 177 /***********************************************************************
alpar@9 178 * NAME
alpar@9 179 *
alpar@9 180 * glp_bf_updated - check if the basis factorization has been updated
alpar@9 181 *
alpar@9 182 * SYNOPSIS
alpar@9 183 *
alpar@9 184 * int glp_bf_updated(glp_prob *lp);
alpar@9 185 *
alpar@9 186 * RETURNS
alpar@9 187 *
alpar@9 188 * If the basis factorization has been just computed from scratch, the
alpar@9 189 * routine glp_bf_updated returns zero. Otherwise, if the factorization
alpar@9 190 * has been updated one or more times, the routine returns non-zero. */
alpar@9 191
alpar@9 192 int glp_bf_updated(glp_prob *lp)
alpar@9 193 { int cnt;
alpar@9 194 if (!(lp->m == 0 || lp->valid))
alpar@9 195 xerror("glp_bf_update: basis factorization does not exist\n");
alpar@9 196 #if 0 /* 15/XI-2009 */
alpar@9 197 cnt = (lp->m == 0 ? 0 : lp->bfd->upd_cnt);
alpar@9 198 #else
alpar@9 199 cnt = (lp->m == 0 ? 0 : bfd_get_count(lp->bfd));
alpar@9 200 #endif
alpar@9 201 return cnt;
alpar@9 202 }
alpar@9 203
alpar@9 204 /***********************************************************************
alpar@9 205 * NAME
alpar@9 206 *
alpar@9 207 * glp_get_bfcp - retrieve basis factorization control parameters
alpar@9 208 *
alpar@9 209 * SYNOPSIS
alpar@9 210 *
alpar@9 211 * void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm);
alpar@9 212 *
alpar@9 213 * DESCRIPTION
alpar@9 214 *
alpar@9 215 * The routine glp_get_bfcp retrieves control parameters, which are
alpar@9 216 * used on computing and updating the basis factorization associated
alpar@9 217 * with the specified problem object.
alpar@9 218 *
alpar@9 219 * Current values of control parameters are stored by the routine in
alpar@9 220 * a glp_bfcp structure, which the parameter parm points to. */
alpar@9 221
alpar@9 222 void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm)
alpar@9 223 { glp_bfcp *bfcp = lp->bfcp;
alpar@9 224 if (bfcp == NULL)
alpar@9 225 { parm->type = GLP_BF_FT;
alpar@9 226 parm->lu_size = 0;
alpar@9 227 parm->piv_tol = 0.10;
alpar@9 228 parm->piv_lim = 4;
alpar@9 229 parm->suhl = GLP_ON;
alpar@9 230 parm->eps_tol = 1e-15;
alpar@9 231 parm->max_gro = 1e+10;
alpar@9 232 parm->nfs_max = 100;
alpar@9 233 parm->upd_tol = 1e-6;
alpar@9 234 parm->nrs_max = 100;
alpar@9 235 parm->rs_size = 0;
alpar@9 236 }
alpar@9 237 else
alpar@9 238 memcpy(parm, bfcp, sizeof(glp_bfcp));
alpar@9 239 return;
alpar@9 240 }
alpar@9 241
alpar@9 242 /***********************************************************************
alpar@9 243 * NAME
alpar@9 244 *
alpar@9 245 * glp_set_bfcp - change basis factorization control parameters
alpar@9 246 *
alpar@9 247 * SYNOPSIS
alpar@9 248 *
alpar@9 249 * void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm);
alpar@9 250 *
alpar@9 251 * DESCRIPTION
alpar@9 252 *
alpar@9 253 * The routine glp_set_bfcp changes control parameters, which are used
alpar@9 254 * by internal GLPK routines in computing and updating the basis
alpar@9 255 * factorization associated with the specified problem object.
alpar@9 256 *
alpar@9 257 * New values of the control parameters should be passed in a structure
alpar@9 258 * glp_bfcp, which the parameter parm points to.
alpar@9 259 *
alpar@9 260 * The parameter parm can be specified as NULL, in which case all
alpar@9 261 * control parameters are reset to their default values. */
alpar@9 262
alpar@9 263 #if 0 /* 15/XI-2009 */
alpar@9 264 static void copy_bfcp(glp_prob *lp)
alpar@9 265 { glp_bfcp _parm, *parm = &_parm;
alpar@9 266 BFD *bfd = lp->bfd;
alpar@9 267 glp_get_bfcp(lp, parm);
alpar@9 268 xassert(bfd != NULL);
alpar@9 269 bfd->type = parm->type;
alpar@9 270 bfd->lu_size = parm->lu_size;
alpar@9 271 bfd->piv_tol = parm->piv_tol;
alpar@9 272 bfd->piv_lim = parm->piv_lim;
alpar@9 273 bfd->suhl = parm->suhl;
alpar@9 274 bfd->eps_tol = parm->eps_tol;
alpar@9 275 bfd->max_gro = parm->max_gro;
alpar@9 276 bfd->nfs_max = parm->nfs_max;
alpar@9 277 bfd->upd_tol = parm->upd_tol;
alpar@9 278 bfd->nrs_max = parm->nrs_max;
alpar@9 279 bfd->rs_size = parm->rs_size;
alpar@9 280 return;
alpar@9 281 }
alpar@9 282 #else
alpar@9 283 static void copy_bfcp(glp_prob *lp)
alpar@9 284 { glp_bfcp _parm, *parm = &_parm;
alpar@9 285 glp_get_bfcp(lp, parm);
alpar@9 286 bfd_set_parm(lp->bfd, parm);
alpar@9 287 return;
alpar@9 288 }
alpar@9 289 #endif
alpar@9 290
alpar@9 291 void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm)
alpar@9 292 { glp_bfcp *bfcp = lp->bfcp;
alpar@9 293 if (parm == NULL)
alpar@9 294 { /* reset to default values */
alpar@9 295 if (bfcp != NULL)
alpar@9 296 xfree(bfcp), lp->bfcp = NULL;
alpar@9 297 }
alpar@9 298 else
alpar@9 299 { /* set to specified values */
alpar@9 300 if (bfcp == NULL)
alpar@9 301 bfcp = lp->bfcp = xmalloc(sizeof(glp_bfcp));
alpar@9 302 memcpy(bfcp, parm, sizeof(glp_bfcp));
alpar@9 303 if (!(bfcp->type == GLP_BF_FT || bfcp->type == GLP_BF_BG ||
alpar@9 304 bfcp->type == GLP_BF_GR))
alpar@9 305 xerror("glp_set_bfcp: type = %d; invalid parameter\n",
alpar@9 306 bfcp->type);
alpar@9 307 if (bfcp->lu_size < 0)
alpar@9 308 xerror("glp_set_bfcp: lu_size = %d; invalid parameter\n",
alpar@9 309 bfcp->lu_size);
alpar@9 310 if (!(0.0 < bfcp->piv_tol && bfcp->piv_tol < 1.0))
alpar@9 311 xerror("glp_set_bfcp: piv_tol = %g; invalid parameter\n",
alpar@9 312 bfcp->piv_tol);
alpar@9 313 if (bfcp->piv_lim < 1)
alpar@9 314 xerror("glp_set_bfcp: piv_lim = %d; invalid parameter\n",
alpar@9 315 bfcp->piv_lim);
alpar@9 316 if (!(bfcp->suhl == GLP_ON || bfcp->suhl == GLP_OFF))
alpar@9 317 xerror("glp_set_bfcp: suhl = %d; invalid parameter\n",
alpar@9 318 bfcp->suhl);
alpar@9 319 if (!(0.0 <= bfcp->eps_tol && bfcp->eps_tol <= 1e-6))
alpar@9 320 xerror("glp_set_bfcp: eps_tol = %g; invalid parameter\n",
alpar@9 321 bfcp->eps_tol);
alpar@9 322 if (bfcp->max_gro < 1.0)
alpar@9 323 xerror("glp_set_bfcp: max_gro = %g; invalid parameter\n",
alpar@9 324 bfcp->max_gro);
alpar@9 325 if (!(1 <= bfcp->nfs_max && bfcp->nfs_max <= 32767))
alpar@9 326 xerror("glp_set_bfcp: nfs_max = %d; invalid parameter\n",
alpar@9 327 bfcp->nfs_max);
alpar@9 328 if (!(0.0 < bfcp->upd_tol && bfcp->upd_tol < 1.0))
alpar@9 329 xerror("glp_set_bfcp: upd_tol = %g; invalid parameter\n",
alpar@9 330 bfcp->upd_tol);
alpar@9 331 if (!(1 <= bfcp->nrs_max && bfcp->nrs_max <= 32767))
alpar@9 332 xerror("glp_set_bfcp: nrs_max = %d; invalid parameter\n",
alpar@9 333 bfcp->nrs_max);
alpar@9 334 if (bfcp->rs_size < 0)
alpar@9 335 xerror("glp_set_bfcp: rs_size = %d; invalid parameter\n",
alpar@9 336 bfcp->nrs_max);
alpar@9 337 if (bfcp->rs_size == 0)
alpar@9 338 bfcp->rs_size = 20 * bfcp->nrs_max;
alpar@9 339 }
alpar@9 340 if (lp->bfd != NULL) copy_bfcp(lp);
alpar@9 341 return;
alpar@9 342 }
alpar@9 343
alpar@9 344 /***********************************************************************
alpar@9 345 * NAME
alpar@9 346 *
alpar@9 347 * glp_get_bhead - retrieve the basis header information
alpar@9 348 *
alpar@9 349 * SYNOPSIS
alpar@9 350 *
alpar@9 351 * int glp_get_bhead(glp_prob *lp, int k);
alpar@9 352 *
alpar@9 353 * DESCRIPTION
alpar@9 354 *
alpar@9 355 * The routine glp_get_bhead returns the basis header information for
alpar@9 356 * the current basis associated with the specified problem object.
alpar@9 357 *
alpar@9 358 * RETURNS
alpar@9 359 *
alpar@9 360 * If xB[k], 1 <= k <= m, is i-th auxiliary variable (1 <= i <= m), the
alpar@9 361 * routine returns i. Otherwise, if xB[k] is j-th structural variable
alpar@9 362 * (1 <= j <= n), the routine returns m+j. Here m is the number of rows
alpar@9 363 * and n is the number of columns in the problem object. */
alpar@9 364
alpar@9 365 int glp_get_bhead(glp_prob *lp, int k)
alpar@9 366 { if (!(lp->m == 0 || lp->valid))
alpar@9 367 xerror("glp_get_bhead: basis factorization does not exist\n");
alpar@9 368 if (!(1 <= k && k <= lp->m))
alpar@9 369 xerror("glp_get_bhead: k = %d; index out of range\n", k);
alpar@9 370 return lp->head[k];
alpar@9 371 }
alpar@9 372
alpar@9 373 /***********************************************************************
alpar@9 374 * NAME
alpar@9 375 *
alpar@9 376 * glp_get_row_bind - retrieve row index in the basis header
alpar@9 377 *
alpar@9 378 * SYNOPSIS
alpar@9 379 *
alpar@9 380 * int glp_get_row_bind(glp_prob *lp, int i);
alpar@9 381 *
alpar@9 382 * RETURNS
alpar@9 383 *
alpar@9 384 * The routine glp_get_row_bind returns the index k of basic variable
alpar@9 385 * xB[k], 1 <= k <= m, which is i-th auxiliary variable, 1 <= i <= m,
alpar@9 386 * in the current basis associated with the specified problem object,
alpar@9 387 * where m is the number of rows. However, if i-th auxiliary variable
alpar@9 388 * is non-basic, the routine returns zero. */
alpar@9 389
alpar@9 390 int glp_get_row_bind(glp_prob *lp, int i)
alpar@9 391 { if (!(lp->m == 0 || lp->valid))
alpar@9 392 xerror("glp_get_row_bind: basis factorization does not exist\n"
alpar@9 393 );
alpar@9 394 if (!(1 <= i && i <= lp->m))
alpar@9 395 xerror("glp_get_row_bind: i = %d; row number out of range\n",
alpar@9 396 i);
alpar@9 397 return lp->row[i]->bind;
alpar@9 398 }
alpar@9 399
alpar@9 400 /***********************************************************************
alpar@9 401 * NAME
alpar@9 402 *
alpar@9 403 * glp_get_col_bind - retrieve column index in the basis header
alpar@9 404 *
alpar@9 405 * SYNOPSIS
alpar@9 406 *
alpar@9 407 * int glp_get_col_bind(glp_prob *lp, int j);
alpar@9 408 *
alpar@9 409 * RETURNS
alpar@9 410 *
alpar@9 411 * The routine glp_get_col_bind returns the index k of basic variable
alpar@9 412 * xB[k], 1 <= k <= m, which is j-th structural variable, 1 <= j <= n,
alpar@9 413 * in the current basis associated with the specified problem object,
alpar@9 414 * where m is the number of rows, n is the number of columns. However,
alpar@9 415 * if j-th structural variable is non-basic, the routine returns zero.*/
alpar@9 416
alpar@9 417 int glp_get_col_bind(glp_prob *lp, int j)
alpar@9 418 { if (!(lp->m == 0 || lp->valid))
alpar@9 419 xerror("glp_get_col_bind: basis factorization does not exist\n"
alpar@9 420 );
alpar@9 421 if (!(1 <= j && j <= lp->n))
alpar@9 422 xerror("glp_get_col_bind: j = %d; column number out of range\n"
alpar@9 423 , j);
alpar@9 424 return lp->col[j]->bind;
alpar@9 425 }
alpar@9 426
alpar@9 427 /***********************************************************************
alpar@9 428 * NAME
alpar@9 429 *
alpar@9 430 * glp_ftran - perform forward transformation (solve system B*x = b)
alpar@9 431 *
alpar@9 432 * SYNOPSIS
alpar@9 433 *
alpar@9 434 * void glp_ftran(glp_prob *lp, double x[]);
alpar@9 435 *
alpar@9 436 * DESCRIPTION
alpar@9 437 *
alpar@9 438 * The routine glp_ftran performs forward transformation, i.e. solves
alpar@9 439 * the system B*x = b, where B is the basis matrix corresponding to the
alpar@9 440 * current basis for the specified problem object, x is the vector of
alpar@9 441 * unknowns to be computed, b is the vector of right-hand sides.
alpar@9 442 *
alpar@9 443 * On entry elements of the vector b should be stored in dense format
alpar@9 444 * in locations x[1], ..., x[m], where m is the number of rows. On exit
alpar@9 445 * the routine stores elements of the vector x in the same locations.
alpar@9 446 *
alpar@9 447 * SCALING/UNSCALING
alpar@9 448 *
alpar@9 449 * Let A~ = (I | -A) is the augmented constraint matrix of the original
alpar@9 450 * (unscaled) problem. In the scaled LP problem instead the matrix A the
alpar@9 451 * scaled matrix A" = R*A*S is actually used, so
alpar@9 452 *
alpar@9 453 * A~" = (I | A") = (I | R*A*S) = (R*I*inv(R) | R*A*S) =
alpar@9 454 * (1)
alpar@9 455 * = R*(I | A)*S~ = R*A~*S~,
alpar@9 456 *
alpar@9 457 * is the scaled augmented constraint matrix, where R and S are diagonal
alpar@9 458 * scaling matrices used to scale rows and columns of the matrix A, and
alpar@9 459 *
alpar@9 460 * S~ = diag(inv(R) | S) (2)
alpar@9 461 *
alpar@9 462 * is an augmented diagonal scaling matrix.
alpar@9 463 *
alpar@9 464 * By definition:
alpar@9 465 *
alpar@9 466 * A~ = (B | N), (3)
alpar@9 467 *
alpar@9 468 * where B is the basic matrix, which consists of basic columns of the
alpar@9 469 * augmented constraint matrix A~, and N is a matrix, which consists of
alpar@9 470 * non-basic columns of A~. From (1) it follows that:
alpar@9 471 *
alpar@9 472 * A~" = (B" | N") = (R*B*SB | R*N*SN), (4)
alpar@9 473 *
alpar@9 474 * where SB and SN are parts of the augmented scaling matrix S~, which
alpar@9 475 * correspond to basic and non-basic variables, respectively. Therefore
alpar@9 476 *
alpar@9 477 * B" = R*B*SB, (5)
alpar@9 478 *
alpar@9 479 * which is the scaled basis matrix. */
alpar@9 480
alpar@9 481 void glp_ftran(glp_prob *lp, double x[])
alpar@9 482 { int m = lp->m;
alpar@9 483 GLPROW **row = lp->row;
alpar@9 484 GLPCOL **col = lp->col;
alpar@9 485 int i, k;
alpar@9 486 /* B*x = b ===> (R*B*SB)*(inv(SB)*x) = R*b ===>
alpar@9 487 B"*x" = b", where b" = R*b, x = SB*x" */
alpar@9 488 if (!(m == 0 || lp->valid))
alpar@9 489 xerror("glp_ftran: basis factorization does not exist\n");
alpar@9 490 /* b" := R*b */
alpar@9 491 for (i = 1; i <= m; i++)
alpar@9 492 x[i] *= row[i]->rii;
alpar@9 493 /* x" := inv(B")*b" */
alpar@9 494 if (m > 0) bfd_ftran(lp->bfd, x);
alpar@9 495 /* x := SB*x" */
alpar@9 496 for (i = 1; i <= m; i++)
alpar@9 497 { k = lp->head[i];
alpar@9 498 if (k <= m)
alpar@9 499 x[i] /= row[k]->rii;
alpar@9 500 else
alpar@9 501 x[i] *= col[k-m]->sjj;
alpar@9 502 }
alpar@9 503 return;
alpar@9 504 }
alpar@9 505
alpar@9 506 /***********************************************************************
alpar@9 507 * NAME
alpar@9 508 *
alpar@9 509 * glp_btran - perform backward transformation (solve system B'*x = b)
alpar@9 510 *
alpar@9 511 * SYNOPSIS
alpar@9 512 *
alpar@9 513 * void glp_btran(glp_prob *lp, double x[]);
alpar@9 514 *
alpar@9 515 * DESCRIPTION
alpar@9 516 *
alpar@9 517 * The routine glp_btran performs backward transformation, i.e. solves
alpar@9 518 * the system B'*x = b, where B' is a matrix transposed to the basis
alpar@9 519 * matrix corresponding to the current basis for the specified problem
alpar@9 520 * problem object, x is the vector of unknowns to be computed, b is the
alpar@9 521 * vector of right-hand sides.
alpar@9 522 *
alpar@9 523 * On entry elements of the vector b should be stored in dense format
alpar@9 524 * in locations x[1], ..., x[m], where m is the number of rows. On exit
alpar@9 525 * the routine stores elements of the vector x in the same locations.
alpar@9 526 *
alpar@9 527 * SCALING/UNSCALING
alpar@9 528 *
alpar@9 529 * See comments to the routine glp_ftran. */
alpar@9 530
alpar@9 531 void glp_btran(glp_prob *lp, double x[])
alpar@9 532 { int m = lp->m;
alpar@9 533 GLPROW **row = lp->row;
alpar@9 534 GLPCOL **col = lp->col;
alpar@9 535 int i, k;
alpar@9 536 /* B'*x = b ===> (SB*B'*R)*(inv(R)*x) = SB*b ===>
alpar@9 537 (B")'*x" = b", where b" = SB*b, x = R*x" */
alpar@9 538 if (!(m == 0 || lp->valid))
alpar@9 539 xerror("glp_btran: basis factorization does not exist\n");
alpar@9 540 /* b" := SB*b */
alpar@9 541 for (i = 1; i <= m; i++)
alpar@9 542 { k = lp->head[i];
alpar@9 543 if (k <= m)
alpar@9 544 x[i] /= row[k]->rii;
alpar@9 545 else
alpar@9 546 x[i] *= col[k-m]->sjj;
alpar@9 547 }
alpar@9 548 /* x" := inv[(B")']*b" */
alpar@9 549 if (m > 0) bfd_btran(lp->bfd, x);
alpar@9 550 /* x := R*x" */
alpar@9 551 for (i = 1; i <= m; i++)
alpar@9 552 x[i] *= row[i]->rii;
alpar@9 553 return;
alpar@9 554 }
alpar@9 555
alpar@9 556 /***********************************************************************
alpar@9 557 * NAME
alpar@9 558 *
alpar@9 559 * glp_warm_up - "warm up" LP basis
alpar@9 560 *
alpar@9 561 * SYNOPSIS
alpar@9 562 *
alpar@9 563 * int glp_warm_up(glp_prob *P);
alpar@9 564 *
alpar@9 565 * DESCRIPTION
alpar@9 566 *
alpar@9 567 * The routine glp_warm_up "warms up" the LP basis for the specified
alpar@9 568 * problem object using current statuses assigned to rows and columns
alpar@9 569 * (that is, to auxiliary and structural variables).
alpar@9 570 *
alpar@9 571 * This operation includes computing factorization of the basis matrix
alpar@9 572 * (if it does not exist), computing primal and dual components of basic
alpar@9 573 * solution, and determining the solution status.
alpar@9 574 *
alpar@9 575 * RETURNS
alpar@9 576 *
alpar@9 577 * 0 The operation has been successfully performed.
alpar@9 578 *
alpar@9 579 * GLP_EBADB
alpar@9 580 * The basis matrix is invalid, i.e. the number of basic (auxiliary
alpar@9 581 * and structural) variables differs from the number of rows in the
alpar@9 582 * problem object.
alpar@9 583 *
alpar@9 584 * GLP_ESING
alpar@9 585 * The basis matrix is singular within the working precision.
alpar@9 586 *
alpar@9 587 * GLP_ECOND
alpar@9 588 * The basis matrix is ill-conditioned. */
alpar@9 589
alpar@9 590 int glp_warm_up(glp_prob *P)
alpar@9 591 { GLPROW *row;
alpar@9 592 GLPCOL *col;
alpar@9 593 GLPAIJ *aij;
alpar@9 594 int i, j, type, ret;
alpar@9 595 double eps, temp, *work;
alpar@9 596 /* invalidate basic solution */
alpar@9 597 P->pbs_stat = P->dbs_stat = GLP_UNDEF;
alpar@9 598 P->obj_val = 0.0;
alpar@9 599 P->some = 0;
alpar@9 600 for (i = 1; i <= P->m; i++)
alpar@9 601 { row = P->row[i];
alpar@9 602 row->prim = row->dual = 0.0;
alpar@9 603 }
alpar@9 604 for (j = 1; j <= P->n; j++)
alpar@9 605 { col = P->col[j];
alpar@9 606 col->prim = col->dual = 0.0;
alpar@9 607 }
alpar@9 608 /* compute the basis factorization, if necessary */
alpar@9 609 if (!glp_bf_exists(P))
alpar@9 610 { ret = glp_factorize(P);
alpar@9 611 if (ret != 0) goto done;
alpar@9 612 }
alpar@9 613 /* allocate working array */
alpar@9 614 work = xcalloc(1+P->m, sizeof(double));
alpar@9 615 /* determine and store values of non-basic variables, compute
alpar@9 616 vector (- N * xN) */
alpar@9 617 for (i = 1; i <= P->m; i++)
alpar@9 618 work[i] = 0.0;
alpar@9 619 for (i = 1; i <= P->m; i++)
alpar@9 620 { row = P->row[i];
alpar@9 621 if (row->stat == GLP_BS)
alpar@9 622 continue;
alpar@9 623 else if (row->stat == GLP_NL)
alpar@9 624 row->prim = row->lb;
alpar@9 625 else if (row->stat == GLP_NU)
alpar@9 626 row->prim = row->ub;
alpar@9 627 else if (row->stat == GLP_NF)
alpar@9 628 row->prim = 0.0;
alpar@9 629 else if (row->stat == GLP_NS)
alpar@9 630 row->prim = row->lb;
alpar@9 631 else
alpar@9 632 xassert(row != row);
alpar@9 633 /* N[j] is i-th column of matrix (I|-A) */
alpar@9 634 work[i] -= row->prim;
alpar@9 635 }
alpar@9 636 for (j = 1; j <= P->n; j++)
alpar@9 637 { col = P->col[j];
alpar@9 638 if (col->stat == GLP_BS)
alpar@9 639 continue;
alpar@9 640 else if (col->stat == GLP_NL)
alpar@9 641 col->prim = col->lb;
alpar@9 642 else if (col->stat == GLP_NU)
alpar@9 643 col->prim = col->ub;
alpar@9 644 else if (col->stat == GLP_NF)
alpar@9 645 col->prim = 0.0;
alpar@9 646 else if (col->stat == GLP_NS)
alpar@9 647 col->prim = col->lb;
alpar@9 648 else
alpar@9 649 xassert(col != col);
alpar@9 650 /* N[j] is (m+j)-th column of matrix (I|-A) */
alpar@9 651 if (col->prim != 0.0)
alpar@9 652 { for (aij = col->ptr; aij != NULL; aij = aij->c_next)
alpar@9 653 work[aij->row->i] += aij->val * col->prim;
alpar@9 654 }
alpar@9 655 }
alpar@9 656 /* compute vector of basic variables xB = - inv(B) * N * xN */
alpar@9 657 glp_ftran(P, work);
alpar@9 658 /* store values of basic variables, check primal feasibility */
alpar@9 659 P->pbs_stat = GLP_FEAS;
alpar@9 660 for (i = 1; i <= P->m; i++)
alpar@9 661 { row = P->row[i];
alpar@9 662 if (row->stat != GLP_BS)
alpar@9 663 continue;
alpar@9 664 row->prim = work[row->bind];
alpar@9 665 type = row->type;
alpar@9 666 if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
alpar@9 667 { eps = 1e-6 + 1e-9 * fabs(row->lb);
alpar@9 668 if (row->prim < row->lb - eps)
alpar@9 669 P->pbs_stat = GLP_INFEAS;
alpar@9 670 }
alpar@9 671 if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
alpar@9 672 { eps = 1e-6 + 1e-9 * fabs(row->ub);
alpar@9 673 if (row->prim > row->ub + eps)
alpar@9 674 P->pbs_stat = GLP_INFEAS;
alpar@9 675 }
alpar@9 676 }
alpar@9 677 for (j = 1; j <= P->n; j++)
alpar@9 678 { col = P->col[j];
alpar@9 679 if (col->stat != GLP_BS)
alpar@9 680 continue;
alpar@9 681 col->prim = work[col->bind];
alpar@9 682 type = col->type;
alpar@9 683 if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
alpar@9 684 { eps = 1e-6 + 1e-9 * fabs(col->lb);
alpar@9 685 if (col->prim < col->lb - eps)
alpar@9 686 P->pbs_stat = GLP_INFEAS;
alpar@9 687 }
alpar@9 688 if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
alpar@9 689 { eps = 1e-6 + 1e-9 * fabs(col->ub);
alpar@9 690 if (col->prim > col->ub + eps)
alpar@9 691 P->pbs_stat = GLP_INFEAS;
alpar@9 692 }
alpar@9 693 }
alpar@9 694 /* compute value of the objective function */
alpar@9 695 P->obj_val = P->c0;
alpar@9 696 for (j = 1; j <= P->n; j++)
alpar@9 697 { col = P->col[j];
alpar@9 698 P->obj_val += col->coef * col->prim;
alpar@9 699 }
alpar@9 700 /* build vector cB of objective coefficients at basic variables */
alpar@9 701 for (i = 1; i <= P->m; i++)
alpar@9 702 work[i] = 0.0;
alpar@9 703 for (j = 1; j <= P->n; j++)
alpar@9 704 { col = P->col[j];
alpar@9 705 if (col->stat == GLP_BS)
alpar@9 706 work[col->bind] = col->coef;
alpar@9 707 }
alpar@9 708 /* compute vector of simplex multipliers pi = inv(B') * cB */
alpar@9 709 glp_btran(P, work);
alpar@9 710 /* compute and store reduced costs of non-basic variables d[j] =
alpar@9 711 c[j] - N'[j] * pi, check dual feasibility */
alpar@9 712 P->dbs_stat = GLP_FEAS;
alpar@9 713 for (i = 1; i <= P->m; i++)
alpar@9 714 { row = P->row[i];
alpar@9 715 if (row->stat == GLP_BS)
alpar@9 716 { row->dual = 0.0;
alpar@9 717 continue;
alpar@9 718 }
alpar@9 719 /* N[j] is i-th column of matrix (I|-A) */
alpar@9 720 row->dual = - work[i];
alpar@9 721 type = row->type;
alpar@9 722 temp = (P->dir == GLP_MIN ? + row->dual : - row->dual);
alpar@9 723 if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
alpar@9 724 (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
alpar@9 725 P->dbs_stat = GLP_INFEAS;
alpar@9 726 }
alpar@9 727 for (j = 1; j <= P->n; j++)
alpar@9 728 { col = P->col[j];
alpar@9 729 if (col->stat == GLP_BS)
alpar@9 730 { col->dual = 0.0;
alpar@9 731 continue;
alpar@9 732 }
alpar@9 733 /* N[j] is (m+j)-th column of matrix (I|-A) */
alpar@9 734 col->dual = col->coef;
alpar@9 735 for (aij = col->ptr; aij != NULL; aij = aij->c_next)
alpar@9 736 col->dual += aij->val * work[aij->row->i];
alpar@9 737 type = col->type;
alpar@9 738 temp = (P->dir == GLP_MIN ? + col->dual : - col->dual);
alpar@9 739 if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
alpar@9 740 (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
alpar@9 741 P->dbs_stat = GLP_INFEAS;
alpar@9 742 }
alpar@9 743 /* free working array */
alpar@9 744 xfree(work);
alpar@9 745 ret = 0;
alpar@9 746 done: return ret;
alpar@9 747 }
alpar@9 748
alpar@9 749 /***********************************************************************
alpar@9 750 * NAME
alpar@9 751 *
alpar@9 752 * glp_eval_tab_row - compute row of the simplex tableau
alpar@9 753 *
alpar@9 754 * SYNOPSIS
alpar@9 755 *
alpar@9 756 * int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[]);
alpar@9 757 *
alpar@9 758 * DESCRIPTION
alpar@9 759 *
alpar@9 760 * The routine glp_eval_tab_row computes a row of the current simplex
alpar@9 761 * tableau for the basic variable, which is specified by the number k:
alpar@9 762 * if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
alpar@9 763 * x[k] is (k-m)-th structural variable, where m is number of rows, and
alpar@9 764 * n is number of columns. The current basis must be available.
alpar@9 765 *
alpar@9 766 * The routine stores column indices and numerical values of non-zero
alpar@9 767 * elements of the computed row using sparse format to the locations
alpar@9 768 * ind[1], ..., ind[len] and val[1], ..., val[len], respectively, where
alpar@9 769 * 0 <= len <= n is number of non-zeros returned on exit.
alpar@9 770 *
alpar@9 771 * Element indices stored in the array ind have the same sense as the
alpar@9 772 * index k, i.e. indices 1 to m denote auxiliary variables and indices
alpar@9 773 * m+1 to m+n denote structural ones (all these variables are obviously
alpar@9 774 * non-basic by definition).
alpar@9 775 *
alpar@9 776 * The computed row shows how the specified basic variable x[k] = xB[i]
alpar@9 777 * depends on non-basic variables:
alpar@9 778 *
alpar@9 779 * xB[i] = alfa[i,1]*xN[1] + alfa[i,2]*xN[2] + ... + alfa[i,n]*xN[n],
alpar@9 780 *
alpar@9 781 * where alfa[i,j] are elements of the simplex table row, xN[j] are
alpar@9 782 * non-basic (auxiliary and structural) variables.
alpar@9 783 *
alpar@9 784 * RETURNS
alpar@9 785 *
alpar@9 786 * The routine returns number of non-zero elements in the simplex table
alpar@9 787 * row stored in the arrays ind and val.
alpar@9 788 *
alpar@9 789 * BACKGROUND
alpar@9 790 *
alpar@9 791 * The system of equality constraints of the LP problem is:
alpar@9 792 *
alpar@9 793 * xR = A * xS, (1)
alpar@9 794 *
alpar@9 795 * where xR is the vector of auxliary variables, xS is the vector of
alpar@9 796 * structural variables, A is the matrix of constraint coefficients.
alpar@9 797 *
alpar@9 798 * The system (1) can be written in homogenous form as follows:
alpar@9 799 *
alpar@9 800 * A~ * x = 0, (2)
alpar@9 801 *
alpar@9 802 * where A~ = (I | -A) is the augmented constraint matrix (has m rows
alpar@9 803 * and m+n columns), x = (xR | xS) is the vector of all (auxiliary and
alpar@9 804 * structural) variables.
alpar@9 805 *
alpar@9 806 * By definition for the current basis we have:
alpar@9 807 *
alpar@9 808 * A~ = (B | N), (3)
alpar@9 809 *
alpar@9 810 * where B is the basis matrix. Thus, the system (2) can be written as:
alpar@9 811 *
alpar@9 812 * B * xB + N * xN = 0. (4)
alpar@9 813 *
alpar@9 814 * From (4) it follows that:
alpar@9 815 *
alpar@9 816 * xB = A^ * xN, (5)
alpar@9 817 *
alpar@9 818 * where the matrix
alpar@9 819 *
alpar@9 820 * A^ = - inv(B) * N (6)
alpar@9 821 *
alpar@9 822 * is called the simplex table.
alpar@9 823 *
alpar@9 824 * It is understood that i-th row of the simplex table is:
alpar@9 825 *
alpar@9 826 * e * A^ = - e * inv(B) * N, (7)
alpar@9 827 *
alpar@9 828 * where e is a unity vector with e[i] = 1.
alpar@9 829 *
alpar@9 830 * To compute i-th row of the simplex table the routine first computes
alpar@9 831 * i-th row of the inverse:
alpar@9 832 *
alpar@9 833 * rho = inv(B') * e, (8)
alpar@9 834 *
alpar@9 835 * where B' is a matrix transposed to B, and then computes elements of
alpar@9 836 * i-th row of the simplex table as scalar products:
alpar@9 837 *
alpar@9 838 * alfa[i,j] = - rho * N[j] for all j, (9)
alpar@9 839 *
alpar@9 840 * where N[j] is a column of the augmented constraint matrix A~, which
alpar@9 841 * corresponds to some non-basic auxiliary or structural variable. */
alpar@9 842
alpar@9 843 int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[])
alpar@9 844 { int m = lp->m;
alpar@9 845 int n = lp->n;
alpar@9 846 int i, t, len, lll, *iii;
alpar@9 847 double alfa, *rho, *vvv;
alpar@9 848 if (!(m == 0 || lp->valid))
alpar@9 849 xerror("glp_eval_tab_row: basis factorization does not exist\n"
alpar@9 850 );
alpar@9 851 if (!(1 <= k && k <= m+n))
alpar@9 852 xerror("glp_eval_tab_row: k = %d; variable number out of range"
alpar@9 853 , k);
alpar@9 854 /* determine xB[i] which corresponds to x[k] */
alpar@9 855 if (k <= m)
alpar@9 856 i = glp_get_row_bind(lp, k);
alpar@9 857 else
alpar@9 858 i = glp_get_col_bind(lp, k-m);
alpar@9 859 if (i == 0)
alpar@9 860 xerror("glp_eval_tab_row: k = %d; variable must be basic", k);
alpar@9 861 xassert(1 <= i && i <= m);
alpar@9 862 /* allocate working arrays */
alpar@9 863 rho = xcalloc(1+m, sizeof(double));
alpar@9 864 iii = xcalloc(1+m, sizeof(int));
alpar@9 865 vvv = xcalloc(1+m, sizeof(double));
alpar@9 866 /* compute i-th row of the inverse; see (8) */
alpar@9 867 for (t = 1; t <= m; t++) rho[t] = 0.0;
alpar@9 868 rho[i] = 1.0;
alpar@9 869 glp_btran(lp, rho);
alpar@9 870 /* compute i-th row of the simplex table */
alpar@9 871 len = 0;
alpar@9 872 for (k = 1; k <= m+n; k++)
alpar@9 873 { if (k <= m)
alpar@9 874 { /* x[k] is auxiliary variable, so N[k] is a unity column */
alpar@9 875 if (glp_get_row_stat(lp, k) == GLP_BS) continue;
alpar@9 876 /* compute alfa[i,j]; see (9) */
alpar@9 877 alfa = - rho[k];
alpar@9 878 }
alpar@9 879 else
alpar@9 880 { /* x[k] is structural variable, so N[k] is a column of the
alpar@9 881 original constraint matrix A with negative sign */
alpar@9 882 if (glp_get_col_stat(lp, k-m) == GLP_BS) continue;
alpar@9 883 /* compute alfa[i,j]; see (9) */
alpar@9 884 lll = glp_get_mat_col(lp, k-m, iii, vvv);
alpar@9 885 alfa = 0.0;
alpar@9 886 for (t = 1; t <= lll; t++) alfa += rho[iii[t]] * vvv[t];
alpar@9 887 }
alpar@9 888 /* store alfa[i,j] */
alpar@9 889 if (alfa != 0.0) len++, ind[len] = k, val[len] = alfa;
alpar@9 890 }
alpar@9 891 xassert(len <= n);
alpar@9 892 /* free working arrays */
alpar@9 893 xfree(rho);
alpar@9 894 xfree(iii);
alpar@9 895 xfree(vvv);
alpar@9 896 /* return to the calling program */
alpar@9 897 return len;
alpar@9 898 }
alpar@9 899
alpar@9 900 /***********************************************************************
alpar@9 901 * NAME
alpar@9 902 *
alpar@9 903 * glp_eval_tab_col - compute column of the simplex tableau
alpar@9 904 *
alpar@9 905 * SYNOPSIS
alpar@9 906 *
alpar@9 907 * int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[]);
alpar@9 908 *
alpar@9 909 * DESCRIPTION
alpar@9 910 *
alpar@9 911 * The routine glp_eval_tab_col computes a column of the current simplex
alpar@9 912 * table for the non-basic variable, which is specified by the number k:
alpar@9 913 * if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
alpar@9 914 * x[k] is (k-m)-th structural variable, where m is number of rows, and
alpar@9 915 * n is number of columns. The current basis must be available.
alpar@9 916 *
alpar@9 917 * The routine stores row indices and numerical values of non-zero
alpar@9 918 * elements of the computed column using sparse format to the locations
alpar@9 919 * ind[1], ..., ind[len] and val[1], ..., val[len] respectively, where
alpar@9 920 * 0 <= len <= m is number of non-zeros returned on exit.
alpar@9 921 *
alpar@9 922 * Element indices stored in the array ind have the same sense as the
alpar@9 923 * index k, i.e. indices 1 to m denote auxiliary variables and indices
alpar@9 924 * m+1 to m+n denote structural ones (all these variables are obviously
alpar@9 925 * basic by the definition).
alpar@9 926 *
alpar@9 927 * The computed column shows how basic variables depend on the specified
alpar@9 928 * non-basic variable x[k] = xN[j]:
alpar@9 929 *
alpar@9 930 * xB[1] = ... + alfa[1,j]*xN[j] + ...
alpar@9 931 * xB[2] = ... + alfa[2,j]*xN[j] + ...
alpar@9 932 * . . . . . .
alpar@9 933 * xB[m] = ... + alfa[m,j]*xN[j] + ...
alpar@9 934 *
alpar@9 935 * where alfa[i,j] are elements of the simplex table column, xB[i] are
alpar@9 936 * basic (auxiliary and structural) variables.
alpar@9 937 *
alpar@9 938 * RETURNS
alpar@9 939 *
alpar@9 940 * The routine returns number of non-zero elements in the simplex table
alpar@9 941 * column stored in the arrays ind and val.
alpar@9 942 *
alpar@9 943 * BACKGROUND
alpar@9 944 *
alpar@9 945 * As it was explained in comments to the routine glp_eval_tab_row (see
alpar@9 946 * above) the simplex table is the following matrix:
alpar@9 947 *
alpar@9 948 * A^ = - inv(B) * N. (1)
alpar@9 949 *
alpar@9 950 * Therefore j-th column of the simplex table is:
alpar@9 951 *
alpar@9 952 * A^ * e = - inv(B) * N * e = - inv(B) * N[j], (2)
alpar@9 953 *
alpar@9 954 * where e is a unity vector with e[j] = 1, B is the basis matrix, N[j]
alpar@9 955 * is a column of the augmented constraint matrix A~, which corresponds
alpar@9 956 * to the given non-basic auxiliary or structural variable. */
alpar@9 957
alpar@9 958 int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[])
alpar@9 959 { int m = lp->m;
alpar@9 960 int n = lp->n;
alpar@9 961 int t, len, stat;
alpar@9 962 double *col;
alpar@9 963 if (!(m == 0 || lp->valid))
alpar@9 964 xerror("glp_eval_tab_col: basis factorization does not exist\n"
alpar@9 965 );
alpar@9 966 if (!(1 <= k && k <= m+n))
alpar@9 967 xerror("glp_eval_tab_col: k = %d; variable number out of range"
alpar@9 968 , k);
alpar@9 969 if (k <= m)
alpar@9 970 stat = glp_get_row_stat(lp, k);
alpar@9 971 else
alpar@9 972 stat = glp_get_col_stat(lp, k-m);
alpar@9 973 if (stat == GLP_BS)
alpar@9 974 xerror("glp_eval_tab_col: k = %d; variable must be non-basic",
alpar@9 975 k);
alpar@9 976 /* obtain column N[k] with negative sign */
alpar@9 977 col = xcalloc(1+m, sizeof(double));
alpar@9 978 for (t = 1; t <= m; t++) col[t] = 0.0;
alpar@9 979 if (k <= m)
alpar@9 980 { /* x[k] is auxiliary variable, so N[k] is a unity column */
alpar@9 981 col[k] = -1.0;
alpar@9 982 }
alpar@9 983 else
alpar@9 984 { /* x[k] is structural variable, so N[k] is a column of the
alpar@9 985 original constraint matrix A with negative sign */
alpar@9 986 len = glp_get_mat_col(lp, k-m, ind, val);
alpar@9 987 for (t = 1; t <= len; t++) col[ind[t]] = val[t];
alpar@9 988 }
alpar@9 989 /* compute column of the simplex table, which corresponds to the
alpar@9 990 specified non-basic variable x[k] */
alpar@9 991 glp_ftran(lp, col);
alpar@9 992 len = 0;
alpar@9 993 for (t = 1; t <= m; t++)
alpar@9 994 { if (col[t] != 0.0)
alpar@9 995 { len++;
alpar@9 996 ind[len] = glp_get_bhead(lp, t);
alpar@9 997 val[len] = col[t];
alpar@9 998 }
alpar@9 999 }
alpar@9 1000 xfree(col);
alpar@9 1001 /* return to the calling program */
alpar@9 1002 return len;
alpar@9 1003 }
alpar@9 1004
alpar@9 1005 /***********************************************************************
alpar@9 1006 * NAME
alpar@9 1007 *
alpar@9 1008 * glp_transform_row - transform explicitly specified row
alpar@9 1009 *
alpar@9 1010 * SYNOPSIS
alpar@9 1011 *
alpar@9 1012 * int glp_transform_row(glp_prob *P, int len, int ind[], double val[]);
alpar@9 1013 *
alpar@9 1014 * DESCRIPTION
alpar@9 1015 *
alpar@9 1016 * The routine glp_transform_row performs the same operation as the
alpar@9 1017 * routine glp_eval_tab_row with exception that the row to be
alpar@9 1018 * transformed is specified explicitly as a sparse vector.
alpar@9 1019 *
alpar@9 1020 * The explicitly specified row may be thought as a linear form:
alpar@9 1021 *
alpar@9 1022 * x = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n], (1)
alpar@9 1023 *
alpar@9 1024 * where x is an auxiliary variable for this row, a[j] are coefficients
alpar@9 1025 * of the linear form, x[m+j] are structural variables.
alpar@9 1026 *
alpar@9 1027 * On entry column indices and numerical values of non-zero elements of
alpar@9 1028 * the row should be stored in locations ind[1], ..., ind[len] and
alpar@9 1029 * val[1], ..., val[len], where len is the number of non-zero elements.
alpar@9 1030 *
alpar@9 1031 * This routine uses the system of equality constraints and the current
alpar@9 1032 * basis in order to express the auxiliary variable x in (1) through the
alpar@9 1033 * current non-basic variables (as if the transformed row were added to
alpar@9 1034 * the problem object and its auxiliary variable were basic), i.e. the
alpar@9 1035 * resultant row has the form:
alpar@9 1036 *
alpar@9 1037 * x = alfa[1]*xN[1] + alfa[2]*xN[2] + ... + alfa[n]*xN[n], (2)
alpar@9 1038 *
alpar@9 1039 * where xN[j] are non-basic (auxiliary or structural) variables, n is
alpar@9 1040 * the number of columns in the LP problem object.
alpar@9 1041 *
alpar@9 1042 * On exit the routine stores indices and numerical values of non-zero
alpar@9 1043 * elements of the resultant row (2) in locations ind[1], ..., ind[len']
alpar@9 1044 * and val[1], ..., val[len'], where 0 <= len' <= n is the number of
alpar@9 1045 * non-zero elements in the resultant row returned by the routine. Note
alpar@9 1046 * that indices (numbers) of non-basic variables stored in the array ind
alpar@9 1047 * correspond to original ordinal numbers of variables: indices 1 to m
alpar@9 1048 * mean auxiliary variables and indices m+1 to m+n mean structural ones.
alpar@9 1049 *
alpar@9 1050 * RETURNS
alpar@9 1051 *
alpar@9 1052 * The routine returns len', which is the number of non-zero elements in
alpar@9 1053 * the resultant row stored in the arrays ind and val.
alpar@9 1054 *
alpar@9 1055 * BACKGROUND
alpar@9 1056 *
alpar@9 1057 * The explicitly specified row (1) is transformed in the same way as it
alpar@9 1058 * were the objective function row.
alpar@9 1059 *
alpar@9 1060 * From (1) it follows that:
alpar@9 1061 *
alpar@9 1062 * x = aB * xB + aN * xN, (3)
alpar@9 1063 *
alpar@9 1064 * where xB is the vector of basic variables, xN is the vector of
alpar@9 1065 * non-basic variables.
alpar@9 1066 *
alpar@9 1067 * The simplex table, which corresponds to the current basis, is:
alpar@9 1068 *
alpar@9 1069 * xB = [-inv(B) * N] * xN. (4)
alpar@9 1070 *
alpar@9 1071 * Therefore substituting xB from (4) to (3) we have:
alpar@9 1072 *
alpar@9 1073 * x = aB * [-inv(B) * N] * xN + aN * xN =
alpar@9 1074 * (5)
alpar@9 1075 * = rho * (-N) * xN + aN * xN = alfa * xN,
alpar@9 1076 *
alpar@9 1077 * where:
alpar@9 1078 *
alpar@9 1079 * rho = inv(B') * aB, (6)
alpar@9 1080 *
alpar@9 1081 * and
alpar@9 1082 *
alpar@9 1083 * alfa = aN + rho * (-N) (7)
alpar@9 1084 *
alpar@9 1085 * is the resultant row computed by the routine. */
alpar@9 1086
alpar@9 1087 int glp_transform_row(glp_prob *P, int len, int ind[], double val[])
alpar@9 1088 { int i, j, k, m, n, t, lll, *iii;
alpar@9 1089 double alfa, *a, *aB, *rho, *vvv;
alpar@9 1090 if (!glp_bf_exists(P))
alpar@9 1091 xerror("glp_transform_row: basis factorization does not exist "
alpar@9 1092 "\n");
alpar@9 1093 m = glp_get_num_rows(P);
alpar@9 1094 n = glp_get_num_cols(P);
alpar@9 1095 /* unpack the row to be transformed to the array a */
alpar@9 1096 a = xcalloc(1+n, sizeof(double));
alpar@9 1097 for (j = 1; j <= n; j++) a[j] = 0.0;
alpar@9 1098 if (!(0 <= len && len <= n))
alpar@9 1099 xerror("glp_transform_row: len = %d; invalid row length\n",
alpar@9 1100 len);
alpar@9 1101 for (t = 1; t <= len; t++)
alpar@9 1102 { j = ind[t];
alpar@9 1103 if (!(1 <= j && j <= n))
alpar@9 1104 xerror("glp_transform_row: ind[%d] = %d; column index out o"
alpar@9 1105 "f range\n", t, j);
alpar@9 1106 if (val[t] == 0.0)
alpar@9 1107 xerror("glp_transform_row: val[%d] = 0; zero coefficient no"
alpar@9 1108 "t allowed\n", t);
alpar@9 1109 if (a[j] != 0.0)
alpar@9 1110 xerror("glp_transform_row: ind[%d] = %d; duplicate column i"
alpar@9 1111 "ndices not allowed\n", t, j);
alpar@9 1112 a[j] = val[t];
alpar@9 1113 }
alpar@9 1114 /* construct the vector aB */
alpar@9 1115 aB = xcalloc(1+m, sizeof(double));
alpar@9 1116 for (i = 1; i <= m; i++)
alpar@9 1117 { k = glp_get_bhead(P, i);
alpar@9 1118 /* xB[i] is k-th original variable */
alpar@9 1119 xassert(1 <= k && k <= m+n);
alpar@9 1120 aB[i] = (k <= m ? 0.0 : a[k-m]);
alpar@9 1121 }
alpar@9 1122 /* solve the system B'*rho = aB to compute the vector rho */
alpar@9 1123 rho = aB, glp_btran(P, rho);
alpar@9 1124 /* compute coefficients at non-basic auxiliary variables */
alpar@9 1125 len = 0;
alpar@9 1126 for (i = 1; i <= m; i++)
alpar@9 1127 { if (glp_get_row_stat(P, i) != GLP_BS)
alpar@9 1128 { alfa = - rho[i];
alpar@9 1129 if (alfa != 0.0)
alpar@9 1130 { len++;
alpar@9 1131 ind[len] = i;
alpar@9 1132 val[len] = alfa;
alpar@9 1133 }
alpar@9 1134 }
alpar@9 1135 }
alpar@9 1136 /* compute coefficients at non-basic structural variables */
alpar@9 1137 iii = xcalloc(1+m, sizeof(int));
alpar@9 1138 vvv = xcalloc(1+m, sizeof(double));
alpar@9 1139 for (j = 1; j <= n; j++)
alpar@9 1140 { if (glp_get_col_stat(P, j) != GLP_BS)
alpar@9 1141 { alfa = a[j];
alpar@9 1142 lll = glp_get_mat_col(P, j, iii, vvv);
alpar@9 1143 for (t = 1; t <= lll; t++) alfa += vvv[t] * rho[iii[t]];
alpar@9 1144 if (alfa != 0.0)
alpar@9 1145 { len++;
alpar@9 1146 ind[len] = m+j;
alpar@9 1147 val[len] = alfa;
alpar@9 1148 }
alpar@9 1149 }
alpar@9 1150 }
alpar@9 1151 xassert(len <= n);
alpar@9 1152 xfree(iii);
alpar@9 1153 xfree(vvv);
alpar@9 1154 xfree(aB);
alpar@9 1155 xfree(a);
alpar@9 1156 return len;
alpar@9 1157 }
alpar@9 1158
alpar@9 1159 /***********************************************************************
alpar@9 1160 * NAME
alpar@9 1161 *
alpar@9 1162 * glp_transform_col - transform explicitly specified column
alpar@9 1163 *
alpar@9 1164 * SYNOPSIS
alpar@9 1165 *
alpar@9 1166 * int glp_transform_col(glp_prob *P, int len, int ind[], double val[]);
alpar@9 1167 *
alpar@9 1168 * DESCRIPTION
alpar@9 1169 *
alpar@9 1170 * The routine glp_transform_col performs the same operation as the
alpar@9 1171 * routine glp_eval_tab_col with exception that the column to be
alpar@9 1172 * transformed is specified explicitly as a sparse vector.
alpar@9 1173 *
alpar@9 1174 * The explicitly specified column may be thought as if it were added
alpar@9 1175 * to the original system of equality constraints:
alpar@9 1176 *
alpar@9 1177 * x[1] = a[1,1]*x[m+1] + ... + a[1,n]*x[m+n] + a[1]*x
alpar@9 1178 * x[2] = a[2,1]*x[m+1] + ... + a[2,n]*x[m+n] + a[2]*x (1)
alpar@9 1179 * . . . . . . . . . . . . . . .
alpar@9 1180 * x[m] = a[m,1]*x[m+1] + ... + a[m,n]*x[m+n] + a[m]*x
alpar@9 1181 *
alpar@9 1182 * where x[i] are auxiliary variables, x[m+j] are structural variables,
alpar@9 1183 * x is a structural variable for the explicitly specified column, a[i]
alpar@9 1184 * are constraint coefficients for x.
alpar@9 1185 *
alpar@9 1186 * On entry row indices and numerical values of non-zero elements of
alpar@9 1187 * the column should be stored in locations ind[1], ..., ind[len] and
alpar@9 1188 * val[1], ..., val[len], where len is the number of non-zero elements.
alpar@9 1189 *
alpar@9 1190 * This routine uses the system of equality constraints and the current
alpar@9 1191 * basis in order to express the current basic variables through the
alpar@9 1192 * structural variable x in (1) (as if the transformed column were added
alpar@9 1193 * to the problem object and the variable x were non-basic), i.e. the
alpar@9 1194 * resultant column has the form:
alpar@9 1195 *
alpar@9 1196 * xB[1] = ... + alfa[1]*x
alpar@9 1197 * xB[2] = ... + alfa[2]*x (2)
alpar@9 1198 * . . . . . .
alpar@9 1199 * xB[m] = ... + alfa[m]*x
alpar@9 1200 *
alpar@9 1201 * where xB are basic (auxiliary and structural) variables, m is the
alpar@9 1202 * number of rows in the problem object.
alpar@9 1203 *
alpar@9 1204 * On exit the routine stores indices and numerical values of non-zero
alpar@9 1205 * elements of the resultant column (2) in locations ind[1], ...,
alpar@9 1206 * ind[len'] and val[1], ..., val[len'], where 0 <= len' <= m is the
alpar@9 1207 * number of non-zero element in the resultant column returned by the
alpar@9 1208 * routine. Note that indices (numbers) of basic variables stored in
alpar@9 1209 * the array ind correspond to original ordinal numbers of variables:
alpar@9 1210 * indices 1 to m mean auxiliary variables and indices m+1 to m+n mean
alpar@9 1211 * structural ones.
alpar@9 1212 *
alpar@9 1213 * RETURNS
alpar@9 1214 *
alpar@9 1215 * The routine returns len', which is the number of non-zero elements
alpar@9 1216 * in the resultant column stored in the arrays ind and val.
alpar@9 1217 *
alpar@9 1218 * BACKGROUND
alpar@9 1219 *
alpar@9 1220 * The explicitly specified column (1) is transformed in the same way
alpar@9 1221 * as any other column of the constraint matrix using the formula:
alpar@9 1222 *
alpar@9 1223 * alfa = inv(B) * a, (3)
alpar@9 1224 *
alpar@9 1225 * where alfa is the resultant column computed by the routine. */
alpar@9 1226
alpar@9 1227 int glp_transform_col(glp_prob *P, int len, int ind[], double val[])
alpar@9 1228 { int i, m, t;
alpar@9 1229 double *a, *alfa;
alpar@9 1230 if (!glp_bf_exists(P))
alpar@9 1231 xerror("glp_transform_col: basis factorization does not exist "
alpar@9 1232 "\n");
alpar@9 1233 m = glp_get_num_rows(P);
alpar@9 1234 /* unpack the column to be transformed to the array a */
alpar@9 1235 a = xcalloc(1+m, sizeof(double));
alpar@9 1236 for (i = 1; i <= m; i++) a[i] = 0.0;
alpar@9 1237 if (!(0 <= len && len <= m))
alpar@9 1238 xerror("glp_transform_col: len = %d; invalid column length\n",
alpar@9 1239 len);
alpar@9 1240 for (t = 1; t <= len; t++)
alpar@9 1241 { i = ind[t];
alpar@9 1242 if (!(1 <= i && i <= m))
alpar@9 1243 xerror("glp_transform_col: ind[%d] = %d; row index out of r"
alpar@9 1244 "ange\n", t, i);
alpar@9 1245 if (val[t] == 0.0)
alpar@9 1246 xerror("glp_transform_col: val[%d] = 0; zero coefficient no"
alpar@9 1247 "t allowed\n", t);
alpar@9 1248 if (a[i] != 0.0)
alpar@9 1249 xerror("glp_transform_col: ind[%d] = %d; duplicate row indi"
alpar@9 1250 "ces not allowed\n", t, i);
alpar@9 1251 a[i] = val[t];
alpar@9 1252 }
alpar@9 1253 /* solve the system B*a = alfa to compute the vector alfa */
alpar@9 1254 alfa = a, glp_ftran(P, alfa);
alpar@9 1255 /* store resultant coefficients */
alpar@9 1256 len = 0;
alpar@9 1257 for (i = 1; i <= m; i++)
alpar@9 1258 { if (alfa[i] != 0.0)
alpar@9 1259 { len++;
alpar@9 1260 ind[len] = glp_get_bhead(P, i);
alpar@9 1261 val[len] = alfa[i];
alpar@9 1262 }
alpar@9 1263 }
alpar@9 1264 xfree(a);
alpar@9 1265 return len;
alpar@9 1266 }
alpar@9 1267
alpar@9 1268 /***********************************************************************
alpar@9 1269 * NAME
alpar@9 1270 *
alpar@9 1271 * glp_prim_rtest - perform primal ratio test
alpar@9 1272 *
alpar@9 1273 * SYNOPSIS
alpar@9 1274 *
alpar@9 1275 * int glp_prim_rtest(glp_prob *P, int len, const int ind[],
alpar@9 1276 * const double val[], int dir, double eps);
alpar@9 1277 *
alpar@9 1278 * DESCRIPTION
alpar@9 1279 *
alpar@9 1280 * The routine glp_prim_rtest performs the primal ratio test using an
alpar@9 1281 * explicitly specified column of the simplex table.
alpar@9 1282 *
alpar@9 1283 * The current basic solution associated with the LP problem object
alpar@9 1284 * must be primal feasible.
alpar@9 1285 *
alpar@9 1286 * The explicitly specified column of the simplex table shows how the
alpar@9 1287 * basic variables xB depend on some non-basic variable x (which is not
alpar@9 1288 * necessarily presented in the problem object):
alpar@9 1289 *
alpar@9 1290 * xB[1] = ... + alfa[1] * x + ...
alpar@9 1291 * xB[2] = ... + alfa[2] * x + ... (*)
alpar@9 1292 * . . . . . . . .
alpar@9 1293 * xB[m] = ... + alfa[m] * x + ...
alpar@9 1294 *
alpar@9 1295 * The column (*) is specifed on entry to the routine using the sparse
alpar@9 1296 * format. Ordinal numbers of basic variables xB[i] should be placed in
alpar@9 1297 * locations ind[1], ..., ind[len], where ordinal number 1 to m denote
alpar@9 1298 * auxiliary variables, and ordinal numbers m+1 to m+n denote structural
alpar@9 1299 * variables. The corresponding non-zero coefficients alfa[i] should be
alpar@9 1300 * placed in locations val[1], ..., val[len]. The arrays ind and val are
alpar@9 1301 * not changed on exit.
alpar@9 1302 *
alpar@9 1303 * The parameter dir specifies direction in which the variable x changes
alpar@9 1304 * on entering the basis: +1 means increasing, -1 means decreasing.
alpar@9 1305 *
alpar@9 1306 * The parameter eps is an absolute tolerance (small positive number)
alpar@9 1307 * used by the routine to skip small alfa[j] of the row (*).
alpar@9 1308 *
alpar@9 1309 * The routine determines which basic variable (among specified in
alpar@9 1310 * ind[1], ..., ind[len]) should leave the basis in order to keep primal
alpar@9 1311 * feasibility.
alpar@9 1312 *
alpar@9 1313 * RETURNS
alpar@9 1314 *
alpar@9 1315 * The routine glp_prim_rtest returns the index piv in the arrays ind
alpar@9 1316 * and val corresponding to the pivot element chosen, 1 <= piv <= len.
alpar@9 1317 * If the adjacent basic solution is primal unbounded and therefore the
alpar@9 1318 * choice cannot be made, the routine returns zero.
alpar@9 1319 *
alpar@9 1320 * COMMENTS
alpar@9 1321 *
alpar@9 1322 * If the non-basic variable x is presented in the LP problem object,
alpar@9 1323 * the column (*) can be computed with the routine glp_eval_tab_col;
alpar@9 1324 * otherwise it can be computed with the routine glp_transform_col. */
alpar@9 1325
alpar@9 1326 int glp_prim_rtest(glp_prob *P, int len, const int ind[],
alpar@9 1327 const double val[], int dir, double eps)
alpar@9 1328 { int k, m, n, piv, t, type, stat;
alpar@9 1329 double alfa, big, beta, lb, ub, temp, teta;
alpar@9 1330 if (glp_get_prim_stat(P) != GLP_FEAS)
alpar@9 1331 xerror("glp_prim_rtest: basic solution is not primal feasible "
alpar@9 1332 "\n");
alpar@9 1333 if (!(dir == +1 || dir == -1))
alpar@9 1334 xerror("glp_prim_rtest: dir = %d; invalid parameter\n", dir);
alpar@9 1335 if (!(0.0 < eps && eps < 1.0))
alpar@9 1336 xerror("glp_prim_rtest: eps = %g; invalid parameter\n", eps);
alpar@9 1337 m = glp_get_num_rows(P);
alpar@9 1338 n = glp_get_num_cols(P);
alpar@9 1339 /* initial settings */
alpar@9 1340 piv = 0, teta = DBL_MAX, big = 0.0;
alpar@9 1341 /* walk through the entries of the specified column */
alpar@9 1342 for (t = 1; t <= len; t++)
alpar@9 1343 { /* get the ordinal number of basic variable */
alpar@9 1344 k = ind[t];
alpar@9 1345 if (!(1 <= k && k <= m+n))
alpar@9 1346 xerror("glp_prim_rtest: ind[%d] = %d; variable number out o"
alpar@9 1347 "f range\n", t, k);
alpar@9 1348 /* determine type, bounds, status and primal value of basic
alpar@9 1349 variable xB[i] = x[k] in the current basic solution */
alpar@9 1350 if (k <= m)
alpar@9 1351 { type = glp_get_row_type(P, k);
alpar@9 1352 lb = glp_get_row_lb(P, k);
alpar@9 1353 ub = glp_get_row_ub(P, k);
alpar@9 1354 stat = glp_get_row_stat(P, k);
alpar@9 1355 beta = glp_get_row_prim(P, k);
alpar@9 1356 }
alpar@9 1357 else
alpar@9 1358 { type = glp_get_col_type(P, k-m);
alpar@9 1359 lb = glp_get_col_lb(P, k-m);
alpar@9 1360 ub = glp_get_col_ub(P, k-m);
alpar@9 1361 stat = glp_get_col_stat(P, k-m);
alpar@9 1362 beta = glp_get_col_prim(P, k-m);
alpar@9 1363 }
alpar@9 1364 if (stat != GLP_BS)
alpar@9 1365 xerror("glp_prim_rtest: ind[%d] = %d; non-basic variable no"
alpar@9 1366 "t allowed\n", t, k);
alpar@9 1367 /* determine influence coefficient at basic variable xB[i]
alpar@9 1368 in the explicitly specified column and turn to the case of
alpar@9 1369 increasing the variable x in order to simplify the program
alpar@9 1370 logic */
alpar@9 1371 alfa = (dir > 0 ? + val[t] : - val[t]);
alpar@9 1372 /* analyze main cases */
alpar@9 1373 if (type == GLP_FR)
alpar@9 1374 { /* xB[i] is free variable */
alpar@9 1375 continue;
alpar@9 1376 }
alpar@9 1377 else if (type == GLP_LO)
alpar@9 1378 lo: { /* xB[i] has an lower bound */
alpar@9 1379 if (alfa > - eps) continue;
alpar@9 1380 temp = (lb - beta) / alfa;
alpar@9 1381 }
alpar@9 1382 else if (type == GLP_UP)
alpar@9 1383 up: { /* xB[i] has an upper bound */
alpar@9 1384 if (alfa < + eps) continue;
alpar@9 1385 temp = (ub - beta) / alfa;
alpar@9 1386 }
alpar@9 1387 else if (type == GLP_DB)
alpar@9 1388 { /* xB[i] has both lower and upper bounds */
alpar@9 1389 if (alfa < 0.0) goto lo; else goto up;
alpar@9 1390 }
alpar@9 1391 else if (type == GLP_FX)
alpar@9 1392 { /* xB[i] is fixed variable */
alpar@9 1393 if (- eps < alfa && alfa < + eps) continue;
alpar@9 1394 temp = 0.0;
alpar@9 1395 }
alpar@9 1396 else
alpar@9 1397 xassert(type != type);
alpar@9 1398 /* if the value of the variable xB[i] violates its lower or
alpar@9 1399 upper bound (slightly, because the current basis is assumed
alpar@9 1400 to be primal feasible), temp is negative; we can think this
alpar@9 1401 happens due to round-off errors and the value is exactly on
alpar@9 1402 the bound; this allows replacing temp by zero */
alpar@9 1403 if (temp < 0.0) temp = 0.0;
alpar@9 1404 /* apply the minimal ratio test */
alpar@9 1405 if (teta > temp || teta == temp && big < fabs(alfa))
alpar@9 1406 piv = t, teta = temp, big = fabs(alfa);
alpar@9 1407 }
alpar@9 1408 /* return index of the pivot element chosen */
alpar@9 1409 return piv;
alpar@9 1410 }
alpar@9 1411
alpar@9 1412 /***********************************************************************
alpar@9 1413 * NAME
alpar@9 1414 *
alpar@9 1415 * glp_dual_rtest - perform dual ratio test
alpar@9 1416 *
alpar@9 1417 * SYNOPSIS
alpar@9 1418 *
alpar@9 1419 * int glp_dual_rtest(glp_prob *P, int len, const int ind[],
alpar@9 1420 * const double val[], int dir, double eps);
alpar@9 1421 *
alpar@9 1422 * DESCRIPTION
alpar@9 1423 *
alpar@9 1424 * The routine glp_dual_rtest performs the dual ratio test using an
alpar@9 1425 * explicitly specified row of the simplex table.
alpar@9 1426 *
alpar@9 1427 * The current basic solution associated with the LP problem object
alpar@9 1428 * must be dual feasible.
alpar@9 1429 *
alpar@9 1430 * The explicitly specified row of the simplex table is a linear form
alpar@9 1431 * that shows how some basic variable x (which is not necessarily
alpar@9 1432 * presented in the problem object) depends on non-basic variables xN:
alpar@9 1433 *
alpar@9 1434 * x = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n]. (*)
alpar@9 1435 *
alpar@9 1436 * The row (*) is specified on entry to the routine using the sparse
alpar@9 1437 * format. Ordinal numbers of non-basic variables xN[j] should be placed
alpar@9 1438 * in locations ind[1], ..., ind[len], where ordinal numbers 1 to m
alpar@9 1439 * denote auxiliary variables, and ordinal numbers m+1 to m+n denote
alpar@9 1440 * structural variables. The corresponding non-zero coefficients alfa[j]
alpar@9 1441 * should be placed in locations val[1], ..., val[len]. The arrays ind
alpar@9 1442 * and val are not changed on exit.
alpar@9 1443 *
alpar@9 1444 * The parameter dir specifies direction in which the variable x changes
alpar@9 1445 * on leaving the basis: +1 means that x goes to its lower bound, and -1
alpar@9 1446 * means that x goes to its upper bound.
alpar@9 1447 *
alpar@9 1448 * The parameter eps is an absolute tolerance (small positive number)
alpar@9 1449 * used by the routine to skip small alfa[j] of the row (*).
alpar@9 1450 *
alpar@9 1451 * The routine determines which non-basic variable (among specified in
alpar@9 1452 * ind[1], ..., ind[len]) should enter the basis in order to keep dual
alpar@9 1453 * feasibility.
alpar@9 1454 *
alpar@9 1455 * RETURNS
alpar@9 1456 *
alpar@9 1457 * The routine glp_dual_rtest returns the index piv in the arrays ind
alpar@9 1458 * and val corresponding to the pivot element chosen, 1 <= piv <= len.
alpar@9 1459 * If the adjacent basic solution is dual unbounded and therefore the
alpar@9 1460 * choice cannot be made, the routine returns zero.
alpar@9 1461 *
alpar@9 1462 * COMMENTS
alpar@9 1463 *
alpar@9 1464 * If the basic variable x is presented in the LP problem object, the
alpar@9 1465 * row (*) can be computed with the routine glp_eval_tab_row; otherwise
alpar@9 1466 * it can be computed with the routine glp_transform_row. */
alpar@9 1467
alpar@9 1468 int glp_dual_rtest(glp_prob *P, int len, const int ind[],
alpar@9 1469 const double val[], int dir, double eps)
alpar@9 1470 { int k, m, n, piv, t, stat;
alpar@9 1471 double alfa, big, cost, obj, temp, teta;
alpar@9 1472 if (glp_get_dual_stat(P) != GLP_FEAS)
alpar@9 1473 xerror("glp_dual_rtest: basic solution is not dual feasible\n")
alpar@9 1474 ;
alpar@9 1475 if (!(dir == +1 || dir == -1))
alpar@9 1476 xerror("glp_dual_rtest: dir = %d; invalid parameter\n", dir);
alpar@9 1477 if (!(0.0 < eps && eps < 1.0))
alpar@9 1478 xerror("glp_dual_rtest: eps = %g; invalid parameter\n", eps);
alpar@9 1479 m = glp_get_num_rows(P);
alpar@9 1480 n = glp_get_num_cols(P);
alpar@9 1481 /* take into account optimization direction */
alpar@9 1482 obj = (glp_get_obj_dir(P) == GLP_MIN ? +1.0 : -1.0);
alpar@9 1483 /* initial settings */
alpar@9 1484 piv = 0, teta = DBL_MAX, big = 0.0;
alpar@9 1485 /* walk through the entries of the specified row */
alpar@9 1486 for (t = 1; t <= len; t++)
alpar@9 1487 { /* get ordinal number of non-basic variable */
alpar@9 1488 k = ind[t];
alpar@9 1489 if (!(1 <= k && k <= m+n))
alpar@9 1490 xerror("glp_dual_rtest: ind[%d] = %d; variable number out o"
alpar@9 1491 "f range\n", t, k);
alpar@9 1492 /* determine status and reduced cost of non-basic variable
alpar@9 1493 x[k] = xN[j] in the current basic solution */
alpar@9 1494 if (k <= m)
alpar@9 1495 { stat = glp_get_row_stat(P, k);
alpar@9 1496 cost = glp_get_row_dual(P, k);
alpar@9 1497 }
alpar@9 1498 else
alpar@9 1499 { stat = glp_get_col_stat(P, k-m);
alpar@9 1500 cost = glp_get_col_dual(P, k-m);
alpar@9 1501 }
alpar@9 1502 if (stat == GLP_BS)
alpar@9 1503 xerror("glp_dual_rtest: ind[%d] = %d; basic variable not al"
alpar@9 1504 "lowed\n", t, k);
alpar@9 1505 /* determine influence coefficient at non-basic variable xN[j]
alpar@9 1506 in the explicitly specified row and turn to the case of
alpar@9 1507 increasing the variable x in order to simplify the program
alpar@9 1508 logic */
alpar@9 1509 alfa = (dir > 0 ? + val[t] : - val[t]);
alpar@9 1510 /* analyze main cases */
alpar@9 1511 if (stat == GLP_NL)
alpar@9 1512 { /* xN[j] is on its lower bound */
alpar@9 1513 if (alfa < + eps) continue;
alpar@9 1514 temp = (obj * cost) / alfa;
alpar@9 1515 }
alpar@9 1516 else if (stat == GLP_NU)
alpar@9 1517 { /* xN[j] is on its upper bound */
alpar@9 1518 if (alfa > - eps) continue;
alpar@9 1519 temp = (obj * cost) / alfa;
alpar@9 1520 }
alpar@9 1521 else if (stat == GLP_NF)
alpar@9 1522 { /* xN[j] is non-basic free variable */
alpar@9 1523 if (- eps < alfa && alfa < + eps) continue;
alpar@9 1524 temp = 0.0;
alpar@9 1525 }
alpar@9 1526 else if (stat == GLP_NS)
alpar@9 1527 { /* xN[j] is non-basic fixed variable */
alpar@9 1528 continue;
alpar@9 1529 }
alpar@9 1530 else
alpar@9 1531 xassert(stat != stat);
alpar@9 1532 /* if the reduced cost of the variable xN[j] violates its zero
alpar@9 1533 bound (slightly, because the current basis is assumed to be
alpar@9 1534 dual feasible), temp is negative; we can think this happens
alpar@9 1535 due to round-off errors and the reduced cost is exact zero;
alpar@9 1536 this allows replacing temp by zero */
alpar@9 1537 if (temp < 0.0) temp = 0.0;
alpar@9 1538 /* apply the minimal ratio test */
alpar@9 1539 if (teta > temp || teta == temp && big < fabs(alfa))
alpar@9 1540 piv = t, teta = temp, big = fabs(alfa);
alpar@9 1541 }
alpar@9 1542 /* return index of the pivot element chosen */
alpar@9 1543 return piv;
alpar@9 1544 }
alpar@9 1545
alpar@9 1546 /***********************************************************************
alpar@9 1547 * NAME
alpar@9 1548 *
alpar@9 1549 * glp_analyze_row - simulate one iteration of dual simplex method
alpar@9 1550 *
alpar@9 1551 * SYNOPSIS
alpar@9 1552 *
alpar@9 1553 * int glp_analyze_row(glp_prob *P, int len, const int ind[],
alpar@9 1554 * const double val[], int type, double rhs, double eps, int *piv,
alpar@9 1555 * double *x, double *dx, double *y, double *dy, double *dz);
alpar@9 1556 *
alpar@9 1557 * DESCRIPTION
alpar@9 1558 *
alpar@9 1559 * Let the current basis be optimal or dual feasible, and there be
alpar@9 1560 * specified a row (constraint), which is violated by the current basic
alpar@9 1561 * solution. The routine glp_analyze_row simulates one iteration of the
alpar@9 1562 * dual simplex method to determine some information on the adjacent
alpar@9 1563 * basis (see below), where the specified row becomes active constraint
alpar@9 1564 * (i.e. its auxiliary variable becomes non-basic).
alpar@9 1565 *
alpar@9 1566 * The current basic solution associated with the problem object passed
alpar@9 1567 * to the routine must be dual feasible, and its primal components must
alpar@9 1568 * be defined.
alpar@9 1569 *
alpar@9 1570 * The row to be analyzed must be previously transformed either with
alpar@9 1571 * the routine glp_eval_tab_row (if the row is in the problem object)
alpar@9 1572 * or with the routine glp_transform_row (if the row is external, i.e.
alpar@9 1573 * not in the problem object). This is needed to express the row only
alpar@9 1574 * through (auxiliary and structural) variables, which are non-basic in
alpar@9 1575 * the current basis:
alpar@9 1576 *
alpar@9 1577 * y = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n],
alpar@9 1578 *
alpar@9 1579 * where y is an auxiliary variable of the row, alfa[j] is an influence
alpar@9 1580 * coefficient, xN[j] is a non-basic variable.
alpar@9 1581 *
alpar@9 1582 * The row is passed to the routine in sparse format. Ordinal numbers
alpar@9 1583 * of non-basic variables are stored in locations ind[1], ..., ind[len],
alpar@9 1584 * where numbers 1 to m denote auxiliary variables while numbers m+1 to
alpar@9 1585 * m+n denote structural variables. Corresponding non-zero coefficients
alpar@9 1586 * alfa[j] are stored in locations val[1], ..., val[len]. The arrays
alpar@9 1587 * ind and val are ot changed on exit.
alpar@9 1588 *
alpar@9 1589 * The parameters type and rhs specify the row type and its right-hand
alpar@9 1590 * side as follows:
alpar@9 1591 *
alpar@9 1592 * type = GLP_LO: y = sum alfa[j] * xN[j] >= rhs
alpar@9 1593 *
alpar@9 1594 * type = GLP_UP: y = sum alfa[j] * xN[j] <= rhs
alpar@9 1595 *
alpar@9 1596 * The parameter eps is an absolute tolerance (small positive number)
alpar@9 1597 * used by the routine to skip small coefficients alfa[j] on performing
alpar@9 1598 * the dual ratio test.
alpar@9 1599 *
alpar@9 1600 * If the operation was successful, the routine stores the following
alpar@9 1601 * information to corresponding location (if some parameter is NULL,
alpar@9 1602 * its value is not stored):
alpar@9 1603 *
alpar@9 1604 * piv index in the array ind and val, 1 <= piv <= len, determining
alpar@9 1605 * the non-basic variable, which would enter the adjacent basis;
alpar@9 1606 *
alpar@9 1607 * x value of the non-basic variable in the current basis;
alpar@9 1608 *
alpar@9 1609 * dx difference between values of the non-basic variable in the
alpar@9 1610 * adjacent and current bases, dx = x.new - x.old;
alpar@9 1611 *
alpar@9 1612 * y value of the row (i.e. of its auxiliary variable) in the
alpar@9 1613 * current basis;
alpar@9 1614 *
alpar@9 1615 * dy difference between values of the row in the adjacent and
alpar@9 1616 * current bases, dy = y.new - y.old;
alpar@9 1617 *
alpar@9 1618 * dz difference between values of the objective function in the
alpar@9 1619 * adjacent and current bases, dz = z.new - z.old. Note that in
alpar@9 1620 * case of minimization dz >= 0, and in case of maximization
alpar@9 1621 * dz <= 0, i.e. in the adjacent basis the objective function
alpar@9 1622 * always gets worse (degrades). */
alpar@9 1623
alpar@9 1624 int _glp_analyze_row(glp_prob *P, int len, const int ind[],
alpar@9 1625 const double val[], int type, double rhs, double eps, int *_piv,
alpar@9 1626 double *_x, double *_dx, double *_y, double *_dy, double *_dz)
alpar@9 1627 { int t, k, dir, piv, ret = 0;
alpar@9 1628 double x, dx, y, dy, dz;
alpar@9 1629 if (P->pbs_stat == GLP_UNDEF)
alpar@9 1630 xerror("glp_analyze_row: primal basic solution components are "
alpar@9 1631 "undefined\n");
alpar@9 1632 if (P->dbs_stat != GLP_FEAS)
alpar@9 1633 xerror("glp_analyze_row: basic solution is not dual feasible\n"
alpar@9 1634 );
alpar@9 1635 /* compute the row value y = sum alfa[j] * xN[j] in the current
alpar@9 1636 basis */
alpar@9 1637 if (!(0 <= len && len <= P->n))
alpar@9 1638 xerror("glp_analyze_row: len = %d; invalid row length\n", len);
alpar@9 1639 y = 0.0;
alpar@9 1640 for (t = 1; t <= len; t++)
alpar@9 1641 { /* determine value of x[k] = xN[j] in the current basis */
alpar@9 1642 k = ind[t];
alpar@9 1643 if (!(1 <= k && k <= P->m+P->n))
alpar@9 1644 xerror("glp_analyze_row: ind[%d] = %d; row/column index out"
alpar@9 1645 " of range\n", t, k);
alpar@9 1646 if (k <= P->m)
alpar@9 1647 { /* x[k] is auxiliary variable */
alpar@9 1648 if (P->row[k]->stat == GLP_BS)
alpar@9 1649 xerror("glp_analyze_row: ind[%d] = %d; basic auxiliary v"
alpar@9 1650 "ariable is not allowed\n", t, k);
alpar@9 1651 x = P->row[k]->prim;
alpar@9 1652 }
alpar@9 1653 else
alpar@9 1654 { /* x[k] is structural variable */
alpar@9 1655 if (P->col[k-P->m]->stat == GLP_BS)
alpar@9 1656 xerror("glp_analyze_row: ind[%d] = %d; basic structural "
alpar@9 1657 "variable is not allowed\n", t, k);
alpar@9 1658 x = P->col[k-P->m]->prim;
alpar@9 1659 }
alpar@9 1660 y += val[t] * x;
alpar@9 1661 }
alpar@9 1662 /* check if the row is primal infeasible in the current basis,
alpar@9 1663 i.e. the constraint is violated at the current point */
alpar@9 1664 if (type == GLP_LO)
alpar@9 1665 { if (y >= rhs)
alpar@9 1666 { /* the constraint is not violated */
alpar@9 1667 ret = 1;
alpar@9 1668 goto done;
alpar@9 1669 }
alpar@9 1670 /* in the adjacent basis y goes to its lower bound */
alpar@9 1671 dir = +1;
alpar@9 1672 }
alpar@9 1673 else if (type == GLP_UP)
alpar@9 1674 { if (y <= rhs)
alpar@9 1675 { /* the constraint is not violated */
alpar@9 1676 ret = 1;
alpar@9 1677 goto done;
alpar@9 1678 }
alpar@9 1679 /* in the adjacent basis y goes to its upper bound */
alpar@9 1680 dir = -1;
alpar@9 1681 }
alpar@9 1682 else
alpar@9 1683 xerror("glp_analyze_row: type = %d; invalid parameter\n",
alpar@9 1684 type);
alpar@9 1685 /* compute dy = y.new - y.old */
alpar@9 1686 dy = rhs - y;
alpar@9 1687 /* perform dual ratio test to determine which non-basic variable
alpar@9 1688 should enter the adjacent basis to keep it dual feasible */
alpar@9 1689 piv = glp_dual_rtest(P, len, ind, val, dir, eps);
alpar@9 1690 if (piv == 0)
alpar@9 1691 { /* no dual feasible adjacent basis exists */
alpar@9 1692 ret = 2;
alpar@9 1693 goto done;
alpar@9 1694 }
alpar@9 1695 /* non-basic variable x[k] = xN[j] should enter the basis */
alpar@9 1696 k = ind[piv];
alpar@9 1697 xassert(1 <= k && k <= P->m+P->n);
alpar@9 1698 /* determine its value in the current basis */
alpar@9 1699 if (k <= P->m)
alpar@9 1700 x = P->row[k]->prim;
alpar@9 1701 else
alpar@9 1702 x = P->col[k-P->m]->prim;
alpar@9 1703 /* compute dx = x.new - x.old = dy / alfa[j] */
alpar@9 1704 xassert(val[piv] != 0.0);
alpar@9 1705 dx = dy / val[piv];
alpar@9 1706 /* compute dz = z.new - z.old = d[j] * dx, where d[j] is reduced
alpar@9 1707 cost of xN[j] in the current basis */
alpar@9 1708 if (k <= P->m)
alpar@9 1709 dz = P->row[k]->dual * dx;
alpar@9 1710 else
alpar@9 1711 dz = P->col[k-P->m]->dual * dx;
alpar@9 1712 /* store the analysis results */
alpar@9 1713 if (_piv != NULL) *_piv = piv;
alpar@9 1714 if (_x != NULL) *_x = x;
alpar@9 1715 if (_dx != NULL) *_dx = dx;
alpar@9 1716 if (_y != NULL) *_y = y;
alpar@9 1717 if (_dy != NULL) *_dy = dy;
alpar@9 1718 if (_dz != NULL) *_dz = dz;
alpar@9 1719 done: return ret;
alpar@9 1720 }
alpar@9 1721
alpar@9 1722 #if 0
alpar@9 1723 int main(void)
alpar@9 1724 { /* example program for the routine glp_analyze_row */
alpar@9 1725 glp_prob *P;
alpar@9 1726 glp_smcp parm;
alpar@9 1727 int i, k, len, piv, ret, ind[1+100];
alpar@9 1728 double rhs, x, dx, y, dy, dz, val[1+100];
alpar@9 1729 P = glp_create_prob();
alpar@9 1730 /* read plan.mps (see glpk/examples) */
alpar@9 1731 ret = glp_read_mps(P, GLP_MPS_DECK, NULL, "plan.mps");
alpar@9 1732 glp_assert(ret == 0);
alpar@9 1733 /* and solve it to optimality */
alpar@9 1734 ret = glp_simplex(P, NULL);
alpar@9 1735 glp_assert(ret == 0);
alpar@9 1736 glp_assert(glp_get_status(P) == GLP_OPT);
alpar@9 1737 /* the optimal objective value is 296.217 */
alpar@9 1738 /* we would like to know what happens if we would add a new row
alpar@9 1739 (constraint) to plan.mps:
alpar@9 1740 .01 * bin1 + .01 * bin2 + .02 * bin4 + .02 * bin5 <= 12 */
alpar@9 1741 /* first, we specify this new row */
alpar@9 1742 glp_create_index(P);
alpar@9 1743 len = 0;
alpar@9 1744 ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
alpar@9 1745 ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
alpar@9 1746 ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
alpar@9 1747 ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
alpar@9 1748 rhs = 12;
alpar@9 1749 /* then we can compute value of the row (i.e. of its auxiliary
alpar@9 1750 variable) in the current basis to see if the constraint is
alpar@9 1751 violated */
alpar@9 1752 y = 0.0;
alpar@9 1753 for (k = 1; k <= len; k++)
alpar@9 1754 y += val[k] * glp_get_col_prim(P, ind[k]);
alpar@9 1755 glp_printf("y = %g\n", y);
alpar@9 1756 /* this prints y = 15.1372, so the constraint is violated, since
alpar@9 1757 we require that y <= rhs = 12 */
alpar@9 1758 /* now we transform the row to express it only through non-basic
alpar@9 1759 (auxiliary and artificial) variables */
alpar@9 1760 len = glp_transform_row(P, len, ind, val);
alpar@9 1761 /* finally, we simulate one step of the dual simplex method to
alpar@9 1762 obtain necessary information for the adjacent basis */
alpar@9 1763 ret = _glp_analyze_row(P, len, ind, val, GLP_UP, rhs, 1e-9, &piv,
alpar@9 1764 &x, &dx, &y, &dy, &dz);
alpar@9 1765 glp_assert(ret == 0);
alpar@9 1766 glp_printf("k = %d, x = %g; dx = %g; y = %g; dy = %g; dz = %g\n",
alpar@9 1767 ind[piv], x, dx, y, dy, dz);
alpar@9 1768 /* this prints dz = 5.64418 and means that in the adjacent basis
alpar@9 1769 the objective function would be 296.217 + 5.64418 = 301.861 */
alpar@9 1770 /* now we actually include the row into the problem object; note
alpar@9 1771 that the arrays ind and val are clobbered, so we need to build
alpar@9 1772 them once again */
alpar@9 1773 len = 0;
alpar@9 1774 ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
alpar@9 1775 ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
alpar@9 1776 ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
alpar@9 1777 ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
alpar@9 1778 rhs = 12;
alpar@9 1779 i = glp_add_rows(P, 1);
alpar@9 1780 glp_set_row_bnds(P, i, GLP_UP, 0, rhs);
alpar@9 1781 glp_set_mat_row(P, i, len, ind, val);
alpar@9 1782 /* and perform one dual simplex iteration */
alpar@9 1783 glp_init_smcp(&parm);
alpar@9 1784 parm.meth = GLP_DUAL;
alpar@9 1785 parm.it_lim = 1;
alpar@9 1786 glp_simplex(P, &parm);
alpar@9 1787 /* the current objective value is 301.861 */
alpar@9 1788 return 0;
alpar@9 1789 }
alpar@9 1790 #endif
alpar@9 1791
alpar@9 1792 /***********************************************************************
alpar@9 1793 * NAME
alpar@9 1794 *
alpar@9 1795 * glp_analyze_bound - analyze active bound of non-basic variable
alpar@9 1796 *
alpar@9 1797 * SYNOPSIS
alpar@9 1798 *
alpar@9 1799 * void glp_analyze_bound(glp_prob *P, int k, double *limit1, int *var1,
alpar@9 1800 * double *limit2, int *var2);
alpar@9 1801 *
alpar@9 1802 * DESCRIPTION
alpar@9 1803 *
alpar@9 1804 * The routine glp_analyze_bound analyzes the effect of varying the
alpar@9 1805 * active bound of specified non-basic variable.
alpar@9 1806 *
alpar@9 1807 * The non-basic variable is specified by the parameter k, where
alpar@9 1808 * 1 <= k <= m means auxiliary variable of corresponding row while
alpar@9 1809 * m+1 <= k <= m+n means structural variable (column).
alpar@9 1810 *
alpar@9 1811 * Note that the current basic solution must be optimal, and the basis
alpar@9 1812 * factorization must exist.
alpar@9 1813 *
alpar@9 1814 * Results of the analysis have the following meaning.
alpar@9 1815 *
alpar@9 1816 * value1 is the minimal value of the active bound, at which the basis
alpar@9 1817 * still remains primal feasible and thus optimal. -DBL_MAX means that
alpar@9 1818 * the active bound has no lower limit.
alpar@9 1819 *
alpar@9 1820 * var1 is the ordinal number of an auxiliary (1 to m) or structural
alpar@9 1821 * (m+1 to n) basic variable, which reaches its bound first and thereby
alpar@9 1822 * limits further decreasing the active bound being analyzed.
alpar@9 1823 * if value1 = -DBL_MAX, var1 is set to 0.
alpar@9 1824 *
alpar@9 1825 * value2 is the maximal value of the active bound, at which the basis
alpar@9 1826 * still remains primal feasible and thus optimal. +DBL_MAX means that
alpar@9 1827 * the active bound has no upper limit.
alpar@9 1828 *
alpar@9 1829 * var2 is the ordinal number of an auxiliary (1 to m) or structural
alpar@9 1830 * (m+1 to n) basic variable, which reaches its bound first and thereby
alpar@9 1831 * limits further increasing the active bound being analyzed.
alpar@9 1832 * if value2 = +DBL_MAX, var2 is set to 0. */
alpar@9 1833
alpar@9 1834 void glp_analyze_bound(glp_prob *P, int k, double *value1, int *var1,
alpar@9 1835 double *value2, int *var2)
alpar@9 1836 { GLPROW *row;
alpar@9 1837 GLPCOL *col;
alpar@9 1838 int m, n, stat, kase, p, len, piv, *ind;
alpar@9 1839 double x, new_x, ll, uu, xx, delta, *val;
alpar@9 1840 /* sanity checks */
alpar@9 1841 if (P == NULL || P->magic != GLP_PROB_MAGIC)
alpar@9 1842 xerror("glp_analyze_bound: P = %p; invalid problem object\n",
alpar@9 1843 P);
alpar@9 1844 m = P->m, n = P->n;
alpar@9 1845 if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
alpar@9 1846 xerror("glp_analyze_bound: optimal basic solution required\n");
alpar@9 1847 if (!(m == 0 || P->valid))
alpar@9 1848 xerror("glp_analyze_bound: basis factorization required\n");
alpar@9 1849 if (!(1 <= k && k <= m+n))
alpar@9 1850 xerror("glp_analyze_bound: k = %d; variable number out of rang"
alpar@9 1851 "e\n", k);
alpar@9 1852 /* retrieve information about the specified non-basic variable
alpar@9 1853 x[k] whose active bound is to be analyzed */
alpar@9 1854 if (k <= m)
alpar@9 1855 { row = P->row[k];
alpar@9 1856 stat = row->stat;
alpar@9 1857 x = row->prim;
alpar@9 1858 }
alpar@9 1859 else
alpar@9 1860 { col = P->col[k-m];
alpar@9 1861 stat = col->stat;
alpar@9 1862 x = col->prim;
alpar@9 1863 }
alpar@9 1864 if (stat == GLP_BS)
alpar@9 1865 xerror("glp_analyze_bound: k = %d; basic variable not allowed "
alpar@9 1866 "\n", k);
alpar@9 1867 /* allocate working arrays */
alpar@9 1868 ind = xcalloc(1+m, sizeof(int));
alpar@9 1869 val = xcalloc(1+m, sizeof(double));
alpar@9 1870 /* compute column of the simplex table corresponding to the
alpar@9 1871 non-basic variable x[k] */
alpar@9 1872 len = glp_eval_tab_col(P, k, ind, val);
alpar@9 1873 xassert(0 <= len && len <= m);
alpar@9 1874 /* perform analysis */
alpar@9 1875 for (kase = -1; kase <= +1; kase += 2)
alpar@9 1876 { /* kase < 0 means active bound of x[k] is decreasing;
alpar@9 1877 kase > 0 means active bound of x[k] is increasing */
alpar@9 1878 /* use the primal ratio test to determine some basic variable
alpar@9 1879 x[p] which reaches its bound first */
alpar@9 1880 piv = glp_prim_rtest(P, len, ind, val, kase, 1e-9);
alpar@9 1881 if (piv == 0)
alpar@9 1882 { /* nothing limits changing the active bound of x[k] */
alpar@9 1883 p = 0;
alpar@9 1884 new_x = (kase < 0 ? -DBL_MAX : +DBL_MAX);
alpar@9 1885 goto store;
alpar@9 1886 }
alpar@9 1887 /* basic variable x[p] limits changing the active bound of
alpar@9 1888 x[k]; determine its value in the current basis */
alpar@9 1889 xassert(1 <= piv && piv <= len);
alpar@9 1890 p = ind[piv];
alpar@9 1891 if (p <= m)
alpar@9 1892 { row = P->row[p];
alpar@9 1893 ll = glp_get_row_lb(P, row->i);
alpar@9 1894 uu = glp_get_row_ub(P, row->i);
alpar@9 1895 stat = row->stat;
alpar@9 1896 xx = row->prim;
alpar@9 1897 }
alpar@9 1898 else
alpar@9 1899 { col = P->col[p-m];
alpar@9 1900 ll = glp_get_col_lb(P, col->j);
alpar@9 1901 uu = glp_get_col_ub(P, col->j);
alpar@9 1902 stat = col->stat;
alpar@9 1903 xx = col->prim;
alpar@9 1904 }
alpar@9 1905 xassert(stat == GLP_BS);
alpar@9 1906 /* determine delta x[p] = bound of x[p] - value of x[p] */
alpar@9 1907 if (kase < 0 && val[piv] > 0.0 ||
alpar@9 1908 kase > 0 && val[piv] < 0.0)
alpar@9 1909 { /* delta x[p] < 0, so x[p] goes toward its lower bound */
alpar@9 1910 xassert(ll != -DBL_MAX);
alpar@9 1911 delta = ll - xx;
alpar@9 1912 }
alpar@9 1913 else
alpar@9 1914 { /* delta x[p] > 0, so x[p] goes toward its upper bound */
alpar@9 1915 xassert(uu != +DBL_MAX);
alpar@9 1916 delta = uu - xx;
alpar@9 1917 }
alpar@9 1918 /* delta x[p] = alfa[p,k] * delta x[k], so new x[k] = x[k] +
alpar@9 1919 delta x[k] = x[k] + delta x[p] / alfa[p,k] is the value of
alpar@9 1920 x[k] in the adjacent basis */
alpar@9 1921 xassert(val[piv] != 0.0);
alpar@9 1922 new_x = x + delta / val[piv];
alpar@9 1923 store: /* store analysis results */
alpar@9 1924 if (kase < 0)
alpar@9 1925 { if (value1 != NULL) *value1 = new_x;
alpar@9 1926 if (var1 != NULL) *var1 = p;
alpar@9 1927 }
alpar@9 1928 else
alpar@9 1929 { if (value2 != NULL) *value2 = new_x;
alpar@9 1930 if (var2 != NULL) *var2 = p;
alpar@9 1931 }
alpar@9 1932 }
alpar@9 1933 /* free working arrays */
alpar@9 1934 xfree(ind);
alpar@9 1935 xfree(val);
alpar@9 1936 return;
alpar@9 1937 }
alpar@9 1938
alpar@9 1939 /***********************************************************************
alpar@9 1940 * NAME
alpar@9 1941 *
alpar@9 1942 * glp_analyze_coef - analyze objective coefficient at basic variable
alpar@9 1943 *
alpar@9 1944 * SYNOPSIS
alpar@9 1945 *
alpar@9 1946 * void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
alpar@9 1947 * double *value1, double *coef2, int *var2, double *value2);
alpar@9 1948 *
alpar@9 1949 * DESCRIPTION
alpar@9 1950 *
alpar@9 1951 * The routine glp_analyze_coef analyzes the effect of varying the
alpar@9 1952 * objective coefficient at specified basic variable.
alpar@9 1953 *
alpar@9 1954 * The basic variable is specified by the parameter k, where
alpar@9 1955 * 1 <= k <= m means auxiliary variable of corresponding row while
alpar@9 1956 * m+1 <= k <= m+n means structural variable (column).
alpar@9 1957 *
alpar@9 1958 * Note that the current basic solution must be optimal, and the basis
alpar@9 1959 * factorization must exist.
alpar@9 1960 *
alpar@9 1961 * Results of the analysis have the following meaning.
alpar@9 1962 *
alpar@9 1963 * coef1 is the minimal value of the objective coefficient, at which
alpar@9 1964 * the basis still remains dual feasible and thus optimal. -DBL_MAX
alpar@9 1965 * means that the objective coefficient has no lower limit.
alpar@9 1966 *
alpar@9 1967 * var1 is the ordinal number of an auxiliary (1 to m) or structural
alpar@9 1968 * (m+1 to n) non-basic variable, whose reduced cost reaches its zero
alpar@9 1969 * bound first and thereby limits further decreasing the objective
alpar@9 1970 * coefficient being analyzed. If coef1 = -DBL_MAX, var1 is set to 0.
alpar@9 1971 *
alpar@9 1972 * value1 is value of the basic variable being analyzed in an adjacent
alpar@9 1973 * basis, which is defined as follows. Let the objective coefficient
alpar@9 1974 * reaches its minimal value (coef1) and continues decreasing. Then the
alpar@9 1975 * reduced cost of the limiting non-basic variable (var1) becomes dual
alpar@9 1976 * infeasible and the current basis becomes non-optimal that forces the
alpar@9 1977 * limiting non-basic variable to enter the basis replacing there some
alpar@9 1978 * basic variable that leaves the basis to keep primal feasibility.
alpar@9 1979 * Should note that on determining the adjacent basis current bounds
alpar@9 1980 * of the basic variable being analyzed are ignored as if it were free
alpar@9 1981 * (unbounded) variable, so it cannot leave the basis. It may happen
alpar@9 1982 * that no dual feasible adjacent basis exists, in which case value1 is
alpar@9 1983 * set to -DBL_MAX or +DBL_MAX.
alpar@9 1984 *
alpar@9 1985 * coef2 is the maximal value of the objective coefficient, at which
alpar@9 1986 * the basis still remains dual feasible and thus optimal. +DBL_MAX
alpar@9 1987 * means that the objective coefficient has no upper limit.
alpar@9 1988 *
alpar@9 1989 * var2 is the ordinal number of an auxiliary (1 to m) or structural
alpar@9 1990 * (m+1 to n) non-basic variable, whose reduced cost reaches its zero
alpar@9 1991 * bound first and thereby limits further increasing the objective
alpar@9 1992 * coefficient being analyzed. If coef2 = +DBL_MAX, var2 is set to 0.
alpar@9 1993 *
alpar@9 1994 * value2 is value of the basic variable being analyzed in an adjacent
alpar@9 1995 * basis, which is defined exactly in the same way as value1 above with
alpar@9 1996 * exception that now the objective coefficient is increasing. */
alpar@9 1997
alpar@9 1998 void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
alpar@9 1999 double *value1, double *coef2, int *var2, double *value2)
alpar@9 2000 { GLPROW *row; GLPCOL *col;
alpar@9 2001 int m, n, type, stat, kase, p, q, dir, clen, cpiv, rlen, rpiv,
alpar@9 2002 *cind, *rind;
alpar@9 2003 double lb, ub, coef, x, lim_coef, new_x, d, delta, ll, uu, xx,
alpar@9 2004 *rval, *cval;
alpar@9 2005 /* sanity checks */
alpar@9 2006 if (P == NULL || P->magic != GLP_PROB_MAGIC)
alpar@9 2007 xerror("glp_analyze_coef: P = %p; invalid problem object\n",
alpar@9 2008 P);
alpar@9 2009 m = P->m, n = P->n;
alpar@9 2010 if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
alpar@9 2011 xerror("glp_analyze_coef: optimal basic solution required\n");
alpar@9 2012 if (!(m == 0 || P->valid))
alpar@9 2013 xerror("glp_analyze_coef: basis factorization required\n");
alpar@9 2014 if (!(1 <= k && k <= m+n))
alpar@9 2015 xerror("glp_analyze_coef: k = %d; variable number out of range"
alpar@9 2016 "\n", k);
alpar@9 2017 /* retrieve information about the specified basic variable x[k]
alpar@9 2018 whose objective coefficient c[k] is to be analyzed */
alpar@9 2019 if (k <= m)
alpar@9 2020 { row = P->row[k];
alpar@9 2021 type = row->type;
alpar@9 2022 lb = row->lb;
alpar@9 2023 ub = row->ub;
alpar@9 2024 coef = 0.0;
alpar@9 2025 stat = row->stat;
alpar@9 2026 x = row->prim;
alpar@9 2027 }
alpar@9 2028 else
alpar@9 2029 { col = P->col[k-m];
alpar@9 2030 type = col->type;
alpar@9 2031 lb = col->lb;
alpar@9 2032 ub = col->ub;
alpar@9 2033 coef = col->coef;
alpar@9 2034 stat = col->stat;
alpar@9 2035 x = col->prim;
alpar@9 2036 }
alpar@9 2037 if (stat != GLP_BS)
alpar@9 2038 xerror("glp_analyze_coef: k = %d; non-basic variable not allow"
alpar@9 2039 "ed\n", k);
alpar@9 2040 /* allocate working arrays */
alpar@9 2041 cind = xcalloc(1+m, sizeof(int));
alpar@9 2042 cval = xcalloc(1+m, sizeof(double));
alpar@9 2043 rind = xcalloc(1+n, sizeof(int));
alpar@9 2044 rval = xcalloc(1+n, sizeof(double));
alpar@9 2045 /* compute row of the simplex table corresponding to the basic
alpar@9 2046 variable x[k] */
alpar@9 2047 rlen = glp_eval_tab_row(P, k, rind, rval);
alpar@9 2048 xassert(0 <= rlen && rlen <= n);
alpar@9 2049 /* perform analysis */
alpar@9 2050 for (kase = -1; kase <= +1; kase += 2)
alpar@9 2051 { /* kase < 0 means objective coefficient c[k] is decreasing;
alpar@9 2052 kase > 0 means objective coefficient c[k] is increasing */
alpar@9 2053 /* note that decreasing c[k] is equivalent to increasing dual
alpar@9 2054 variable lambda[k] and vice versa; we need to correctly set
alpar@9 2055 the dir flag as required by the routine glp_dual_rtest */
alpar@9 2056 if (P->dir == GLP_MIN)
alpar@9 2057 dir = - kase;
alpar@9 2058 else if (P->dir == GLP_MAX)
alpar@9 2059 dir = + kase;
alpar@9 2060 else
alpar@9 2061 xassert(P != P);
alpar@9 2062 /* use the dual ratio test to determine non-basic variable
alpar@9 2063 x[q] whose reduced cost d[q] reaches zero bound first */
alpar@9 2064 rpiv = glp_dual_rtest(P, rlen, rind, rval, dir, 1e-9);
alpar@9 2065 if (rpiv == 0)
alpar@9 2066 { /* nothing limits changing c[k] */
alpar@9 2067 lim_coef = (kase < 0 ? -DBL_MAX : +DBL_MAX);
alpar@9 2068 q = 0;
alpar@9 2069 /* x[k] keeps its current value */
alpar@9 2070 new_x = x;
alpar@9 2071 goto store;
alpar@9 2072 }
alpar@9 2073 /* non-basic variable x[q] limits changing coefficient c[k];
alpar@9 2074 determine its status and reduced cost d[k] in the current
alpar@9 2075 basis */
alpar@9 2076 xassert(1 <= rpiv && rpiv <= rlen);
alpar@9 2077 q = rind[rpiv];
alpar@9 2078 xassert(1 <= q && q <= m+n);
alpar@9 2079 if (q <= m)
alpar@9 2080 { row = P->row[q];
alpar@9 2081 stat = row->stat;
alpar@9 2082 d = row->dual;
alpar@9 2083 }
alpar@9 2084 else
alpar@9 2085 { col = P->col[q-m];
alpar@9 2086 stat = col->stat;
alpar@9 2087 d = col->dual;
alpar@9 2088 }
alpar@9 2089 /* note that delta d[q] = new d[q] - d[q] = - d[q], because
alpar@9 2090 new d[q] = 0; delta d[q] = alfa[k,q] * delta c[k], so
alpar@9 2091 delta c[k] = delta d[q] / alfa[k,q] = - d[q] / alfa[k,q] */
alpar@9 2092 xassert(rval[rpiv] != 0.0);
alpar@9 2093 delta = - d / rval[rpiv];
alpar@9 2094 /* compute new c[k] = c[k] + delta c[k], which is the limiting
alpar@9 2095 value of the objective coefficient c[k] */
alpar@9 2096 lim_coef = coef + delta;
alpar@9 2097 /* let c[k] continue decreasing/increasing that makes d[q]
alpar@9 2098 dual infeasible and forces x[q] to enter the basis;
alpar@9 2099 to perform the primal ratio test we need to know in which
alpar@9 2100 direction x[q] changes on entering the basis; we determine
alpar@9 2101 that analyzing the sign of delta d[q] (see above), since
alpar@9 2102 d[q] may be close to zero having wrong sign */
alpar@9 2103 /* let, for simplicity, the problem is minimization */
alpar@9 2104 if (kase < 0 && rval[rpiv] > 0.0 ||
alpar@9 2105 kase > 0 && rval[rpiv] < 0.0)
alpar@9 2106 { /* delta d[q] < 0, so d[q] being non-negative will become
alpar@9 2107 negative, so x[q] will increase */
alpar@9 2108 dir = +1;
alpar@9 2109 }
alpar@9 2110 else
alpar@9 2111 { /* delta d[q] > 0, so d[q] being non-positive will become
alpar@9 2112 positive, so x[q] will decrease */
alpar@9 2113 dir = -1;
alpar@9 2114 }
alpar@9 2115 /* if the problem is maximization, correct the direction */
alpar@9 2116 if (P->dir == GLP_MAX) dir = - dir;
alpar@9 2117 /* check that we didn't make a silly mistake */
alpar@9 2118 if (dir > 0)
alpar@9 2119 xassert(stat == GLP_NL || stat == GLP_NF);
alpar@9 2120 else
alpar@9 2121 xassert(stat == GLP_NU || stat == GLP_NF);
alpar@9 2122 /* compute column of the simplex table corresponding to the
alpar@9 2123 non-basic variable x[q] */
alpar@9 2124 clen = glp_eval_tab_col(P, q, cind, cval);
alpar@9 2125 /* make x[k] temporarily free (unbounded) */
alpar@9 2126 if (k <= m)
alpar@9 2127 { row = P->row[k];
alpar@9 2128 row->type = GLP_FR;
alpar@9 2129 row->lb = row->ub = 0.0;
alpar@9 2130 }
alpar@9 2131 else
alpar@9 2132 { col = P->col[k-m];
alpar@9 2133 col->type = GLP_FR;
alpar@9 2134 col->lb = col->ub = 0.0;
alpar@9 2135 }
alpar@9 2136 /* use the primal ratio test to determine some basic variable
alpar@9 2137 which leaves the basis */
alpar@9 2138 cpiv = glp_prim_rtest(P, clen, cind, cval, dir, 1e-9);
alpar@9 2139 /* restore original bounds of the basic variable x[k] */
alpar@9 2140 if (k <= m)
alpar@9 2141 { row = P->row[k];
alpar@9 2142 row->type = type;
alpar@9 2143 row->lb = lb, row->ub = ub;
alpar@9 2144 }
alpar@9 2145 else
alpar@9 2146 { col = P->col[k-m];
alpar@9 2147 col->type = type;
alpar@9 2148 col->lb = lb, col->ub = ub;
alpar@9 2149 }
alpar@9 2150 if (cpiv == 0)
alpar@9 2151 { /* non-basic variable x[q] can change unlimitedly */
alpar@9 2152 if (dir < 0 && rval[rpiv] > 0.0 ||
alpar@9 2153 dir > 0 && rval[rpiv] < 0.0)
alpar@9 2154 { /* delta x[k] = alfa[k,q] * delta x[q] < 0 */
alpar@9 2155 new_x = -DBL_MAX;
alpar@9 2156 }
alpar@9 2157 else
alpar@9 2158 { /* delta x[k] = alfa[k,q] * delta x[q] > 0 */
alpar@9 2159 new_x = +DBL_MAX;
alpar@9 2160 }
alpar@9 2161 goto store;
alpar@9 2162 }
alpar@9 2163 /* some basic variable x[p] limits changing non-basic variable
alpar@9 2164 x[q] in the adjacent basis */
alpar@9 2165 xassert(1 <= cpiv && cpiv <= clen);
alpar@9 2166 p = cind[cpiv];
alpar@9 2167 xassert(1 <= p && p <= m+n);
alpar@9 2168 xassert(p != k);
alpar@9 2169 if (p <= m)
alpar@9 2170 { row = P->row[p];
alpar@9 2171 xassert(row->stat == GLP_BS);
alpar@9 2172 ll = glp_get_row_lb(P, row->i);
alpar@9 2173 uu = glp_get_row_ub(P, row->i);
alpar@9 2174 xx = row->prim;
alpar@9 2175 }
alpar@9 2176 else
alpar@9 2177 { col = P->col[p-m];
alpar@9 2178 xassert(col->stat == GLP_BS);
alpar@9 2179 ll = glp_get_col_lb(P, col->j);
alpar@9 2180 uu = glp_get_col_ub(P, col->j);
alpar@9 2181 xx = col->prim;
alpar@9 2182 }
alpar@9 2183 /* determine delta x[p] = new x[p] - x[p] */
alpar@9 2184 if (dir < 0 && cval[cpiv] > 0.0 ||
alpar@9 2185 dir > 0 && cval[cpiv] < 0.0)
alpar@9 2186 { /* delta x[p] < 0, so x[p] goes toward its lower bound */
alpar@9 2187 xassert(ll != -DBL_MAX);
alpar@9 2188 delta = ll - xx;
alpar@9 2189 }
alpar@9 2190 else
alpar@9 2191 { /* delta x[p] > 0, so x[p] goes toward its upper bound */
alpar@9 2192 xassert(uu != +DBL_MAX);
alpar@9 2193 delta = uu - xx;
alpar@9 2194 }
alpar@9 2195 /* compute new x[k] = x[k] + alfa[k,q] * delta x[q], where
alpar@9 2196 delta x[q] = delta x[p] / alfa[p,q] */
alpar@9 2197 xassert(cval[cpiv] != 0.0);
alpar@9 2198 new_x = x + (rval[rpiv] / cval[cpiv]) * delta;
alpar@9 2199 store: /* store analysis results */
alpar@9 2200 if (kase < 0)
alpar@9 2201 { if (coef1 != NULL) *coef1 = lim_coef;
alpar@9 2202 if (var1 != NULL) *var1 = q;
alpar@9 2203 if (value1 != NULL) *value1 = new_x;
alpar@9 2204 }
alpar@9 2205 else
alpar@9 2206 { if (coef2 != NULL) *coef2 = lim_coef;
alpar@9 2207 if (var2 != NULL) *var2 = q;
alpar@9 2208 if (value2 != NULL) *value2 = new_x;
alpar@9 2209 }
alpar@9 2210 }
alpar@9 2211 /* free working arrays */
alpar@9 2212 xfree(cind);
alpar@9 2213 xfree(cval);
alpar@9 2214 xfree(rind);
alpar@9 2215 xfree(rval);
alpar@9 2216 return;
alpar@9 2217 }
alpar@9 2218
alpar@9 2219 /* eof */