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 */ |
---|