[9] | 1 | /* glpini01.c */ |
---|
| 2 | |
---|
| 3 | /*********************************************************************** |
---|
| 4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
| 5 | * |
---|
| 6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
| 7 | * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, |
---|
| 8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
| 9 | * E-mail: <mao@gnu.org>. |
---|
| 10 | * |
---|
| 11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
| 12 | * under the terms of the GNU General Public License as published by |
---|
| 13 | * the Free Software Foundation, either version 3 of the License, or |
---|
| 14 | * (at your option) any later version. |
---|
| 15 | * |
---|
| 16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
| 17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
| 18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
| 19 | * License for more details. |
---|
| 20 | * |
---|
| 21 | * You should have received a copy of the GNU General Public License |
---|
| 22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
| 23 | ***********************************************************************/ |
---|
| 24 | |
---|
| 25 | #include "glpapi.h" |
---|
| 26 | |
---|
| 27 | /*---------------------------------------------------------------------- |
---|
| 28 | -- triang - find maximal triangular part of a rectangular matrix. |
---|
| 29 | -- |
---|
| 30 | -- *Synopsis* |
---|
| 31 | -- |
---|
| 32 | -- int triang(int m, int n, |
---|
| 33 | -- void *info, int (*mat)(void *info, int k, int ndx[]), |
---|
| 34 | -- int rn[], int cn[]); |
---|
| 35 | -- |
---|
| 36 | -- *Description* |
---|
| 37 | -- |
---|
| 38 | -- For a given rectangular (sparse) matrix A with m rows and n columns |
---|
| 39 | -- the routine triang tries to find such permutation matrices P and Q |
---|
| 40 | -- that the first rows and columns of the matrix B = P*A*Q form a lower |
---|
| 41 | -- triangular submatrix of as greatest size as possible: |
---|
| 42 | -- |
---|
| 43 | -- 1 n |
---|
| 44 | -- 1 * . . . . . . x x x x x x |
---|
| 45 | -- * * . . . . . x x x x x x |
---|
| 46 | -- * * * . . . . x x x x x x |
---|
| 47 | -- * * * * . . . x x x x x x |
---|
| 48 | -- B = P*A*Q = * * * * * . . x x x x x x |
---|
| 49 | -- * * * * * * . x x x x x x |
---|
| 50 | -- * * * * * * * x x x x x x |
---|
| 51 | -- x x x x x x x x x x x x x |
---|
| 52 | -- x x x x x x x x x x x x x |
---|
| 53 | -- m x x x x x x x x x x x x x |
---|
| 54 | -- |
---|
| 55 | -- where: '*' - elements of the lower triangular part, '.' - structural |
---|
| 56 | -- zeros, 'x' - other (either non-zero or zero) elements. |
---|
| 57 | -- |
---|
| 58 | -- The parameter info is a transit pointer passed to the formal routine |
---|
| 59 | -- mat (see below). |
---|
| 60 | -- |
---|
| 61 | -- The formal routine mat specifies the given matrix A in both row- and |
---|
| 62 | -- column-wise formats. In order to obtain an i-th row of the matrix A |
---|
| 63 | -- the routine triang calls the routine mat with the parameter k = +i, |
---|
| 64 | -- 1 <= i <= m. In response the routine mat should store column indices |
---|
| 65 | -- of (non-zero) elements of the i-th row to the locations ndx[1], ..., |
---|
| 66 | -- ndx[len], where len is number of non-zeros in the i-th row returned |
---|
| 67 | -- on exit. Analogously, in order to obtain a j-th column of the matrix |
---|
| 68 | -- A, the routine mat is called with the parameter k = -j, 1 <= j <= n, |
---|
| 69 | -- and should return pattern of the j-th column in the same way as for |
---|
| 70 | -- row patterns. Note that the routine mat may be called more than once |
---|
| 71 | -- for the same rows and columns. |
---|
| 72 | -- |
---|
| 73 | -- On exit the routine computes two resultant arrays rn and cn, which |
---|
| 74 | -- define the permutation matrices P and Q, respectively. The array rn |
---|
| 75 | -- should have at least 1+m locations, where rn[i] = i' (1 <= i <= m) |
---|
| 76 | -- means that i-th row of the original matrix A corresponds to i'-th row |
---|
| 77 | -- of the matrix B = P*A*Q. Similarly, the array cn should have at least |
---|
| 78 | -- 1+n locations, where cn[j] = j' (1 <= j <= n) means that j-th column |
---|
| 79 | -- of the matrix A corresponds to j'-th column of the matrix B. |
---|
| 80 | -- |
---|
| 81 | -- *Returns* |
---|
| 82 | -- |
---|
| 83 | -- The routine triang returns the size of the lower tringular part of |
---|
| 84 | -- the matrix B = P*A*Q (see the figure above). |
---|
| 85 | -- |
---|
| 86 | -- *Complexity* |
---|
| 87 | -- |
---|
| 88 | -- The time complexity of the routine triang is O(nnz), where nnz is |
---|
| 89 | -- number of non-zeros in the given matrix A. |
---|
| 90 | -- |
---|
| 91 | -- *Algorithm* |
---|
| 92 | -- |
---|
| 93 | -- The routine triang starts from the matrix B = P*Q*A, where P and Q |
---|
| 94 | -- are unity matrices, so initially B = A. |
---|
| 95 | -- |
---|
| 96 | -- Before the next iteration B = (B1 | B2 | B3), where B1 is partially |
---|
| 97 | -- built a lower triangular submatrix, B2 is the active submatrix, and |
---|
| 98 | -- B3 is a submatrix that contains rejected columns. Thus, the current |
---|
| 99 | -- matrix B looks like follows (initially k1 = 1 and k2 = n): |
---|
| 100 | -- |
---|
| 101 | -- 1 k1 k2 n |
---|
| 102 | -- 1 x . . . . . . . . . . . . . # # # |
---|
| 103 | -- x x . . . . . . . . . . . . # # # |
---|
| 104 | -- x x x . . . . . . . . . . # # # # |
---|
| 105 | -- x x x x . . . . . . . . . # # # # |
---|
| 106 | -- x x x x x . . . . . . . # # # # # |
---|
| 107 | -- k1 x x x x x * * * * * * * # # # # # |
---|
| 108 | -- x x x x x * * * * * * * # # # # # |
---|
| 109 | -- x x x x x * * * * * * * # # # # # |
---|
| 110 | -- x x x x x * * * * * * * # # # # # |
---|
| 111 | -- m x x x x x * * * * * * * # # # # # |
---|
| 112 | -- <--B1---> <----B2-----> <---B3--> |
---|
| 113 | -- |
---|
| 114 | -- On each iteartion the routine looks for a singleton row, i.e. some |
---|
| 115 | -- row that has the only non-zero in the active submatrix B2. If such |
---|
| 116 | -- row exists and the corresponding non-zero is b[i,j], where (by the |
---|
| 117 | -- definition) k1 <= i <= m and k1 <= j <= k2, the routine permutes |
---|
| 118 | -- k1-th and i-th rows and k1-th and j-th columns of the matrix B (in |
---|
| 119 | -- order to place the element in the position b[k1,k1]), removes the |
---|
| 120 | -- k1-th column from the active submatrix B2, and adds this column to |
---|
| 121 | -- the submatrix B1. If no row singletons exist, but B2 is not empty |
---|
| 122 | -- yet, the routine chooses a j-th column, which has maximal number of |
---|
| 123 | -- non-zeros among other columns of B2, removes this column from B2 and |
---|
| 124 | -- adds it to the submatrix B3 in the hope that new row singletons will |
---|
| 125 | -- appear in the active submatrix. */ |
---|
| 126 | |
---|
| 127 | static int triang(int m, int n, |
---|
| 128 | void *info, int (*mat)(void *info, int k, int ndx[]), |
---|
| 129 | int rn[], int cn[]) |
---|
| 130 | { int *ndx; /* int ndx[1+max(m,n)]; */ |
---|
| 131 | /* this array is used for querying row and column patterns of the |
---|
| 132 | given matrix A (the third parameter to the routine mat) */ |
---|
| 133 | int *rs_len; /* int rs_len[1+m]; */ |
---|
| 134 | /* rs_len[0] is not used; |
---|
| 135 | rs_len[i], 1 <= i <= m, is number of non-zeros in the i-th row |
---|
| 136 | of the matrix A, which (non-zeros) belong to the current active |
---|
| 137 | submatrix */ |
---|
| 138 | int *rs_head; /* int rs_head[1+n]; */ |
---|
| 139 | /* rs_head[len], 0 <= len <= n, is the number i of the first row |
---|
| 140 | of the matrix A, for which rs_len[i] = len */ |
---|
| 141 | int *rs_prev; /* int rs_prev[1+m]; */ |
---|
| 142 | /* rs_prev[0] is not used; |
---|
| 143 | rs_prev[i], 1 <= i <= m, is a number i' of the previous row of |
---|
| 144 | the matrix A, for which rs_len[i] = rs_len[i'] (zero marks the |
---|
| 145 | end of this linked list) */ |
---|
| 146 | int *rs_next; /* int rs_next[1+m]; */ |
---|
| 147 | /* rs_next[0] is not used; |
---|
| 148 | rs_next[i], 1 <= i <= m, is a number i' of the next row of the |
---|
| 149 | matrix A, for which rs_len[i] = rs_len[i'] (zero marks the end |
---|
| 150 | this linked list) */ |
---|
| 151 | int cs_head; |
---|
| 152 | /* is a number j of the first column of the matrix A, which has |
---|
| 153 | maximal number of non-zeros among other columns */ |
---|
| 154 | int *cs_prev; /* cs_prev[1+n]; */ |
---|
| 155 | /* cs_prev[0] is not used; |
---|
| 156 | cs_prev[j], 1 <= j <= n, is a number of the previous column of |
---|
| 157 | the matrix A with the same or greater number of non-zeros than |
---|
| 158 | in the j-th column (zero marks the end of this linked list) */ |
---|
| 159 | int *cs_next; /* cs_next[1+n]; */ |
---|
| 160 | /* cs_next[0] is not used; |
---|
| 161 | cs_next[j], 1 <= j <= n, is a number of the next column of |
---|
| 162 | the matrix A with the same or lesser number of non-zeros than |
---|
| 163 | in the j-th column (zero marks the end of this linked list) */ |
---|
| 164 | int i, j, ii, jj, k1, k2, len, t, size = 0; |
---|
| 165 | int *head, *rn_inv, *cn_inv; |
---|
| 166 | if (!(m > 0 && n > 0)) |
---|
| 167 | xerror("triang: m = %d; n = %d; invalid dimension\n", m, n); |
---|
| 168 | /* allocate working arrays */ |
---|
| 169 | ndx = xcalloc(1+(m >= n ? m : n), sizeof(int)); |
---|
| 170 | rs_len = xcalloc(1+m, sizeof(int)); |
---|
| 171 | rs_head = xcalloc(1+n, sizeof(int)); |
---|
| 172 | rs_prev = xcalloc(1+m, sizeof(int)); |
---|
| 173 | rs_next = xcalloc(1+m, sizeof(int)); |
---|
| 174 | cs_prev = xcalloc(1+n, sizeof(int)); |
---|
| 175 | cs_next = xcalloc(1+n, sizeof(int)); |
---|
| 176 | /* build linked lists of columns of the matrix A with the same |
---|
| 177 | number of non-zeros */ |
---|
| 178 | head = rs_len; /* currently rs_len is used as working array */ |
---|
| 179 | for (len = 0; len <= m; len ++) head[len] = 0; |
---|
| 180 | for (j = 1; j <= n; j++) |
---|
| 181 | { /* obtain length of the j-th column */ |
---|
| 182 | len = mat(info, -j, ndx); |
---|
| 183 | xassert(0 <= len && len <= m); |
---|
| 184 | /* include the j-th column in the corresponding linked list */ |
---|
| 185 | cs_prev[j] = head[len]; |
---|
| 186 | head[len] = j; |
---|
| 187 | } |
---|
| 188 | /* merge all linked lists of columns in one linked list, where |
---|
| 189 | columns are ordered by descending of their lengths */ |
---|
| 190 | cs_head = 0; |
---|
| 191 | for (len = 0; len <= m; len++) |
---|
| 192 | { for (j = head[len]; j != 0; j = cs_prev[j]) |
---|
| 193 | { cs_next[j] = cs_head; |
---|
| 194 | cs_head = j; |
---|
| 195 | } |
---|
| 196 | } |
---|
| 197 | jj = 0; |
---|
| 198 | for (j = cs_head; j != 0; j = cs_next[j]) |
---|
| 199 | { cs_prev[j] = jj; |
---|
| 200 | jj = j; |
---|
| 201 | } |
---|
| 202 | /* build initial doubly linked lists of rows of the matrix A with |
---|
| 203 | the same number of non-zeros */ |
---|
| 204 | for (len = 0; len <= n; len++) rs_head[len] = 0; |
---|
| 205 | for (i = 1; i <= m; i++) |
---|
| 206 | { /* obtain length of the i-th row */ |
---|
| 207 | rs_len[i] = len = mat(info, +i, ndx); |
---|
| 208 | xassert(0 <= len && len <= n); |
---|
| 209 | /* include the i-th row in the correspondng linked list */ |
---|
| 210 | rs_prev[i] = 0; |
---|
| 211 | rs_next[i] = rs_head[len]; |
---|
| 212 | if (rs_next[i] != 0) rs_prev[rs_next[i]] = i; |
---|
| 213 | rs_head[len] = i; |
---|
| 214 | } |
---|
| 215 | /* initially all rows and columns of the matrix A are active */ |
---|
| 216 | for (i = 1; i <= m; i++) rn[i] = 0; |
---|
| 217 | for (j = 1; j <= n; j++) cn[j] = 0; |
---|
| 218 | /* set initial bounds of the active submatrix */ |
---|
| 219 | k1 = 1, k2 = n; |
---|
| 220 | /* main loop starts here */ |
---|
| 221 | while (k1 <= k2) |
---|
| 222 | { i = rs_head[1]; |
---|
| 223 | if (i != 0) |
---|
| 224 | { /* the i-th row of the matrix A is a row singleton, since |
---|
| 225 | it has the only non-zero in the active submatrix */ |
---|
| 226 | xassert(rs_len[i] == 1); |
---|
| 227 | /* determine the number j of an active column of the matrix |
---|
| 228 | A, in which this non-zero is placed */ |
---|
| 229 | j = 0; |
---|
| 230 | t = mat(info, +i, ndx); |
---|
| 231 | xassert(0 <= t && t <= n); |
---|
| 232 | for (t = t; t >= 1; t--) |
---|
| 233 | { jj = ndx[t]; |
---|
| 234 | xassert(1 <= jj && jj <= n); |
---|
| 235 | if (cn[jj] == 0) |
---|
| 236 | { xassert(j == 0); |
---|
| 237 | j = jj; |
---|
| 238 | } |
---|
| 239 | } |
---|
| 240 | xassert(j != 0); |
---|
| 241 | /* the singleton is a[i,j]; move a[i,j] to the position |
---|
| 242 | b[k1,k1] of the matrix B */ |
---|
| 243 | rn[i] = cn[j] = k1; |
---|
| 244 | /* shift the left bound of the active submatrix */ |
---|
| 245 | k1++; |
---|
| 246 | /* increase the size of the lower triangular part */ |
---|
| 247 | size++; |
---|
| 248 | } |
---|
| 249 | else |
---|
| 250 | { /* the current active submatrix has no row singletons */ |
---|
| 251 | /* remove an active column with maximal number of non-zeros |
---|
| 252 | from the active submatrix */ |
---|
| 253 | j = cs_head; |
---|
| 254 | xassert(j != 0); |
---|
| 255 | cn[j] = k2; |
---|
| 256 | /* shift the right bound of the active submatrix */ |
---|
| 257 | k2--; |
---|
| 258 | } |
---|
| 259 | /* the j-th column of the matrix A has been removed from the |
---|
| 260 | active submatrix */ |
---|
| 261 | /* remove the j-th column from the linked list */ |
---|
| 262 | if (cs_prev[j] == 0) |
---|
| 263 | cs_head = cs_next[j]; |
---|
| 264 | else |
---|
| 265 | cs_next[cs_prev[j]] = cs_next[j]; |
---|
| 266 | if (cs_next[j] == 0) |
---|
| 267 | /* nop */; |
---|
| 268 | else |
---|
| 269 | cs_prev[cs_next[j]] = cs_prev[j]; |
---|
| 270 | /* go through non-zeros of the j-th columns and update active |
---|
| 271 | lengths of the corresponding rows */ |
---|
| 272 | t = mat(info, -j, ndx); |
---|
| 273 | xassert(0 <= t && t <= m); |
---|
| 274 | for (t = t; t >= 1; t--) |
---|
| 275 | { i = ndx[t]; |
---|
| 276 | xassert(1 <= i && i <= m); |
---|
| 277 | /* the non-zero a[i,j] has left the active submatrix */ |
---|
| 278 | len = rs_len[i]; |
---|
| 279 | xassert(len >= 1); |
---|
| 280 | /* remove the i-th row from the linked list of rows with |
---|
| 281 | active length len */ |
---|
| 282 | if (rs_prev[i] == 0) |
---|
| 283 | rs_head[len] = rs_next[i]; |
---|
| 284 | else |
---|
| 285 | rs_next[rs_prev[i]] = rs_next[i]; |
---|
| 286 | if (rs_next[i] == 0) |
---|
| 287 | /* nop */; |
---|
| 288 | else |
---|
| 289 | rs_prev[rs_next[i]] = rs_prev[i]; |
---|
| 290 | /* decrease the active length of the i-th row */ |
---|
| 291 | rs_len[i] = --len; |
---|
| 292 | /* return the i-th row to the corresponding linked list */ |
---|
| 293 | rs_prev[i] = 0; |
---|
| 294 | rs_next[i] = rs_head[len]; |
---|
| 295 | if (rs_next[i] != 0) rs_prev[rs_next[i]] = i; |
---|
| 296 | rs_head[len] = i; |
---|
| 297 | } |
---|
| 298 | } |
---|
| 299 | /* other rows of the matrix A, which are still active, correspond |
---|
| 300 | to rows k1, ..., m of the matrix B (in arbitrary order) */ |
---|
| 301 | for (i = 1; i <= m; i++) if (rn[i] == 0) rn[i] = k1++; |
---|
| 302 | /* but for columns this is not needed, because now the submatrix |
---|
| 303 | B2 has no columns */ |
---|
| 304 | for (j = 1; j <= n; j++) xassert(cn[j] != 0); |
---|
| 305 | /* perform some optional checks */ |
---|
| 306 | /* make sure that rn is a permutation of {1, ..., m} and cn is a |
---|
| 307 | permutation of {1, ..., n} */ |
---|
| 308 | rn_inv = rs_len; /* used as working array */ |
---|
| 309 | for (ii = 1; ii <= m; ii++) rn_inv[ii] = 0; |
---|
| 310 | for (i = 1; i <= m; i++) |
---|
| 311 | { ii = rn[i]; |
---|
| 312 | xassert(1 <= ii && ii <= m); |
---|
| 313 | xassert(rn_inv[ii] == 0); |
---|
| 314 | rn_inv[ii] = i; |
---|
| 315 | } |
---|
| 316 | cn_inv = rs_head; /* used as working array */ |
---|
| 317 | for (jj = 1; jj <= n; jj++) cn_inv[jj] = 0; |
---|
| 318 | for (j = 1; j <= n; j++) |
---|
| 319 | { jj = cn[j]; |
---|
| 320 | xassert(1 <= jj && jj <= n); |
---|
| 321 | xassert(cn_inv[jj] == 0); |
---|
| 322 | cn_inv[jj] = j; |
---|
| 323 | } |
---|
| 324 | /* make sure that the matrix B = P*A*Q really has the form, which |
---|
| 325 | was declared */ |
---|
| 326 | for (ii = 1; ii <= size; ii++) |
---|
| 327 | { int diag = 0; |
---|
| 328 | i = rn_inv[ii]; |
---|
| 329 | t = mat(info, +i, ndx); |
---|
| 330 | xassert(0 <= t && t <= n); |
---|
| 331 | for (t = t; t >= 1; t--) |
---|
| 332 | { j = ndx[t]; |
---|
| 333 | xassert(1 <= j && j <= n); |
---|
| 334 | jj = cn[j]; |
---|
| 335 | if (jj <= size) xassert(jj <= ii); |
---|
| 336 | if (jj == ii) |
---|
| 337 | { xassert(!diag); |
---|
| 338 | diag = 1; |
---|
| 339 | } |
---|
| 340 | } |
---|
| 341 | xassert(diag); |
---|
| 342 | } |
---|
| 343 | /* free working arrays */ |
---|
| 344 | xfree(ndx); |
---|
| 345 | xfree(rs_len); |
---|
| 346 | xfree(rs_head); |
---|
| 347 | xfree(rs_prev); |
---|
| 348 | xfree(rs_next); |
---|
| 349 | xfree(cs_prev); |
---|
| 350 | xfree(cs_next); |
---|
| 351 | /* return to the calling program */ |
---|
| 352 | return size; |
---|
| 353 | } |
---|
| 354 | |
---|
| 355 | /*---------------------------------------------------------------------- |
---|
| 356 | -- adv_basis - construct advanced initial LP basis. |
---|
| 357 | -- |
---|
| 358 | -- *Synopsis* |
---|
| 359 | -- |
---|
| 360 | -- #include "glpini.h" |
---|
| 361 | -- void adv_basis(glp_prob *lp); |
---|
| 362 | -- |
---|
| 363 | -- *Description* |
---|
| 364 | -- |
---|
| 365 | -- The routine adv_basis constructs an advanced initial basis for an LP |
---|
| 366 | -- problem object, which the parameter lp points to. |
---|
| 367 | -- |
---|
| 368 | -- In order to build the initial basis the routine does the following: |
---|
| 369 | -- |
---|
| 370 | -- 1) includes in the basis all non-fixed auxiliary variables; |
---|
| 371 | -- |
---|
| 372 | -- 2) includes in the basis as many as possible non-fixed structural |
---|
| 373 | -- variables preserving triangular form of the basis matrix; |
---|
| 374 | -- |
---|
| 375 | -- 3) includes in the basis appropriate (fixed) auxiliary variables |
---|
| 376 | -- in order to complete the basis. |
---|
| 377 | -- |
---|
| 378 | -- As a result the initial basis has minimum of fixed variables and the |
---|
| 379 | -- corresponding basis matrix is triangular. */ |
---|
| 380 | |
---|
| 381 | static int mat(void *info, int k, int ndx[]) |
---|
| 382 | { /* this auxiliary routine returns the pattern of a given row or |
---|
| 383 | a given column of the augmented constraint matrix A~ = (I|-A), |
---|
| 384 | in which columns of fixed variables are implicitly cleared */ |
---|
| 385 | LPX *lp = info; |
---|
| 386 | int m = lpx_get_num_rows(lp); |
---|
| 387 | int n = lpx_get_num_cols(lp); |
---|
| 388 | int typx, i, j, lll, len = 0; |
---|
| 389 | if (k > 0) |
---|
| 390 | { /* the pattern of the i-th row is required */ |
---|
| 391 | i = +k; |
---|
| 392 | xassert(1 <= i && i <= m); |
---|
| 393 | #if 0 /* 22/XII-2003 */ |
---|
| 394 | /* if the auxiliary variable x[i] is non-fixed, include its |
---|
| 395 | element (placed in the i-th column) in the pattern */ |
---|
| 396 | lpx_get_row_bnds(lp, i, &typx, NULL, NULL); |
---|
| 397 | if (typx != LPX_FX) ndx[++len] = i; |
---|
| 398 | /* include in the pattern elements placed in columns, which |
---|
| 399 | correspond to non-fixed structural varables */ |
---|
| 400 | i_beg = aa_ptr[i]; |
---|
| 401 | i_end = i_beg + aa_len[i] - 1; |
---|
| 402 | for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) |
---|
| 403 | { j = m + sv_ndx[i_ptr]; |
---|
| 404 | lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL); |
---|
| 405 | if (typx != LPX_FX) ndx[++len] = j; |
---|
| 406 | } |
---|
| 407 | #else |
---|
| 408 | lll = lpx_get_mat_row(lp, i, ndx, NULL); |
---|
| 409 | for (k = 1; k <= lll; k++) |
---|
| 410 | { lpx_get_col_bnds(lp, ndx[k], &typx, NULL, NULL); |
---|
| 411 | if (typx != LPX_FX) ndx[++len] = m + ndx[k]; |
---|
| 412 | } |
---|
| 413 | lpx_get_row_bnds(lp, i, &typx, NULL, NULL); |
---|
| 414 | if (typx != LPX_FX) ndx[++len] = i; |
---|
| 415 | #endif |
---|
| 416 | } |
---|
| 417 | else |
---|
| 418 | { /* the pattern of the j-th column is required */ |
---|
| 419 | j = -k; |
---|
| 420 | xassert(1 <= j && j <= m+n); |
---|
| 421 | /* if the (auxiliary or structural) variable x[j] is fixed, |
---|
| 422 | the pattern of its column is empty */ |
---|
| 423 | if (j <= m) |
---|
| 424 | lpx_get_row_bnds(lp, j, &typx, NULL, NULL); |
---|
| 425 | else |
---|
| 426 | lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL); |
---|
| 427 | if (typx != LPX_FX) |
---|
| 428 | { if (j <= m) |
---|
| 429 | { /* x[j] is non-fixed auxiliary variable */ |
---|
| 430 | ndx[++len] = j; |
---|
| 431 | } |
---|
| 432 | else |
---|
| 433 | { /* x[j] is non-fixed structural variables */ |
---|
| 434 | #if 0 /* 22/XII-2003 */ |
---|
| 435 | j_beg = aa_ptr[j]; |
---|
| 436 | j_end = j_beg + aa_len[j] - 1; |
---|
| 437 | for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) |
---|
| 438 | ndx[++len] = sv_ndx[j_ptr]; |
---|
| 439 | #else |
---|
| 440 | len = lpx_get_mat_col(lp, j-m, ndx, NULL); |
---|
| 441 | #endif |
---|
| 442 | } |
---|
| 443 | } |
---|
| 444 | } |
---|
| 445 | /* return the length of the row/column pattern */ |
---|
| 446 | return len; |
---|
| 447 | } |
---|
| 448 | |
---|
| 449 | static void adv_basis(glp_prob *lp) |
---|
| 450 | { int m = lpx_get_num_rows(lp); |
---|
| 451 | int n = lpx_get_num_cols(lp); |
---|
| 452 | int i, j, jj, k, size; |
---|
| 453 | int *rn, *cn, *rn_inv, *cn_inv; |
---|
| 454 | int typx, *tagx = xcalloc(1+m+n, sizeof(int)); |
---|
| 455 | double lb, ub; |
---|
| 456 | xprintf("Constructing initial basis...\n"); |
---|
| 457 | #if 0 /* 13/V-2009 */ |
---|
| 458 | if (m == 0) |
---|
| 459 | xerror("glp_adv_basis: problem has no rows\n"); |
---|
| 460 | if (n == 0) |
---|
| 461 | xerror("glp_adv_basis: problem has no columns\n"); |
---|
| 462 | #else |
---|
| 463 | if (m == 0 || n == 0) |
---|
| 464 | { glp_std_basis(lp); |
---|
| 465 | return; |
---|
| 466 | } |
---|
| 467 | #endif |
---|
| 468 | /* use the routine triang (see above) to find maximal triangular |
---|
| 469 | part of the augmented constraint matrix A~ = (I|-A); in order |
---|
| 470 | to prevent columns of fixed variables to be included in the |
---|
| 471 | triangular part, such columns are implictly removed from the |
---|
| 472 | matrix A~ by the routine adv_mat */ |
---|
| 473 | rn = xcalloc(1+m, sizeof(int)); |
---|
| 474 | cn = xcalloc(1+m+n, sizeof(int)); |
---|
| 475 | size = triang(m, m+n, lp, mat, rn, cn); |
---|
| 476 | if (lpx_get_int_parm(lp, LPX_K_MSGLEV) >= 3) |
---|
| 477 | xprintf("Size of triangular part = %d\n", size); |
---|
| 478 | /* the first size rows and columns of the matrix P*A~*Q (where |
---|
| 479 | P and Q are permutation matrices defined by the arrays rn and |
---|
| 480 | cn) form a lower triangular matrix; build the arrays (rn_inv |
---|
| 481 | and cn_inv), which define the matrices inv(P) and inv(Q) */ |
---|
| 482 | rn_inv = xcalloc(1+m, sizeof(int)); |
---|
| 483 | cn_inv = xcalloc(1+m+n, sizeof(int)); |
---|
| 484 | for (i = 1; i <= m; i++) rn_inv[rn[i]] = i; |
---|
| 485 | for (j = 1; j <= m+n; j++) cn_inv[cn[j]] = j; |
---|
| 486 | /* include the columns of the matrix A~, which correspond to the |
---|
| 487 | first size columns of the matrix P*A~*Q, in the basis */ |
---|
| 488 | for (k = 1; k <= m+n; k++) tagx[k] = -1; |
---|
| 489 | for (jj = 1; jj <= size; jj++) |
---|
| 490 | { j = cn_inv[jj]; |
---|
| 491 | /* the j-th column of A~ is the jj-th column of P*A~*Q */ |
---|
| 492 | tagx[j] = LPX_BS; |
---|
| 493 | } |
---|
| 494 | /* if size < m, we need to add appropriate columns of auxiliary |
---|
| 495 | variables to the basis */ |
---|
| 496 | for (jj = size + 1; jj <= m; jj++) |
---|
| 497 | { /* the jj-th column of P*A~*Q should be replaced by the column |
---|
| 498 | of the auxiliary variable, for which the only unity element |
---|
| 499 | is placed in the position [jj,jj] */ |
---|
| 500 | i = rn_inv[jj]; |
---|
| 501 | /* the jj-th row of P*A~*Q is the i-th row of A~, but in the |
---|
| 502 | i-th row of A~ the unity element belongs to the i-th column |
---|
| 503 | of A~; therefore the disired column corresponds to the i-th |
---|
| 504 | auxiliary variable (note that this column doesn't belong to |
---|
| 505 | the triangular part found by the routine triang) */ |
---|
| 506 | xassert(1 <= i && i <= m); |
---|
| 507 | xassert(cn[i] > size); |
---|
| 508 | tagx[i] = LPX_BS; |
---|
| 509 | } |
---|
| 510 | /* free working arrays */ |
---|
| 511 | xfree(rn); |
---|
| 512 | xfree(cn); |
---|
| 513 | xfree(rn_inv); |
---|
| 514 | xfree(cn_inv); |
---|
| 515 | /* build tags of non-basic variables */ |
---|
| 516 | for (k = 1; k <= m+n; k++) |
---|
| 517 | { if (tagx[k] != LPX_BS) |
---|
| 518 | { if (k <= m) |
---|
| 519 | lpx_get_row_bnds(lp, k, &typx, &lb, &ub); |
---|
| 520 | else |
---|
| 521 | lpx_get_col_bnds(lp, k-m, &typx, &lb, &ub); |
---|
| 522 | switch (typx) |
---|
| 523 | { case LPX_FR: |
---|
| 524 | tagx[k] = LPX_NF; break; |
---|
| 525 | case LPX_LO: |
---|
| 526 | tagx[k] = LPX_NL; break; |
---|
| 527 | case LPX_UP: |
---|
| 528 | tagx[k] = LPX_NU; break; |
---|
| 529 | case LPX_DB: |
---|
| 530 | tagx[k] = |
---|
| 531 | (fabs(lb) <= fabs(ub) ? LPX_NL : LPX_NU); |
---|
| 532 | break; |
---|
| 533 | case LPX_FX: |
---|
| 534 | tagx[k] = LPX_NS; break; |
---|
| 535 | default: |
---|
| 536 | xassert(typx != typx); |
---|
| 537 | } |
---|
| 538 | } |
---|
| 539 | } |
---|
| 540 | for (k = 1; k <= m+n; k++) |
---|
| 541 | { if (k <= m) |
---|
| 542 | lpx_set_row_stat(lp, k, tagx[k]); |
---|
| 543 | else |
---|
| 544 | lpx_set_col_stat(lp, k-m, tagx[k]); |
---|
| 545 | } |
---|
| 546 | xfree(tagx); |
---|
| 547 | return; |
---|
| 548 | } |
---|
| 549 | |
---|
| 550 | /*********************************************************************** |
---|
| 551 | * NAME |
---|
| 552 | * |
---|
| 553 | * glp_adv_basis - construct advanced initial LP basis |
---|
| 554 | * |
---|
| 555 | * SYNOPSIS |
---|
| 556 | * |
---|
| 557 | * void glp_adv_basis(glp_prob *lp, int flags); |
---|
| 558 | * |
---|
| 559 | * DESCRIPTION |
---|
| 560 | * |
---|
| 561 | * The routine glp_adv_basis constructs an advanced initial basis for |
---|
| 562 | * the specified problem object. |
---|
| 563 | * |
---|
| 564 | * The parameter flags is reserved for use in the future and must be |
---|
| 565 | * specified as zero. */ |
---|
| 566 | |
---|
| 567 | void glp_adv_basis(glp_prob *lp, int flags) |
---|
| 568 | { if (flags != 0) |
---|
| 569 | xerror("glp_adv_basis: flags = %d; invalid flags\n", flags); |
---|
| 570 | if (lp->m == 0 || lp->n == 0) |
---|
| 571 | glp_std_basis(lp); |
---|
| 572 | else |
---|
| 573 | adv_basis(lp); |
---|
| 574 | return; |
---|
| 575 | } |
---|
| 576 | |
---|
| 577 | /* eof */ |
---|