lemon-project-template-glpk
comparison deps/glpk/src/glpnet02.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 |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:325e7da0cc6f |
---|---|
1 /* glpnet02.c (permutations to block triangular form) */ | |
2 | |
3 /*********************************************************************** | |
4 * This code is part of GLPK (GNU Linear Programming Kit). | |
5 * | |
6 * This code is the result of translation of the Fortran subroutines | |
7 * MC13D and MC13E associated with the following paper: | |
8 * | |
9 * I.S.Duff, J.K.Reid, Algorithm 529: Permutations to block triangular | |
10 * form, ACM Trans. on Math. Softw. 4 (1978), 189-192. | |
11 * | |
12 * Use of ACM Algorithms is subject to the ACM Software Copyright and | |
13 * License Agreement. See <http://www.acm.org/publications/policies>. | |
14 * | |
15 * The translation was made by Andrew Makhorin <mao@gnu.org>. | |
16 * | |
17 * GLPK is free software: you can redistribute it and/or modify it | |
18 * under the terms of the GNU General Public License as published by | |
19 * the Free Software Foundation, either version 3 of the License, or | |
20 * (at your option) any later version. | |
21 * | |
22 * GLPK is distributed in the hope that it will be useful, but WITHOUT | |
23 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY | |
24 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public | |
25 * License for more details. | |
26 * | |
27 * You should have received a copy of the GNU General Public License | |
28 * along with GLPK. If not, see <http://www.gnu.org/licenses/>. | |
29 ***********************************************************************/ | |
30 | |
31 #include "glpnet.h" | |
32 | |
33 /*********************************************************************** | |
34 * NAME | |
35 * | |
36 * mc13d - permutations to block triangular form | |
37 * | |
38 * SYNOPSIS | |
39 * | |
40 * #include "glpnet.h" | |
41 * int mc13d(int n, const int icn[], const int ip[], const int lenr[], | |
42 * int ior[], int ib[], int lowl[], int numb[], int prev[]); | |
43 * | |
44 * DESCRIPTION | |
45 * | |
46 * Given the column numbers of the nonzeros in each row of the sparse | |
47 * matrix, the routine mc13d finds a symmetric permutation that makes | |
48 * the matrix block lower triangular. | |
49 * | |
50 * INPUT PARAMETERS | |
51 * | |
52 * n order of the matrix. | |
53 * | |
54 * icn array containing the column indices of the non-zeros. Those | |
55 * belonging to a single row must be contiguous but the ordering | |
56 * of column indices within each row is unimportant and wasted | |
57 * space between rows is permitted. | |
58 * | |
59 * ip ip[i], i = 1,2,...,n, is the position in array icn of the | |
60 * first column index of a non-zero in row i. | |
61 * | |
62 * lenr lenr[i], i = 1,2,...,n, is the number of non-zeros in row i. | |
63 * | |
64 * OUTPUT PARAMETERS | |
65 * | |
66 * ior ior[i], i = 1,2,...,n, gives the position on the original | |
67 * ordering of the row or column which is in position i in the | |
68 * permuted form. | |
69 * | |
70 * ib ib[i], i = 1,2,...,num, is the row number in the permuted | |
71 * matrix of the beginning of block i, 1 <= num <= n. | |
72 * | |
73 * WORKING ARRAYS | |
74 * | |
75 * arp working array of length [1+n], where arp[0] is not used. | |
76 * arp[i] is one less than the number of unsearched edges leaving | |
77 * node i. At the end of the algorithm it is set to a permutation | |
78 * which puts the matrix in block lower triangular form. | |
79 * | |
80 * ib working array of length [1+n], where ib[0] is not used. | |
81 * ib[i] is the position in the ordering of the start of the ith | |
82 * block. ib[n+1-i] holds the node number of the ith node on the | |
83 * stack. | |
84 * | |
85 * lowl working array of length [1+n], where lowl[0] is not used. | |
86 * lowl[i] is the smallest stack position of any node to which a | |
87 * path from node i has been found. It is set to n+1 when node i | |
88 * is removed from the stack. | |
89 * | |
90 * numb working array of length [1+n], where numb[0] is not used. | |
91 * numb[i] is the position of node i in the stack if it is on it, | |
92 * is the permuted order of node i for those nodes whose final | |
93 * position has been found and is otherwise zero. | |
94 * | |
95 * prev working array of length [1+n], where prev[0] is not used. | |
96 * prev[i] is the node at the end of the path when node i was | |
97 * placed on the stack. | |
98 * | |
99 * RETURNS | |
100 * | |
101 * The routine mc13d returns num, the number of blocks found. */ | |
102 | |
103 int mc13d(int n, const int icn[], const int ip[], const int lenr[], | |
104 int ior[], int ib[], int lowl[], int numb[], int prev[]) | |
105 { int *arp = ior; | |
106 int dummy, i, i1, i2, icnt, ii, isn, ist, ist1, iv, iw, j, lcnt, | |
107 nnm1, num, stp; | |
108 /* icnt is the number of nodes whose positions in final ordering | |
109 have been found. */ | |
110 icnt = 0; | |
111 /* num is the number of blocks that have been found. */ | |
112 num = 0; | |
113 nnm1 = n + n - 1; | |
114 /* Initialization of arrays. */ | |
115 for (j = 1; j <= n; j++) | |
116 { numb[j] = 0; | |
117 arp[j] = lenr[j] - 1; | |
118 } | |
119 for (isn = 1; isn <= n; isn++) | |
120 { /* Look for a starting node. */ | |
121 if (numb[isn] != 0) continue; | |
122 iv = isn; | |
123 /* ist is the number of nodes on the stack ... it is the stack | |
124 pointer. */ | |
125 ist = 1; | |
126 /* Put node iv at beginning of stack. */ | |
127 lowl[iv] = numb[iv] = 1; | |
128 ib[n] = iv; | |
129 /* The body of this loop puts a new node on the stack or | |
130 backtracks. */ | |
131 for (dummy = 1; dummy <= nnm1; dummy++) | |
132 { i1 = arp[iv]; | |
133 /* Have all edges leaving node iv been searched? */ | |
134 if (i1 >= 0) | |
135 { i2 = ip[iv] + lenr[iv] - 1; | |
136 i1 = i2 - i1; | |
137 /* Look at edges leaving node iv until one enters a new | |
138 node or all edges are exhausted. */ | |
139 for (ii = i1; ii <= i2; ii++) | |
140 { iw = icn[ii]; | |
141 /* Has node iw been on stack already? */ | |
142 if (numb[iw] == 0) goto L70; | |
143 /* Update value of lowl[iv] if necessary. */ | |
144 if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw]; | |
145 } | |
146 /* There are no more edges leaving node iv. */ | |
147 arp[iv] = -1; | |
148 } | |
149 /* Is node iv the root of a block? */ | |
150 if (lowl[iv] < numb[iv]) goto L60; | |
151 /* Order nodes in a block. */ | |
152 num++; | |
153 ist1 = n + 1 - ist; | |
154 lcnt = icnt + 1; | |
155 /* Peel block off the top of the stack starting at the top | |
156 and working down to the root of the block. */ | |
157 for (stp = ist1; stp <= n; stp++) | |
158 { iw = ib[stp]; | |
159 lowl[iw] = n + 1; | |
160 numb[iw] = ++icnt; | |
161 if (iw == iv) break; | |
162 } | |
163 ist = n - stp; | |
164 ib[num] = lcnt; | |
165 /* Are there any nodes left on the stack? */ | |
166 if (ist != 0) goto L60; | |
167 /* Have all the nodes been ordered? */ | |
168 if (icnt < n) break; | |
169 goto L100; | |
170 L60: /* Backtrack to previous node on path. */ | |
171 iw = iv; | |
172 iv = prev[iv]; | |
173 /* Update value of lowl[iv] if necessary. */ | |
174 if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw]; | |
175 continue; | |
176 L70: /* Put new node on the stack. */ | |
177 arp[iv] = i2 - ii - 1; | |
178 prev[iw] = iv; | |
179 iv = iw; | |
180 lowl[iv] = numb[iv] = ++ist; | |
181 ib[n+1-ist] = iv; | |
182 } | |
183 } | |
184 L100: /* Put permutation in the required form. */ | |
185 for (i = 1; i <= n; i++) | |
186 arp[numb[i]] = i; | |
187 return num; | |
188 } | |
189 | |
190 /**********************************************************************/ | |
191 | |
192 #if 0 | |
193 #include "glplib.h" | |
194 | |
195 void test(int n, int ipp); | |
196 | |
197 int main(void) | |
198 { /* test program for routine mc13d */ | |
199 test( 1, 0); | |
200 test( 2, 1); | |
201 test( 2, 2); | |
202 test( 3, 3); | |
203 test( 4, 4); | |
204 test( 5, 10); | |
205 test(10, 10); | |
206 test(10, 20); | |
207 test(20, 20); | |
208 test(20, 50); | |
209 test(50, 50); | |
210 test(50, 200); | |
211 return 0; | |
212 } | |
213 | |
214 void fa01bs(int max, int *nrand); | |
215 | |
216 void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[]); | |
217 | |
218 void test(int n, int ipp) | |
219 { int ip[1+50], icn[1+1000], ior[1+50], ib[1+51], iw[1+150], | |
220 lenr[1+50]; | |
221 char a[1+50][1+50], hold[1+100]; | |
222 int i, ii, iblock, ij, index, j, jblock, jj, k9, num; | |
223 xprintf("\n\n\nMatrix is of order %d and has %d off-diagonal non-" | |
224 "zeros\n", n, ipp); | |
225 for (j = 1; j <= n; j++) | |
226 { for (i = 1; i <= n; i++) | |
227 a[i][j] = 0; | |
228 a[j][j] = 1; | |
229 } | |
230 for (k9 = 1; k9 <= ipp; k9++) | |
231 { /* these statements should be replaced by calls to your | |
232 favorite random number generator to place two pseudo-random | |
233 numbers between 1 and n in the variables i and j */ | |
234 for (;;) | |
235 { fa01bs(n, &i); | |
236 fa01bs(n, &j); | |
237 if (!a[i][j]) break; | |
238 } | |
239 a[i][j] = 1; | |
240 } | |
241 /* setup converts matrix a[i,j] to required sparsity-oriented | |
242 storage format */ | |
243 setup(n, a, ip, icn, lenr); | |
244 num = mc13d(n, icn, ip, lenr, ior, ib, &iw[0], &iw[n], &iw[n+n]); | |
245 /* output reordered matrix with blocking to improve clarity */ | |
246 xprintf("\nThe reordered matrix which has %d block%s is of the fo" | |
247 "rm\n", num, num == 1 ? "" : "s"); | |
248 ib[num+1] = n + 1; | |
249 index = 100; | |
250 iblock = 1; | |
251 for (i = 1; i <= n; i++) | |
252 { for (ij = 1; ij <= index; ij++) | |
253 hold[ij] = ' '; | |
254 if (i == ib[iblock]) | |
255 { xprintf("\n"); | |
256 iblock++; | |
257 } | |
258 jblock = 1; | |
259 index = 0; | |
260 for (j = 1; j <= n; j++) | |
261 { if (j == ib[jblock]) | |
262 { hold[++index] = ' '; | |
263 jblock++; | |
264 } | |
265 ii = ior[i]; | |
266 jj = ior[j]; | |
267 hold[++index] = (char)(a[ii][jj] ? 'X' : '0'); | |
268 } | |
269 xprintf("%.*s\n", index, &hold[1]); | |
270 } | |
271 xprintf("\nThe starting point for each block is given by\n"); | |
272 for (i = 1; i <= num; i++) | |
273 { if ((i - 1) % 12 == 0) xprintf("\n"); | |
274 xprintf(" %4d", ib[i]); | |
275 } | |
276 xprintf("\n"); | |
277 return; | |
278 } | |
279 | |
280 void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[]) | |
281 { int i, j, ind; | |
282 for (i = 1; i <= n; i++) | |
283 lenr[i] = 0; | |
284 ind = 1; | |
285 for (i = 1; i <= n; i++) | |
286 { ip[i] = ind; | |
287 for (j = 1; j <= n; j++) | |
288 { if (a[i][j]) | |
289 { lenr[i]++; | |
290 icn[ind++] = j; | |
291 } | |
292 } | |
293 } | |
294 return; | |
295 } | |
296 | |
297 double g = 1431655765.0; | |
298 | |
299 double fa01as(int i) | |
300 { /* random number generator */ | |
301 g = fmod(g * 9228907.0, 4294967296.0); | |
302 if (i >= 0) | |
303 return g / 4294967296.0; | |
304 else | |
305 return 2.0 * g / 4294967296.0 - 1.0; | |
306 } | |
307 | |
308 void fa01bs(int max, int *nrand) | |
309 { *nrand = (int)(fa01as(1) * (double)max) + 1; | |
310 return; | |
311 } | |
312 #endif | |
313 | |
314 /* eof */ |