1 /* glpssx.h (simplex method, rational arithmetic) */ 2 3 /*********************************************************************** 4 * This code is part of GLPK (GNU Linear Programming Kit). 5 * Copyright (C) 2003-2013 Free Software Foundation, Inc. 6 * Written by Andrew Makhorin <mao@gnu.org>. 7 * 8 * GLPK is free software: you can redistribute it and/or modify it 9 * under the terms of the GNU General Public License as published by 10 * the Free Software Foundation, either version 3 of the License, or 11 * (at your option) any later version. 12 * 13 * GLPK is distributed in the hope that it will be useful, but WITHOUT 14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public 16 * License for more details. 17 * 18 * You should have received a copy of the GNU General Public License 19 * along with GLPK. If not, see <http://www.gnu.org/licenses/>. 20 ***********************************************************************/ 21 22 #ifndef GLPSSX_H 23 #define GLPSSX_H 24 25 #include "bfx.h" 26 #include "env.h" 27 #if 1 /* 25/XI-2017 */ 28 #include "glpk.h" 29 #endif 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 #if 1 /* 25/XI-2017 */ 309 int msg_lev; 310 /* verbosity level: 311 GLP_MSG_OFF no output 312 GLP_MSG_ERR report errors and warnings 313 GLP_MSG_ON normal output 314 GLP_MSG_ALL highest verbosity */ 315 #endif 316 int it_lim; 317 /* simplex iterations limit; if this value is positive, it is 318 decreased by one each time when one simplex iteration has been 319 performed, and reaching zero value signals the solver to stop 320 the search; negative value means no iterations limit */ 321 int it_cnt; 322 /* simplex iterations count; this count is increased by one each 323 time when one simplex iteration has been performed */ 324 double tm_lim; 325 /* searching time limit, in seconds; if this value is positive, 326 it is decreased each time when one simplex iteration has been 327 performed by the amount of time spent for the iteration, and 328 reaching zero value signals the solver to stop the search; 329 negative value means no time limit */ 330 double out_frq; 331 /* output frequency, in seconds; this parameter specifies how 332 frequently the solver sends information about the progress of 333 the search to the standard output */ 334 #if 0 /* 10/VI-2013 */ 335 glp_long tm_beg; 336 #else 337 double tm_beg; 338 #endif 339 /* starting time of the search, in seconds; the total time of the 340 search is the difference between xtime() and tm_beg */ 341 #if 0 /* 10/VI-2013 */ 342 glp_long tm_lag; 343 #else 344 double tm_lag; 345 #endif 346 /* the most recent time, in seconds, at which the progress of the 347 the search was displayed */ 348 }; 349 350 #define ssx_create _glp_ssx_create 351 #define ssx_factorize _glp_ssx_factorize 352 #define ssx_get_xNj _glp_ssx_get_xNj 353 #define ssx_eval_bbar _glp_ssx_eval_bbar 354 #define ssx_eval_pi _glp_ssx_eval_pi 355 #define ssx_eval_dj _glp_ssx_eval_dj 356 #define ssx_eval_cbar _glp_ssx_eval_cbar 357 #define ssx_eval_rho _glp_ssx_eval_rho 358 #define ssx_eval_row _glp_ssx_eval_row 359 #define ssx_eval_col _glp_ssx_eval_col 360 #define ssx_chuzc _glp_ssx_chuzc 361 #define ssx_chuzr _glp_ssx_chuzr 362 #define ssx_update_bbar _glp_ssx_update_bbar 363 #define ssx_update_pi _glp_ssx_update_pi 364 #define ssx_update_cbar _glp_ssx_update_cbar 365 #define ssx_change_basis _glp_ssx_change_basis 366 #define ssx_delete _glp_ssx_delete 367 368 #define ssx_phase_I _glp_ssx_phase_I 369 #define ssx_phase_II _glp_ssx_phase_II 370 #define ssx_driver _glp_ssx_driver 371 372 SSX *ssx_create(int m, int n, int nnz); 373 /* create simplex solver workspace */ 374 375 int ssx_factorize(SSX *ssx); 376 /* factorize the current basis matrix */ 377 378 void ssx_get_xNj(SSX *ssx, int j, mpq_t x); 379 /* determine value of non-basic variable */ 380 381 void ssx_eval_bbar(SSX *ssx); 382 /* compute values of basic variables */ 383 384 void ssx_eval_pi(SSX *ssx); 385 /* compute values of simplex multipliers */ 386 387 void ssx_eval_dj(SSX *ssx, int j, mpq_t dj); 388 /* compute reduced cost of non-basic variable */ 389 390 void ssx_eval_cbar(SSX *ssx); 391 /* compute reduced costs of all non-basic variables */ 392 393 void ssx_eval_rho(SSX *ssx); 394 /* compute p-th row of the inverse */ 395 396 void ssx_eval_row(SSX *ssx); 397 /* compute pivot row of the simplex table */ 398 399 void ssx_eval_col(SSX *ssx); 400 /* compute pivot column of the simplex table */ 401 402 void ssx_chuzc(SSX *ssx); 403 /* choose pivot column */ 404 405 void ssx_chuzr(SSX *ssx); 406 /* choose pivot row */ 407 408 void ssx_update_bbar(SSX *ssx); 409 /* update values of basic variables */ 410 411 void ssx_update_pi(SSX *ssx); 412 /* update simplex multipliers */ 413 414 void ssx_update_cbar(SSX *ssx); 415 /* update reduced costs of non-basic variables */ 416 417 void ssx_change_basis(SSX *ssx); 418 /* change current basis to adjacent one */ 419 420 void ssx_delete(SSX *ssx); 421 /* delete simplex solver workspace */ 422 423 int ssx_phase_I(SSX *ssx); 424 /* find primal feasible solution */ 425 426 int ssx_phase_II(SSX *ssx); 427 /* find optimal solution */ 428 429 int ssx_driver(SSX *ssx); 430 /* base driver to exact simplex method */ 431 432 #endif 433 434 /* eof */ 435