lemon-project-template-glpk
view 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 |
line source
1 /* glpnet02.c (permutations to block triangular form) */
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 ***********************************************************************/
31 #include "glpnet.h"
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. */
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 }
190 /**********************************************************************/
192 #if 0
193 #include "glplib.h"
195 void test(int n, int ipp);
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 }
214 void fa01bs(int max, int *nrand);
216 void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[]);
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 }
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 }
297 double g = 1431655765.0;
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 }
308 void fa01bs(int max, int *nrand)
309 { *nrand = (int)(fa01as(1) * (double)max) + 1;
310 return;
311 }
312 #endif
314 /* eof */