1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpini01.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,577 @@
1.4 +/* glpini01.c */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpapi.h"
1.29 +
1.30 +/*----------------------------------------------------------------------
1.31 +-- triang - find maximal triangular part of a rectangular matrix.
1.32 +--
1.33 +-- *Synopsis*
1.34 +--
1.35 +-- int triang(int m, int n,
1.36 +-- void *info, int (*mat)(void *info, int k, int ndx[]),
1.37 +-- int rn[], int cn[]);
1.38 +--
1.39 +-- *Description*
1.40 +--
1.41 +-- For a given rectangular (sparse) matrix A with m rows and n columns
1.42 +-- the routine triang tries to find such permutation matrices P and Q
1.43 +-- that the first rows and columns of the matrix B = P*A*Q form a lower
1.44 +-- triangular submatrix of as greatest size as possible:
1.45 +--
1.46 +-- 1 n
1.47 +-- 1 * . . . . . . x x x x x x
1.48 +-- * * . . . . . x x x x x x
1.49 +-- * * * . . . . x x x x x x
1.50 +-- * * * * . . . x x x x x x
1.51 +-- B = P*A*Q = * * * * * . . x x x x x x
1.52 +-- * * * * * * . x x x x x x
1.53 +-- * * * * * * * x x x x x x
1.54 +-- x x x x x x x x x x x x x
1.55 +-- x x x x x x x x x x x x x
1.56 +-- m x x x x x x x x x x x x x
1.57 +--
1.58 +-- where: '*' - elements of the lower triangular part, '.' - structural
1.59 +-- zeros, 'x' - other (either non-zero or zero) elements.
1.60 +--
1.61 +-- The parameter info is a transit pointer passed to the formal routine
1.62 +-- mat (see below).
1.63 +--
1.64 +-- The formal routine mat specifies the given matrix A in both row- and
1.65 +-- column-wise formats. In order to obtain an i-th row of the matrix A
1.66 +-- the routine triang calls the routine mat with the parameter k = +i,
1.67 +-- 1 <= i <= m. In response the routine mat should store column indices
1.68 +-- of (non-zero) elements of the i-th row to the locations ndx[1], ...,
1.69 +-- ndx[len], where len is number of non-zeros in the i-th row returned
1.70 +-- on exit. Analogously, in order to obtain a j-th column of the matrix
1.71 +-- A, the routine mat is called with the parameter k = -j, 1 <= j <= n,
1.72 +-- and should return pattern of the j-th column in the same way as for
1.73 +-- row patterns. Note that the routine mat may be called more than once
1.74 +-- for the same rows and columns.
1.75 +--
1.76 +-- On exit the routine computes two resultant arrays rn and cn, which
1.77 +-- define the permutation matrices P and Q, respectively. The array rn
1.78 +-- should have at least 1+m locations, where rn[i] = i' (1 <= i <= m)
1.79 +-- means that i-th row of the original matrix A corresponds to i'-th row
1.80 +-- of the matrix B = P*A*Q. Similarly, the array cn should have at least
1.81 +-- 1+n locations, where cn[j] = j' (1 <= j <= n) means that j-th column
1.82 +-- of the matrix A corresponds to j'-th column of the matrix B.
1.83 +--
1.84 +-- *Returns*
1.85 +--
1.86 +-- The routine triang returns the size of the lower tringular part of
1.87 +-- the matrix B = P*A*Q (see the figure above).
1.88 +--
1.89 +-- *Complexity*
1.90 +--
1.91 +-- The time complexity of the routine triang is O(nnz), where nnz is
1.92 +-- number of non-zeros in the given matrix A.
1.93 +--
1.94 +-- *Algorithm*
1.95 +--
1.96 +-- The routine triang starts from the matrix B = P*Q*A, where P and Q
1.97 +-- are unity matrices, so initially B = A.
1.98 +--
1.99 +-- Before the next iteration B = (B1 | B2 | B3), where B1 is partially
1.100 +-- built a lower triangular submatrix, B2 is the active submatrix, and
1.101 +-- B3 is a submatrix that contains rejected columns. Thus, the current
1.102 +-- matrix B looks like follows (initially k1 = 1 and k2 = n):
1.103 +--
1.104 +-- 1 k1 k2 n
1.105 +-- 1 x . . . . . . . . . . . . . # # #
1.106 +-- x x . . . . . . . . . . . . # # #
1.107 +-- x x x . . . . . . . . . . # # # #
1.108 +-- x x x x . . . . . . . . . # # # #
1.109 +-- x x x x x . . . . . . . # # # # #
1.110 +-- k1 x x x x x * * * * * * * # # # # #
1.111 +-- x x x x x * * * * * * * # # # # #
1.112 +-- x x x x x * * * * * * * # # # # #
1.113 +-- x x x x x * * * * * * * # # # # #
1.114 +-- m x x x x x * * * * * * * # # # # #
1.115 +-- <--B1---> <----B2-----> <---B3-->
1.116 +--
1.117 +-- On each iteartion the routine looks for a singleton row, i.e. some
1.118 +-- row that has the only non-zero in the active submatrix B2. If such
1.119 +-- row exists and the corresponding non-zero is b[i,j], where (by the
1.120 +-- definition) k1 <= i <= m and k1 <= j <= k2, the routine permutes
1.121 +-- k1-th and i-th rows and k1-th and j-th columns of the matrix B (in
1.122 +-- order to place the element in the position b[k1,k1]), removes the
1.123 +-- k1-th column from the active submatrix B2, and adds this column to
1.124 +-- the submatrix B1. If no row singletons exist, but B2 is not empty
1.125 +-- yet, the routine chooses a j-th column, which has maximal number of
1.126 +-- non-zeros among other columns of B2, removes this column from B2 and
1.127 +-- adds it to the submatrix B3 in the hope that new row singletons will
1.128 +-- appear in the active submatrix. */
1.129 +
1.130 +static int triang(int m, int n,
1.131 + void *info, int (*mat)(void *info, int k, int ndx[]),
1.132 + int rn[], int cn[])
1.133 +{ int *ndx; /* int ndx[1+max(m,n)]; */
1.134 + /* this array is used for querying row and column patterns of the
1.135 + given matrix A (the third parameter to the routine mat) */
1.136 + int *rs_len; /* int rs_len[1+m]; */
1.137 + /* rs_len[0] is not used;
1.138 + rs_len[i], 1 <= i <= m, is number of non-zeros in the i-th row
1.139 + of the matrix A, which (non-zeros) belong to the current active
1.140 + submatrix */
1.141 + int *rs_head; /* int rs_head[1+n]; */
1.142 + /* rs_head[len], 0 <= len <= n, is the number i of the first row
1.143 + of the matrix A, for which rs_len[i] = len */
1.144 + int *rs_prev; /* int rs_prev[1+m]; */
1.145 + /* rs_prev[0] is not used;
1.146 + rs_prev[i], 1 <= i <= m, is a number i' of the previous row of
1.147 + the matrix A, for which rs_len[i] = rs_len[i'] (zero marks the
1.148 + end of this linked list) */
1.149 + int *rs_next; /* int rs_next[1+m]; */
1.150 + /* rs_next[0] is not used;
1.151 + rs_next[i], 1 <= i <= m, is a number i' of the next row of the
1.152 + matrix A, for which rs_len[i] = rs_len[i'] (zero marks the end
1.153 + this linked list) */
1.154 + int cs_head;
1.155 + /* is a number j of the first column of the matrix A, which has
1.156 + maximal number of non-zeros among other columns */
1.157 + int *cs_prev; /* cs_prev[1+n]; */
1.158 + /* cs_prev[0] is not used;
1.159 + cs_prev[j], 1 <= j <= n, is a number of the previous column of
1.160 + the matrix A with the same or greater number of non-zeros than
1.161 + in the j-th column (zero marks the end of this linked list) */
1.162 + int *cs_next; /* cs_next[1+n]; */
1.163 + /* cs_next[0] is not used;
1.164 + cs_next[j], 1 <= j <= n, is a number of the next column of
1.165 + the matrix A with the same or lesser number of non-zeros than
1.166 + in the j-th column (zero marks the end of this linked list) */
1.167 + int i, j, ii, jj, k1, k2, len, t, size = 0;
1.168 + int *head, *rn_inv, *cn_inv;
1.169 + if (!(m > 0 && n > 0))
1.170 + xerror("triang: m = %d; n = %d; invalid dimension\n", m, n);
1.171 + /* allocate working arrays */
1.172 + ndx = xcalloc(1+(m >= n ? m : n), sizeof(int));
1.173 + rs_len = xcalloc(1+m, sizeof(int));
1.174 + rs_head = xcalloc(1+n, sizeof(int));
1.175 + rs_prev = xcalloc(1+m, sizeof(int));
1.176 + rs_next = xcalloc(1+m, sizeof(int));
1.177 + cs_prev = xcalloc(1+n, sizeof(int));
1.178 + cs_next = xcalloc(1+n, sizeof(int));
1.179 + /* build linked lists of columns of the matrix A with the same
1.180 + number of non-zeros */
1.181 + head = rs_len; /* currently rs_len is used as working array */
1.182 + for (len = 0; len <= m; len ++) head[len] = 0;
1.183 + for (j = 1; j <= n; j++)
1.184 + { /* obtain length of the j-th column */
1.185 + len = mat(info, -j, ndx);
1.186 + xassert(0 <= len && len <= m);
1.187 + /* include the j-th column in the corresponding linked list */
1.188 + cs_prev[j] = head[len];
1.189 + head[len] = j;
1.190 + }
1.191 + /* merge all linked lists of columns in one linked list, where
1.192 + columns are ordered by descending of their lengths */
1.193 + cs_head = 0;
1.194 + for (len = 0; len <= m; len++)
1.195 + { for (j = head[len]; j != 0; j = cs_prev[j])
1.196 + { cs_next[j] = cs_head;
1.197 + cs_head = j;
1.198 + }
1.199 + }
1.200 + jj = 0;
1.201 + for (j = cs_head; j != 0; j = cs_next[j])
1.202 + { cs_prev[j] = jj;
1.203 + jj = j;
1.204 + }
1.205 + /* build initial doubly linked lists of rows of the matrix A with
1.206 + the same number of non-zeros */
1.207 + for (len = 0; len <= n; len++) rs_head[len] = 0;
1.208 + for (i = 1; i <= m; i++)
1.209 + { /* obtain length of the i-th row */
1.210 + rs_len[i] = len = mat(info, +i, ndx);
1.211 + xassert(0 <= len && len <= n);
1.212 + /* include the i-th row in the correspondng linked list */
1.213 + rs_prev[i] = 0;
1.214 + rs_next[i] = rs_head[len];
1.215 + if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
1.216 + rs_head[len] = i;
1.217 + }
1.218 + /* initially all rows and columns of the matrix A are active */
1.219 + for (i = 1; i <= m; i++) rn[i] = 0;
1.220 + for (j = 1; j <= n; j++) cn[j] = 0;
1.221 + /* set initial bounds of the active submatrix */
1.222 + k1 = 1, k2 = n;
1.223 + /* main loop starts here */
1.224 + while (k1 <= k2)
1.225 + { i = rs_head[1];
1.226 + if (i != 0)
1.227 + { /* the i-th row of the matrix A is a row singleton, since
1.228 + it has the only non-zero in the active submatrix */
1.229 + xassert(rs_len[i] == 1);
1.230 + /* determine the number j of an active column of the matrix
1.231 + A, in which this non-zero is placed */
1.232 + j = 0;
1.233 + t = mat(info, +i, ndx);
1.234 + xassert(0 <= t && t <= n);
1.235 + for (t = t; t >= 1; t--)
1.236 + { jj = ndx[t];
1.237 + xassert(1 <= jj && jj <= n);
1.238 + if (cn[jj] == 0)
1.239 + { xassert(j == 0);
1.240 + j = jj;
1.241 + }
1.242 + }
1.243 + xassert(j != 0);
1.244 + /* the singleton is a[i,j]; move a[i,j] to the position
1.245 + b[k1,k1] of the matrix B */
1.246 + rn[i] = cn[j] = k1;
1.247 + /* shift the left bound of the active submatrix */
1.248 + k1++;
1.249 + /* increase the size of the lower triangular part */
1.250 + size++;
1.251 + }
1.252 + else
1.253 + { /* the current active submatrix has no row singletons */
1.254 + /* remove an active column with maximal number of non-zeros
1.255 + from the active submatrix */
1.256 + j = cs_head;
1.257 + xassert(j != 0);
1.258 + cn[j] = k2;
1.259 + /* shift the right bound of the active submatrix */
1.260 + k2--;
1.261 + }
1.262 + /* the j-th column of the matrix A has been removed from the
1.263 + active submatrix */
1.264 + /* remove the j-th column from the linked list */
1.265 + if (cs_prev[j] == 0)
1.266 + cs_head = cs_next[j];
1.267 + else
1.268 + cs_next[cs_prev[j]] = cs_next[j];
1.269 + if (cs_next[j] == 0)
1.270 + /* nop */;
1.271 + else
1.272 + cs_prev[cs_next[j]] = cs_prev[j];
1.273 + /* go through non-zeros of the j-th columns and update active
1.274 + lengths of the corresponding rows */
1.275 + t = mat(info, -j, ndx);
1.276 + xassert(0 <= t && t <= m);
1.277 + for (t = t; t >= 1; t--)
1.278 + { i = ndx[t];
1.279 + xassert(1 <= i && i <= m);
1.280 + /* the non-zero a[i,j] has left the active submatrix */
1.281 + len = rs_len[i];
1.282 + xassert(len >= 1);
1.283 + /* remove the i-th row from the linked list of rows with
1.284 + active length len */
1.285 + if (rs_prev[i] == 0)
1.286 + rs_head[len] = rs_next[i];
1.287 + else
1.288 + rs_next[rs_prev[i]] = rs_next[i];
1.289 + if (rs_next[i] == 0)
1.290 + /* nop */;
1.291 + else
1.292 + rs_prev[rs_next[i]] = rs_prev[i];
1.293 + /* decrease the active length of the i-th row */
1.294 + rs_len[i] = --len;
1.295 + /* return the i-th row to the corresponding linked list */
1.296 + rs_prev[i] = 0;
1.297 + rs_next[i] = rs_head[len];
1.298 + if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
1.299 + rs_head[len] = i;
1.300 + }
1.301 + }
1.302 + /* other rows of the matrix A, which are still active, correspond
1.303 + to rows k1, ..., m of the matrix B (in arbitrary order) */
1.304 + for (i = 1; i <= m; i++) if (rn[i] == 0) rn[i] = k1++;
1.305 + /* but for columns this is not needed, because now the submatrix
1.306 + B2 has no columns */
1.307 + for (j = 1; j <= n; j++) xassert(cn[j] != 0);
1.308 + /* perform some optional checks */
1.309 + /* make sure that rn is a permutation of {1, ..., m} and cn is a
1.310 + permutation of {1, ..., n} */
1.311 + rn_inv = rs_len; /* used as working array */
1.312 + for (ii = 1; ii <= m; ii++) rn_inv[ii] = 0;
1.313 + for (i = 1; i <= m; i++)
1.314 + { ii = rn[i];
1.315 + xassert(1 <= ii && ii <= m);
1.316 + xassert(rn_inv[ii] == 0);
1.317 + rn_inv[ii] = i;
1.318 + }
1.319 + cn_inv = rs_head; /* used as working array */
1.320 + for (jj = 1; jj <= n; jj++) cn_inv[jj] = 0;
1.321 + for (j = 1; j <= n; j++)
1.322 + { jj = cn[j];
1.323 + xassert(1 <= jj && jj <= n);
1.324 + xassert(cn_inv[jj] == 0);
1.325 + cn_inv[jj] = j;
1.326 + }
1.327 + /* make sure that the matrix B = P*A*Q really has the form, which
1.328 + was declared */
1.329 + for (ii = 1; ii <= size; ii++)
1.330 + { int diag = 0;
1.331 + i = rn_inv[ii];
1.332 + t = mat(info, +i, ndx);
1.333 + xassert(0 <= t && t <= n);
1.334 + for (t = t; t >= 1; t--)
1.335 + { j = ndx[t];
1.336 + xassert(1 <= j && j <= n);
1.337 + jj = cn[j];
1.338 + if (jj <= size) xassert(jj <= ii);
1.339 + if (jj == ii)
1.340 + { xassert(!diag);
1.341 + diag = 1;
1.342 + }
1.343 + }
1.344 + xassert(diag);
1.345 + }
1.346 + /* free working arrays */
1.347 + xfree(ndx);
1.348 + xfree(rs_len);
1.349 + xfree(rs_head);
1.350 + xfree(rs_prev);
1.351 + xfree(rs_next);
1.352 + xfree(cs_prev);
1.353 + xfree(cs_next);
1.354 + /* return to the calling program */
1.355 + return size;
1.356 +}
1.357 +
1.358 +/*----------------------------------------------------------------------
1.359 +-- adv_basis - construct advanced initial LP basis.
1.360 +--
1.361 +-- *Synopsis*
1.362 +--
1.363 +-- #include "glpini.h"
1.364 +-- void adv_basis(glp_prob *lp);
1.365 +--
1.366 +-- *Description*
1.367 +--
1.368 +-- The routine adv_basis constructs an advanced initial basis for an LP
1.369 +-- problem object, which the parameter lp points to.
1.370 +--
1.371 +-- In order to build the initial basis the routine does the following:
1.372 +--
1.373 +-- 1) includes in the basis all non-fixed auxiliary variables;
1.374 +--
1.375 +-- 2) includes in the basis as many as possible non-fixed structural
1.376 +-- variables preserving triangular form of the basis matrix;
1.377 +--
1.378 +-- 3) includes in the basis appropriate (fixed) auxiliary variables
1.379 +-- in order to complete the basis.
1.380 +--
1.381 +-- As a result the initial basis has minimum of fixed variables and the
1.382 +-- corresponding basis matrix is triangular. */
1.383 +
1.384 +static int mat(void *info, int k, int ndx[])
1.385 +{ /* this auxiliary routine returns the pattern of a given row or
1.386 + a given column of the augmented constraint matrix A~ = (I|-A),
1.387 + in which columns of fixed variables are implicitly cleared */
1.388 + LPX *lp = info;
1.389 + int m = lpx_get_num_rows(lp);
1.390 + int n = lpx_get_num_cols(lp);
1.391 + int typx, i, j, lll, len = 0;
1.392 + if (k > 0)
1.393 + { /* the pattern of the i-th row is required */
1.394 + i = +k;
1.395 + xassert(1 <= i && i <= m);
1.396 +#if 0 /* 22/XII-2003 */
1.397 + /* if the auxiliary variable x[i] is non-fixed, include its
1.398 + element (placed in the i-th column) in the pattern */
1.399 + lpx_get_row_bnds(lp, i, &typx, NULL, NULL);
1.400 + if (typx != LPX_FX) ndx[++len] = i;
1.401 + /* include in the pattern elements placed in columns, which
1.402 + correspond to non-fixed structural varables */
1.403 + i_beg = aa_ptr[i];
1.404 + i_end = i_beg + aa_len[i] - 1;
1.405 + for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
1.406 + { j = m + sv_ndx[i_ptr];
1.407 + lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL);
1.408 + if (typx != LPX_FX) ndx[++len] = j;
1.409 + }
1.410 +#else
1.411 + lll = lpx_get_mat_row(lp, i, ndx, NULL);
1.412 + for (k = 1; k <= lll; k++)
1.413 + { lpx_get_col_bnds(lp, ndx[k], &typx, NULL, NULL);
1.414 + if (typx != LPX_FX) ndx[++len] = m + ndx[k];
1.415 + }
1.416 + lpx_get_row_bnds(lp, i, &typx, NULL, NULL);
1.417 + if (typx != LPX_FX) ndx[++len] = i;
1.418 +#endif
1.419 + }
1.420 + else
1.421 + { /* the pattern of the j-th column is required */
1.422 + j = -k;
1.423 + xassert(1 <= j && j <= m+n);
1.424 + /* if the (auxiliary or structural) variable x[j] is fixed,
1.425 + the pattern of its column is empty */
1.426 + if (j <= m)
1.427 + lpx_get_row_bnds(lp, j, &typx, NULL, NULL);
1.428 + else
1.429 + lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL);
1.430 + if (typx != LPX_FX)
1.431 + { if (j <= m)
1.432 + { /* x[j] is non-fixed auxiliary variable */
1.433 + ndx[++len] = j;
1.434 + }
1.435 + else
1.436 + { /* x[j] is non-fixed structural variables */
1.437 +#if 0 /* 22/XII-2003 */
1.438 + j_beg = aa_ptr[j];
1.439 + j_end = j_beg + aa_len[j] - 1;
1.440 + for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
1.441 + ndx[++len] = sv_ndx[j_ptr];
1.442 +#else
1.443 + len = lpx_get_mat_col(lp, j-m, ndx, NULL);
1.444 +#endif
1.445 + }
1.446 + }
1.447 + }
1.448 + /* return the length of the row/column pattern */
1.449 + return len;
1.450 +}
1.451 +
1.452 +static void adv_basis(glp_prob *lp)
1.453 +{ int m = lpx_get_num_rows(lp);
1.454 + int n = lpx_get_num_cols(lp);
1.455 + int i, j, jj, k, size;
1.456 + int *rn, *cn, *rn_inv, *cn_inv;
1.457 + int typx, *tagx = xcalloc(1+m+n, sizeof(int));
1.458 + double lb, ub;
1.459 + xprintf("Constructing initial basis...\n");
1.460 +#if 0 /* 13/V-2009 */
1.461 + if (m == 0)
1.462 + xerror("glp_adv_basis: problem has no rows\n");
1.463 + if (n == 0)
1.464 + xerror("glp_adv_basis: problem has no columns\n");
1.465 +#else
1.466 + if (m == 0 || n == 0)
1.467 + { glp_std_basis(lp);
1.468 + return;
1.469 + }
1.470 +#endif
1.471 + /* use the routine triang (see above) to find maximal triangular
1.472 + part of the augmented constraint matrix A~ = (I|-A); in order
1.473 + to prevent columns of fixed variables to be included in the
1.474 + triangular part, such columns are implictly removed from the
1.475 + matrix A~ by the routine adv_mat */
1.476 + rn = xcalloc(1+m, sizeof(int));
1.477 + cn = xcalloc(1+m+n, sizeof(int));
1.478 + size = triang(m, m+n, lp, mat, rn, cn);
1.479 + if (lpx_get_int_parm(lp, LPX_K_MSGLEV) >= 3)
1.480 + xprintf("Size of triangular part = %d\n", size);
1.481 + /* the first size rows and columns of the matrix P*A~*Q (where
1.482 + P and Q are permutation matrices defined by the arrays rn and
1.483 + cn) form a lower triangular matrix; build the arrays (rn_inv
1.484 + and cn_inv), which define the matrices inv(P) and inv(Q) */
1.485 + rn_inv = xcalloc(1+m, sizeof(int));
1.486 + cn_inv = xcalloc(1+m+n, sizeof(int));
1.487 + for (i = 1; i <= m; i++) rn_inv[rn[i]] = i;
1.488 + for (j = 1; j <= m+n; j++) cn_inv[cn[j]] = j;
1.489 + /* include the columns of the matrix A~, which correspond to the
1.490 + first size columns of the matrix P*A~*Q, in the basis */
1.491 + for (k = 1; k <= m+n; k++) tagx[k] = -1;
1.492 + for (jj = 1; jj <= size; jj++)
1.493 + { j = cn_inv[jj];
1.494 + /* the j-th column of A~ is the jj-th column of P*A~*Q */
1.495 + tagx[j] = LPX_BS;
1.496 + }
1.497 + /* if size < m, we need to add appropriate columns of auxiliary
1.498 + variables to the basis */
1.499 + for (jj = size + 1; jj <= m; jj++)
1.500 + { /* the jj-th column of P*A~*Q should be replaced by the column
1.501 + of the auxiliary variable, for which the only unity element
1.502 + is placed in the position [jj,jj] */
1.503 + i = rn_inv[jj];
1.504 + /* the jj-th row of P*A~*Q is the i-th row of A~, but in the
1.505 + i-th row of A~ the unity element belongs to the i-th column
1.506 + of A~; therefore the disired column corresponds to the i-th
1.507 + auxiliary variable (note that this column doesn't belong to
1.508 + the triangular part found by the routine triang) */
1.509 + xassert(1 <= i && i <= m);
1.510 + xassert(cn[i] > size);
1.511 + tagx[i] = LPX_BS;
1.512 + }
1.513 + /* free working arrays */
1.514 + xfree(rn);
1.515 + xfree(cn);
1.516 + xfree(rn_inv);
1.517 + xfree(cn_inv);
1.518 + /* build tags of non-basic variables */
1.519 + for (k = 1; k <= m+n; k++)
1.520 + { if (tagx[k] != LPX_BS)
1.521 + { if (k <= m)
1.522 + lpx_get_row_bnds(lp, k, &typx, &lb, &ub);
1.523 + else
1.524 + lpx_get_col_bnds(lp, k-m, &typx, &lb, &ub);
1.525 + switch (typx)
1.526 + { case LPX_FR:
1.527 + tagx[k] = LPX_NF; break;
1.528 + case LPX_LO:
1.529 + tagx[k] = LPX_NL; break;
1.530 + case LPX_UP:
1.531 + tagx[k] = LPX_NU; break;
1.532 + case LPX_DB:
1.533 + tagx[k] =
1.534 + (fabs(lb) <= fabs(ub) ? LPX_NL : LPX_NU);
1.535 + break;
1.536 + case LPX_FX:
1.537 + tagx[k] = LPX_NS; break;
1.538 + default:
1.539 + xassert(typx != typx);
1.540 + }
1.541 + }
1.542 + }
1.543 + for (k = 1; k <= m+n; k++)
1.544 + { if (k <= m)
1.545 + lpx_set_row_stat(lp, k, tagx[k]);
1.546 + else
1.547 + lpx_set_col_stat(lp, k-m, tagx[k]);
1.548 + }
1.549 + xfree(tagx);
1.550 + return;
1.551 +}
1.552 +
1.553 +/***********************************************************************
1.554 +* NAME
1.555 +*
1.556 +* glp_adv_basis - construct advanced initial LP basis
1.557 +*
1.558 +* SYNOPSIS
1.559 +*
1.560 +* void glp_adv_basis(glp_prob *lp, int flags);
1.561 +*
1.562 +* DESCRIPTION
1.563 +*
1.564 +* The routine glp_adv_basis constructs an advanced initial basis for
1.565 +* the specified problem object.
1.566 +*
1.567 +* The parameter flags is reserved for use in the future and must be
1.568 +* specified as zero. */
1.569 +
1.570 +void glp_adv_basis(glp_prob *lp, int flags)
1.571 +{ if (flags != 0)
1.572 + xerror("glp_adv_basis: flags = %d; invalid flags\n", flags);
1.573 + if (lp->m == 0 || lp->n == 0)
1.574 + glp_std_basis(lp);
1.575 + else
1.576 + adv_basis(lp);
1.577 + return;
1.578 +}
1.579 +
1.580 +/* eof */