1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpssx.h Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,418 @@
1.4 +/* glpssx.h (simplex method, bignum arithmetic) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#ifndef GLPSSX_H
1.29 +#define GLPSSX_H
1.30 +
1.31 +#include "glpbfx.h"
1.32 +#include "glpenv.h"
1.33 +
1.34 +typedef struct SSX SSX;
1.35 +
1.36 +struct SSX
1.37 +{ /* simplex solver workspace */
1.38 +/*----------------------------------------------------------------------
1.39 +// LP PROBLEM DATA
1.40 +//
1.41 +// It is assumed that LP problem has the following statement:
1.42 +//
1.43 +// minimize (or maximize)
1.44 +//
1.45 +// z = c[1]*x[1] + ... + c[m+n]*x[m+n] + c[0] (1)
1.46 +//
1.47 +// subject to equality constraints
1.48 +//
1.49 +// x[1] - a[1,1]*x[m+1] - ... - a[1,n]*x[m+n] = 0
1.50 +//
1.51 +// . . . . . . . (2)
1.52 +//
1.53 +// x[m] - a[m,1]*x[m+1] + ... - a[m,n]*x[m+n] = 0
1.54 +//
1.55 +// and bounds of variables
1.56 +//
1.57 +// l[1] <= x[1] <= u[1]
1.58 +//
1.59 +// . . . . . . . (3)
1.60 +//
1.61 +// l[m+n] <= x[m+n] <= u[m+n]
1.62 +//
1.63 +// where:
1.64 +// x[1], ..., x[m] - auxiliary variables;
1.65 +// x[m+1], ..., x[m+n] - structural variables;
1.66 +// z - objective function;
1.67 +// c[1], ..., c[m+n] - coefficients of the objective function;
1.68 +// c[0] - constant term of the objective function;
1.69 +// a[1,1], ..., a[m,n] - constraint coefficients;
1.70 +// l[1], ..., l[m+n] - lower bounds of variables;
1.71 +// u[1], ..., u[m+n] - upper bounds of variables.
1.72 +//
1.73 +// Bounds of variables can be finite as well as inifinite. Besides,
1.74 +// lower and upper bounds can be equal to each other. So the following
1.75 +// five types of variables are possible:
1.76 +//
1.77 +// Bounds of variable Type of variable
1.78 +// -------------------------------------------------
1.79 +// -inf < x[k] < +inf Free (unbounded) variable
1.80 +// l[k] <= x[k] < +inf Variable with lower bound
1.81 +// -inf < x[k] <= u[k] Variable with upper bound
1.82 +// l[k] <= x[k] <= u[k] Double-bounded variable
1.83 +// l[k] = x[k] = u[k] Fixed variable
1.84 +//
1.85 +// Using vector-matrix notations the LP problem (1)-(3) can be written
1.86 +// as follows:
1.87 +//
1.88 +// minimize (or maximize)
1.89 +//
1.90 +// z = c * x + c[0] (4)
1.91 +//
1.92 +// subject to equality constraints
1.93 +//
1.94 +// xR - A * xS = 0 (5)
1.95 +//
1.96 +// and bounds of variables
1.97 +//
1.98 +// l <= x <= u (6)
1.99 +//
1.100 +// where:
1.101 +// xR - vector of auxiliary variables;
1.102 +// xS - vector of structural variables;
1.103 +// x = (xR, xS) - vector of all variables;
1.104 +// z - objective function;
1.105 +// c - vector of objective coefficients;
1.106 +// c[0] - constant term of the objective function;
1.107 +// A - matrix of constraint coefficients (has m rows
1.108 +// and n columns);
1.109 +// l - vector of lower bounds of variables;
1.110 +// u - vector of upper bounds of variables.
1.111 +//
1.112 +// The simplex method makes no difference between auxiliary and
1.113 +// structural variables, so it is convenient to think the system of
1.114 +// equality constraints (5) written in a homogeneous form:
1.115 +//
1.116 +// (I | -A) * x = 0, (7)
1.117 +//
1.118 +// where (I | -A) is an augmented (m+n)xm constraint matrix, I is mxm
1.119 +// unity matrix whose columns correspond to auxiliary variables, and A
1.120 +// is the original mxn constraint matrix whose columns correspond to
1.121 +// structural variables. Note that only the matrix A is stored.
1.122 +----------------------------------------------------------------------*/
1.123 + int m;
1.124 + /* number of rows (auxiliary variables), m > 0 */
1.125 + int n;
1.126 + /* number of columns (structural variables), n > 0 */
1.127 + int *type; /* int type[1+m+n]; */
1.128 + /* type[0] is not used;
1.129 + type[k], 1 <= k <= m+n, is the type of variable x[k]: */
1.130 +#define SSX_FR 0 /* free (unbounded) variable */
1.131 +#define SSX_LO 1 /* variable with lower bound */
1.132 +#define SSX_UP 2 /* variable with upper bound */
1.133 +#define SSX_DB 3 /* double-bounded variable */
1.134 +#define SSX_FX 4 /* fixed variable */
1.135 + mpq_t *lb; /* mpq_t lb[1+m+n]; alias: l */
1.136 + /* lb[0] is not used;
1.137 + lb[k], 1 <= k <= m+n, is an lower bound of variable x[k];
1.138 + if x[k] has no lower bound, lb[k] is zero */
1.139 + mpq_t *ub; /* mpq_t ub[1+m+n]; alias: u */
1.140 + /* ub[0] is not used;
1.141 + ub[k], 1 <= k <= m+n, is an upper bound of variable x[k];
1.142 + if x[k] has no upper bound, ub[k] is zero;
1.143 + if x[k] is of fixed type, ub[k] is equal to lb[k] */
1.144 + int dir;
1.145 + /* optimization direction (sense of the objective function): */
1.146 +#define SSX_MIN 0 /* minimization */
1.147 +#define SSX_MAX 1 /* maximization */
1.148 + mpq_t *coef; /* mpq_t coef[1+m+n]; alias: c */
1.149 + /* coef[0] is a constant term of the objective function;
1.150 + coef[k], 1 <= k <= m+n, is a coefficient of the objective
1.151 + function at variable x[k];
1.152 + note that auxiliary variables also may have non-zero objective
1.153 + coefficients */
1.154 + int *A_ptr; /* int A_ptr[1+n+1]; */
1.155 + int *A_ind; /* int A_ind[A_ptr[n+1]]; */
1.156 + mpq_t *A_val; /* mpq_t A_val[A_ptr[n+1]]; */
1.157 + /* constraint matrix A (see (5)) in storage-by-columns format */
1.158 +/*----------------------------------------------------------------------
1.159 +// LP BASIS AND CURRENT BASIC SOLUTION
1.160 +//
1.161 +// The LP basis is defined by the following partition of the augmented
1.162 +// constraint matrix (7):
1.163 +//
1.164 +// (B | N) = (I | -A) * Q, (8)
1.165 +//
1.166 +// where B is a mxm non-singular basis matrix whose columns correspond
1.167 +// to basic variables xB, N is a mxn matrix whose columns correspond to
1.168 +// non-basic variables xN, and Q is a permutation (m+n)x(m+n) matrix.
1.169 +//
1.170 +// From (7) and (8) it follows that
1.171 +//
1.172 +// (I | -A) * x = (I | -A) * Q * Q' * x = (B | N) * (xB, xN),
1.173 +//
1.174 +// therefore
1.175 +//
1.176 +// (xB, xN) = Q' * x, (9)
1.177 +//
1.178 +// where x is the vector of all variables in the original order, xB is
1.179 +// a vector of basic variables, xN is a vector of non-basic variables,
1.180 +// Q' = inv(Q) is a matrix transposed to Q.
1.181 +//
1.182 +// Current values of non-basic variables xN[j], j = 1, ..., n, are not
1.183 +// stored; they are defined implicitly by their statuses as follows:
1.184 +//
1.185 +// 0, if xN[j] is free variable
1.186 +// lN[j], if xN[j] is on its lower bound (10)
1.187 +// uN[j], if xN[j] is on its upper bound
1.188 +// lN[j] = uN[j], if xN[j] is fixed variable
1.189 +//
1.190 +// where lN[j] and uN[j] are lower and upper bounds of xN[j].
1.191 +//
1.192 +// Current values of basic variables xB[i], i = 1, ..., m, are computed
1.193 +// as follows:
1.194 +//
1.195 +// beta = - inv(B) * N * xN, (11)
1.196 +//
1.197 +// where current values of xN are defined by (10).
1.198 +//
1.199 +// Current values of simplex multipliers pi[i], i = 1, ..., m (which
1.200 +// are values of Lagrange multipliers for equality constraints (7) also
1.201 +// called shadow prices) are computed as follows:
1.202 +//
1.203 +// pi = inv(B') * cB, (12)
1.204 +//
1.205 +// where B' is a matrix transposed to B, cB is a vector of objective
1.206 +// coefficients at basic variables xB.
1.207 +//
1.208 +// Current values of reduced costs d[j], j = 1, ..., n, (which are
1.209 +// values of Langrange multipliers for active inequality constraints
1.210 +// corresponding to non-basic variables) are computed as follows:
1.211 +//
1.212 +// d = cN - N' * pi, (13)
1.213 +//
1.214 +// where N' is a matrix transposed to N, cN is a vector of objective
1.215 +// coefficients at non-basic variables xN.
1.216 +----------------------------------------------------------------------*/
1.217 + int *stat; /* int stat[1+m+n]; */
1.218 + /* stat[0] is not used;
1.219 + stat[k], 1 <= k <= m+n, is the status of variable x[k]: */
1.220 +#define SSX_BS 0 /* basic variable */
1.221 +#define SSX_NL 1 /* non-basic variable on lower bound */
1.222 +#define SSX_NU 2 /* non-basic variable on upper bound */
1.223 +#define SSX_NF 3 /* non-basic free variable */
1.224 +#define SSX_NS 4 /* non-basic fixed variable */
1.225 + int *Q_row; /* int Q_row[1+m+n]; */
1.226 + /* matrix Q in row-like format;
1.227 + Q_row[0] is not used;
1.228 + Q_row[i] = j means that q[i,j] = 1 */
1.229 + int *Q_col; /* int Q_col[1+m+n]; */
1.230 + /* matrix Q in column-like format;
1.231 + Q_col[0] is not used;
1.232 + Q_col[j] = i means that q[i,j] = 1 */
1.233 + /* if k-th column of the matrix (I | A) is k'-th column of the
1.234 + matrix (B | N), then Q_row[k] = k' and Q_col[k'] = k;
1.235 + if x[k] is xB[i], then Q_row[k] = i and Q_col[i] = k;
1.236 + if x[k] is xN[j], then Q_row[k] = m+j and Q_col[m+j] = k */
1.237 + BFX *binv;
1.238 + /* invertable form of the basis matrix B */
1.239 + mpq_t *bbar; /* mpq_t bbar[1+m]; alias: beta */
1.240 + /* bbar[0] is a value of the objective function;
1.241 + bbar[i], 1 <= i <= m, is a value of basic variable xB[i] */
1.242 + mpq_t *pi; /* mpq_t pi[1+m]; */
1.243 + /* pi[0] is not used;
1.244 + pi[i], 1 <= i <= m, is a simplex multiplier corresponding to
1.245 + i-th row (equality constraint) */
1.246 + mpq_t *cbar; /* mpq_t cbar[1+n]; alias: d */
1.247 + /* cbar[0] is not used;
1.248 + cbar[j], 1 <= j <= n, is a reduced cost of non-basic variable
1.249 + xN[j] */
1.250 +/*----------------------------------------------------------------------
1.251 +// SIMPLEX TABLE
1.252 +//
1.253 +// Due to (8) and (9) the system of equality constraints (7) for the
1.254 +// current basis can be written as follows:
1.255 +//
1.256 +// xB = A~ * xN, (14)
1.257 +//
1.258 +// where
1.259 +//
1.260 +// A~ = - inv(B) * N (15)
1.261 +//
1.262 +// is a mxn matrix called the simplex table.
1.263 +//
1.264 +// The revised simplex method uses only two components of A~, namely,
1.265 +// pivot column corresponding to non-basic variable xN[q] chosen to
1.266 +// enter the basis, and pivot row corresponding to basic variable xB[p]
1.267 +// chosen to leave the basis.
1.268 +//
1.269 +// Pivot column alfa_q is q-th column of A~, so
1.270 +//
1.271 +// alfa_q = A~ * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q], (16)
1.272 +//
1.273 +// where N[q] is q-th column of the matrix N.
1.274 +//
1.275 +// Pivot row alfa_p is p-th row of A~ or, equivalently, p-th column of
1.276 +// A~', a matrix transposed to A~, so
1.277 +//
1.278 +// alfa_p = A~' * e[p] = - N' * inv(B') * e[p] = - N' * rho_p, (17)
1.279 +//
1.280 +// where (*)' means transposition, and
1.281 +//
1.282 +// rho_p = inv(B') * e[p], (18)
1.283 +//
1.284 +// is p-th column of inv(B') or, that is the same, p-th row of inv(B).
1.285 +----------------------------------------------------------------------*/
1.286 + int p;
1.287 + /* number of basic variable xB[p], 1 <= p <= m, chosen to leave
1.288 + the basis */
1.289 + mpq_t *rho; /* mpq_t rho[1+m]; */
1.290 + /* p-th row of the inverse inv(B); see (18) */
1.291 + mpq_t *ap; /* mpq_t ap[1+n]; */
1.292 + /* p-th row of the simplex table; see (17) */
1.293 + int q;
1.294 + /* number of non-basic variable xN[q], 1 <= q <= n, chosen to
1.295 + enter the basis */
1.296 + mpq_t *aq; /* mpq_t aq[1+m]; */
1.297 + /* q-th column of the simplex table; see (16) */
1.298 +/*--------------------------------------------------------------------*/
1.299 + int q_dir;
1.300 + /* direction in which non-basic variable xN[q] should change on
1.301 + moving to the adjacent vertex of the polyhedron:
1.302 + +1 means that xN[q] increases
1.303 + -1 means that xN[q] decreases */
1.304 + int p_stat;
1.305 + /* non-basic status which should be assigned to basic variable
1.306 + xB[p] when it has left the basis and become xN[q] */
1.307 + mpq_t delta;
1.308 + /* actual change of xN[q] in the adjacent basis (it has the same
1.309 + sign as q_dir) */
1.310 +/*--------------------------------------------------------------------*/
1.311 + int it_lim;
1.312 + /* simplex iterations limit; if this value is positive, it is
1.313 + decreased by one each time when one simplex iteration has been
1.314 + performed, and reaching zero value signals the solver to stop
1.315 + the search; negative value means no iterations limit */
1.316 + int it_cnt;
1.317 + /* simplex iterations count; this count is increased by one each
1.318 + time when one simplex iteration has been performed */
1.319 + double tm_lim;
1.320 + /* searching time limit, in seconds; if this value is positive,
1.321 + it is decreased each time when one simplex iteration has been
1.322 + performed by the amount of time spent for the iteration, and
1.323 + reaching zero value signals the solver to stop the search;
1.324 + negative value means no time limit */
1.325 + double out_frq;
1.326 + /* output frequency, in seconds; this parameter specifies how
1.327 + frequently the solver sends information about the progress of
1.328 + the search to the standard output */
1.329 + glp_long tm_beg;
1.330 + /* starting time of the search, in seconds; the total time of the
1.331 + search is the difference between xtime() and tm_beg */
1.332 + glp_long tm_lag;
1.333 + /* the most recent time, in seconds, at which the progress of the
1.334 + the search was displayed */
1.335 +};
1.336 +
1.337 +#define ssx_create _glp_ssx_create
1.338 +#define ssx_factorize _glp_ssx_factorize
1.339 +#define ssx_get_xNj _glp_ssx_get_xNj
1.340 +#define ssx_eval_bbar _glp_ssx_eval_bbar
1.341 +#define ssx_eval_pi _glp_ssx_eval_pi
1.342 +#define ssx_eval_dj _glp_ssx_eval_dj
1.343 +#define ssx_eval_cbar _glp_ssx_eval_cbar
1.344 +#define ssx_eval_rho _glp_ssx_eval_rho
1.345 +#define ssx_eval_row _glp_ssx_eval_row
1.346 +#define ssx_eval_col _glp_ssx_eval_col
1.347 +#define ssx_chuzc _glp_ssx_chuzc
1.348 +#define ssx_chuzr _glp_ssx_chuzr
1.349 +#define ssx_update_bbar _glp_ssx_update_bbar
1.350 +#define ssx_update_pi _glp_ssx_update_pi
1.351 +#define ssx_update_cbar _glp_ssx_update_cbar
1.352 +#define ssx_change_basis _glp_ssx_change_basis
1.353 +#define ssx_delete _glp_ssx_delete
1.354 +
1.355 +#define ssx_phase_I _glp_ssx_phase_I
1.356 +#define ssx_phase_II _glp_ssx_phase_II
1.357 +#define ssx_driver _glp_ssx_driver
1.358 +
1.359 +SSX *ssx_create(int m, int n, int nnz);
1.360 +/* create simplex solver workspace */
1.361 +
1.362 +int ssx_factorize(SSX *ssx);
1.363 +/* factorize the current basis matrix */
1.364 +
1.365 +void ssx_get_xNj(SSX *ssx, int j, mpq_t x);
1.366 +/* determine value of non-basic variable */
1.367 +
1.368 +void ssx_eval_bbar(SSX *ssx);
1.369 +/* compute values of basic variables */
1.370 +
1.371 +void ssx_eval_pi(SSX *ssx);
1.372 +/* compute values of simplex multipliers */
1.373 +
1.374 +void ssx_eval_dj(SSX *ssx, int j, mpq_t dj);
1.375 +/* compute reduced cost of non-basic variable */
1.376 +
1.377 +void ssx_eval_cbar(SSX *ssx);
1.378 +/* compute reduced costs of all non-basic variables */
1.379 +
1.380 +void ssx_eval_rho(SSX *ssx);
1.381 +/* compute p-th row of the inverse */
1.382 +
1.383 +void ssx_eval_row(SSX *ssx);
1.384 +/* compute pivot row of the simplex table */
1.385 +
1.386 +void ssx_eval_col(SSX *ssx);
1.387 +/* compute pivot column of the simplex table */
1.388 +
1.389 +void ssx_chuzc(SSX *ssx);
1.390 +/* choose pivot column */
1.391 +
1.392 +void ssx_chuzr(SSX *ssx);
1.393 +/* choose pivot row */
1.394 +
1.395 +void ssx_update_bbar(SSX *ssx);
1.396 +/* update values of basic variables */
1.397 +
1.398 +void ssx_update_pi(SSX *ssx);
1.399 +/* update simplex multipliers */
1.400 +
1.401 +void ssx_update_cbar(SSX *ssx);
1.402 +/* update reduced costs of non-basic variables */
1.403 +
1.404 +void ssx_change_basis(SSX *ssx);
1.405 +/* change current basis to adjacent one */
1.406 +
1.407 +void ssx_delete(SSX *ssx);
1.408 +/* delete simplex solver workspace */
1.409 +
1.410 +int ssx_phase_I(SSX *ssx);
1.411 +/* find primal feasible solution */
1.412 +
1.413 +int ssx_phase_II(SSX *ssx);
1.414 +/* find optimal solution */
1.415 +
1.416 +int ssx_driver(SSX *ssx);
1.417 +/* base driver to exact simplex method */
1.418 +
1.419 +#endif
1.420 +
1.421 +/* eof */