lemon-project-template-glpk
view deps/glpk/src/glpssx.h @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 21:43:29 +0100 |
| parents | |
| children |
line source
1 /* glpssx.h (simplex method, bignum arithmetic) */
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, 2011 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 ***********************************************************************/
25 #ifndef GLPSSX_H
26 #define GLPSSX_H
28 #include "glpbfx.h"
29 #include "glpenv.h"
31 typedef struct SSX SSX;
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 };
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
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
356 SSX *ssx_create(int m, int n, int nnz);
357 /* create simplex solver workspace */
359 int ssx_factorize(SSX *ssx);
360 /* factorize the current basis matrix */
362 void ssx_get_xNj(SSX *ssx, int j, mpq_t x);
363 /* determine value of non-basic variable */
365 void ssx_eval_bbar(SSX *ssx);
366 /* compute values of basic variables */
368 void ssx_eval_pi(SSX *ssx);
369 /* compute values of simplex multipliers */
371 void ssx_eval_dj(SSX *ssx, int j, mpq_t dj);
372 /* compute reduced cost of non-basic variable */
374 void ssx_eval_cbar(SSX *ssx);
375 /* compute reduced costs of all non-basic variables */
377 void ssx_eval_rho(SSX *ssx);
378 /* compute p-th row of the inverse */
380 void ssx_eval_row(SSX *ssx);
381 /* compute pivot row of the simplex table */
383 void ssx_eval_col(SSX *ssx);
384 /* compute pivot column of the simplex table */
386 void ssx_chuzc(SSX *ssx);
387 /* choose pivot column */
389 void ssx_chuzr(SSX *ssx);
390 /* choose pivot row */
392 void ssx_update_bbar(SSX *ssx);
393 /* update values of basic variables */
395 void ssx_update_pi(SSX *ssx);
396 /* update simplex multipliers */
398 void ssx_update_cbar(SSX *ssx);
399 /* update reduced costs of non-basic variables */
401 void ssx_change_basis(SSX *ssx);
402 /* change current basis to adjacent one */
404 void ssx_delete(SSX *ssx);
405 /* delete simplex solver workspace */
407 int ssx_phase_I(SSX *ssx);
408 /* find primal feasible solution */
410 int ssx_phase_II(SSX *ssx);
411 /* find optimal solution */
413 int ssx_driver(SSX *ssx);
414 /* base driver to exact simplex method */
416 #endif
418 /* eof */