1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/doc/glpk04.tex	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,1411 @@
     1.4 +%* glpk04.tex *%
     1.5 +
     1.6 +\chapter{Advanced API Routines}
     1.7 +
     1.8 +\section{Background}
     1.9 +\label{basbgd}
    1.10 +
    1.11 +Using vector and matrix notations LP problem (1.1)---(1.3) (see Section
    1.12 +\ref{seclp}, page \pageref{seclp}) can be stated as follows:
    1.13 +
    1.14 +\medskip
    1.15 +
    1.16 +\noindent
    1.17 +\hspace{.5in} minimize (or maximize)
    1.18 +$$z=c^Tx_S+c_0\eqno(3.1)$$
    1.19 +\hspace{.5in} subject to linear constraints
    1.20 +$$x_R=Ax_S\eqno(3.2)$$
    1.21 +\hspace{.5in} and bounds of variables
    1.22 +$$
    1.23 +\begin{array}{l@{\ }c@{\ }l@{\ }c@{\ }l}
    1.24 +l_R&\leq&x_R&\leq&u_R\\
    1.25 +l_S&\leq&x_S&\leq&u_S\\
    1.26 +\end{array}\eqno(3.3)
    1.27 +$$
    1.28 +where:
    1.29 +
    1.30 +\noindent
    1.31 +$x_R=(x_1,\dots,x_m)$ is the vector of auxiliary variables;
    1.32 +
    1.33 +\noindent
    1.34 +$x_S=(x_{m+1},\dots,x_{m+n})$ is the vector of structural
    1.35 +variables;
    1.36 +
    1.37 +\noindent
    1.38 +$z$ is the objective function;
    1.39 +
    1.40 +\noindent
    1.41 +$c=(c_1,\dots,c_n)$ is the vector of objective coefficients;
    1.42 +
    1.43 +\noindent
    1.44 +$c_0$ is the constant term (``shift'') of the objective function;
    1.45 +
    1.46 +\noindent
    1.47 +$A=(a_{11},\dots,a_{mn})$ is the constraint matrix;
    1.48 +
    1.49 +\noindent
    1.50 +$l_R=(l_1,\dots,l_m)$ is the vector of lower bounds of auxiliary
    1.51 +variables;
    1.52 +
    1.53 +\noindent
    1.54 +$u_R=(u_1,\dots,u_m)$ is the vector of upper bounds of auxiliary
    1.55 +variables;
    1.56 +
    1.57 +\noindent
    1.58 +$l_S=(l_{m+1},\dots,l_{m+n})$ is the vector of lower bounds of
    1.59 +structural variables;
    1.60 +
    1.61 +\noindent
    1.62 +$u_S=(u_{m+1},\dots,u_{m+n})$ is the vector of upper bounds of
    1.63 +structural variables.
    1.64 +
    1.65 +\medskip
    1.66 +
    1.67 +From the simplex method's standpoint there is no difference between
    1.68 +auxiliary and structural variables. This allows combining all these
    1.69 +variables into one vector that leads to the following problem statement:
    1.70 +
    1.71 +\medskip
    1.72 +
    1.73 +\noindent
    1.74 +\hspace{.5in} minimize (or maximize)
    1.75 +$$z=(0\ |\ c)^Tx+c_0\eqno(3.4)$$
    1.76 +\hspace{.5in} subject to linear constraints
    1.77 +$$(I\ |-\!A)x=0\eqno(3.5)$$
    1.78 +\hspace{.5in} and bounds of variables
    1.79 +$$l\leq x\leq u\eqno(3.6)$$
    1.80 +where:
    1.81 +
    1.82 +\noindent
    1.83 +$x=(x_R\ |\ x_S)$ is the $(m+n)$-vector of (all) variables;
    1.84 +
    1.85 +\noindent
    1.86 +$(0\ |\ c)$ is the $(m+n)$-vector of objective
    1.87 +coefficients;\footnote{Subvector 0 corresponds to objective coefficients
    1.88 +at auxiliary variables.}
    1.89 +
    1.90 +\noindent
    1.91 +$(I\ |-\!A)$ is the {\it augmented} constraint
    1.92 +$m\times(m+n)$-matrix;\footnote{Note that due to auxiliary variables
    1.93 +matrix $(I\ |-\!A)$ contains the unity submatrix and therefore has full
    1.94 +rank. This means, in particular, that the system (3.5) has no linearly
    1.95 +dependent constraints.}
    1.96 +
    1.97 +\noindent
    1.98 +$l=(l_R\ |\ l_S)$ is the $(m+n)$-vector of lower bounds of (all)
    1.99 +variables;
   1.100 +
   1.101 +\noindent
   1.102 +$u=(u_R\ |\ u_S)$ is the $(m+n)$-vector of upper bounds of (all)
   1.103 +variables.
   1.104 +
   1.105 +\medskip
   1.106 +
   1.107 +By definition an {\it LP basic solution} geometrically is a point in
   1.108 +the space of all variables, which is the intersection of planes
   1.109 +corresponding to active constraints\footnote{A constraint is called
   1.110 +{\it active} if in a given point it is satisfied as equality, otherwise
   1.111 +it is called {\it inactive}.}. The space of all variables has the
   1.112 +dimension $m+n$, therefore, to define some basic solution we have to
   1.113 +define $m+n$ active constraints. Note that $m$ constraints (3.5) being
   1.114 +linearly independent equalities are always active, so remaining $n$
   1.115 +active constraints can be chosen only from bound constraints (3.6).
   1.116 +
   1.117 +A variable is called {\it non-basic}, if its (lower or upper) bound is
   1.118 +active, otherwise it is called {\it basic}. Since, as was said above,
   1.119 +exactly $n$ bound constraints must be active, in any basic solution
   1.120 +there are always $n$ non-basic variables and $m$ basic variables.
   1.121 +(Note that a free variable also can be non-basic. Although such
   1.122 +variable has no bounds, we can think it as the difference between two
   1.123 +non-negative variables, which both are non-basic in this case.)
   1.124 +
   1.125 +Now consider how to determine numeric values of all variables for a
   1.126 +given basic solution.
   1.127 +
   1.128 +Let $\Pi$ be an appropriate permutation matrix of the order $(m+n)$.
   1.129 +Then we can write:
   1.130 +$$\left(\begin{array}{@{}c@{}}x_B\\x_N\\\end{array}\right)=
   1.131 +\Pi\left(\begin{array}{@{}c@{}}x_R\\x_S\\\end{array}\right)=\Pi x,
   1.132 +\eqno(3.7)$$
   1.133 +where $x_B$ is the vector of basic variables, $x_N$ is the vector of
   1.134 +non-basic variables, $x=(x_R\ |\ x_S)$ is the vector of all variables
   1.135 +in the original order. In this case the system of linear constraints
   1.136 +(3.5) can be rewritten as follows:
   1.137 +$$(I\ |-\!A)\Pi^T\Pi x=0\ \ \ \Rightarrow\ \ \ (B\ |\ N)
   1.138 +\left(\begin{array}{@{}c@{}}x_B\\x_N\\\end{array}\right)=0,\eqno(3.8)$$
   1.139 +where
   1.140 +$$(B\ |\ N)=(I\ |-\!A)\Pi^T.\eqno(3.9)$$
   1.141 +Matrix $B$ is a square non-singular $m\times m$-matrix, which is
   1.142 +composed from columns of the augmented constraint matrix corresponding
   1.143 +to basic variables. It is called the {\it basis matrix} or simply the
   1.144 +{\it basis}. Matrix $N$ is a rectangular $m\times n$-matrix, which is
   1.145 +composed from columns of the augmented constraint matrix corresponding
   1.146 +to non-basic variables.
   1.147 +
   1.148 +From (3.8) it follows that:
   1.149 +$$Bx_B+Nx_N=0,\eqno(3.10)$$
   1.150 +therefore,
   1.151 +$$x_B=-B^{-1}Nx_N.\eqno(3.11)$$
   1.152 +Thus, the formula (3.11) shows how to determine numeric values of basic
   1.153 +variables $x_B$ assuming that non-basic variables $x_N$ are fixed on
   1.154 +their active bounds.
   1.155 +
   1.156 +The $m\times n$-matrix
   1.157 +$$\Xi=-B^{-1}N,\eqno(3.12)$$
   1.158 +which appears in (3.11), is called the {\it simplex
   1.159 +tableau}.\footnote{This definition corresponds to the GLPK
   1.160 +implementation.} It shows how basic variables depend on non-basic
   1.161 +variables:
   1.162 +$$x_B=\Xi x_N.\eqno(3.13)$$
   1.163 +
   1.164 +The system (3.13) is equivalent to the system (3.5) in the sense that
   1.165 +they both define the same set of points in the space of (primal)
   1.166 +variables, which satisfy to these systems. If, moreover, values of all
   1.167 +basic variables satisfy to their bound constraints (3.3), the
   1.168 +corresponding basic solution is called {\it (primal) feasible},
   1.169 +otherwise {\it (primal) infeasible}. It is understood that any (primal)
   1.170 +feasible basic solution satisfy to all constraints (3.2) and (3.3).
   1.171 +
   1.172 +The LP theory says that if LP has optimal solution, it has (at least
   1.173 +one) basic feasible solution, which corresponds to the optimum. And the
   1.174 +most natural way to determine whether a given basic solution is optimal
   1.175 +or not is to use the Karush---Kuhn---Tucker optimality conditions.
   1.176 +
   1.177 +\def\arraystretch{1.5}
   1.178 +
   1.179 +For the problem statement (3.4)---(3.6) the optimality conditions are
   1.180 +the following:\footnote{These conditions can be appiled to any solution,
   1.181 +not only to a basic solution.}
   1.182 +$$(I\ |-\!A)x=0\eqno(3.14)$$
   1.183 +$$(I\ |-\!A)^T\pi+\lambda_l+\lambda_u=\nabla z=(0\ |\ c)^T\eqno(3.15)$$
   1.184 +$$l\leq x\leq u\eqno(3.16)$$
   1.185 +$$\lambda_l\geq 0,\ \ \lambda_u\leq 0\ \ \mbox{(minimization)}
   1.186 +\eqno(3.17)$$
   1.187 +$$\lambda_l\leq 0,\ \ \lambda_u\geq 0\ \ \mbox{(maximization)}
   1.188 +\eqno(3.18)$$
   1.189 +$$(\lambda_l)_k(x_k-l_k)=0,\ \ (\lambda_u)_k(x_k-u_k)=0,\ \ k=1,2,\dots,
   1.190 +m+n\eqno(3.19)$$
   1.191 +where:
   1.192 +$\pi=(\pi_1,\pi_2,\dots,\pi_m)$ is a $m$-vector of Lagrange
   1.193 +multipliers for equality constraints (3.5);
   1.194 +$\lambda_l=[(\lambda_l)_1,(\lambda_l)_2,\dots,(\lambda_l)_n]$ is a
   1.195 +$n$-vector of Lagrange multipliers for lower bound constraints (3.6);
   1.196 +$\lambda_u=[(\lambda_u)_1,(\lambda_u)_2,\dots,(\lambda_u)_n]$ is a
   1.197 +$n$-vector of Lagrange multipliers for upper bound constraints (3.6).
   1.198 +
   1.199 +Condition (3.14) is the {\it primal} (original) system of equality
   1.200 +constraints (3.5).
   1.201 +
   1.202 +Condition (3.15) is the {\it dual} system of equality constraints.
   1.203 +It requires the gradient of the objective function to be a linear
   1.204 +combination of normals to the planes defined by constraints of the
   1.205 +original problem.
   1.206 +
   1.207 +Condition (3.16) is the primal (original) system of bound constraints
   1.208 +(3.6).
   1.209 +
   1.210 +Condition (3.17) (or (3.18) in case of maximization) is the dual system
   1.211 +of bound constraints.
   1.212 +
   1.213 +Condition (3.19) is the {\it complementary slackness condition}. It
   1.214 +requires, for each original (auxiliary or structural) variable $x_k$,
   1.215 +that either its (lower or upper) bound must be active, or zero bound of
   1.216 +the corresponding Lagrange multiplier ($(\lambda_l)_k$ or
   1.217 +$(\lambda_u)_k$) must be active.
   1.218 +
   1.219 +In GLPK two multipliers $(\lambda_l)_k$ and $(\lambda_u)_k$ for each
   1.220 +primal (original) variable $x_k$, $k=1,2,\dots,m+n$, are combined into
   1.221 +one multiplier:
   1.222 +$$\lambda_k=(\lambda_l)_k+(\lambda_u)_k,\eqno(3.20)$$
   1.223 +which is called a {\it dual variable} for $x_k$. This {\it cannot} lead
   1.224 +to the ambiguity, because both lower and upper bounds of $x_k$ cannot be
   1.225 +active at the same time,\footnote{If $x_k$ is a fixed variable, we can
   1.226 +think it as double-bounded variable $l_k\leq x_k\leq u_k$, where
   1.227 +$l_k=u_k.$} so at least one of $(\lambda_l)_k$ and $(\lambda_u)_k$ must
   1.228 +be equal to zero, and because these multipliers have different signs,
   1.229 +the combined multiplier, which is their sum, uniquely defines each of
   1.230 +them.
   1.231 +
   1.232 +\def\arraystretch{1}
   1.233 +
   1.234 +Using dual variables $\lambda_k$ the dual system of bound constraints
   1.235 +(3.17) and (3.18) can be written in the form of so called {\it ``rule of
   1.236 +signs''} as follows:
   1.237 +
   1.238 +\begin{center}
   1.239 +\begin{tabular}{|@{\,}c@{$\,$}|@{$\,$}c@{$\,$}|@{$\,$}c@{$\,$}|
   1.240 +@{$\,$}c|c@{$\,$}|@{$\,$}c@{$\,$}|@{$\,$}c@{$\,$}|}
   1.241 +\hline
   1.242 +Original bound&\multicolumn{3}{c|}{Minimization}&\multicolumn{3}{c|}
   1.243 +{Maximization}\\
   1.244 +\cline{2-7}
   1.245 +constraint&$(\lambda_l)_k$&$(\lambda_u)_k$&$(\lambda_l)_k+
   1.246 +(\lambda_u)_k$&$(\lambda_l)_k$&$(\lambda_u)_k$&$(\lambda_l)_k+
   1.247 +(\lambda_u)_k$\\
   1.248 +\hline
   1.249 +$-\infty<x_k<+\infty$&$=0$&$=0$&$\lambda_k=0$&$=0$&$=0$&$\lambda_k=0$\\
   1.250 +$x_k\geq l_k$&$\geq 0$&$=0$&$\lambda_k\geq 0$&$\leq 0$&$=0$&$\lambda_k
   1.251 +\leq0$\\
   1.252 +$x_k\leq u_k$&$=0$&$\leq 0$&$\lambda_k\leq 0$&$=0$&$\geq 0$&$\lambda_k
   1.253 +\geq0$\\
   1.254 +$l_k\leq x_k\leq u_k$&$\geq 0$& $\leq 0$& $-\infty\!<\!\lambda_k\!<
   1.255 +\!+\infty$
   1.256 +&$\leq 0$& $\geq 0$& $-\infty\!<\!\lambda_k\!<\!+\infty$\\
   1.257 +$x_k=l_k=u_k$&$\geq 0$& $\leq 0$& $-\infty\!<\!\lambda_k\!<\!+\infty$&
   1.258 +$\leq 0$&
   1.259 +$\geq 0$& $-\infty\!<\!\lambda_k\!<\!+\infty$\\
   1.260 +\hline
   1.261 +\end{tabular}
   1.262 +\end{center}
   1.263 +
   1.264 +May note that each primal variable $x_k$ has its dual counterpart
   1.265 +$\lambda_k$ and vice versa. This allows applying the same partition for
   1.266 +the vector of dual variables as (3.7):
   1.267 +$$\left(\begin{array}{@{}c@{}}\lambda_B\\\lambda_N\\\end{array}\right)=
   1.268 +\Pi\lambda,\eqno(3.21)$$
   1.269 +where $\lambda_B$ is a vector of dual variables for basic variables
   1.270 +$x_B$, $\lambda_N$ is a vector of dual variables for non-basic variables
   1.271 +$x_N$.
   1.272 +
   1.273 +By definition, bounds of basic variables are inactive constraints, so in
   1.274 +any basic solution $\lambda_B=0$. Corresponding values of dual variables
   1.275 +$\lambda_N$ for non-basic variables $x_N$ can be determined in the
   1.276 +following way. From the dual system (3.15) we have:
   1.277 +$$(I\ |-\!A)^T\pi+\lambda=(0\ |\ c)^T,\eqno(3.22)$$
   1.278 +so multiplying both sides of (3.22) by matrix $\Pi$ gives:
   1.279 +$$\Pi(I\ |-\!A)^T\pi+\Pi\lambda=\Pi(0\ |\ c)^T.\eqno(3.23)$$
   1.280 +From (3.9) it follows that
   1.281 +$$\Pi(I\ |-\!A)^T=[(I\ |-\!A)\Pi^T]^T=(B\ |\ N)^T.\eqno(3.24)$$
   1.282 +Further, we can apply the partition (3.7) also to the vector of
   1.283 +objective coefficients (see (3.4)):
   1.284 +$$\left(\begin{array}{@{}c@{}}c_B\\c_N\\\end{array}\right)=
   1.285 +\Pi\left(\begin{array}{@{}c@{}}0\\c\\\end{array}\right),\eqno(3.25)$$
   1.286 +where $c_B$ is a vector of objective coefficients at basic variables,
   1.287 +$c_N$ is a vector of objective coefficients at non-basic variables.
   1.288 +Now, substituting (3.24), (3.21), and (3.25) into (3.23), leads to:
   1.289 +$$(B\ |\ N)^T\pi+(\lambda_B\ |\ \lambda_N)^T=(c_B\ |\ c_N)^T,
   1.290 +\eqno(3.26)$$
   1.291 +and transposing both sides of (3.26) gives the system:
   1.292 +$$\left(\begin{array}{@{}c@{}}B^T\\N^T\\\end{array}\right)\pi+
   1.293 +\left(\begin{array}{@{}c@{}}\lambda_B\\\lambda_N\\\end{array}\right)=
   1.294 +\left(\begin{array}{@{}c@{}}c_B\\c_T\\\end{array}\right),\eqno(3.27)$$
   1.295 +which can be written as follows:
   1.296 +$$\left\{
   1.297 +\begin{array}{@{\ }r@{\ }c@{\ }r@{\ }c@{\ }l@{\ }}
   1.298 +B^T\pi&+&\lambda_B&=&c_B\\
   1.299 +N^T\pi&+&\lambda_N&=&c_N\\
   1.300 +\end{array}
   1.301 +\right.\eqno(3.28)
   1.302 +$$
   1.303 +Lagrange multipliers $\pi=(\pi_i)$ correspond to equality constraints
   1.304 +(3.5) and therefore can have any sign. This allows resolving the first
   1.305 +subsystem of (3.28) as follows:\footnote{$B^{-T}$ means $(B^T)^{-1}=
   1.306 +(B^{-1})^T$.}
   1.307 +$$\pi=B^{-T}(c_B-\lambda_B)=-B^{-T}\lambda_B+B^{-T}c_B,\eqno(3.29)$$
   1.308 +and substitution of $\pi$ from (3.29) into the second subsystem of
   1.309 +(3.28) gives:
   1.310 +$$\lambda_N=-N^T\pi+c_N=N^TB^{-T}\lambda_B+(c_N-N^TB^{-T}c_B).
   1.311 +\eqno(3.30)$$
   1.312 +The latter system can be written in the following final form:
   1.313 +$$\lambda_N=-\Xi^T\lambda_B+d,\eqno(3.31)$$
   1.314 +where $\Xi$ is the simplex tableau (see (3.12)), and
   1.315 +$$d=c_N-N^TB^{-T}c_B=c_N+\Xi^Tc_B\eqno(3.32)$$
   1.316 +is the vector of so called {\it reduced costs} of non-basic variables.
   1.317 +
   1.318 +\pagebreak
   1.319 +
   1.320 +Above it was said that in any basic solution $\lambda_B=0$, so
   1.321 +$\lambda_N=d$ as it follows from (3.31).
   1.322 +
   1.323 +The system (3.31) is equivalent to the system (3.15) in the sense that
   1.324 +they both define the same set of points in the space of dual variables
   1.325 +$\lambda$, which satisfy to these systems. If, moreover, values of all
   1.326 +dual variables $\lambda_N$ (i.e. reduced costs $d$) satisfy to their
   1.327 +bound constraints (i.e. to the ``rule of signs''; see the table above),
   1.328 +the corresponding basic solution is called {\it dual feasible},
   1.329 +otherwise {\it dual infeasible}. It is understood that any dual feasible
   1.330 +solution satisfy to all constraints (3.15) and (3.17) (or (3.18) in case
   1.331 +of maximization).
   1.332 +
   1.333 +It can be easily shown that the complementary slackness condition
   1.334 +(3.19) is always satisfied for {\it any} basic solution. Therefore,
   1.335 +a basic solution\footnote{It is assumed that a complete basic solution
   1.336 +has the form $(x,\lambda)$, i.e. it includes primal as well as dual
   1.337 +variables.} is {\it optimal} if and only if it is primal and dual
   1.338 +feasible, because in this case it satifies to all the optimality
   1.339 +conditions (3.14)---(3.19).
   1.340 +
   1.341 +\def\arraystretch{1.5}
   1.342 +
   1.343 +The meaning of reduced costs $d=(d_j)$ of non-basic variables can be
   1.344 +explained in the following way. From (3.4), (3.7), and (3.25) it follows
   1.345 +that:
   1.346 +$$z=c_B^Tx_B+c_N^Tx_N+c_0.\eqno(3.33)$$
   1.347 +Substituting $x_B$ from (3.11) into (3.33) we can eliminate basic
   1.348 +variables and express the objective only through non-basic variables:
   1.349 +$$
   1.350 +\begin{array}{r@{\ }c@{\ }l}
   1.351 +z&=&c_B^T(-B^{-1}Nx_N)+c_N^Tx_N+c_0=\\
   1.352 +&=&(c_N^T-c_B^TB^{-1}N)x_N+c_0=\\
   1.353 +&=&(c_N-N^TB^{-T}c_B)^Tx_N+c_0=\\
   1.354 +&=&d^Tx_N+c_0.
   1.355 +\end{array}\eqno(3.34)
   1.356 +$$
   1.357 +From (3.34) it is seen that reduced cost $d_j$ shows how the objective
   1.358 +function $z$ depends on non-basic variable $(x_N)_j$ in the neighborhood
   1.359 +of the current basic solution, i.e. while the current basis remains
   1.360 +unchanged.
   1.361 +
   1.362 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   1.363 +
   1.364 +\newpage
   1.365 +
   1.366 +\section{LP basis routines}
   1.367 +\label{lpbasis}
   1.368 +
   1.369 +\subsection{glp\_bf\_exists---check if the basis factorization exists}
   1.370 +
   1.371 +\subsubsection*{Synopsis}
   1.372 +
   1.373 +\begin{verbatim}
   1.374 +int glp_bf_exists(glp_prob *lp);
   1.375 +\end{verbatim}
   1.376 +
   1.377 +\subsubsection*{Returns}
   1.378 +
   1.379 +If the basis factorization for the current basis associated with the
   1.380 +specified problem object exists and therefore is available for
   1.381 +computations, the routine \verb|glp_bf_exists| returns non-zero.
   1.382 +Otherwise the routine returns zero.
   1.383 +
   1.384 +\subsubsection*{Comments}
   1.385 +
   1.386 +Let the problem object have $m$ rows and $n$ columns. In GLPK the
   1.387 +{\it basis matrix} $B$ is a square non-singular matrix of the order $m$,
   1.388 +whose columns correspond to basic (auxiliary and/or structural)
   1.389 +variables. It is defined by the following main
   1.390 +equality:\footnote{For more details see Subsection \ref{basbgd},
   1.391 +page \pageref{basbgd}.}
   1.392 +$$(B\ |\ N)=(I\ |-\!A)\Pi^T,$$
   1.393 +where $I$ is the unity matrix of the order $m$, whose columns correspond
   1.394 +to auxiliary variables; $A$ is the original constraint
   1.395 +$m\times n$-matrix, whose columns correspond to structural variables;
   1.396 +$(I\ |-\!A)$ is the augmented constraint\linebreak
   1.397 +$m\times(m+n)$-matrix, whose columns correspond to all (auxiliary and
   1.398 +structural) variables following in the original order; $\Pi$ is a
   1.399 +permutation matrix of the order $m+n$; and $N$ is a rectangular
   1.400 +$m\times n$-matrix, whose columns correspond to non-basic (auxiliary
   1.401 +and/or structural) variables.
   1.402 +
   1.403 +For various reasons it may be necessary to solve linear systems with
   1.404 +matrix $B$. To provide this possibility the GLPK implementation
   1.405 +maintains an invertable form of $B$ (that is, some representation of
   1.406 +$B^{-1}$) called the {\it basis factorization}, which is an internal
   1.407 +component of the problem object. Typically, the basis factorization is
   1.408 +computed by the simplex solver, which keeps it in the problem object
   1.409 +to be available for other computations.
   1.410 +
   1.411 +Should note that any changes in the problem object, which affects the
   1.412 +basis matrix (e.g. changing the status of a row or column, changing
   1.413 +a basic column of the constraint matrix, removing an active constraint,
   1.414 +etc.), invalidates the basis factorization. So before calling any API
   1.415 +routine, which uses the basis factorization, the application program
   1.416 +must make sure (using the routine \verb|glp_bf_exists|) that the
   1.417 +factorization exists and therefore available for computations.
   1.418 +
   1.419 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   1.420 +
   1.421 +\subsection{glp\_factorize---compute the basis factorization}
   1.422 +
   1.423 +\subsubsection*{Synopsis}
   1.424 +
   1.425 +\begin{verbatim}
   1.426 +int glp_factorize(glp_prob *lp);
   1.427 +\end{verbatim}
   1.428 +
   1.429 +\subsubsection*{Description}
   1.430 +
   1.431 +The routine \verb|glp_factorize| computes the basis factorization for
   1.432 +the current basis associated with the specified problem
   1.433 +object.\footnote{The current basis is defined by the current statuses
   1.434 +of rows (auxiliary variables) and columns (structural variables).}
   1.435 +
   1.436 +The basis factorization is computed from ``scratch'' even if it exists,
   1.437 +so the application program may use the routine \verb|glp_bf_exists|,
   1.438 +and, if the basis factorization already exists, not to call the routine
   1.439 +\verb|glp_factorize| to prevent an extra work.
   1.440 +
   1.441 +The routine \verb|glp_factorize| {\it does not} compute components of
   1.442 +the basic solution (i.e. primal and dual values).
   1.443 +
   1.444 +\subsubsection*{Returns}
   1.445 +
   1.446 +\begin{tabular}{@{}p{25mm}p{97.3mm}@{}}
   1.447 +0 & The basis factorization has been successfully computed.\\
   1.448 +\verb|GLP_EBADB| & The basis matrix is invalid, because the number of
   1.449 +basic (auxiliary and structural) variables is not the same as the number
   1.450 +of rows in the problem object.\\
   1.451 +\verb|GLP_ESING| & The basis matrix is singular within the working
   1.452 +precision.\\
   1.453 +\verb|GLP_ECOND| & The basis matrix is ill-conditioned, i.e. its
   1.454 +condition number is too large.\\
   1.455 +\end{tabular}
   1.456 +
   1.457 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   1.458 +
   1.459 +\newpage
   1.460 +
   1.461 +\subsection{glp\_bf\_updated---check if the basis factorization has\\
   1.462 +been updated}
   1.463 +
   1.464 +\subsubsection*{Synopsis}
   1.465 +
   1.466 +\begin{verbatim}
   1.467 +int glp_bf_updated(glp_prob *lp);
   1.468 +\end{verbatim}
   1.469 +
   1.470 +\subsubsection*{Returns}
   1.471 +
   1.472 +If the basis factorization has been just computed from ``scratch'', the
   1.473 +routine \verb|glp_bf_updated| returns zero. Otherwise, if the
   1.474 +factorization has been updated at least once, the routine returns
   1.475 +non-zero.
   1.476 +
   1.477 +\subsubsection*{Comments}
   1.478 +
   1.479 +{\it Updating} the basis factorization means recomputing it to reflect
   1.480 +changes in the basis matrix. For example, on every iteration of the
   1.481 +simplex method some column of the current basis matrix is replaced by a
   1.482 +new column that gives a new basis matrix corresponding to the adjacent
   1.483 +basis. In this case computing the basis factorization for the adjacent
   1.484 +basis from ``scratch'' (as the routine \verb|glp_factorize| does) would
   1.485 +be too time-consuming.
   1.486 +
   1.487 +On the other hand, since the basis factorization update is a numeric
   1.488 +computational procedure, applying it many times may lead to accumulating
   1.489 +round-off errors. Therefore the basis is periodically refactorized
   1.490 +(reinverted) from ``scratch'' (with the routine \verb|glp_factorize|)
   1.491 +that allows improving its numerical properties.
   1.492 +
   1.493 +The routine \verb|glp_bf_updated| allows determining if the basis
   1.494 +factorization has been updated at least once since it was computed from
   1.495 +``scratch''.
   1.496 +
   1.497 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   1.498 +
   1.499 +\newpage
   1.500 +
   1.501 +\subsection{glp\_get\_bfcp---retrieve basis factorization control
   1.502 +parameters}
   1.503 +
   1.504 +\subsubsection*{Synopsis}
   1.505 +
   1.506 +\begin{verbatim}
   1.507 +void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm);
   1.508 +\end{verbatim}
   1.509 +
   1.510 +\subsubsection*{Description}
   1.511 +
   1.512 +The routine \verb|glp_get_bfcp| retrieves control parameters, which are
   1.513 +used on computing and updating the basis factorization associated with
   1.514 +the specified problem object.
   1.515 +
   1.516 +Current values of the control parameters are stored in a \verb|glp_bfcp|
   1.517 +structure, which the parameter \verb|parm| points to. For a detailed
   1.518 +description of the structure \verb|glp_bfcp| see comments to the routine
   1.519 +\verb|glp_set_bfcp| in the next subsection.
   1.520 +
   1.521 +\subsubsection*{Comments}
   1.522 +
   1.523 +The purpose of the routine \verb|glp_get_bfcp| is two-fold. First, it
   1.524 +allows the application program obtaining current values of control
   1.525 +parameters used by internal GLPK routines, which compute and update the
   1.526 +basis factorization.
   1.527 +
   1.528 +The second purpose of this routine is to provide proper values for all
   1.529 +fields of the structure \verb|glp_bfcp| in the case when the application
   1.530 +program needs to change some control parameters.
   1.531 +
   1.532 +\subsection{glp\_set\_bfcp---change basis factorization control
   1.533 +parameters}
   1.534 +
   1.535 +\subsubsection*{Synopsis}
   1.536 +
   1.537 +\begin{verbatim}
   1.538 +void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm);
   1.539 +\end{verbatim}
   1.540 +
   1.541 +\subsubsection*{Description}
   1.542 +
   1.543 +The routine \verb|glp_set_bfcp| changes control parameters, which are
   1.544 +used by internal GLPK routines on computing and updating the basis
   1.545 +factorization associated with the specified problem object.
   1.546 +
   1.547 +New values of the control parameters should be passed in a structure
   1.548 +\verb|glp_bfcp|, which the parameter \verb|parm| points to. For a
   1.549 +detailed description of the structure \verb|glp_bfcp| see paragraph
   1.550 +``Control parameters'' below.
   1.551 +
   1.552 +The parameter \verb|parm| can be specified as \verb|NULL|, in which case
   1.553 +all control parameters are reset to their default values.
   1.554 +
   1.555 +\subsubsection*{Comments}
   1.556 +
   1.557 +Before changing some control parameters with the routine
   1.558 +\verb|glp_set_bfcp| the application program should retrieve current
   1.559 +values of all control parameters with the routine \verb|glp_get_bfcp|.
   1.560 +This is needed for backward compatibility, because in the future there
   1.561 +may appear new members in the structure \verb|glp_bfcp|.
   1.562 +
   1.563 +Note that new values of control parameters come into effect on a next
   1.564 +computation of the basis factorization, not immediately.
   1.565 +
   1.566 +\subsubsection*{Example}
   1.567 +
   1.568 +\begin{verbatim}
   1.569 +glp_prob *lp;
   1.570 +glp_bfcp parm;
   1.571 +. . .
   1.572 +/* retrieve current values of control parameters */
   1.573 +glp_get_bfcp(lp, &parm);
   1.574 +/* change the threshold pivoting tolerance */
   1.575 +parm.piv_tol = 0.05;
   1.576 +/* set new values of control parameters */
   1.577 +glp_set_bfcp(lp, &parm);
   1.578 +. . .
   1.579 +\end{verbatim}
   1.580 +
   1.581 +\subsubsection*{Control parameters}
   1.582 +
   1.583 +This paragraph describes all basis factorization control parameters
   1.584 +currently used in the package. Symbolic names of control parameters are
   1.585 +names of corresponding members in the structure \verb|glp_bfcp|.
   1.586 +
   1.587 +\def\arraystretch{1}
   1.588 +
   1.589 +\medskip
   1.590 +
   1.591 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.592 +\multicolumn{2}{@{}l}{{\tt int type} (default: {\tt GLP\_BF\_FT})} \\
   1.593 +&Basis factorization type:\\
   1.594 +&\verb|GLP_BF_FT|---$LU$ + Forrest--Tomlin update;\\
   1.595 +&\verb|GLP_BF_BG|---$LU$ + Schur complement + Bartels--Golub update;\\
   1.596 +&\verb|GLP_BF_GR|---$LU$ + Schur complement + Givens rotation update.
   1.597 +\\
   1.598 +&In case of \verb|GLP_BF_FT| the update is applied to matrix $U$, while
   1.599 +in cases of \verb|GLP_BF_BG| and \verb|GLP_BF_GR| the update is applied
   1.600 +to the Schur complement.
   1.601 +\end{tabular}
   1.602 +
   1.603 +\medskip
   1.604 +
   1.605 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.606 +\multicolumn{2}{@{}l}{{\tt int lu\_size} (default: {\tt 0})} \\
   1.607 +&The initial size of the Sparse Vector Area, in non-zeros, used on
   1.608 +computing $LU$-factorization of the basis matrix for the first time.
   1.609 +If this parameter is set to 0, the initial SVA size is determined
   1.610 +automatically.\\
   1.611 +\end{tabular}
   1.612 +
   1.613 +\medskip
   1.614 +
   1.615 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.616 +\multicolumn{2}{@{}l}{{\tt double piv\_tol} (default: {\tt 0.10})} \\
   1.617 +&Threshold pivoting (Markowitz) tolerance, 0 $<$ \verb|piv_tol| $<$ 1,
   1.618 +used on computing $LU$-factorization of the basis matrix. Element
   1.619 +$u_{ij}$ of the active submatrix of factor $U$ fits to be pivot if it
   1.620 +satisfies to the stability criterion
   1.621 +$|u_{ij}| >= {\tt piv\_tol}\cdot\max|u_{i*}|$, i.e. if it is not very
   1.622 +small in the magnitude among other elements in the same row. Decreasing
   1.623 +this parameter may lead to better sparsity at the expense of numerical
   1.624 +accuracy, and vice versa.\\
   1.625 +\end{tabular}
   1.626 +
   1.627 +\medskip
   1.628 +
   1.629 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.630 +\multicolumn{2}{@{}l}{{\tt int piv\_lim} (default: {\tt 4})} \\
   1.631 +&This parameter is used on computing $LU$-factorization of the basis
   1.632 +matrix and specifies how many pivot candidates needs to be considered
   1.633 +on choosing a pivot element, \verb|piv_lim| $\geq$ 1. If \verb|piv_lim|
   1.634 +candidates have been considered, the pivoting routine prematurely
   1.635 +terminates the search with the best candidate found.\\
   1.636 +\end{tabular}
   1.637 +
   1.638 +\medskip
   1.639 +
   1.640 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.641 +\multicolumn{2}{@{}l}{{\tt int suhl} (default: {\tt GLP\_ON})} \\
   1.642 +&This parameter is used on computing $LU$-factorization of the basis
   1.643 +matrix. Being set to {\tt GLP\_ON} it enables applying the following
   1.644 +heuristic proposed by Uwe Suhl: if a column of the active submatrix has
   1.645 +no eligible pivot candidates, it is no more considered until it becomes
   1.646 +a column singleton. In many cases this allows reducing the time needed
   1.647 +for pivot searching. To disable this heuristic the parameter \verb|suhl|
   1.648 +should be set to {\tt GLP\_OFF}.\\
   1.649 +\end{tabular}
   1.650 +
   1.651 +\medskip
   1.652 +
   1.653 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.654 +\multicolumn{2}{@{}l}{{\tt double eps\_tol} (default: {\tt 1e-15})} \\
   1.655 +&Epsilon tolerance, \verb|eps_tol| $\geq$ 0, used on computing
   1.656 +$LU$-factorization of the basis matrix. If an element of the active
   1.657 +submatrix of factor $U$ is less than \verb|eps_tol| in the magnitude,
   1.658 +it is replaced by exact zero.\\
   1.659 +\end{tabular}
   1.660 +
   1.661 +\medskip
   1.662 +
   1.663 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.664 +\multicolumn{2}{@{}l}{{\tt double max\_gro} (default: {\tt 1e+10})} \\
   1.665 +&Maximal growth of elements of factor $U$, \verb|max_gro| $\geq$ 1,
   1.666 +allowable on computing $LU$-factorization of the basis matrix. If on
   1.667 +some elimination step the ratio $u_{big}/b_{max}$ (where $u_{big}$ is
   1.668 +the largest magnitude of elements of factor $U$ appeared in its active
   1.669 +submatrix during all the factorization process, $b_{max}$ is the largest
   1.670 +magnitude of elements of the basis matrix to be factorized), the basis
   1.671 +matrix is considered as ill-conditioned.\\
   1.672 +\end{tabular}
   1.673 +
   1.674 +\medskip
   1.675 +
   1.676 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.677 +\multicolumn{2}{@{}l}{{\tt int nfs\_max} (default: {\tt 100})} \\
   1.678 +&Maximal number of additional row-like factors (entries of the eta
   1.679 +file), \verb|nfs_max| $\geq$ 1, which can be added to $LU$-factorization
   1.680 +of the basis matrix on updating it with the Forrest--Tomlin technique.
   1.681 +This parameter is used only once, before $LU$-factorization is computed
   1.682 +for the first time, to allocate working arrays. As a rule, each update
   1.683 +adds one new factor (however, some updates may need no addition), so
   1.684 +this parameter limits the number of updates between refactorizations.\\
   1.685 +\end{tabular}
   1.686 +
   1.687 +\medskip
   1.688 +
   1.689 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.690 +\multicolumn{2}{@{}l}{{\tt double upd\_tol} (default: {\tt 1e-6})} \\
   1.691 +&Update tolerance, 0 $<$ \verb|upd_tol| $<$ 1, used on updating
   1.692 +$LU$-factorization of the basis matrix with the Forrest--Tomlin
   1.693 +technique. If after updating the magnitude of some diagonal element
   1.694 +$u_{kk}$ of factor $U$ becomes less than
   1.695 +${\tt upd\_tol}\cdot\max(|u_{k*}|, |u_{*k}|)$, the factorization is
   1.696 +considered as inaccurate.\\
   1.697 +\end{tabular}
   1.698 +
   1.699 +\medskip
   1.700 +
   1.701 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.702 +\multicolumn{2}{@{}l}{{\tt int nrs\_max} (default: {\tt 100})} \\
   1.703 +&Maximal number of additional rows and columns, \verb|nrs_max| $\geq$ 1,
   1.704 +which can be added to $LU$-factorization of the basis matrix on updating
   1.705 +it with the Schur complement technique. This parameter is used only
   1.706 +once, before $LU$-factorization is computed for the first time, to
   1.707 +allocate working arrays. As a rule, each update adds one new row and
   1.708 +column (however, some updates may need no addition), so this parameter
   1.709 +limits the number of updates between refactorizations.\\
   1.710 +\end{tabular}
   1.711 +
   1.712 +\medskip
   1.713 +
   1.714 +\noindent\begin{tabular}{@{}p{17pt}@{}p{120.5mm}@{}}
   1.715 +\multicolumn{2}{@{}l}{{\tt int rs\_size} (default: {\tt 0})} \\
   1.716 +&The initial size of the Sparse Vector Area, in non-zeros, used to
   1.717 +store non-zero elements of additional rows and columns introduced on
   1.718 +updating $LU$-factorization of the basis matrix with the Schur
   1.719 +complement technique. If this parameter is set to 0, the initial SVA
   1.720 +size is determined automatically.\\
   1.721 +\end{tabular}
   1.722 +
   1.723 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   1.724 +
   1.725 +\newpage
   1.726 +
   1.727 +\subsection{glp\_get\_bhead---retrieve the basis header information}
   1.728 +
   1.729 +\subsubsection*{Synopsis}
   1.730 +
   1.731 +\begin{verbatim}
   1.732 +int glp_get_bhead(glp_prob *lp, int k);
   1.733 +\end{verbatim}
   1.734 +
   1.735 +\subsubsection*{Description}
   1.736 +
   1.737 +The routine \verb|glp_get_bhead| returns the basis header information
   1.738 +for the current basis associated with the specified problem object.
   1.739 +
   1.740 +\subsubsection*{Returns}
   1.741 +
   1.742 +If basic variable $(x_B)_k$, $1\leq k\leq m$, is $i$-th auxiliary
   1.743 +variable ($1\leq i\leq m$), the routine returns $i$. Otherwise, if
   1.744 +$(x_B)_k$ is $j$-th structural variable ($1\leq j\leq n$), the routine
   1.745 +returns $m+j$. Here $m$ is the number of rows and $n$ is the number of
   1.746 +columns in the problem object.
   1.747 +
   1.748 +\subsubsection*{Comments}
   1.749 +
   1.750 +Sometimes the application program may need to know which original
   1.751 +(auxiliary and structural) variable correspond to a given basic
   1.752 +variable, or, that is the same, which column of the augmented constraint
   1.753 +matrix $(I\ |-\!A)$ correspond to a given column of the basis matrix
   1.754 +$B$.
   1.755 +
   1.756 +\def\arraystretch{1}
   1.757 +
   1.758 +The correspondence is defined as follows:\footnote{For more details see
   1.759 +Subsection \ref{basbgd}, page \pageref{basbgd}.}
   1.760 +$$\left(\begin{array}{@{}c@{}}x_B\\x_N\\\end{array}\right)=
   1.761 +\Pi\left(\begin{array}{@{}c@{}}x_R\\x_S\\\end{array}\right)
   1.762 +\ \ \Leftrightarrow
   1.763 +\ \ \left(\begin{array}{@{}c@{}}x_R\\x_S\\\end{array}\right)=
   1.764 +\Pi^T\left(\begin{array}{@{}c@{}}x_B\\x_N\\\end{array}\right),$$
   1.765 +where $x_B$ is the vector of basic variables, $x_N$ is the vector of
   1.766 +non-basic variables, $x_R$ is the vector of auxiliary variables
   1.767 +following in their original order,\footnote{The original order of
   1.768 +auxiliary and structural variables is defined by the ordinal numbers
   1.769 +of corresponding rows and columns in the problem object.} $x_S$ is the
   1.770 +vector of structural variables following in their original order, $\Pi$
   1.771 +is a permutation matrix (which is a component of the basis
   1.772 +factorization).
   1.773 +
   1.774 +Thus, if $(x_B)_k=(x_R)_i$ is $i$-th auxiliary variable, the routine
   1.775 +returns $i$, and if $(x_B)_k=(x_S)_j$ is $j$-th structural variable,
   1.776 +the routine returns $m+j$, where $m$ is the number of rows in the
   1.777 +problem object.
   1.778 +
   1.779 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   1.780 +
   1.781 +\newpage
   1.782 +
   1.783 +\subsection{glp\_get\_row\_bind---retrieve row index in the basis\\
   1.784 +header}
   1.785 +
   1.786 +\subsubsection*{Synopsis}
   1.787 +
   1.788 +\begin{verbatim}
   1.789 +int glp_get_row_bind(glp_prob *lp, int i);
   1.790 +\end{verbatim}
   1.791 +
   1.792 +\subsubsection*{Returns}
   1.793 +
   1.794 +The routine \verb|glp_get_row_bind| returns the index $k$ of basic
   1.795 +variable $(x_B)_k$, $1\leq k\leq m$, which is $i$-th auxiliary variable
   1.796 +(that is, the auxiliary variable corresponding to $i$-th row),
   1.797 +$1\leq i\leq m$, in the current basis associated with the specified
   1.798 +problem object, where $m$ is the number of rows. However, if $i$-th
   1.799 +auxiliary variable is non-basic, the routine returns zero.
   1.800 +
   1.801 +\subsubsection*{Comments}
   1.802 +
   1.803 +The routine \verb|glp_get_row_bind| is an inverse to the routine
   1.804 +\verb|glp_get_bhead|: if \verb|glp_get_bhead|$(lp,k)$ returns $i$,
   1.805 +\verb|glp_get_row_bind|$(lp,i)$ returns $k$, and vice versa.
   1.806 +
   1.807 +\subsection{glp\_get\_col\_bind---retrieve column index in the basis
   1.808 +header}
   1.809 +
   1.810 +\subsubsection*{Synopsis}
   1.811 +
   1.812 +\begin{verbatim}
   1.813 +int glp_get_col_bind(glp_prob *lp, int j);
   1.814 +\end{verbatim}
   1.815 +
   1.816 +\subsubsection*{Returns}
   1.817 +
   1.818 +The routine \verb|glp_get_col_bind| returns the index $k$ of basic
   1.819 +variable $(x_B)_k$, $1\leq k\leq m$, which is $j$-th structural
   1.820 +variable (that is, the structural variable corresponding to $j$-th
   1.821 +column), $1\leq j\leq n$, in the current basis associated with the
   1.822 +specified problem object, where $m$ is the number of rows, $n$ is the
   1.823 +number of columns. However, if $j$-th structural variable is non-basic,
   1.824 +the routine returns zero.
   1.825 +
   1.826 +\subsubsection*{Comments}
   1.827 +
   1.828 +The routine \verb|glp_get_col_bind| is an inverse to the routine
   1.829 +\verb|glp_get_bhead|: if \verb|glp_get_bhead|$(lp,k)$ returns $m+j$,
   1.830 +\verb|glp_get_col_bind|$(lp,j)$ returns $k$, and vice versa.
   1.831 +
   1.832 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   1.833 +
   1.834 +\newpage
   1.835 +
   1.836 +\subsection{glp\_ftran---perform forward transformation}
   1.837 +
   1.838 +\subsubsection*{Synopsis}
   1.839 +
   1.840 +\begin{verbatim}
   1.841 +void glp_ftran(glp_prob *lp, double x[]);
   1.842 +\end{verbatim}
   1.843 +
   1.844 +\subsubsection*{Description}
   1.845 +
   1.846 +The routine \verb|glp_ftran| performs forward transformation (FTRAN),
   1.847 +i.e. it solves the system $Bx=b$, where $B$ is the basis matrix
   1.848 +associated with the specified problem object, $x$ is the vector of
   1.849 +unknowns to be computed, $b$ is the vector of right-hand sides.
   1.850 +
   1.851 +On entry to the routine elements of the vector $b$ should be stored in
   1.852 +locations \verb|x[1]|, \dots, \verb|x[m]|, where $m$ is the number of
   1.853 +rows. On exit the routine stores elements of the vector $x$ in the same
   1.854 +locations.
   1.855 +
   1.856 +\subsection{glp\_btran---perform backward transformation}
   1.857 +
   1.858 +\subsubsection*{Synopsis}
   1.859 +
   1.860 +\begin{verbatim}
   1.861 +void glp_btran(glp_prob *lp, double x[]);
   1.862 +\end{verbatim}
   1.863 +
   1.864 +\subsubsection*{Description}
   1.865 +
   1.866 +The routine \verb|glp_btran| performs backward transformation (BTRAN),
   1.867 +i.e. it solves the system $B^Tx=b$, where $B^T$ is a matrix transposed
   1.868 +to the basis matrix $B$ associated with the specified problem object,
   1.869 +$x$ is the vector of unknowns to be computed, $b$ is the vector of
   1.870 +right-hand sides.
   1.871 +
   1.872 +On entry to the routine elements of the vector $b$ should be stored in
   1.873 +locations \verb|x[1]|, \dots, \verb|x[m]|, where $m$ is the number of
   1.874 +rows. On exit the routine stores elements of the vector $x$ in the same
   1.875 +locations.
   1.876 +
   1.877 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   1.878 +
   1.879 +\newpage
   1.880 +
   1.881 +\subsection{glp\_warm\_up---``warm up'' LP basis}
   1.882 +
   1.883 +\subsubsection*{Synopsis}
   1.884 +
   1.885 +\begin{verbatim}
   1.886 +int glp_warm_up(glp_prob *P);
   1.887 +\end{verbatim}
   1.888 +
   1.889 +\subsubsection*{Description}
   1.890 +
   1.891 +The routine \verb|glp_warm_up| ``warms up'' the LP basis for the
   1.892 +specified problem object using current statuses assigned to rows and
   1.893 +columns (that is, to auxiliary and structural variables).
   1.894 +
   1.895 +This operation includes computing factorization of the basis matrix
   1.896 +(if it does not exist), computing primal and dual components of basic
   1.897 +solution, and determining the solution status.
   1.898 +
   1.899 +\subsubsection*{Returns}
   1.900 +
   1.901 +\begin{tabular}{@{}p{25mm}p{97.3mm}@{}}
   1.902 +0 & The operation has been successfully performed.\\
   1.903 +\verb|GLP_EBADB| & The basis matrix is invalid, because the number of
   1.904 +basic (auxiliary and structural) variables is not the same as the number
   1.905 +of rows in the problem object.\\
   1.906 +\verb|GLP_ESING| & The basis matrix is singular within the working
   1.907 +precision.\\
   1.908 +\verb|GLP_ECOND| & The basis matrix is ill-conditioned, i.e. its
   1.909 +condition number is too large.\\
   1.910 +\end{tabular}
   1.911 +
   1.912 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   1.913 +
   1.914 +\newpage
   1.915 +
   1.916 +\section{Simplex tableau routines}
   1.917 +
   1.918 +\subsection{glp\_eval\_tab\_row---compute row of the tableau}
   1.919 +
   1.920 +\subsubsection*{Synopsis}
   1.921 +
   1.922 +\begin{verbatim}
   1.923 +int glp_eval_tab_row(glp_prob *lp, int k, int ind[],
   1.924 +      double val[]);
   1.925 +\end{verbatim}
   1.926 +
   1.927 +\subsubsection*{Description}
   1.928 +
   1.929 +The routine \verb|glp_eval_tab_row| computes a row of the current
   1.930 +simplex tableau (see Subsection 3.1.1, formula (3.12)), which (row)
   1.931 +corresponds to some basic variable specified by the parameter $k$ as
   1.932 +follows: if $1\leq k\leq m$, the basic variable is $k$-th auxiliary
   1.933 +variable, and if $m+1\leq k\leq m+n$, the basic variable is $(k-m)$-th
   1.934 +structural variable, where $m$ is the number of rows and $n$ is the
   1.935 +number of columns in the specified problem object. The basis
   1.936 +factorization must exist.
   1.937 +
   1.938 +The computed row shows how the specified basic variable depends on
   1.939 +non-basic variables:
   1.940 +$$x_k=(x_B)_i=\xi_{i1}(x_N)_1+\xi_{i2}(x_N)_2+\dots+\xi_{in}(x_N)_n,$$
   1.941 +where $\xi_{i1}$, $\xi_{i2}$, \dots, $\xi_{in}$ are elements of the
   1.942 +simplex table row, $(x_N)_1$, $(x_N)_2$, \dots, $(x_N)_n$ are non-basic
   1.943 +(auxiliary and structural) variables.
   1.944 +
   1.945 +The routine stores column indices and corresponding numeric values of
   1.946 +non-zero elements of the computed row in unordered sparse format in
   1.947 +locations \verb|ind[1]|, \dots, \verb|ind[len]| and \verb|val[1]|,
   1.948 +\dots, \verb|val[len]|, respectively, where $0\leq{\tt len}\leq n$ is
   1.949 +the number of non-zero elements in the row returned on exit.
   1.950 +
   1.951 +Element indices stored in the array \verb|ind| have the same sense as
   1.952 +index $k$, i.e. indices 1 to $m$ denote auxiliary variables while
   1.953 +indices $m+1$ to $m+n$ denote structural variables (all these variables
   1.954 +are obviously non-basic by definition).
   1.955 +
   1.956 +\subsubsection*{Returns}
   1.957 +
   1.958 +The routine \verb|glp_eval_tab_row| returns \verb|len|, which is the
   1.959 +number of non-zero elements in the simplex table row stored in the
   1.960 +arrays \verb|ind| and \verb|val|.
   1.961 +
   1.962 +\subsubsection*{Comments}
   1.963 +
   1.964 +A row of the simplex table is computed as follows. At first, the
   1.965 +routine checks that the specified variable $x_k$ is basic and uses the
   1.966 +permutation matrix $\Pi$ (3.7) to determine index $i$ of basic variable
   1.967 +$(x_B)_i$, which corresponds to $x_k$.
   1.968 +
   1.969 +The row to be computed is $i$-th row of the matrix $\Xi$ (3.12),
   1.970 +therefore:
   1.971 +$$\xi_i=e_i^T\Xi=-e_i^TB^{-1}N=-(B^{-T}e_i)^TN,$$
   1.972 +where $e_i$ is $i$-th unity vector. So the routine performs BTRAN to
   1.973 +obtain $i$-th row of the inverse $B^{-1}$:
   1.974 +$$\varrho_i=B^{-T}e_i,$$
   1.975 +and then computes elements of the simplex table row as inner products:
   1.976 +$$\xi_{ij}=-\varrho_i^TN_j,\ \ j=1,2,\dots,n,$$
   1.977 +where $N_j$ is $j$-th column of matrix $N$ (3.9), which (column)
   1.978 +corresponds to non-basic variable $(x_N)_j$. The permutation matrix
   1.979 +$\Pi$ is used again to convert indices $j$ of non-basic columns to
   1.980 +original ordinal numbers of auxiliary and structural variables.
   1.981 +
   1.982 +\subsection{glp\_eval\_tab\_col---compute column of the tableau}
   1.983 +
   1.984 +\subsubsection*{Synopsis}
   1.985 +
   1.986 +\begin{verbatim}
   1.987 +int glp_eval_tab_col(glp_prob *lp, int k, int ind[],
   1.988 +      double val[]);
   1.989 +\end{verbatim}
   1.990 +
   1.991 +\subsubsection*{Description}
   1.992 +
   1.993 +The routine \verb|glp_eval_tab_col| computes a column of the current
   1.994 +simplex tableau (see Subsection 3.1.1, formula (3.12)), which (column)
   1.995 +corresponds to some non-basic variable specified by the parameter $k$:
   1.996 +if $1\leq k\leq m$, the non-basic variable is $k$-th auxiliary variable,
   1.997 +and if $m+1\leq k\leq m+n$, the non-basic variable is $(k-m)$-th
   1.998 +structural variable, where $m$ is the number of rows and $n$ is the
   1.999 +number of columns in the specified problem object. The basis
  1.1000 +factorization must exist.
  1.1001 +
  1.1002 +The computed column shows how basic variables depends on the specified
  1.1003 +non-basic variable $x_k=(x_N)_j$:
  1.1004 +$$
  1.1005 +\begin{array}{r@{\ }c@{\ }l@{\ }l}
  1.1006 +(x_B)_1&=&\dots+\xi_{1j}(x_N)_j&+\dots\\
  1.1007 +(x_B)_2&=&\dots+\xi_{2j}(x_N)_j&+\dots\\
  1.1008 +.\ \ .&.&.\ \ .\ \ .\ \ .\ \ .\ \ .\ \ .\\
  1.1009 +(x_B)_m&=&\dots+\xi_{mj}(x_N)_j&+\dots\\
  1.1010 +\end{array}
  1.1011 +$$
  1.1012 +where $\xi_{1j}$, $\xi_{2j}$, \dots, $\xi_{mj}$ are elements of the
  1.1013 +simplex table column, $(x_B)_1$, $(x_B)_2$, \dots, $(x_B)_m$ are basic
  1.1014 +(auxiliary and structural) variables.
  1.1015 +
  1.1016 +The routine stores row indices and corresponding numeric values of
  1.1017 +non-zero elements of the computed column in unordered sparse format in
  1.1018 +locations \verb|ind[1]|, \dots, \verb|ind[len]| and \verb|val[1]|,
  1.1019 +\dots, \verb|val[len]|, respectively, where $0\leq{\tt len}\leq m$ is
  1.1020 +the number of non-zero elements in the column returned on exit.
  1.1021 +
  1.1022 +Element indices stored in the array \verb|ind| have the same sense as
  1.1023 +index $k$, i.e. indices 1 to $m$ denote auxiliary variables while
  1.1024 +indices $m+1$ to $m+n$ denote structural variables (all these variables
  1.1025 +are obviously basic by definition).
  1.1026 +
  1.1027 +\subsubsection*{Returns}
  1.1028 +
  1.1029 +The routine \verb|glp_eval_tab_col| returns \verb|len|, which is the
  1.1030 +number of non-zero elements in the simplex table column stored in the
  1.1031 +arrays \verb|ind| and \verb|val|.
  1.1032 +
  1.1033 +\subsubsection*{Comments}
  1.1034 +
  1.1035 +A column of the simplex table is computed as follows. At first, the
  1.1036 +routine checks that the specified variable $x_k$ is non-basic and uses
  1.1037 +the permutation matrix $\Pi$ (3.7) to determine index $j$ of non-basic
  1.1038 +variable $(x_N)_j$, which corresponds to $x_k$.
  1.1039 +
  1.1040 +The column to be computed is $j$-th column of the matrix $\Xi$ (3.12),
  1.1041 +therefore:
  1.1042 +$$\Xi_j=\Xi e_j=-B^{-1}Ne_j=-B^{-1}N_j,$$
  1.1043 +where $e_j$ is $j$-th unity vector, $N_j$ is $j$-th column of matrix
  1.1044 +$N$ (3.9). So the routine performs FTRAN to transform $N_j$ to the
  1.1045 +simplex table column $\Xi_j=(\xi_{ij})$ and uses the permutation matrix
  1.1046 +$\Pi$ to convert row indices $i$ to original ordinal numbers of
  1.1047 +auxiliary and structural variables.
  1.1048 +
  1.1049 +\newpage
  1.1050 +
  1.1051 +\subsection{glp\_transform\_row---transform explicitly specified\\
  1.1052 +row}
  1.1053 +
  1.1054 +\subsubsection*{Synopsis}
  1.1055 +
  1.1056 +\begin{verbatim}
  1.1057 +int glp_transform_row(glp_prob *P, int len, int ind[],
  1.1058 +      double val[]);
  1.1059 +\end{verbatim}
  1.1060 +
  1.1061 +\subsubsection*{Description}
  1.1062 +
  1.1063 +The routine \verb|glp_transform_row| performs the same operation as the
  1.1064 +routine \verb|glp_eval_tab_row| with exception that the row to be
  1.1065 +transformed is specified explicitly as a sparse vector.
  1.1066 +
  1.1067 +The explicitly specified row may be thought as a linear form:
  1.1068 +$$x=a_1x_{m+1}+a_2x_{m+2}+\dots+a_nx_{m+n},$$
  1.1069 +where $x$ is an auxiliary variable for this row, $a_j$ are coefficients
  1.1070 +of the linear form, $x_{m+j}$ are structural variables.
  1.1071 +
  1.1072 +On entry column indices and numerical values of non-zero coefficients
  1.1073 +$a_j$ of the specified row should be placed in locations \verb|ind[1]|,
  1.1074 +\dots, \verb|ind[len]| and \verb|val[1]|, \dots, \verb|val[len]|, where
  1.1075 +\verb|len| is number of non-zero coefficients.
  1.1076 +
  1.1077 +This routine uses the system of equality constraints and the current
  1.1078 +basis in order to express the auxiliary variable $x$ through the current
  1.1079 +non-basic variables (as if the transformed row were added to the problem
  1.1080 +object and the auxiliary variable $x$ were basic), i.e. the resultant
  1.1081 +row has the form:
  1.1082 +$$x=\xi_1(x_N)_1+\xi_2(x_N)_2+\dots+\xi_n(x_N)_n,$$
  1.1083 +where $\xi_j$ are influence coefficients, $(x_N)_j$ are non-basic
  1.1084 +(auxiliary and structural) variables, $n$ is the number of columns in
  1.1085 +the problem object.
  1.1086 +
  1.1087 +On exit the routine stores indices and numerical values of non-zero
  1.1088 +coefficients $\xi_j$ of the resultant row in locations \verb|ind[1]|,
  1.1089 +\dots, \verb|ind[len']| and \verb|val[1]|, \dots, \verb|val[len']|,
  1.1090 +where $0\leq{\tt len'}\leq n$ is the number of non-zero coefficients in
  1.1091 +the resultant row returned by the routine. Note that indices of
  1.1092 +non-basic variables stored in the array \verb|ind| correspond to
  1.1093 +original ordinal numbers of variables: indices 1 to $m$ mean auxiliary
  1.1094 +variables and indices $m+1$ to $m+n$ mean structural ones.
  1.1095 +
  1.1096 +\subsubsection*{Returns}
  1.1097 +
  1.1098 +The routine \verb|glp_transform_row| returns \verb|len'|, the number of
  1.1099 +non-zero coefficients in the resultant row stored in the arrays
  1.1100 +\verb|ind| and \verb|val|.
  1.1101 +
  1.1102 +\subsection{glp\_transform\_col---transform explicitly specified\\
  1.1103 +column}
  1.1104 +
  1.1105 +\subsubsection*{Synopsis}
  1.1106 +
  1.1107 +\begin{verbatim}
  1.1108 +int glp_transform_col(glp_prob *P, int len, int ind[],
  1.1109 +      double val[]);
  1.1110 +\end{verbatim}
  1.1111 +
  1.1112 +\subsubsection*{Description}
  1.1113 +
  1.1114 +The routine \verb|glp_transform_col| performs the same operation as the
  1.1115 +routine \verb|glp_eval_tab_col| with exception that the column to be
  1.1116 +transformed is specified explicitly as a sparse vector.
  1.1117 +
  1.1118 +The explicitly specified column may be thought as it were added to
  1.1119 +the original system of equality constraints:
  1.1120 +$$
  1.1121 +\begin{array}{l@{\ }c@{\ }r@{\ }c@{\ }r@{\ }c@{\ }r}
  1.1122 +x_1&=&a_{11}x_{m+1}&+\dots+&a_{1n}x_{m+n}&+&a_1x \\
  1.1123 +x_2&=&a_{21}x_{m+1}&+\dots+&a_{2n}x_{m+n}&+&a_2x \\
  1.1124 +\multicolumn{7}{c}
  1.1125 +{.\ \ .\ \ .\ \ .\ \ .\ \ .\ \ .\ \ .\ \ .\ \ .\ \ .\ \ .\ \ .\ \ .}\\
  1.1126 +x_m&=&a_{m1}x_{m+1}&+\dots+&a_{mn}x_{m+n}&+&a_mx \\
  1.1127 +\end{array}
  1.1128 +$$
  1.1129 +where $x_i$ are auxiliary variables, $x_{m+j}$ are structural variables
  1.1130 +(presented in the problem object), $x$ is a structural variable for the
  1.1131 +explicitly specified column, $a_i$ are constraint coefficients at $x$.
  1.1132 +
  1.1133 +On entry row indices and numerical values of non-zero coefficients
  1.1134 +$a_i$ of the specified column should be placed in locations
  1.1135 +\verb|ind[1]|, \dots, \verb|ind[len]| and \verb|val[1]|, \dots,
  1.1136 +\verb|val[len]|, where \verb|len| is number of non-zero coefficients.
  1.1137 +
  1.1138 +This routine uses the system of equality constraints and the current
  1.1139 +basis in order to express the current basic variables through the
  1.1140 +structural variable $x$ (as if the transformed column were added to the
  1.1141 +problem object and the variable $x$ were non-basic):
  1.1142 +$$
  1.1143 +\begin{array}{l@{\ }c@{\ }r}
  1.1144 +(x_B)_1&=\dots+&\xi_{1}x\\
  1.1145 +(x_B)_2&=\dots+&\xi_{2}x\\
  1.1146 +\multicolumn{3}{c}{.\ \ .\ \ .\ \ .\ \ .\ \ .}\\
  1.1147 +(x_B)_m&=\dots+&\xi_{m}x\\
  1.1148 +\end{array}
  1.1149 +$$
  1.1150 +where $\xi_i$ are influence coefficients, $x_B$ are basic (auxiliary
  1.1151 +and structural) variables, $m$ is the number of rows in the problem
  1.1152 +object.
  1.1153 +
  1.1154 +On exit the routine stores indices and numerical values of non-zero
  1.1155 +coefficients $\xi_i$ of the resultant column in locations \verb|ind[1]|,
  1.1156 +\dots, \verb|ind[len']| and \verb|val[1]|, \dots, \verb|val[len']|,
  1.1157 +where $0\leq{\tt len'}\leq m$ is the number of non-zero coefficients in
  1.1158 +the resultant column returned by the routine. Note that indices of basic
  1.1159 +variables stored in the array \verb|ind| correspond to original ordinal
  1.1160 +numbers of variables, i.e. indices 1 to $m$ mean auxiliary variables,
  1.1161 +indices $m+1$ to $m+n$ mean structural ones.
  1.1162 +
  1.1163 +\subsubsection*{Returns}
  1.1164 +
  1.1165 +The routine \verb|glp_transform_col| returns \verb|len'|, the number of
  1.1166 +non-zero coefficients in the resultant column stored in the arrays
  1.1167 +\verb|ind| and \verb|val|.
  1.1168 +
  1.1169 +\subsection{glp\_prim\_rtest---perform primal ratio test}
  1.1170 +
  1.1171 +\subsubsection*{Synopsis}
  1.1172 +
  1.1173 +\begin{verbatim}
  1.1174 +int glp_prim_rtest(glp_prob *P, int len, const int ind[],
  1.1175 +      const double val[], int dir, double eps);
  1.1176 +\end{verbatim}
  1.1177 +
  1.1178 +\subsubsection*{Description}
  1.1179 +
  1.1180 +The routine \verb|glp_prim_rtest| performs the primal ratio test using
  1.1181 +an explicitly specified column of the simplex table.
  1.1182 +
  1.1183 +The current basic solution associated with the LP problem object must be
  1.1184 +primal feasible.
  1.1185 +
  1.1186 +The explicitly specified column of the simplex table shows how the basic
  1.1187 +variables $x_B$ depend on some non-basic variable $x$ (which is not
  1.1188 +necessarily presented in the problem object):
  1.1189 +$$
  1.1190 +\begin{array}{l@{\ }c@{\ }r}
  1.1191 +(x_B)_1&=\dots+&\xi_{1}x\\
  1.1192 +(x_B)_2&=\dots+&\xi_{2}x\\
  1.1193 +\multicolumn{3}{c}{.\ \ .\ \ .\ \ .\ \ .\ \ .}\\
  1.1194 +(x_B)_m&=\dots+&\xi_{m}x\\
  1.1195 +\end{array}
  1.1196 +$$
  1.1197 +
  1.1198 +The column is specifed on entry to the routine in sparse format. Ordinal
  1.1199 +numbers of basic variables $(x_B)_i$ should be placed in locations
  1.1200 +\verb|ind[1]|, \dots, \verb|ind[len]|, where ordinal number 1 to $m$
  1.1201 +denote auxiliary variables, and ordinal numbers $m+1$ to $m+n$ denote
  1.1202 +structural variables. The corresponding non-zero coefficients $\xi_i$
  1.1203 +should be placed in locations \verb|val[1]|, \dots, \verb|val[len]|. The
  1.1204 +arrays \verb|ind| and \verb|val| are not changed by the routine.
  1.1205 +
  1.1206 +The parameter \verb|dir| specifies direction in which the variable $x$
  1.1207 +changes on entering the basis: $+1$ means increasing, $-1$ means
  1.1208 +decreasing.
  1.1209 +
  1.1210 +The parameter \verb|eps| is an absolute tolerance (small positive
  1.1211 +number, say, $10^{-9}$) used by the routine to skip $\xi_i$'s whose
  1.1212 +magnitude is less than \verb|eps|.
  1.1213 +
  1.1214 +The routine determines which basic variable (among those specified in
  1.1215 +\verb|ind[1]|, \dots, \verb|ind[len]|) reaches its (lower or upper)
  1.1216 +bound first before any other basic variables do, and which, therefore,
  1.1217 +should leave the basis in order to keep primal feasibility.
  1.1218 +
  1.1219 +\subsubsection*{Returns}
  1.1220 +
  1.1221 +The routine \verb|glp_prim_rtest| returns the index, \verb|piv|, in the
  1.1222 +arrays \verb|ind| and \verb|val| corresponding to the pivot element
  1.1223 +chosen, $1\leq$ \verb|piv| $\leq$ \verb|len|. If the adjacent basic
  1.1224 +solution is primal unbounded, and therefore the choice cannot be made,
  1.1225 +the routine returns zero.
  1.1226 +
  1.1227 +\subsubsection*{Comments}
  1.1228 +
  1.1229 +If the non-basic variable $x$ is presented in the LP problem object, the
  1.1230 +input column can be computed with the routine \verb|glp_eval_tab_col|;
  1.1231 +otherwise, it can be computed with the routine \verb|glp_transform_col|.
  1.1232 +
  1.1233 +\subsection{glp\_dual\_rtest---perform dual ratio test}
  1.1234 +
  1.1235 +\subsubsection*{Synopsis}
  1.1236 +
  1.1237 +\begin{verbatim}
  1.1238 +int glp_dual_rtest(glp_prob *P, int len, const int ind[],
  1.1239 +      const double val[], int dir, double eps);
  1.1240 +\end{verbatim}
  1.1241 +
  1.1242 +\subsubsection*{Description}
  1.1243 +
  1.1244 +The routine \verb|glp_dual_rtest| performs the dual ratio test using
  1.1245 +an explicitly specified row of the simplex table.
  1.1246 +
  1.1247 +The current basic solution associated with the LP problem object must be
  1.1248 +dual feasible.
  1.1249 +
  1.1250 +The explicitly specified row of the simplex table is a linear form
  1.1251 +that shows how some basic variable $x$ (which is not necessarily
  1.1252 +presented in the problem object) depends on non-basic variables $x_N$:
  1.1253 +$$x=\xi_1(x_N)_1+\xi_2(x_N)_2+\dots+\xi_n(x_N)_n.$$
  1.1254 +
  1.1255 +The row is specified on entry to the routine in sparse format. Ordinal
  1.1256 +numbers of non-basic variables $(x_N)_j$ should be placed in locations
  1.1257 +\verb|ind[1]|, \dots, \verb|ind[len]|, where ordinal numbers 1 to $m$
  1.1258 +denote auxiliary variables, and ordinal numbers $m+1$ to $m+n$ denote
  1.1259 +structural variables. The corresponding non-zero coefficients $\xi_j$
  1.1260 +should be placed in locations \verb|val[1]|, \dots, \verb|val[len]|.
  1.1261 +The arrays \verb|ind| and \verb|val| are not changed by the routine.
  1.1262 +
  1.1263 +The parameter \verb|dir| specifies direction in which the variable $x$
  1.1264 +changes on leaving the basis: $+1$ means that $x$ goes on its lower
  1.1265 +bound, so its reduced cost (dual variable) is increasing (minimization)
  1.1266 +or decreasing (maximization); $-1$ means that $x$ goes on its upper
  1.1267 +bound, so its reduced cost is decreasing (minimization) or increasing
  1.1268 +(maximization).
  1.1269 +
  1.1270 +The parameter \verb|eps| is an absolute tolerance (small positive
  1.1271 +number, say, $10^{-9}$) used by the routine to skip $\xi_j$'s whose
  1.1272 +magnitude is less than \verb|eps|.
  1.1273 +
  1.1274 +The routine determines which non-basic variable (among those specified
  1.1275 +in \verb|ind[1]|, \dots, \verb|ind[len]|) should enter the basis in
  1.1276 +order to keep dual feasibility, because its reduced cost reaches the
  1.1277 +(zero) bound first before this occurs for any other non-basic variables.
  1.1278 +
  1.1279 +\subsubsection*{Returns}
  1.1280 +
  1.1281 +The routine \verb|glp_dual_rtest| returns the index, \verb|piv|, in the
  1.1282 +arrays \verb|ind| and \verb|val| corresponding to the pivot element
  1.1283 +chosen, $1\leq$ \verb|piv| $\leq$ \verb|len|. If the adjacent basic
  1.1284 +solution is dual unbounded, and therefore the choice cannot be made,
  1.1285 +the routine returns zero.
  1.1286 +
  1.1287 +\subsubsection*{Comments}
  1.1288 +
  1.1289 +If the basic variable $x$ is presented in the LP problem object, the
  1.1290 +input row can be computed with the routine \verb|glp_eval_tab_row|;
  1.1291 +otherwise, it can be computed with the routine \verb|glp_transform_row|.
  1.1292 +
  1.1293 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  1.1294 +
  1.1295 +\newpage
  1.1296 +
  1.1297 +\section{Post-optimal analysis routines}
  1.1298 +
  1.1299 +\subsection{glp\_analyze\_bound---analyze active bound of non-basic
  1.1300 +variable}
  1.1301 +
  1.1302 +\subsubsection*{Synopsis}
  1.1303 +
  1.1304 +\begin{verbatim}
  1.1305 +void glp_analyze_bound(glp_prob *P, int k, double *limit1,
  1.1306 +      int *var1, double *limit2, int *var2);
  1.1307 +\end{verbatim}
  1.1308 +
  1.1309 +\subsubsection*{Description}
  1.1310 +
  1.1311 +The routine \verb|glp_analyze_bound| analyzes the effect of varying the
  1.1312 +active bound of specified non-basic variable.
  1.1313 +
  1.1314 +The non-basic variable is specified by the parameter $k$, where
  1.1315 +$1\leq k\leq m$ means auxiliary variable of corresponding row, and
  1.1316 +$m+1\leq k\leq m+n$ means structural variable (column).
  1.1317 +
  1.1318 +Note that the current basic solution must be optimal, and the basis
  1.1319 +factorization must exist.
  1.1320 +
  1.1321 +Results of the analysis have the following meaning.
  1.1322 +
  1.1323 +\verb|value1| is the minimal value of the active bound, at which the
  1.1324 +basis still remains primal feasible and thus optimal. \verb|-DBL_MAX|
  1.1325 +means that the active bound has no lower limit.
  1.1326 +
  1.1327 +\verb|var1| is the ordinal number of an auxiliary (1 to $m$) or
  1.1328 +structural ($m+1$ to $m+n$) basic variable, which reaches its bound
  1.1329 +first and thereby limits further decreasing the active bound being
  1.1330 +analyzed. if \verb|value1| = \verb|-DBL_MAX|, \verb|var1| is set to 0.
  1.1331 +
  1.1332 +\verb|value2| is the maximal value of the active bound, at which the
  1.1333 +basis still remains primal feasible and thus optimal. \verb|+DBL_MAX|
  1.1334 +means that the active bound has no upper limit.
  1.1335 +
  1.1336 +\verb|var2| is the ordinal number of an auxiliary (1 to $m$) or
  1.1337 +structural ($m+1$ to $m+n$) basic variable, which reaches its bound
  1.1338 +first and thereby limits further increasing the active bound being
  1.1339 +analyzed. if \verb|value2| = \verb|+DBL_MAX|, \verb|var2| is set to 0.
  1.1340 +
  1.1341 +The parameters \verb|value1|, \verb|var1|, \verb|value2|, \verb|var2|
  1.1342 +can be specified as \verb|NULL|, in which case corresponding information
  1.1343 +is not stored.
  1.1344 +
  1.1345 +\newpage
  1.1346 +
  1.1347 +\subsection{glp\_analyze\_coef---analyze objective coefficient at basic
  1.1348 +variable}
  1.1349 +
  1.1350 +\subsubsection*{Synopsis}
  1.1351 +
  1.1352 +\begin{verbatim}
  1.1353 +void glp_analyze_coef(glp_prob *P, int k, double *coef1,
  1.1354 +      int *var1, double *value1, double *coef2, int *var2,
  1.1355 +      double *value2);
  1.1356 +\end{verbatim}
  1.1357 +
  1.1358 +\subsubsection*{Description}
  1.1359 +
  1.1360 +The routine \verb|glp_analyze_coef| analyzes the effect of varying the
  1.1361 +objective coefficient at specified basic variable.
  1.1362 +
  1.1363 +The basic variable is specified by the parameter $k$, where
  1.1364 +$1\leq k\leq m$ means auxiliary variable of corresponding row, and
  1.1365 +$m+1\leq k\leq m+n$ means structural variable (column).
  1.1366 +
  1.1367 +Note that the current basic solution must be optimal, and the basis
  1.1368 +factorization must exist.
  1.1369 +
  1.1370 +Results of the analysis have the following meaning.
  1.1371 +
  1.1372 +\verb|coef1| is the minimal value of the objective coefficient, at
  1.1373 +which the basis still remains dual feasible and thus optimal.
  1.1374 +\verb|-DBL_MAX| means that the objective coefficient has no lower limit.
  1.1375 +
  1.1376 +\verb|var1| is the ordinal number of an auxiliary (1 to $m$) or
  1.1377 +structural ($m+1$ to $m+n$) non-basic variable, whose reduced cost
  1.1378 +reaches its zero bound first and thereby limits further decreasing the
  1.1379 +objective coefficient being analyzed. If \verb|coef1| = \verb|-DBL_MAX|,
  1.1380 +\verb|var1| is set to 0.
  1.1381 +
  1.1382 +\verb|value1| is value of the basic variable being analyzed in an
  1.1383 +adjacent basis, which is defined as follows. Let the objective
  1.1384 +coefficient reaches its minimal value (\verb|coef1|) and continues
  1.1385 +decreasing. Then the reduced cost of the limiting non-basic variable
  1.1386 +(\verb|var1|) becomes dual infeasible and the current basis becomes
  1.1387 +non-optimal that forces the limiting non-basic variable to enter the
  1.1388 +basis replacing there some basic variable that leaves the basis to keep
  1.1389 +primal feasibility. Should note that on determining the adjacent basis
  1.1390 +current bounds of the basic variable being analyzed are ignored as if
  1.1391 +it were free (unbounded) variable, so it cannot leave the basis. It may
  1.1392 +happen that no dual feasible adjacent basis exists, in which case
  1.1393 +\verb|value1| is set to \verb|-DBL_MAX| or \verb|+DBL_MAX|.
  1.1394 +
  1.1395 +\verb|coef2| is the maximal value of the objective coefficient, at
  1.1396 +which the basis still remains dual feasible and thus optimal.
  1.1397 +\verb|+DBL_MAX| means that the objective coefficient has no upper limit.
  1.1398 +
  1.1399 +\verb|var2| is the ordinal number of an auxiliary (1 to $m$) or
  1.1400 +structural ($m+1$ to $m+n$) non-basic variable, whose reduced cost
  1.1401 +reaches its zero bound first and thereby limits further increasing the
  1.1402 +objective coefficient being analyzed. If \verb|coef2| = \verb|+DBL_MAX|,
  1.1403 +\verb|var2| is set to 0.
  1.1404 +
  1.1405 +\verb|value2| is value of the basic variable being analyzed in an
  1.1406 +adjacent basis, which is defined exactly in the same way as
  1.1407 +\verb|value1| above with exception that now the objective coefficient
  1.1408 +is increasing.
  1.1409 +
  1.1410 +The parameters \verb|coef1|, \verb|var1|, \verb|value1|, \verb|coef2|,
  1.1411 +\verb|var2|, \verb|value2| can be specified as \verb|NULL|, in which
  1.1412 +case corresponding information is not stored.
  1.1413 +
  1.1414 +%* eof *%