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
302
303
304
305
306
307
308
309
310
311
|
/* triang.c (find maximal triangular part of rectangular matrix) */
/***********************************************************************
* This code is part of GLPK (GNU Linear Programming Kit).
*
* Copyright (C) 2012-2013 Andrew Makhorin, Department for Applied
* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
* reserved. E-mail: <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 "env.h"
#include "triang.h"
/***********************************************************************
* triang - find maximal triangular part of rectangular matrix
*
* Given a mxn sparse matrix A this routine finds permutation matrices
* P and Q such that matrix A' = P * A * Q has the following structure:
*
* 1 s n
* 1 * . . . . . x x x x x
* * * . . . . x x x x x
* * * * . . . x x x x x
* * * * * . . x x x x x
* * * * * * . x x x x x
* s * * * * * * x x x x x
* x x x x x x x x x x x
* x x x x x x x x x x x
* m x x x x x x x x x x x
*
* where '*' are elements of the triangular part, '.' are structural
* zeros, 'x' are other elements.
*
* The formal routine mat specifies the original matrix A in both row-
* and column-wise format. If the routine mat is called with k = +i,
* 1 <= i <= m, it should store column indices and values of non-zero
* elements of i-th row of A in locations ind[1], ..., ind[len] and
* val[1], ..., val[len], resp., where len is the returned number of
* non-zeros in the row, 0 <= len <= n. Similarly, if the routine mat
* is called with k = -j, 1 <= j <= n, it should store row indices and
* values of non-zero elements of j-th column of A and return len, the
* number of non-zeros in the column, 0 <= len <= m. Should note that
* duplicate indices are not allowed.
*
* The parameter info is a transit pointer passed to the routine mat.
*
* The parameter tol is a tolerance. The routine triang guarantees that
* each diagonal element in the triangular part of matrix A' is not
* less in magnitude than tol * max, where max is the maximal magnitude
* of elements in corresponding column.
*
* On exit the routine triang stores information on the triangular part
* found in the arrays rn and cn. Elements rn[1], ..., rn[s] specify
* row numbers and elements cn[1], ..., cn[s] specify column numbers
* of the original matrix A, which correspond to rows/columns 1, ..., s
* of matrix A', where s is the size of the triangular part returned by
* the routine, 0 <= s <= min(m, n). The order of rows and columns that
* are not included in the triangular part remains unspecified.
*
* ALGORITHM
*
* The routine triang uses a simple greedy heuristic.
*
* At some step the matrix A' = P * A * Q has the following structure:
*
* 1 n
* 1 * . . . . . . . x x x
* * * . . . . . . x x x
* * * * . . . . . x x x
* * * * * . . . . x x x
* x x x x # # # # x x x
* x x x x # # # # x x x
* x x x x # # # # x x x
* x x x x # # # # x x x
* m x x x x # # # # x x x
*
* where '#' are elements of active submatrix. Initially P = Q = I, so
* the active submatrix is the original matrix A = A'.
*
* If some row has exactly one non-zero in the active submatrix (row
* singleton), the routine includes this row and corresponding column
* in the triangular part, and removes the column from the active
* submatrix. Otherwise, the routine simply removes a column having
* maximal number of non-zeros from the active submatrix in the hope
* that new row singleton(s) will appear.
*
* COMPLEXITY
*
* The time complexity of the routine triang is O(nnz), where nnz is
* number of non-zeros in the original matrix A. */
int triang(int m, int n, int (*mat)(void *info, int k, int ind[],
double val[]), void *info, double tol, int rn[], int cn[])
{ int head, i, j, jj, k, kk, ks, len, len2, next_j, ns, size;
int *cind, *rind, *cnt, *ptr, *list, *prev, *next;
double *cval, *rval, *big;
char *flag;
/* allocate working arrays */
cind = talloc(1+m, int);
cval = talloc(1+m, double);
rind = talloc(1+n, int);
rval = talloc(1+n, double);
cnt = ptr = talloc(1+m, int);
list = talloc(1+n, int);
prev = talloc(1+n, int);
next = talloc(1+n, int);
big = talloc(1+n, double);
flag = talloc(1+n, char);
/*--------------------------------------------------------------*/
/* build linked lists of columns having equal lengths */
/*--------------------------------------------------------------*/
/* ptr[len], 0 <= len <= m, is number of first column of length
* len;
* next[j], 1 <= j <= n, is number of next column having the same
* length as column j;
* big[j], 1 <= j <= n, is maximal magnitude of elements in j-th
* column */
for (len = 0; len <= m; len++)
ptr[len] = 0;
for (j = 1; j <= n; j++)
{ /* get j-th column */
len = mat(info, -j, cind, cval);
xassert(0 <= len && len <= m);
/* add this column to beginning of list ptr[len] */
next[j] = ptr[len];
ptr[len] = j;
/* determine maximal magnitude of elements in this column */
big[j] = 0.0;
for (k = 1; k <= len; k++)
{ if (big[j] < fabs(cval[k]))
big[j] = fabs(cval[k]);
}
}
/*--------------------------------------------------------------*/
/* build doubly linked list of columns ordered by decreasing */
/* column lengths */
/*--------------------------------------------------------------*/
/* head is number of first column in the list;
* prev[j], 1 <= j <= n, is number of column that precedes j-th
* column in the list;
* next[j], 1 <= j <= n, is number of column that follows j-th
* column in the list */
head = 0;
for (len = 0; len <= m; len++)
{ /* walk thru list of columns of length len */
for (j = ptr[len]; j != 0; j = next_j)
{ next_j = next[j];
/* add j-th column to beginning of the column list */
prev[j] = 0;
next[j] = head;
if (head != 0)
prev[head] = j;
head = j;
}
}
/*--------------------------------------------------------------*/
/* build initial singleton list */
/*--------------------------------------------------------------*/
/* there are used two list of columns:
* 1) doubly linked list of active columns, in which all columns
* are ordered by decreasing column lengths;
* 2) singleton list; an active column is included in this list
* if it has at least one row singleton in active submatrix */
/* flag[j], 1 <= j <= n, is a flag of j-th column:
* 0 j-th column is inactive;
* 1 j-th column is active;
* 2 j-th column is active and has row singleton(s) */
/* initially all columns are active */
for (j = 1; j <= n; j++)
flag[j] = 1;
/* initialize row counts and build initial singleton list */
/* cnt[i], 1 <= i <= m, is number of non-zeros, which i-th row
* has in active submatrix;
* ns is size of singleton list;
* list[1], ..., list[ns] are numbers of active columns included
* in the singleton list */
ns = 0;
for (i = 1; i <= m; i++)
{ /* get i-th row */
len = cnt[i] = mat(info, +i, rind, rval);
xassert(0 <= len && len <= n);
if (len == 1)
{ /* a[i,j] is row singleton */
j = rind[1];
xassert(1 <= j && j <= n);
if (flag[j] != 2)
{ /* include j-th column in singleton list */
flag[j] = 2;
list[++ns] = j;
}
}
}
/*--------------------------------------------------------------*/
/* main loop */
/*--------------------------------------------------------------*/
size = 0; /* size of triangular part */
/* loop until active column list is non-empty, i.e. until the
* active submatrix has at least one column */
while (head != 0)
{ if (ns == 0)
{ /* singleton list is empty */
/* remove from the active submatrix a column of maximal
* length in the hope that some row singletons appear */
j = head;
len = mat(info, -j, cind, cval);
xassert(0 <= len && len <= m);
goto drop;
}
/* take column j from the singleton list */
j = list[ns--];
xassert(flag[j] == 2);
/* j-th column has at least one row singleton in the active
* submatrix; choose one having maximal magnitude */
len = mat(info, -j, cind, cval);
xassert(0 <= len && len <= m);
kk = 0;
for (k = 1; k <= len; k++)
{ i = cind[k];
xassert(1 <= i && i <= m);
if (cnt[i] == 1)
{ /* a[i,j] is row singleton */
if (kk == 0 || fabs(cval[kk]) < fabs(cval[k]))
kk = k;
}
}
xassert(kk > 0);
/* check magnitude of the row singleton chosen */
if (fabs(cval[kk]) < tol * big[j])
{ /* all row singletons are too small in magnitude; drop j-th
* column */
goto drop;
}
/* row singleton a[i,j] is ok; add i-th row and j-th column to
* the triangular part */
size++;
rn[size] = cind[kk];
cn[size] = j;
drop: /* remove j-th column from the active submatrix */
xassert(flag[j]);
flag[j] = 0;
if (prev[j] == 0)
head = next[j];
else
next[prev[j]] = next[j];
if (next[j] == 0)
;
else
prev[next[j]] = prev[j];
/* decrease row counts */
for (k = 1; k <= len; k++)
{ i = cind[k];
xassert(1 <= i && i <= m);
xassert(cnt[i] > 0);
cnt[i]--;
if (cnt[i] == 1)
{ /* new singleton appeared in i-th row; determine number
* of corresponding column (it is the only active column
* in this row) */
len2 = mat(info, +i, rind, rval);
xassert(0 <= len2 && len2 <= n);
ks = 0;
for (kk = 1; kk <= len2; kk++)
{ jj = rind[kk];
xassert(1 <= jj && jj <= n);
if (flag[jj])
{ xassert(ks == 0);
ks = kk;
}
}
xassert(ks > 0);
/* a[i,jj] is new row singleton */
jj = rind[ks];
if (flag[jj] != 2)
{ /* include jj-th column in the singleton list */
flag[jj] = 2;
list[++ns] = jj;
}
}
}
}
/* now all row counts should be zero */
for (i = 1; i <= m; i++)
xassert(cnt[i] == 0);
/* deallocate working arrays */
tfree(cind);
tfree(cval);
tfree(rind);
tfree(rval);
tfree(ptr);
tfree(list);
tfree(prev);
tfree(next);
tfree(big);
tfree(flag);
return size;
}
/* eof */
|
