lemon-project-template-glpk
view deps/glpk/src/glpnet01.c @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 21:43:29 +0100 |
| parents | |
| children |
line source
1 /* glpnet01.c (permutations for zero-free diagonal) */
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 ***********************************************************************/
31 #include "glpnet.h"
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. */
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 }
186 /**********************************************************************/
188 #if 0
189 #include "glplib.h"
191 int sing;
193 void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum,
194 int iw[]);
196 void fa01bs(int max, int *nrand);
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 }
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 }
284 double g = 1431655765.0;
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 }
295 void fa01bs(int max, int *nrand)
296 { *nrand = (int)(fa01as(1) * (double)max) + 1;
297 return;
298 }
299 #endif
301 /* eof */