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