1 /*
2 * Copyright 1997, Regents of the University of Minnesota
3 *
4 * smbfactor.c
5 *
6 * This file performs the symbolic factorization of a matrix
7 *
8 * Started 8/1/97
9 * George
10 *
11 * $Id: smbfactor.c 10154 2011-06-09 21:27:35Z karypis $
12 *
13 */
14
15 #include "metisbin.h"
16
17
18 /*************************************************************************/
19 /*! This function sets up data structures for fill-in computations */
20 /*************************************************************************/
ComputeFillIn(graph_t * graph,idx_t * perm,idx_t * iperm,size_t * r_maxlnz,size_t * r_opc)21 void ComputeFillIn(graph_t *graph, idx_t *perm, idx_t *iperm,
22 size_t *r_maxlnz, size_t *r_opc)
23 {
24 idx_t i, j, k, nvtxs, maxlnz, maxsub;
25 idx_t *xadj, *adjncy;
26 idx_t *xlnz, *xnzsub, *nzsub;
27 size_t opc;
28
29 /*
30 printf("\nSymbolic factorization... --------------------------------------------\n");
31 */
32
33 nvtxs = graph->nvtxs;
34 xadj = graph->xadj;
35 adjncy = graph->adjncy;
36
37 maxsub = 8*(nvtxs+xadj[nvtxs]);
38
39 /* Relabel the vertices so that it starts from 1 */
40 for (i=0; i<xadj[nvtxs]; i++)
41 adjncy[i]++;
42 for (i=0; i<nvtxs+1; i++)
43 xadj[i]++;
44 for (i=0; i<nvtxs; i++) {
45 iperm[i]++;
46 perm[i]++;
47 }
48
49 /* Allocate the required memory */
50 xlnz = imalloc(nvtxs+2, "ComputeFillIn: xlnz");
51 xnzsub = imalloc(nvtxs+2, "ComputeFillIn: xnzsub");
52 nzsub = imalloc(maxsub+1, "ComputeFillIn: nzsub");
53
54
55 /* Call sparspak's routine. */
56 if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) {
57 printf("Realocating nzsub...\n");
58 gk_free((void **)&nzsub, LTERM);
59
60 maxsub *= 2;
61 nzsub = imalloc(maxsub+1, "ComputeFillIn: nzsub");
62 if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub))
63 errexit("MAXSUB is too small!");
64 }
65
66 for (i=0; i<nvtxs; i++)
67 xlnz[i]--;
68 for (opc=0, i=0; i<nvtxs; i++)
69 opc += (xlnz[i+1]-xlnz[i])*(xlnz[i+1]-xlnz[i]) - (xlnz[i+1]-xlnz[i]);
70
71 *r_maxlnz = maxlnz;
72 *r_opc = opc;
73
74 gk_free((void **)&xlnz, &xnzsub, &nzsub, LTERM);
75
76 /* Relabel the vertices so that it starts from 0 */
77 for (i=0; i<nvtxs; i++) {
78 iperm[i]--;
79 perm[i]--;
80 }
81 for (i=0; i<nvtxs+1; i++)
82 xadj[i]--;
83 for (i=0; i<xadj[nvtxs]; i++)
84 adjncy[i]--;
85
86 }
87
88
89
90 /*************************************************************************/
91 /*!
92 PURPOSE - THIS ROUTINE PERFORMS SYMBOLIC FACTORIZATION
93 ON A PERMUTED LINEAR SYSTEM AND IT ALSO SETS UP THE
94 COMPRESSED DATA STRUCTURE FOR THE SYSTEM.
95
96 INPUT PARAMETERS -
97 NEQNS - NUMBER OF EQUATIONS.
98 (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE.
99 (PERM, INVP) - THE PERMUTATION VECTOR AND ITS INVERSE.
100
101 UPDATED PARAMETERS -
102 MAXSUB - SIZE OF THE SUBSCRIPT ARRAY NZSUB. ON RETURN,
103 IT CONTAINS THE NUMBER OF SUBSCRIPTS USED
104
105 OUTPUT PARAMETERS -
106 XLNZ - INDEX INTO THE NONZERO STORAGE VECTOR LNZ.
107 (XNZSUB, NZSUB) - THE COMPRESSED SUBSCRIPT VECTORS.
108 MAXLNZ - THE NUMBER OF NONZEROS FOUND.
109 */
110 /*************************************************************************/
smbfct(idx_t neqns,idx_t * xadj,idx_t * adjncy,idx_t * perm,idx_t * invp,idx_t * xlnz,idx_t * maxlnz,idx_t * xnzsub,idx_t * nzsub,idx_t * maxsub)111 idx_t smbfct(idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t *perm, idx_t *invp,
112 idx_t *xlnz, idx_t *maxlnz, idx_t *xnzsub, idx_t *nzsub,
113 idx_t *maxsub)
114 {
115 /* Local variables */
116 idx_t node, rchm, mrgk, lmax, i, j, k, m, nabor, nzbeg, nzend;
117 idx_t kxsub, jstop, jstrt, mrkflg, inz, knz, flag;
118 idx_t *mrglnk, *marker, *rchlnk;
119
120 rchlnk = ismalloc(neqns+1, 0, "smbfct: rchlnk");
121 marker = ismalloc(neqns+1, 0, "smbfct: marker");
122 mrglnk = ismalloc(neqns+1, 0, "smbfct: mgrlnk");
123
124 /* Parameter adjustments */
125 --marker;
126 --mrglnk;
127 --rchlnk;
128 --nzsub;
129 --xnzsub;
130 --xlnz;
131 --invp;
132 --perm;
133 --adjncy;
134 --xadj;
135
136 /* Function Body */
137 flag = 0;
138 nzbeg = 1;
139 nzend = 0;
140 xlnz[1] = 1;
141
142 /* FOR EACH COLUMN KNZ COUNTS THE NUMBER OF NONZEROS IN COLUMN K ACCUMULATED IN RCHLNK. */
143 for (k=1; k<=neqns; k++) {
144 xnzsub[k] = nzend;
145 node = perm[k];
146 knz = 0;
147 mrgk = mrglnk[k];
148 mrkflg = 0;
149 marker[k] = k;
150 if (mrgk != 0) {
151 assert(mrgk > 0 && mrgk <= neqns);
152 marker[k] = marker[mrgk];
153 }
154
155 if (xadj[node] >= xadj[node+1]) {
156 xlnz[k+1] = xlnz[k];
157 continue;
158 }
159
160 /* USE RCHLNK TO LINK THROUGH THE STRUCTURE OF A(*,K) BELOW DIAGONAL */
161 assert(k <= neqns && k > 0);
162 rchlnk[k] = neqns+1;
163 for (j=xadj[node]; j<xadj[node+1]; j++) {
164 nabor = invp[adjncy[j]];
165 if (nabor <= k)
166 continue;
167 rchm = k;
168
169 do {
170 m = rchm;
171 assert(m > 0 && m <= neqns);
172 rchm = rchlnk[m];
173 } while (rchm <= nabor);
174
175 knz++;
176 assert(m > 0 && m <= neqns);
177 rchlnk[m] = nabor;
178 assert(nabor > 0 && nabor <= neqns);
179 rchlnk[nabor] = rchm;
180 assert(k > 0 && k <= neqns);
181 if (marker[nabor] != marker[k])
182 mrkflg = 1;
183 }
184
185
186 /* TEST FOR MASS SYMBOLIC ELIMINATION */
187 lmax = 0;
188 assert(mrgk >= 0 && mrgk <= neqns);
189 if (mrkflg != 0 || mrgk == 0 || mrglnk[mrgk] != 0)
190 goto L350;
191 xnzsub[k] = xnzsub[mrgk] + 1;
192 knz = xlnz[mrgk + 1] - (xlnz[mrgk] + 1);
193 goto L1400;
194
195
196 L350:
197 /* LINK THROUGH EACH COLUMN I THAT AFFECTS L(*,K) */
198 i = k;
199 assert(i > 0 && i <= neqns);
200 while ((i = mrglnk[i]) != 0) {
201 assert(i > 0 && i <= neqns);
202 inz = xlnz[i+1] - (xlnz[i]+1);
203 jstrt = xnzsub[i] + 1;
204 jstop = xnzsub[i] + inz;
205
206 if (inz > lmax) {
207 lmax = inz;
208 xnzsub[k] = jstrt;
209 }
210
211 /* MERGE STRUCTURE OF L(*,I) IN NZSUB INTO RCHLNK. */
212 rchm = k;
213 for (j=jstrt; j<=jstop; j++) {
214 nabor = nzsub[j];
215 do {
216 m = rchm;
217 assert(m > 0 && m <= neqns);
218 rchm = rchlnk[m];
219 } while (rchm < nabor);
220
221 if (rchm != nabor) {
222 knz++;
223 assert(m > 0 && m <= neqns);
224 rchlnk[m] = nabor;
225 assert(nabor > 0 && nabor <= neqns);
226 rchlnk[nabor] = rchm;
227 rchm = nabor;
228 }
229 }
230 }
231
232
233 /* CHECK IF SUBSCRIPTS DUPLICATE THOSE OF ANOTHER COLUMN */
234 if (knz == lmax)
235 goto L1400;
236
237 /* OR IF TAIL OF K-1ST COLUMN MATCHES HEAD OF KTH */
238 if (nzbeg > nzend)
239 goto L1200;
240
241 assert(k > 0 && k <= neqns);
242 i = rchlnk[k];
243 for (jstrt = nzbeg; jstrt <= nzend; ++jstrt) {
244 if (nzsub[jstrt] < i)
245 continue;
246
247 if (nzsub[jstrt] == i)
248 goto L1000;
249 else
250 goto L1200;
251 }
252 goto L1200;
253
254
255 L1000:
256 xnzsub[k] = jstrt;
257 for (j = jstrt; j <= nzend; ++j) {
258 if (nzsub[j] != i)
259 goto L1200;
260
261 assert(i > 0 && i <= neqns);
262 i = rchlnk[i];
263 if (i > neqns)
264 goto L1400;
265 }
266 nzend = jstrt - 1;
267
268
269 /* COPY THE STRUCTURE OF L(*,K) FROM RCHLNK TO THE DATA STRUCTURE (XNZSUB, NZSUB) */
270 L1200:
271 nzbeg = nzend + 1;
272 nzend += knz;
273
274 if (nzend >= *maxsub) {
275 flag = 1; /* Out of memory */
276 break;
277 }
278
279 i = k;
280 for (j=nzbeg; j<=nzend; j++) {
281 assert(i > 0 && i <= neqns);
282 i = rchlnk[i];
283 nzsub[j] = i;
284 assert(i > 0 && i <= neqns);
285 marker[i] = k;
286 }
287 xnzsub[k] = nzbeg;
288 assert(k > 0 && k <= neqns);
289 marker[k] = k;
290
291 /*
292 * UPDATE THE VECTOR MRGLNK. NOTE COLUMN L(*,K) JUST FOUND
293 * IS REQUIRED TO DETERMINE COLUMN L(*,J), WHERE
294 * L(J,K) IS THE FIRST NONZERO IN L(*,K) BELOW DIAGONAL.
295 */
296 L1400:
297 if (knz > 1) {
298 kxsub = xnzsub[k];
299 i = nzsub[kxsub];
300 assert(i > 0 && i <= neqns);
301 assert(k > 0 && k <= neqns);
302 mrglnk[k] = mrglnk[i];
303 mrglnk[i] = k;
304 }
305
306 xlnz[k + 1] = xlnz[k] + knz;
307 }
308
309 if (flag == 0) {
310 *maxlnz = xlnz[neqns] - 1;
311 *maxsub = xnzsub[neqns];
312 xnzsub[neqns + 1] = xnzsub[neqns];
313 }
314
315
316 marker++;
317 mrglnk++;
318 rchlnk++;
319 nzsub++;
320 xnzsub++;
321 xlnz++;
322 invp++;
323 perm++;
324 adjncy++;
325 xadj++;
326
327 gk_free((void **)&rchlnk, &mrglnk, &marker, LTERM);
328
329 return flag;
330
331 }
332
333