lemon-project-template-glpk
comparison deps/glpk/src/glpnet01.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:2522d5fcbfb8 |
---|---|
1 /* glpnet01.c (permutations for zero-free diagonal) */ | |
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 * MC21A and MC21B associated with the following paper: | |
8 * | |
9 * I.S.Duff, Algorithm 575: Permutations for zero-free diagonal, ACM | |
10 * Trans. on Math. Softw. 7 (1981), 387-390. | |
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 * mc21a - permutations for zero-free diagonal | |
37 * | |
38 * SYNOPSIS | |
39 * | |
40 * #include "glpnet.h" | |
41 * int mc21a(int n, const int icn[], const int ip[], const int lenr[], | |
42 * int iperm[], int pr[], int arp[], int cv[], int out[]); | |
43 * | |
44 * DESCRIPTION | |
45 * | |
46 * Given the pattern of nonzeros of a sparse matrix, the routine mc21a | |
47 * attempts to find a permutation of its rows that makes the matrix have | |
48 * no zeros on its diagonal. | |
49 * | |
50 * INPUT PARAMETERS | |
51 * | |
52 * n order of 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 PARAMETER | |
65 * | |
66 * iperm contains permutation to make diagonal have the smallest | |
67 * number of zeros on it. Elements (iperm[i], i), i = 1,2,...,n, | |
68 * are non-zero at the end of the algorithm unless the matrix is | |
69 * structurally singular. In this case, (iperm[i], i) will be | |
70 * zero for n - numnz entries. | |
71 * | |
72 * WORKING ARRAYS | |
73 * | |
74 * pr working array of length [1+n], where pr[0] is not used. | |
75 * pr[i] is the previous row to i in the depth first search. | |
76 * | |
77 * arp working array of length [1+n], where arp[0] is not used. | |
78 * arp[i] is one less than the number of non-zeros in row i which | |
79 * have not been scanned when looking for a cheap assignment. | |
80 * | |
81 * cv working array of length [1+n], where cv[0] is not used. | |
82 * cv[i] is the most recent row extension at which column i was | |
83 * visited. | |
84 * | |
85 * out working array of length [1+n], where out[0] is not used. | |
86 * out[i] is one less than the number of non-zeros in row i | |
87 * which have not been scanned during one pass through the main | |
88 * loop. | |
89 * | |
90 * RETURNS | |
91 * | |
92 * The routine mc21a returns numnz, the number of non-zeros on diagonal | |
93 * of permuted matrix. */ | |
94 | |
95 int mc21a(int n, const int icn[], const int ip[], const int lenr[], | |
96 int iperm[], int pr[], int arp[], int cv[], int out[]) | |
97 { int i, ii, in1, in2, j, j1, jord, k, kk, numnz; | |
98 /* Initialization of arrays. */ | |
99 for (i = 1; i <= n; i++) | |
100 { arp[i] = lenr[i] - 1; | |
101 cv[i] = iperm[i] = 0; | |
102 } | |
103 numnz = 0; | |
104 /* Main loop. */ | |
105 /* Each pass round this loop either results in a new assignment | |
106 or gives a row with no assignment. */ | |
107 for (jord = 1; jord <= n; jord++) | |
108 { j = jord; | |
109 pr[j] = -1; | |
110 for (k = 1; k <= jord; k++) | |
111 { /* Look for a cheap assignment. */ | |
112 in1 = arp[j]; | |
113 if (in1 >= 0) | |
114 { in2 = ip[j] + lenr[j] - 1; | |
115 in1 = in2 - in1; | |
116 for (ii = in1; ii <= in2; ii++) | |
117 { i = icn[ii]; | |
118 if (iperm[i] == 0) goto L110; | |
119 } | |
120 /* No cheap assignment in row. */ | |
121 arp[j] = -1; | |
122 } | |
123 /* Begin looking for assignment chain starting with row j.*/ | |
124 out[j] = lenr[j] - 1; | |
125 /* Inner loop. Extends chain by one or backtracks. */ | |
126 for (kk = 1; kk <= jord; kk++) | |
127 { in1 = out[j]; | |
128 if (in1 >= 0) | |
129 { in2 = ip[j] + lenr[j] - 1; | |
130 in1 = in2 - in1; | |
131 /* Forward scan. */ | |
132 for (ii = in1; ii <= in2; ii++) | |
133 { i = icn[ii]; | |
134 if (cv[i] != jord) | |
135 { /* Column i has not yet been accessed during | |
136 this pass. */ | |
137 j1 = j; | |
138 j = iperm[i]; | |
139 cv[i] = jord; | |
140 pr[j] = j1; | |
141 out[j1] = in2 - ii - 1; | |
142 goto L100; | |
143 } | |
144 } | |
145 } | |
146 /* Backtracking step. */ | |
147 j = pr[j]; | |
148 if (j == -1) goto L130; | |
149 } | |
150 L100: ; | |
151 } | |
152 L110: /* New assignment is made. */ | |
153 iperm[i] = j; | |
154 arp[j] = in2 - ii - 1; | |
155 numnz++; | |
156 for (k = 1; k <= jord; k++) | |
157 { j = pr[j]; | |
158 if (j == -1) break; | |
159 ii = ip[j] + lenr[j] - out[j] - 2; | |
160 i = icn[ii]; | |
161 iperm[i] = j; | |
162 } | |
163 L130: ; | |
164 } | |
165 /* If matrix is structurally singular, we now complete the | |
166 permutation iperm. */ | |
167 if (numnz < n) | |
168 { for (i = 1; i <= n; i++) | |
169 arp[i] = 0; | |
170 k = 0; | |
171 for (i = 1; i <= n; i++) | |
172 { if (iperm[i] == 0) | |
173 out[++k] = i; | |
174 else | |
175 arp[iperm[i]] = i; | |
176 } | |
177 k = 0; | |
178 for (i = 1; i <= n; i++) | |
179 { if (arp[i] == 0) | |
180 iperm[out[++k]] = i; | |
181 } | |
182 } | |
183 return numnz; | |
184 } | |
185 | |
186 /**********************************************************************/ | |
187 | |
188 #if 0 | |
189 #include "glplib.h" | |
190 | |
191 int sing; | |
192 | |
193 void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum, | |
194 int iw[]); | |
195 | |
196 void fa01bs(int max, int *nrand); | |
197 | |
198 int main(void) | |
199 { /* test program for the routine mc21a */ | |
200 /* these runs on random matrices cause all possible statements in | |
201 mc21a to be executed */ | |
202 int i, iold, j, j1, j2, jj, knum, l, licn, n, nov4, num, numnz; | |
203 int ip[1+21], icn[1+1000], iperm[1+20], lenr[1+20], iw1[1+80]; | |
204 licn = 1000; | |
205 /* run on random matrices of orders 1 through 20 */ | |
206 for (n = 1; n <= 20; n++) | |
207 { nov4 = n / 4; | |
208 if (nov4 < 1) nov4 = 1; | |
209 L10: fa01bs(nov4, &l); | |
210 knum = l * n; | |
211 /* knum is requested number of non-zeros in random matrix */ | |
212 if (knum > licn) goto L10; | |
213 /* if sing is false, matrix is guaranteed structurally | |
214 non-singular */ | |
215 sing = ((n / 2) * 2 == n); | |
216 /* call to subroutine to generate random matrix */ | |
217 ranmat(n, n, icn, ip, n+1, &knum, iw1); | |
218 /* knum is now actual number of non-zeros in random matrix */ | |
219 if (knum > licn) goto L10; | |
220 xprintf("n = %2d; nz = %4d; sing = %d\n", n, knum, sing); | |
221 /* set up array of row lengths */ | |
222 for (i = 1; i <= n; i++) | |
223 lenr[i] = ip[i+1] - ip[i]; | |
224 /* call to mc21a */ | |
225 numnz = mc21a(n, icn, ip, lenr, iperm, &iw1[0], &iw1[n], | |
226 &iw1[n+n], &iw1[n+n+n]); | |
227 /* testing to see if there are numnz non-zeros on the diagonal | |
228 of the permuted matrix. */ | |
229 num = 0; | |
230 for (i = 1; i <= n; i++) | |
231 { iold = iperm[i]; | |
232 j1 = ip[iold]; | |
233 j2 = j1 + lenr[iold] - 1; | |
234 if (j2 < j1) continue; | |
235 for (jj = j1; jj <= j2; jj++) | |
236 { j = icn[jj]; | |
237 if (j == i) | |
238 { num++; | |
239 break; | |
240 } | |
241 } | |
242 } | |
243 if (num != numnz) | |
244 xprintf("Failure in mc21a, numnz = %d instead of %d\n", | |
245 numnz, num); | |
246 } | |
247 return 0; | |
248 } | |
249 | |
250 void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum, | |
251 int iw[]) | |
252 { /* subroutine to generate random matrix */ | |
253 int i, ii, inum, j, lrow, matnum; | |
254 inum = (*knum / n) * 2; | |
255 if (inum > n-1) inum = n-1; | |
256 matnum = 1; | |
257 /* each pass through this loop generates a row of the matrix */ | |
258 for (j = 1; j <= m; j++) | |
259 { iptr[j] = matnum; | |
260 if (!(sing || j > n)) | |
261 icn[matnum++] = j; | |
262 if (n == 1) continue; | |
263 for (i = 1; i <= n; i++) iw[i] = 0; | |
264 if (!sing) iw[j] = 1; | |
265 fa01bs(inum, &lrow); | |
266 lrow--; | |
267 if (lrow == 0) continue; | |
268 /* lrow off-diagonal non-zeros in row j of the matrix */ | |
269 for (ii = 1; ii <= lrow; ii++) | |
270 { for (;;) | |
271 { fa01bs(n, &i); | |
272 if (iw[i] != 1) break; | |
273 } | |
274 iw[i] = 1; | |
275 icn[matnum++] = i; | |
276 } | |
277 } | |
278 for (i = m+1; i <= nnnp1; i++) | |
279 iptr[i] = matnum; | |
280 *knum = matnum - 1; | |
281 return; | |
282 } | |
283 | |
284 double g = 1431655765.0; | |
285 | |
286 double fa01as(int i) | |
287 { /* random number generator */ | |
288 g = fmod(g * 9228907.0, 4294967296.0); | |
289 if (i >= 0) | |
290 return g / 4294967296.0; | |
291 else | |
292 return 2.0 * g / 4294967296.0 - 1.0; | |
293 } | |
294 | |
295 void fa01bs(int max, int *nrand) | |
296 { *nrand = (int)(fa01as(1) * (double)max) + 1; | |
297 return; | |
298 } | |
299 #endif | |
300 | |
301 /* eof */ |