1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
|
/* mc21a.c (permutations for zero-free diagonal) */
/***********************************************************************
* This code is part of GLPK (GNU Linear Programming Kit).
*
* This code is the result of translation of the Fortran subroutines
* MC21A and MC21B associated with the following paper:
*
* I.S.Duff, Algorithm 575: Permutations for zero-free diagonal, ACM
* Trans. on Math. Softw. 7 (1981), 387-390.
*
* Use of ACM Algorithms is subject to the ACM Software Copyright and
* License Agreement. See <http://www.acm.org/publications/policies>.
*
* The translation was made by Andrew Makhorin <mao@gnu.org>.
*
* GLPK is free software: you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* GLPK is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
* License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#include "mc21a.h"
/***********************************************************************
* NAME
*
* mc21a - permutations for zero-free diagonal
*
* SYNOPSIS
*
* #include "mc21a.h"
* int mc21a(int n, const int icn[], const int ip[], const int lenr[],
* int iperm[], int pr[], int arp[], int cv[], int out[]);
*
* DESCRIPTION
*
* Given the pattern of nonzeros of a sparse matrix, the routine mc21a
* attempts to find a permutation of its rows that makes the matrix have
* no zeros on its diagonal.
*
* INPUT PARAMETERS
*
* n order of matrix.
*
* icn array containing the column indices of the non-zeros. Those
* belonging to a single row must be contiguous but the ordering
* of column indices within each row is unimportant and wasted
* space between rows is permitted.
*
* ip ip[i], i = 1,2,...,n, is the position in array icn of the
* first column index of a non-zero in row i.
*
* lenr lenr[i], i = 1,2,...,n, is the number of non-zeros in row i.
*
* OUTPUT PARAMETER
*
* iperm contains permutation to make diagonal have the smallest
* number of zeros on it. Elements (iperm[i], i), i = 1,2,...,n,
* are non-zero at the end of the algorithm unless the matrix is
* structurally singular. In this case, (iperm[i], i) will be
* zero for n - numnz entries.
*
* WORKING ARRAYS
*
* pr working array of length [1+n], where pr[0] is not used.
* pr[i] is the previous row to i in the depth first search.
*
* arp working array of length [1+n], where arp[0] is not used.
* arp[i] is one less than the number of non-zeros in row i which
* have not been scanned when looking for a cheap assignment.
*
* cv working array of length [1+n], where cv[0] is not used.
* cv[i] is the most recent row extension at which column i was
* visited.
*
* out working array of length [1+n], where out[0] is not used.
* out[i] is one less than the number of non-zeros in row i
* which have not been scanned during one pass through the main
* loop.
*
* RETURNS
*
* The routine mc21a returns numnz, the number of non-zeros on diagonal
* of permuted matrix. */
int mc21a(int n, const int icn[], const int ip[], const int lenr[],
int iperm[], int pr[], int arp[], int cv[], int out[])
{ int i, ii, in1, in2, j, j1, jord, k, kk, numnz;
/* Initialization of arrays. */
for (i = 1; i <= n; i++)
{ arp[i] = lenr[i] - 1;
cv[i] = iperm[i] = 0;
}
numnz = 0;
/* Main loop. */
/* Each pass round this loop either results in a new assignment
* or gives a row with no assignment. */
for (jord = 1; jord <= n; jord++)
{ j = jord;
pr[j] = -1;
for (k = 1; k <= jord; k++)
{ /* Look for a cheap assignment. */
in1 = arp[j];
if (in1 >= 0)
{ in2 = ip[j] + lenr[j] - 1;
in1 = in2 - in1;
for (ii = in1; ii <= in2; ii++)
{ i = icn[ii];
if (iperm[i] == 0) goto L110;
}
/* No cheap assignment in row. */
arp[j] = -1;
}
/* Begin looking for assignment chain starting with row j.*/
out[j] = lenr[j] - 1;
/* Inner loop. Extends chain by one or backtracks. */
for (kk = 1; kk <= jord; kk++)
{ in1 = out[j];
if (in1 >= 0)
{ in2 = ip[j] + lenr[j] - 1;
in1 = in2 - in1;
/* Forward scan. */
for (ii = in1; ii <= in2; ii++)
{ i = icn[ii];
if (cv[i] != jord)
{ /* Column i has not yet been accessed during
* this pass. */
j1 = j;
j = iperm[i];
cv[i] = jord;
pr[j] = j1;
out[j1] = in2 - ii - 1;
goto L100;
}
}
}
/* Backtracking step. */
j = pr[j];
if (j == -1) goto L130;
}
L100: ;
}
L110: /* New assignment is made. */
iperm[i] = j;
arp[j] = in2 - ii - 1;
numnz++;
for (k = 1; k <= jord; k++)
{ j = pr[j];
if (j == -1) break;
ii = ip[j] + lenr[j] - out[j] - 2;
i = icn[ii];
iperm[i] = j;
}
L130: ;
}
/* If matrix is structurally singular, we now complete the
* permutation iperm. */
if (numnz < n)
{ for (i = 1; i <= n; i++)
arp[i] = 0;
k = 0;
for (i = 1; i <= n; i++)
{ if (iperm[i] == 0)
out[++k] = i;
else
arp[iperm[i]] = i;
}
k = 0;
for (i = 1; i <= n; i++)
{ if (arp[i] == 0)
iperm[out[++k]] = i;
}
}
return numnz;
}
/**********************************************************************/
#ifdef GLP_TEST
#include "env.h"
int sing;
void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum,
int iw[]);
void fa01bs(int max, int *nrand);
int main(void)
{ /* test program for the routine mc21a */
/* these runs on random matrices cause all possible statements in
* mc21a to be executed */
int i, iold, j, j1, j2, jj, knum, l, licn, n, nov4, num, numnz;
int ip[1+21], icn[1+1000], iperm[1+20], lenr[1+20], iw1[1+80];
licn = 1000;
/* run on random matrices of orders 1 through 20 */
for (n = 1; n <= 20; n++)
{ nov4 = n / 4;
if (nov4 < 1) nov4 = 1;
L10: fa01bs(nov4, &l);
knum = l * n;
/* knum is requested number of non-zeros in random matrix */
if (knum > licn) goto L10;
/* if sing is false, matrix is guaranteed structurally
* non-singular */
sing = ((n / 2) * 2 == n);
/* call to subroutine to generate random matrix */
ranmat(n, n, icn, ip, n+1, &knum, iw1);
/* knum is now actual number of non-zeros in random matrix */
if (knum > licn) goto L10;
xprintf("n = %2d; nz = %4d; sing = %d\n", n, knum, sing);
/* set up array of row lengths */
for (i = 1; i <= n; i++)
lenr[i] = ip[i+1] - ip[i];
/* call to mc21a */
numnz = mc21a(n, icn, ip, lenr, iperm, &iw1[0], &iw1[n],
&iw1[n+n], &iw1[n+n+n]);
/* testing to see if there are numnz non-zeros on the diagonal
* of the permuted matrix. */
num = 0;
for (i = 1; i <= n; i++)
{ iold = iperm[i];
j1 = ip[iold];
j2 = j1 + lenr[iold] - 1;
if (j2 < j1) continue;
for (jj = j1; jj <= j2; jj++)
{ j = icn[jj];
if (j == i)
{ num++;
break;
}
}
}
if (num != numnz)
xprintf("Failure in mc21a, numnz = %d instead of %d\n",
numnz, num);
}
return 0;
}
void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum,
int iw[])
{ /* subroutine to generate random matrix */
int i, ii, inum, j, lrow, matnum;
inum = (*knum / n) * 2;
if (inum > n-1) inum = n-1;
matnum = 1;
/* each pass through this loop generates a row of the matrix */
for (j = 1; j <= m; j++)
{ iptr[j] = matnum;
if (!(sing || j > n))
icn[matnum++] = j;
if (n == 1) continue;
for (i = 1; i <= n; i++) iw[i] = 0;
if (!sing) iw[j] = 1;
fa01bs(inum, &lrow);
lrow--;
if (lrow == 0) continue;
/* lrow off-diagonal non-zeros in row j of the matrix */
for (ii = 1; ii <= lrow; ii++)
{ for (;;)
{ fa01bs(n, &i);
if (iw[i] != 1) break;
}
iw[i] = 1;
icn[matnum++] = i;
}
}
for (i = m+1; i <= nnnp1; i++)
iptr[i] = matnum;
*knum = matnum - 1;
return;
}
double g = 1431655765.0;
double fa01as(int i)
{ /* random number generator */
g = fmod(g * 9228907.0, 4294967296.0);
if (i >= 0)
return g / 4294967296.0;
else
return 2.0 * g / 4294967296.0 - 1.0;
}
void fa01bs(int max, int *nrand)
{ *nrand = (int)(fa01as(1) * (double)max) + 1;
return;
}
#endif
/* eof */
|
