COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpssx.h @ 10:5545663ca997

subpack-glpk
Last change on this file since 10:5545663ca997 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import GLPK 4.47

File size: 16.1 KB
Line 
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, 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***********************************************************************/
24
25#ifndef GLPSSX_H
26#define GLPSSX_H
27
28#include "glpbfx.h"
29#include "glpenv.h"
30
31typedef struct SSX SSX;
32
33struct 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
356SSX *ssx_create(int m, int n, int nnz);
357/* create simplex solver workspace */
358
359int ssx_factorize(SSX *ssx);
360/* factorize the current basis matrix */
361
362void ssx_get_xNj(SSX *ssx, int j, mpq_t x);
363/* determine value of non-basic variable */
364
365void ssx_eval_bbar(SSX *ssx);
366/* compute values of basic variables */
367
368void ssx_eval_pi(SSX *ssx);
369/* compute values of simplex multipliers */
370
371void ssx_eval_dj(SSX *ssx, int j, mpq_t dj);
372/* compute reduced cost of non-basic variable */
373
374void ssx_eval_cbar(SSX *ssx);
375/* compute reduced costs of all non-basic variables */
376
377void ssx_eval_rho(SSX *ssx);
378/* compute p-th row of the inverse */
379
380void ssx_eval_row(SSX *ssx);
381/* compute pivot row of the simplex table */
382
383void ssx_eval_col(SSX *ssx);
384/* compute pivot column of the simplex table */
385
386void ssx_chuzc(SSX *ssx);
387/* choose pivot column */
388
389void ssx_chuzr(SSX *ssx);
390/* choose pivot row */
391
392void ssx_update_bbar(SSX *ssx);
393/* update values of basic variables */
394
395void ssx_update_pi(SSX *ssx);
396/* update simplex multipliers */
397
398void ssx_update_cbar(SSX *ssx);
399/* update reduced costs of non-basic variables */
400
401void ssx_change_basis(SSX *ssx);
402/* change current basis to adjacent one */
403
404void ssx_delete(SSX *ssx);
405/* delete simplex solver workspace */
406
407int ssx_phase_I(SSX *ssx);
408/* find primal feasible solution */
409
410int ssx_phase_II(SSX *ssx);
411/* find optimal solution */
412
413int ssx_driver(SSX *ssx);
414/* base driver to exact simplex method */
415
416#endif
417
418/* eof */
Note: See TracBrowser for help on using the repository browser.