lemon-project-template-glpk: 33de93886c88 deps/glpk/src/glpini01.c

lemon-project-template-glpk

view deps/glpk/src/glpini01.c @ 9:33de93886c88

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

line source

1 /* glpini01.c */

3 /***********************************************************************

4 * This code is part of GLPK (GNU Linear Programming Kit).

5 *

7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,

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 ***********************************************************************/

25 #include "glpapi.h"

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 */