alpar@1: /* glpssx.h (simplex method, bignum arithmetic) */ alpar@1: alpar@1: /*********************************************************************** alpar@1: * This code is part of GLPK (GNU Linear Programming Kit). alpar@1: * alpar@1: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@1: * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, alpar@1: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@1: * E-mail: . alpar@1: * alpar@1: * GLPK is free software: you can redistribute it and/or modify it alpar@1: * under the terms of the GNU General Public License as published by alpar@1: * the Free Software Foundation, either version 3 of the License, or alpar@1: * (at your option) any later version. alpar@1: * alpar@1: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@1: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@1: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@1: * License for more details. alpar@1: * alpar@1: * You should have received a copy of the GNU General Public License alpar@1: * along with GLPK. If not, see . alpar@1: ***********************************************************************/ alpar@1: alpar@1: #ifndef GLPSSX_H alpar@1: #define GLPSSX_H alpar@1: alpar@1: #include "glpbfx.h" alpar@1: #include "glpenv.h" alpar@1: alpar@1: typedef struct SSX SSX; alpar@1: alpar@1: struct SSX alpar@1: { /* simplex solver workspace */ alpar@1: /*---------------------------------------------------------------------- alpar@1: // LP PROBLEM DATA alpar@1: // alpar@1: // It is assumed that LP problem has the following statement: alpar@1: // alpar@1: // minimize (or maximize) alpar@1: // alpar@1: // z = c[1]*x[1] + ... + c[m+n]*x[m+n] + c[0] (1) alpar@1: // alpar@1: // subject to equality constraints alpar@1: // alpar@1: // x[1] - a[1,1]*x[m+1] - ... - a[1,n]*x[m+n] = 0 alpar@1: // alpar@1: // . . . . . . . (2) alpar@1: // alpar@1: // x[m] - a[m,1]*x[m+1] + ... - a[m,n]*x[m+n] = 0 alpar@1: // alpar@1: // and bounds of variables alpar@1: // alpar@1: // l[1] <= x[1] <= u[1] alpar@1: // alpar@1: // . . . . . . . (3) alpar@1: // alpar@1: // l[m+n] <= x[m+n] <= u[m+n] alpar@1: // alpar@1: // where: alpar@1: // x[1], ..., x[m] - auxiliary variables; alpar@1: // x[m+1], ..., x[m+n] - structural variables; alpar@1: // z - objective function; alpar@1: // c[1], ..., c[m+n] - coefficients of the objective function; alpar@1: // c[0] - constant term of the objective function; alpar@1: // a[1,1], ..., a[m,n] - constraint coefficients; alpar@1: // l[1], ..., l[m+n] - lower bounds of variables; alpar@1: // u[1], ..., u[m+n] - upper bounds of variables. alpar@1: // alpar@1: // Bounds of variables can be finite as well as inifinite. Besides, alpar@1: // lower and upper bounds can be equal to each other. So the following alpar@1: // five types of variables are possible: alpar@1: // alpar@1: // Bounds of variable Type of variable alpar@1: // ------------------------------------------------- alpar@1: // -inf < x[k] < +inf Free (unbounded) variable alpar@1: // l[k] <= x[k] < +inf Variable with lower bound alpar@1: // -inf < x[k] <= u[k] Variable with upper bound alpar@1: // l[k] <= x[k] <= u[k] Double-bounded variable alpar@1: // l[k] = x[k] = u[k] Fixed variable alpar@1: // alpar@1: // Using vector-matrix notations the LP problem (1)-(3) can be written alpar@1: // as follows: alpar@1: // alpar@1: // minimize (or maximize) alpar@1: // alpar@1: // z = c * x + c[0] (4) alpar@1: // alpar@1: // subject to equality constraints alpar@1: // alpar@1: // xR - A * xS = 0 (5) alpar@1: // alpar@1: // and bounds of variables alpar@1: // alpar@1: // l <= x <= u (6) alpar@1: // alpar@1: // where: alpar@1: // xR - vector of auxiliary variables; alpar@1: // xS - vector of structural variables; alpar@1: // x = (xR, xS) - vector of all variables; alpar@1: // z - objective function; alpar@1: // c - vector of objective coefficients; alpar@1: // c[0] - constant term of the objective function; alpar@1: // A - matrix of constraint coefficients (has m rows alpar@1: // and n columns); alpar@1: // l - vector of lower bounds of variables; alpar@1: // u - vector of upper bounds of variables. alpar@1: // alpar@1: // The simplex method makes no difference between auxiliary and alpar@1: // structural variables, so it is convenient to think the system of alpar@1: // equality constraints (5) written in a homogeneous form: alpar@1: // alpar@1: // (I | -A) * x = 0, (7) alpar@1: // alpar@1: // where (I | -A) is an augmented (m+n)xm constraint matrix, I is mxm alpar@1: // unity matrix whose columns correspond to auxiliary variables, and A alpar@1: // is the original mxn constraint matrix whose columns correspond to alpar@1: // structural variables. Note that only the matrix A is stored. alpar@1: ----------------------------------------------------------------------*/ alpar@1: int m; alpar@1: /* number of rows (auxiliary variables), m > 0 */ alpar@1: int n; alpar@1: /* number of columns (structural variables), n > 0 */ alpar@1: int *type; /* int type[1+m+n]; */ alpar@1: /* type[0] is not used; alpar@1: type[k], 1 <= k <= m+n, is the type of variable x[k]: */ alpar@1: #define SSX_FR 0 /* free (unbounded) variable */ alpar@1: #define SSX_LO 1 /* variable with lower bound */ alpar@1: #define SSX_UP 2 /* variable with upper bound */ alpar@1: #define SSX_DB 3 /* double-bounded variable */ alpar@1: #define SSX_FX 4 /* fixed variable */ alpar@1: mpq_t *lb; /* mpq_t lb[1+m+n]; alias: l */ alpar@1: /* lb[0] is not used; alpar@1: lb[k], 1 <= k <= m+n, is an lower bound of variable x[k]; alpar@1: if x[k] has no lower bound, lb[k] is zero */ alpar@1: mpq_t *ub; /* mpq_t ub[1+m+n]; alias: u */ alpar@1: /* ub[0] is not used; alpar@1: ub[k], 1 <= k <= m+n, is an upper bound of variable x[k]; alpar@1: if x[k] has no upper bound, ub[k] is zero; alpar@1: if x[k] is of fixed type, ub[k] is equal to lb[k] */ alpar@1: int dir; alpar@1: /* optimization direction (sense of the objective function): */ alpar@1: #define SSX_MIN 0 /* minimization */ alpar@1: #define SSX_MAX 1 /* maximization */ alpar@1: mpq_t *coef; /* mpq_t coef[1+m+n]; alias: c */ alpar@1: /* coef[0] is a constant term of the objective function; alpar@1: coef[k], 1 <= k <= m+n, is a coefficient of the objective alpar@1: function at variable x[k]; alpar@1: note that auxiliary variables also may have non-zero objective alpar@1: coefficients */ alpar@1: int *A_ptr; /* int A_ptr[1+n+1]; */ alpar@1: int *A_ind; /* int A_ind[A_ptr[n+1]]; */ alpar@1: mpq_t *A_val; /* mpq_t A_val[A_ptr[n+1]]; */ alpar@1: /* constraint matrix A (see (5)) in storage-by-columns format */ alpar@1: /*---------------------------------------------------------------------- alpar@1: // LP BASIS AND CURRENT BASIC SOLUTION alpar@1: // alpar@1: // The LP basis is defined by the following partition of the augmented alpar@1: // constraint matrix (7): alpar@1: // alpar@1: // (B | N) = (I | -A) * Q, (8) alpar@1: // alpar@1: // where B is a mxm non-singular basis matrix whose columns correspond alpar@1: // to basic variables xB, N is a mxn matrix whose columns correspond to alpar@1: // non-basic variables xN, and Q is a permutation (m+n)x(m+n) matrix. alpar@1: // alpar@1: // From (7) and (8) it follows that alpar@1: // alpar@1: // (I | -A) * x = (I | -A) * Q * Q' * x = (B | N) * (xB, xN), alpar@1: // alpar@1: // therefore alpar@1: // alpar@1: // (xB, xN) = Q' * x, (9) alpar@1: // alpar@1: // where x is the vector of all variables in the original order, xB is alpar@1: // a vector of basic variables, xN is a vector of non-basic variables, alpar@1: // Q' = inv(Q) is a matrix transposed to Q. alpar@1: // alpar@1: // Current values of non-basic variables xN[j], j = 1, ..., n, are not alpar@1: // stored; they are defined implicitly by their statuses as follows: alpar@1: // alpar@1: // 0, if xN[j] is free variable alpar@1: // lN[j], if xN[j] is on its lower bound (10) alpar@1: // uN[j], if xN[j] is on its upper bound alpar@1: // lN[j] = uN[j], if xN[j] is fixed variable alpar@1: // alpar@1: // where lN[j] and uN[j] are lower and upper bounds of xN[j]. alpar@1: // alpar@1: // Current values of basic variables xB[i], i = 1, ..., m, are computed alpar@1: // as follows: alpar@1: // alpar@1: // beta = - inv(B) * N * xN, (11) alpar@1: // alpar@1: // where current values of xN are defined by (10). alpar@1: // alpar@1: // Current values of simplex multipliers pi[i], i = 1, ..., m (which alpar@1: // are values of Lagrange multipliers for equality constraints (7) also alpar@1: // called shadow prices) are computed as follows: alpar@1: // alpar@1: // pi = inv(B') * cB, (12) alpar@1: // alpar@1: // where B' is a matrix transposed to B, cB is a vector of objective alpar@1: // coefficients at basic variables xB. alpar@1: // alpar@1: // Current values of reduced costs d[j], j = 1, ..., n, (which are alpar@1: // values of Langrange multipliers for active inequality constraints alpar@1: // corresponding to non-basic variables) are computed as follows: alpar@1: // alpar@1: // d = cN - N' * pi, (13) alpar@1: // alpar@1: // where N' is a matrix transposed to N, cN is a vector of objective alpar@1: // coefficients at non-basic variables xN. alpar@1: ----------------------------------------------------------------------*/ alpar@1: int *stat; /* int stat[1+m+n]; */ alpar@1: /* stat[0] is not used; alpar@1: stat[k], 1 <= k <= m+n, is the status of variable x[k]: */ alpar@1: #define SSX_BS 0 /* basic variable */ alpar@1: #define SSX_NL 1 /* non-basic variable on lower bound */ alpar@1: #define SSX_NU 2 /* non-basic variable on upper bound */ alpar@1: #define SSX_NF 3 /* non-basic free variable */ alpar@1: #define SSX_NS 4 /* non-basic fixed variable */ alpar@1: int *Q_row; /* int Q_row[1+m+n]; */ alpar@1: /* matrix Q in row-like format; alpar@1: Q_row[0] is not used; alpar@1: Q_row[i] = j means that q[i,j] = 1 */ alpar@1: int *Q_col; /* int Q_col[1+m+n]; */ alpar@1: /* matrix Q in column-like format; alpar@1: Q_col[0] is not used; alpar@1: Q_col[j] = i means that q[i,j] = 1 */ alpar@1: /* if k-th column of the matrix (I | A) is k'-th column of the alpar@1: matrix (B | N), then Q_row[k] = k' and Q_col[k'] = k; alpar@1: if x[k] is xB[i], then Q_row[k] = i and Q_col[i] = k; alpar@1: if x[k] is xN[j], then Q_row[k] = m+j and Q_col[m+j] = k */ alpar@1: BFX *binv; alpar@1: /* invertable form of the basis matrix B */ alpar@1: mpq_t *bbar; /* mpq_t bbar[1+m]; alias: beta */ alpar@1: /* bbar[0] is a value of the objective function; alpar@1: bbar[i], 1 <= i <= m, is a value of basic variable xB[i] */ alpar@1: mpq_t *pi; /* mpq_t pi[1+m]; */ alpar@1: /* pi[0] is not used; alpar@1: pi[i], 1 <= i <= m, is a simplex multiplier corresponding to alpar@1: i-th row (equality constraint) */ alpar@1: mpq_t *cbar; /* mpq_t cbar[1+n]; alias: d */ alpar@1: /* cbar[0] is not used; alpar@1: cbar[j], 1 <= j <= n, is a reduced cost of non-basic variable alpar@1: xN[j] */ alpar@1: /*---------------------------------------------------------------------- alpar@1: // SIMPLEX TABLE alpar@1: // alpar@1: // Due to (8) and (9) the system of equality constraints (7) for the alpar@1: // current basis can be written as follows: alpar@1: // alpar@1: // xB = A~ * xN, (14) alpar@1: // alpar@1: // where alpar@1: // alpar@1: // A~ = - inv(B) * N (15) alpar@1: // alpar@1: // is a mxn matrix called the simplex table. alpar@1: // alpar@1: // The revised simplex method uses only two components of A~, namely, alpar@1: // pivot column corresponding to non-basic variable xN[q] chosen to alpar@1: // enter the basis, and pivot row corresponding to basic variable xB[p] alpar@1: // chosen to leave the basis. alpar@1: // alpar@1: // Pivot column alfa_q is q-th column of A~, so alpar@1: // alpar@1: // alfa_q = A~ * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q], (16) alpar@1: // alpar@1: // where N[q] is q-th column of the matrix N. alpar@1: // alpar@1: // Pivot row alfa_p is p-th row of A~ or, equivalently, p-th column of alpar@1: // A~', a matrix transposed to A~, so alpar@1: // alpar@1: // alfa_p = A~' * e[p] = - N' * inv(B') * e[p] = - N' * rho_p, (17) alpar@1: // alpar@1: // where (*)' means transposition, and alpar@1: // alpar@1: // rho_p = inv(B') * e[p], (18) alpar@1: // alpar@1: // is p-th column of inv(B') or, that is the same, p-th row of inv(B). alpar@1: ----------------------------------------------------------------------*/ alpar@1: int p; alpar@1: /* number of basic variable xB[p], 1 <= p <= m, chosen to leave alpar@1: the basis */ alpar@1: mpq_t *rho; /* mpq_t rho[1+m]; */ alpar@1: /* p-th row of the inverse inv(B); see (18) */ alpar@1: mpq_t *ap; /* mpq_t ap[1+n]; */ alpar@1: /* p-th row of the simplex table; see (17) */ alpar@1: int q; alpar@1: /* number of non-basic variable xN[q], 1 <= q <= n, chosen to alpar@1: enter the basis */ alpar@1: mpq_t *aq; /* mpq_t aq[1+m]; */ alpar@1: /* q-th column of the simplex table; see (16) */ alpar@1: /*--------------------------------------------------------------------*/ alpar@1: int q_dir; alpar@1: /* direction in which non-basic variable xN[q] should change on alpar@1: moving to the adjacent vertex of the polyhedron: alpar@1: +1 means that xN[q] increases alpar@1: -1 means that xN[q] decreases */ alpar@1: int p_stat; alpar@1: /* non-basic status which should be assigned to basic variable alpar@1: xB[p] when it has left the basis and become xN[q] */ alpar@1: mpq_t delta; alpar@1: /* actual change of xN[q] in the adjacent basis (it has the same alpar@1: sign as q_dir) */ alpar@1: /*--------------------------------------------------------------------*/ alpar@1: int it_lim; alpar@1: /* simplex iterations limit; if this value is positive, it is alpar@1: decreased by one each time when one simplex iteration has been alpar@1: performed, and reaching zero value signals the solver to stop alpar@1: the search; negative value means no iterations limit */ alpar@1: int it_cnt; alpar@1: /* simplex iterations count; this count is increased by one each alpar@1: time when one simplex iteration has been performed */ alpar@1: double tm_lim; alpar@1: /* searching time limit, in seconds; if this value is positive, alpar@1: it is decreased each time when one simplex iteration has been alpar@1: performed by the amount of time spent for the iteration, and alpar@1: reaching zero value signals the solver to stop the search; alpar@1: negative value means no time limit */ alpar@1: double out_frq; alpar@1: /* output frequency, in seconds; this parameter specifies how alpar@1: frequently the solver sends information about the progress of alpar@1: the search to the standard output */ alpar@1: glp_long tm_beg; alpar@1: /* starting time of the search, in seconds; the total time of the alpar@1: search is the difference between xtime() and tm_beg */ alpar@1: glp_long tm_lag; alpar@1: /* the most recent time, in seconds, at which the progress of the alpar@1: the search was displayed */ alpar@1: }; alpar@1: alpar@1: #define ssx_create _glp_ssx_create alpar@1: #define ssx_factorize _glp_ssx_factorize alpar@1: #define ssx_get_xNj _glp_ssx_get_xNj alpar@1: #define ssx_eval_bbar _glp_ssx_eval_bbar alpar@1: #define ssx_eval_pi _glp_ssx_eval_pi alpar@1: #define ssx_eval_dj _glp_ssx_eval_dj alpar@1: #define ssx_eval_cbar _glp_ssx_eval_cbar alpar@1: #define ssx_eval_rho _glp_ssx_eval_rho alpar@1: #define ssx_eval_row _glp_ssx_eval_row alpar@1: #define ssx_eval_col _glp_ssx_eval_col alpar@1: #define ssx_chuzc _glp_ssx_chuzc alpar@1: #define ssx_chuzr _glp_ssx_chuzr alpar@1: #define ssx_update_bbar _glp_ssx_update_bbar alpar@1: #define ssx_update_pi _glp_ssx_update_pi alpar@1: #define ssx_update_cbar _glp_ssx_update_cbar alpar@1: #define ssx_change_basis _glp_ssx_change_basis alpar@1: #define ssx_delete _glp_ssx_delete alpar@1: alpar@1: #define ssx_phase_I _glp_ssx_phase_I alpar@1: #define ssx_phase_II _glp_ssx_phase_II alpar@1: #define ssx_driver _glp_ssx_driver alpar@1: alpar@1: SSX *ssx_create(int m, int n, int nnz); alpar@1: /* create simplex solver workspace */ alpar@1: alpar@1: int ssx_factorize(SSX *ssx); alpar@1: /* factorize the current basis matrix */ alpar@1: alpar@1: void ssx_get_xNj(SSX *ssx, int j, mpq_t x); alpar@1: /* determine value of non-basic variable */ alpar@1: alpar@1: void ssx_eval_bbar(SSX *ssx); alpar@1: /* compute values of basic variables */ alpar@1: alpar@1: void ssx_eval_pi(SSX *ssx); alpar@1: /* compute values of simplex multipliers */ alpar@1: alpar@1: void ssx_eval_dj(SSX *ssx, int j, mpq_t dj); alpar@1: /* compute reduced cost of non-basic variable */ alpar@1: alpar@1: void ssx_eval_cbar(SSX *ssx); alpar@1: /* compute reduced costs of all non-basic variables */ alpar@1: alpar@1: void ssx_eval_rho(SSX *ssx); alpar@1: /* compute p-th row of the inverse */ alpar@1: alpar@1: void ssx_eval_row(SSX *ssx); alpar@1: /* compute pivot row of the simplex table */ alpar@1: alpar@1: void ssx_eval_col(SSX *ssx); alpar@1: /* compute pivot column of the simplex table */ alpar@1: alpar@1: void ssx_chuzc(SSX *ssx); alpar@1: /* choose pivot column */ alpar@1: alpar@1: void ssx_chuzr(SSX *ssx); alpar@1: /* choose pivot row */ alpar@1: alpar@1: void ssx_update_bbar(SSX *ssx); alpar@1: /* update values of basic variables */ alpar@1: alpar@1: void ssx_update_pi(SSX *ssx); alpar@1: /* update simplex multipliers */ alpar@1: alpar@1: void ssx_update_cbar(SSX *ssx); alpar@1: /* update reduced costs of non-basic variables */ alpar@1: alpar@1: void ssx_change_basis(SSX *ssx); alpar@1: /* change current basis to adjacent one */ alpar@1: alpar@1: void ssx_delete(SSX *ssx); alpar@1: /* delete simplex solver workspace */ alpar@1: alpar@1: int ssx_phase_I(SSX *ssx); alpar@1: /* find primal feasible solution */ alpar@1: alpar@1: int ssx_phase_II(SSX *ssx); alpar@1: /* find optimal solution */ alpar@1: alpar@1: int ssx_driver(SSX *ssx); alpar@1: /* base driver to exact simplex method */ alpar@1: alpar@1: #endif alpar@1: alpar@1: /* eof */