1 /* glpssx.h (simplex method, bignum arithmetic) */ 2 3 /*********************************************************************** 4 * This code is part of GLPK (GNU Linear Programming Kit). 5 * 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 7 * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved. 9 * E-mail: <mao@gnu.org>. 10 * 11 * GLPK is free software: you can redistribute it and/or modify it 12 * under the terms of the GNU General Public License as published by 13 * the Free Software Foundation, either version 3 of the License, or 14 * (at your option) any later version. 15 * 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public 19 * License for more details. 20 * 21 * You should have received a copy of the GNU General Public License 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>. 23 ***********************************************************************/ 24 25 #ifndef GLPSSX_H 26 #define GLPSSX_H 27 28 #include "glpbfx.h" 29 #include "glpenv.h" 30 31 typedef struct SSX SSX; 32 33 struct SSX 34 { /* simplex solver workspace */ 35 /*---------------------------------------------------------------------- 36 // LP PROBLEM DATA 37 // 38 // It is assumed that LP problem has the following statement: 39 // 40 // minimize (or maximize) 41 // 42 // z = c[1]*x[1] + ... + c[m+n]*x[m+n] + c[0] (1) 43 // 44 // subject to equality constraints 45 // 46 // x[1] - a[1,1]*x[m+1] - ... - a[1,n]*x[m+n] = 0 47 // 48 // . . . . . . . (2) 49 // 50 // x[m] - a[m,1]*x[m+1] + ... - a[m,n]*x[m+n] = 0 51 // 52 // and bounds of variables 53 // 54 // l[1] <= x[1] <= u[1] 55 // 56 // . . . . . . . (3) 57 // 58 // l[m+n] <= x[m+n] <= u[m+n] 59 // 60 // where: 61 // x[1], ..., x[m] - auxiliary variables; 62 // x[m+1], ..., x[m+n] - structural variables; 63 // z - objective function; 64 // c[1], ..., c[m+n] - coefficients of the objective function; 65 // c[0] - constant term of the objective function; 66 // a[1,1], ..., a[m,n] - constraint coefficients; 67 // l[1], ..., l[m+n] - lower bounds of variables; 68 // u[1], ..., u[m+n] - upper bounds of variables. 69 // 70 // Bounds of variables can be finite as well as inifinite. Besides, 71 // lower and upper bounds can be equal to each other. So the following 72 // five types of variables are possible: 73 // 74 // Bounds of variable Type of variable 75 // ------------------------------------------------- 76 // -inf < x[k] < +inf Free (unbounded) variable 77 // l[k] <= x[k] < +inf Variable with lower bound 78 // -inf < x[k] <= u[k] Variable with upper bound 79 // l[k] <= x[k] <= u[k] Double-bounded variable 80 // l[k] = x[k] = u[k] Fixed variable 81 // 82 // Using vector-matrix notations the LP problem (1)-(3) can be written 83 // as follows: 84 // 85 // minimize (or maximize) 86 // 87 // z = c * x + c[0] (4) 88 // 89 // subject to equality constraints 90 // 91 // xR - A * xS = 0 (5) 92 // 93 // and bounds of variables 94 // 95 // l <= x <= u (6) 96 // 97 // where: 98 // xR - vector of auxiliary variables; 99 // xS - vector of structural variables; 100 // x = (xR, xS) - vector of all variables; 101 // z - objective function; 102 // c - vector of objective coefficients; 103 // c[0] - constant term of the objective function; 104 // A - matrix of constraint coefficients (has m rows 105 // and n columns); 106 // l - vector of lower bounds of variables; 107 // u - vector of upper bounds of variables. 108 // 109 // The simplex method makes no difference between auxiliary and 110 // structural variables, so it is convenient to think the system of 111 // equality constraints (5) written in a homogeneous form: 112 // 113 // (I | -A) * x = 0, (7) 114 // 115 // where (I | -A) is an augmented (m+n)xm constraint matrix, I is mxm 116 // unity matrix whose columns correspond to auxiliary variables, and A 117 // is the original mxn constraint matrix whose columns correspond to 118 // structural variables. Note that only the matrix A is stored. 119 ----------------------------------------------------------------------*/ 120 int m; 121 /* number of rows (auxiliary variables), m > 0 */ 122 int n; 123 /* number of columns (structural variables), n > 0 */ 124 int *type; /* int type[1+m+n]; */ 125 /* type[0] is not used; 126 type[k], 1 <= k <= m+n, is the type of variable x[k]: */ 127 #define SSX_FR 0 /* free (unbounded) variable */ 128 #define SSX_LO 1 /* variable with lower bound */ 129 #define SSX_UP 2 /* variable with upper bound */ 130 #define SSX_DB 3 /* double-bounded variable */ 131 #define SSX_FX 4 /* fixed variable */ 132 mpq_t *lb; /* mpq_t lb[1+m+n]; alias: l */ 133 /* lb[0] is not used; 134 lb[k], 1 <= k <= m+n, is an lower bound of variable x[k]; 135 if x[k] has no lower bound, lb[k] is zero */ 136 mpq_t *ub; /* mpq_t ub[1+m+n]; alias: u */ 137 /* ub[0] is not used; 138 ub[k], 1 <= k <= m+n, is an upper bound of variable x[k]; 139 if x[k] has no upper bound, ub[k] is zero; 140 if x[k] is of fixed type, ub[k] is equal to lb[k] */ 141 int dir; 142 /* optimization direction (sense of the objective function): */ 143 #define SSX_MIN 0 /* minimization */ 144 #define SSX_MAX 1 /* maximization */ 145 mpq_t *coef; /* mpq_t coef[1+m+n]; alias: c */ 146 /* coef[0] is a constant term of the objective function; 147 coef[k], 1 <= k <= m+n, is a coefficient of the objective 148 function at variable x[k]; 149 note that auxiliary variables also may have non-zero objective 150 coefficients */ 151 int *A_ptr; /* int A_ptr[1+n+1]; */ 152 int *A_ind; /* int A_ind[A_ptr[n+1]]; */ 153 mpq_t *A_val; /* mpq_t A_val[A_ptr[n+1]]; */ 154 /* constraint matrix A (see (5)) in storage-by-columns format */ 155 /*---------------------------------------------------------------------- 156 // LP BASIS AND CURRENT BASIC SOLUTION 157 // 158 // The LP basis is defined by the following partition of the augmented 159 // constraint matrix (7): 160 // 161 // (B | N) = (I | -A) * Q, (8) 162 // 163 // where B is a mxm non-singular basis matrix whose columns correspond 164 // to basic variables xB, N is a mxn matrix whose columns correspond to 165 // non-basic variables xN, and Q is a permutation (m+n)x(m+n) matrix. 166 // 167 // From (7) and (8) it follows that 168 // 169 // (I | -A) * x = (I | -A) * Q * Q' * x = (B | N) * (xB, xN), 170 // 171 // therefore 172 // 173 // (xB, xN) = Q' * x, (9) 174 // 175 // where x is the vector of all variables in the original order, xB is 176 // a vector of basic variables, xN is a vector of non-basic variables, 177 // Q' = inv(Q) is a matrix transposed to Q. 178 // 179 // Current values of non-basic variables xN[j], j = 1, ..., n, are not 180 // stored; they are defined implicitly by their statuses as follows: 181 // 182 // 0, if xN[j] is free variable 183 // lN[j], if xN[j] is on its lower bound (10) 184 // uN[j], if xN[j] is on its upper bound 185 // lN[j] = uN[j], if xN[j] is fixed variable 186 // 187 // where lN[j] and uN[j] are lower and upper bounds of xN[j]. 188 // 189 // Current values of basic variables xB[i], i = 1, ..., m, are computed 190 // as follows: 191 // 192 // beta = - inv(B) * N * xN, (11) 193 // 194 // where current values of xN are defined by (10). 195 // 196 // Current values of simplex multipliers pi[i], i = 1, ..., m (which 197 // are values of Lagrange multipliers for equality constraints (7) also 198 // called shadow prices) are computed as follows: 199 // 200 // pi = inv(B') * cB, (12) 201 // 202 // where B' is a matrix transposed to B, cB is a vector of objective 203 // coefficients at basic variables xB. 204 // 205 // Current values of reduced costs d[j], j = 1, ..., n, (which are 206 // values of Langrange multipliers for active inequality constraints 207 // corresponding to non-basic variables) are computed as follows: 208 // 209 // d = cN - N' * pi, (13) 210 // 211 // where N' is a matrix transposed to N, cN is a vector of objective 212 // coefficients at non-basic variables xN. 213 ----------------------------------------------------------------------*/ 214 int *stat; /* int stat[1+m+n]; */ 215 /* stat[0] is not used; 216 stat[k], 1 <= k <= m+n, is the status of variable x[k]: */ 217 #define SSX_BS 0 /* basic variable */ 218 #define SSX_NL 1 /* non-basic variable on lower bound */ 219 #define SSX_NU 2 /* non-basic variable on upper bound */ 220 #define SSX_NF 3 /* non-basic free variable */ 221 #define SSX_NS 4 /* non-basic fixed variable */ 222 int *Q_row; /* int Q_row[1+m+n]; */ 223 /* matrix Q in row-like format; 224 Q_row[0] is not used; 225 Q_row[i] = j means that q[i,j] = 1 */ 226 int *Q_col; /* int Q_col[1+m+n]; */ 227 /* matrix Q in column-like format; 228 Q_col[0] is not used; 229 Q_col[j] = i means that q[i,j] = 1 */ 230 /* if k-th column of the matrix (I | A) is k'-th column of the 231 matrix (B | N), then Q_row[k] = k' and Q_col[k'] = k; 232 if x[k] is xB[i], then Q_row[k] = i and Q_col[i] = k; 233 if x[k] is xN[j], then Q_row[k] = m+j and Q_col[m+j] = k */ 234 BFX *binv; 235 /* invertable form of the basis matrix B */ 236 mpq_t *bbar; /* mpq_t bbar[1+m]; alias: beta */ 237 /* bbar[0] is a value of the objective function; 238 bbar[i], 1 <= i <= m, is a value of basic variable xB[i] */ 239 mpq_t *pi; /* mpq_t pi[1+m]; */ 240 /* pi[0] is not used; 241 pi[i], 1 <= i <= m, is a simplex multiplier corresponding to 242 i-th row (equality constraint) */ 243 mpq_t *cbar; /* mpq_t cbar[1+n]; alias: d */ 244 /* cbar[0] is not used; 245 cbar[j], 1 <= j <= n, is a reduced cost of non-basic variable 246 xN[j] */ 247 /*---------------------------------------------------------------------- 248 // SIMPLEX TABLE 249 // 250 // Due to (8) and (9) the system of equality constraints (7) for the 251 // current basis can be written as follows: 252 // 253 // xB = A~ * xN, (14) 254 // 255 // where 256 // 257 // A~ = - inv(B) * N (15) 258 // 259 // is a mxn matrix called the simplex table. 260 // 261 // The revised simplex method uses only two components of A~, namely, 262 // pivot column corresponding to non-basic variable xN[q] chosen to 263 // enter the basis, and pivot row corresponding to basic variable xB[p] 264 // chosen to leave the basis. 265 // 266 // Pivot column alfa_q is q-th column of A~, so 267 // 268 // alfa_q = A~ * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q], (16) 269 // 270 // where N[q] is q-th column of the matrix N. 271 // 272 // Pivot row alfa_p is p-th row of A~ or, equivalently, p-th column of 273 // A~', a matrix transposed to A~, so 274 // 275 // alfa_p = A~' * e[p] = - N' * inv(B') * e[p] = - N' * rho_p, (17) 276 // 277 // where (*)' means transposition, and 278 // 279 // rho_p = inv(B') * e[p], (18) 280 // 281 // is p-th column of inv(B') or, that is the same, p-th row of inv(B). 282 ----------------------------------------------------------------------*/ 283 int p; 284 /* number of basic variable xB[p], 1 <= p <= m, chosen to leave 285 the basis */ 286 mpq_t *rho; /* mpq_t rho[1+m]; */ 287 /* p-th row of the inverse inv(B); see (18) */ 288 mpq_t *ap; /* mpq_t ap[1+n]; */ 289 /* p-th row of the simplex table; see (17) */ 290 int q; 291 /* number of non-basic variable xN[q], 1 <= q <= n, chosen to 292 enter the basis */ 293 mpq_t *aq; /* mpq_t aq[1+m]; */ 294 /* q-th column of the simplex table; see (16) */ 295 /*--------------------------------------------------------------------*/ 296 int q_dir; 297 /* direction in which non-basic variable xN[q] should change on 298 moving to the adjacent vertex of the polyhedron: 299 +1 means that xN[q] increases 300 -1 means that xN[q] decreases */ 301 int p_stat; 302 /* non-basic status which should be assigned to basic variable 303 xB[p] when it has left the basis and become xN[q] */ 304 mpq_t delta; 305 /* actual change of xN[q] in the adjacent basis (it has the same 306 sign as q_dir) */ 307 /*--------------------------------------------------------------------*/ 308 int it_lim; 309 /* simplex iterations limit; if this value is positive, it is 310 decreased by one each time when one simplex iteration has been 311 performed, and reaching zero value signals the solver to stop 312 the search; negative value means no iterations limit */ 313 int it_cnt; 314 /* simplex iterations count; this count is increased by one each 315 time when one simplex iteration has been performed */ 316 double tm_lim; 317 /* searching time limit, in seconds; if this value is positive, 318 it is decreased each time when one simplex iteration has been 319 performed by the amount of time spent for the iteration, and 320 reaching zero value signals the solver to stop the search; 321 negative value means no time limit */ 322 double out_frq; 323 /* output frequency, in seconds; this parameter specifies how 324 frequently the solver sends information about the progress of 325 the search to the standard output */ 326 glp_long tm_beg; 327 /* starting time of the search, in seconds; the total time of the 328 search is the difference between xtime() and tm_beg */ 329 glp_long tm_lag; 330 /* the most recent time, in seconds, at which the progress of the 331 the search was displayed */ 332 }; 333 334 #define ssx_create _glp_ssx_create 335 #define ssx_factorize _glp_ssx_factorize 336 #define ssx_get_xNj _glp_ssx_get_xNj 337 #define ssx_eval_bbar _glp_ssx_eval_bbar 338 #define ssx_eval_pi _glp_ssx_eval_pi 339 #define ssx_eval_dj _glp_ssx_eval_dj 340 #define ssx_eval_cbar _glp_ssx_eval_cbar 341 #define ssx_eval_rho _glp_ssx_eval_rho 342 #define ssx_eval_row _glp_ssx_eval_row 343 #define ssx_eval_col _glp_ssx_eval_col 344 #define ssx_chuzc _glp_ssx_chuzc 345 #define ssx_chuzr _glp_ssx_chuzr 346 #define ssx_update_bbar _glp_ssx_update_bbar 347 #define ssx_update_pi _glp_ssx_update_pi 348 #define ssx_update_cbar _glp_ssx_update_cbar 349 #define ssx_change_basis _glp_ssx_change_basis 350 #define ssx_delete _glp_ssx_delete 351 352 #define ssx_phase_I _glp_ssx_phase_I 353 #define ssx_phase_II _glp_ssx_phase_II 354 #define ssx_driver _glp_ssx_driver 355 356 SSX *ssx_create(int m, int n, int nnz); 357 /* create simplex solver workspace */ 358 359 int ssx_factorize(SSX *ssx); 360 /* factorize the current basis matrix */ 361 362 void ssx_get_xNj(SSX *ssx, int j, mpq_t x); 363 /* determine value of non-basic variable */ 364 365 void ssx_eval_bbar(SSX *ssx); 366 /* compute values of basic variables */ 367 368 void ssx_eval_pi(SSX *ssx); 369 /* compute values of simplex multipliers */ 370 371 void ssx_eval_dj(SSX *ssx, int j, mpq_t dj); 372 /* compute reduced cost of non-basic variable */ 373 374 void ssx_eval_cbar(SSX *ssx); 375 /* compute reduced costs of all non-basic variables */ 376 377 void ssx_eval_rho(SSX *ssx); 378 /* compute p-th row of the inverse */ 379 380 void ssx_eval_row(SSX *ssx); 381 /* compute pivot row of the simplex table */ 382 383 void ssx_eval_col(SSX *ssx); 384 /* compute pivot column of the simplex table */ 385 386 void ssx_chuzc(SSX *ssx); 387 /* choose pivot column */ 388 389 void ssx_chuzr(SSX *ssx); 390 /* choose pivot row */ 391 392 void ssx_update_bbar(SSX *ssx); 393 /* update values of basic variables */ 394 395 void ssx_update_pi(SSX *ssx); 396 /* update simplex multipliers */ 397 398 void ssx_update_cbar(SSX *ssx); 399 /* update reduced costs of non-basic variables */ 400 401 void ssx_change_basis(SSX *ssx); 402 /* change current basis to adjacent one */ 403 404 void ssx_delete(SSX *ssx); 405 /* delete simplex solver workspace */ 406 407 int ssx_phase_I(SSX *ssx); 408 /* find primal feasible solution */ 409 410 int ssx_phase_II(SSX *ssx); 411 /* find optimal solution */ 412 413 int ssx_driver(SSX *ssx); 414 /* base driver to exact simplex method */ 415 416 #endif 417 418 /* eof */ 419