lemon-project-template-glpk
comparison deps/glpk/src/amd/amd_2.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 |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:434f0c381400 |
---|---|
1 /* ========================================================================= */ | |
2 /* === AMD_2 =============================================================== */ | |
3 /* ========================================================================= */ | |
4 | |
5 /* ------------------------------------------------------------------------- */ | |
6 /* AMD, Copyright (c) Timothy A. Davis, */ | |
7 /* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */ | |
8 /* email: davis at cise.ufl.edu CISE Department, Univ. of Florida. */ | |
9 /* web: http://www.cise.ufl.edu/research/sparse/amd */ | |
10 /* ------------------------------------------------------------------------- */ | |
11 | |
12 /* AMD_2: performs the AMD ordering on a symmetric sparse matrix A, followed | |
13 * by a postordering (via depth-first search) of the assembly tree using the | |
14 * AMD_postorder routine. | |
15 */ | |
16 | |
17 #include "amd_internal.h" | |
18 | |
19 /* ========================================================================= */ | |
20 /* === clear_flag ========================================================== */ | |
21 /* ========================================================================= */ | |
22 | |
23 static Int clear_flag (Int wflg, Int wbig, Int W [ ], Int n) | |
24 { | |
25 Int x ; | |
26 if (wflg < 2 || wflg >= wbig) | |
27 { | |
28 for (x = 0 ; x < n ; x++) | |
29 { | |
30 if (W [x] != 0) W [x] = 1 ; | |
31 } | |
32 wflg = 2 ; | |
33 } | |
34 /* at this point, W [0..n-1] < wflg holds */ | |
35 return (wflg) ; | |
36 } | |
37 | |
38 | |
39 /* ========================================================================= */ | |
40 /* === AMD_2 =============================================================== */ | |
41 /* ========================================================================= */ | |
42 | |
43 GLOBAL void AMD_2 | |
44 ( | |
45 Int n, /* A is n-by-n, where n > 0 */ | |
46 Int Pe [ ], /* Pe [0..n-1]: index in Iw of row i on input */ | |
47 Int Iw [ ], /* workspace of size iwlen. Iw [0..pfree-1] | |
48 * holds the matrix on input */ | |
49 Int Len [ ], /* Len [0..n-1]: length for row/column i on input */ | |
50 Int iwlen, /* length of Iw. iwlen >= pfree + n */ | |
51 Int pfree, /* Iw [pfree ... iwlen-1] is empty on input */ | |
52 | |
53 /* 7 size-n workspaces, not defined on input: */ | |
54 Int Nv [ ], /* the size of each supernode on output */ | |
55 Int Next [ ], /* the output inverse permutation */ | |
56 Int Last [ ], /* the output permutation */ | |
57 Int Head [ ], | |
58 Int Elen [ ], /* the size columns of L for each supernode */ | |
59 Int Degree [ ], | |
60 Int W [ ], | |
61 | |
62 /* control parameters and output statistics */ | |
63 double Control [ ], /* array of size AMD_CONTROL */ | |
64 double Info [ ] /* array of size AMD_INFO */ | |
65 ) | |
66 { | |
67 | |
68 /* | |
69 * Given a representation of the nonzero pattern of a symmetric matrix, A, | |
70 * (excluding the diagonal) perform an approximate minimum (UMFPACK/MA38-style) | |
71 * degree ordering to compute a pivot order such that the introduction of | |
72 * nonzeros (fill-in) in the Cholesky factors A = LL' is kept low. At each | |
73 * step, the pivot selected is the one with the minimum UMFAPACK/MA38-style | |
74 * upper-bound on the external degree. This routine can optionally perform | |
75 * aggresive absorption (as done by MC47B in the Harwell Subroutine | |
76 * Library). | |
77 * | |
78 * The approximate degree algorithm implemented here is the symmetric analog of | |
79 * the degree update algorithm in MA38 and UMFPACK (the Unsymmetric-pattern | |
80 * MultiFrontal PACKage, both by Davis and Duff). The routine is based on the | |
81 * MA27 minimum degree ordering algorithm by Iain Duff and John Reid. | |
82 * | |
83 * This routine is a translation of the original AMDBAR and MC47B routines, | |
84 * in Fortran, with the following modifications: | |
85 * | |
86 * (1) dense rows/columns are removed prior to ordering the matrix, and placed | |
87 * last in the output order. The presence of a dense row/column can | |
88 * increase the ordering time by up to O(n^2), unless they are removed | |
89 * prior to ordering. | |
90 * | |
91 * (2) the minimum degree ordering is followed by a postordering (depth-first | |
92 * search) of the assembly tree. Note that mass elimination (discussed | |
93 * below) combined with the approximate degree update can lead to the mass | |
94 * elimination of nodes with lower exact degree than the current pivot | |
95 * element. No additional fill-in is caused in the representation of the | |
96 * Schur complement. The mass-eliminated nodes merge with the current | |
97 * pivot element. They are ordered prior to the current pivot element. | |
98 * Because they can have lower exact degree than the current element, the | |
99 * merger of two or more of these nodes in the current pivot element can | |
100 * lead to a single element that is not a "fundamental supernode". The | |
101 * diagonal block can have zeros in it. Thus, the assembly tree used here | |
102 * is not guaranteed to be the precise supernodal elemination tree (with | |
103 * "funadmental" supernodes), and the postordering performed by this | |
104 * routine is not guaranteed to be a precise postordering of the | |
105 * elimination tree. | |
106 * | |
107 * (3) input parameters are added, to control aggressive absorption and the | |
108 * detection of "dense" rows/columns of A. | |
109 * | |
110 * (4) additional statistical information is returned, such as the number of | |
111 * nonzeros in L, and the flop counts for subsequent LDL' and LU | |
112 * factorizations. These are slight upper bounds, because of the mass | |
113 * elimination issue discussed above. | |
114 * | |
115 * (5) additional routines are added to interface this routine to MATLAB | |
116 * to provide a simple C-callable user-interface, to check inputs for | |
117 * errors, compute the symmetry of the pattern of A and the number of | |
118 * nonzeros in each row/column of A+A', to compute the pattern of A+A', | |
119 * to perform the assembly tree postordering, and to provide debugging | |
120 * ouput. Many of these functions are also provided by the Fortran | |
121 * Harwell Subroutine Library routine MC47A. | |
122 * | |
123 * (6) both int and UF_long versions are provided. In the descriptions below | |
124 * and integer is and int or UF_long depending on which version is | |
125 * being used. | |
126 | |
127 ********************************************************************** | |
128 ***** CAUTION: ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT. ****** | |
129 ********************************************************************** | |
130 ** If you want error checking, a more versatile input format, and a ** | |
131 ** simpler user interface, use amd_order or amd_l_order instead. ** | |
132 ** This routine is not meant to be user-callable. ** | |
133 ********************************************************************** | |
134 | |
135 * ---------------------------------------------------------------------------- | |
136 * References: | |
137 * ---------------------------------------------------------------------------- | |
138 * | |
139 * [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern multifrontal | |
140 * method for sparse LU factorization", SIAM J. Matrix Analysis and | |
141 * Applications, vol. 18, no. 1, pp. 140-158. Discusses UMFPACK / MA38, | |
142 * which first introduced the approximate minimum degree used by this | |
143 * routine. | |
144 * | |
145 * [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, "An approximate | |
146 * minimum degree ordering algorithm," SIAM J. Matrix Analysis and | |
147 * Applications, vol. 17, no. 4, pp. 886-905, 1996. Discusses AMDBAR and | |
148 * MC47B, which are the Fortran versions of this routine. | |
149 * | |
150 * [3] Alan George and Joseph Liu, "The evolution of the minimum degree | |
151 * ordering algorithm," SIAM Review, vol. 31, no. 1, pp. 1-19, 1989. | |
152 * We list below the features mentioned in that paper that this code | |
153 * includes: | |
154 * | |
155 * mass elimination: | |
156 * Yes. MA27 relied on supervariable detection for mass elimination. | |
157 * | |
158 * indistinguishable nodes: | |
159 * Yes (we call these "supervariables"). This was also in the MA27 | |
160 * code - although we modified the method of detecting them (the | |
161 * previous hash was the true degree, which we no longer keep track | |
162 * of). A supervariable is a set of rows with identical nonzero | |
163 * pattern. All variables in a supervariable are eliminated together. | |
164 * Each supervariable has as its numerical name that of one of its | |
165 * variables (its principal variable). | |
166 * | |
167 * quotient graph representation: | |
168 * Yes. We use the term "element" for the cliques formed during | |
169 * elimination. This was also in the MA27 code. The algorithm can | |
170 * operate in place, but it will work more efficiently if given some | |
171 * "elbow room." | |
172 * | |
173 * element absorption: | |
174 * Yes. This was also in the MA27 code. | |
175 * | |
176 * external degree: | |
177 * Yes. The MA27 code was based on the true degree. | |
178 * | |
179 * incomplete degree update and multiple elimination: | |
180 * No. This was not in MA27, either. Our method of degree update | |
181 * within MC47B is element-based, not variable-based. It is thus | |
182 * not well-suited for use with incomplete degree update or multiple | |
183 * elimination. | |
184 * | |
185 * Authors, and Copyright (C) 2004 by: | |
186 * Timothy A. Davis, Patrick Amestoy, Iain S. Duff, John K. Reid. | |
187 * | |
188 * Acknowledgements: This work (and the UMFPACK package) was supported by the | |
189 * National Science Foundation (ASC-9111263, DMS-9223088, and CCR-0203270). | |
190 * The UMFPACK/MA38 approximate degree update algorithm, the unsymmetric analog | |
191 * which forms the basis of AMD, was developed while Tim Davis was supported by | |
192 * CERFACS (Toulouse, France) in a post-doctoral position. This C version, and | |
193 * the etree postorder, were written while Tim Davis was on sabbatical at | |
194 * Stanford University and Lawrence Berkeley National Laboratory. | |
195 | |
196 * ---------------------------------------------------------------------------- | |
197 * INPUT ARGUMENTS (unaltered): | |
198 * ---------------------------------------------------------------------------- | |
199 | |
200 * n: The matrix order. Restriction: n >= 1. | |
201 * | |
202 * iwlen: The size of the Iw array. On input, the matrix is stored in | |
203 * Iw [0..pfree-1]. However, Iw [0..iwlen-1] should be slightly larger | |
204 * than what is required to hold the matrix, at least iwlen >= pfree + n. | |
205 * Otherwise, excessive compressions will take place. The recommended | |
206 * value of iwlen is 1.2 * pfree + n, which is the value used in the | |
207 * user-callable interface to this routine (amd_order.c). The algorithm | |
208 * will not run at all if iwlen < pfree. Restriction: iwlen >= pfree + n. | |
209 * Note that this is slightly more restrictive than the actual minimum | |
210 * (iwlen >= pfree), but AMD_2 will be very slow with no elbow room. | |
211 * Thus, this routine enforces a bare minimum elbow room of size n. | |
212 * | |
213 * pfree: On input the tail end of the array, Iw [pfree..iwlen-1], is empty, | |
214 * and the matrix is stored in Iw [0..pfree-1]. During execution, | |
215 * additional data is placed in Iw, and pfree is modified so that | |
216 * Iw [pfree..iwlen-1] is always the unused part of Iw. | |
217 * | |
218 * Control: A double array of size AMD_CONTROL containing input parameters | |
219 * that affect how the ordering is computed. If NULL, then default | |
220 * settings are used. | |
221 * | |
222 * Control [AMD_DENSE] is used to determine whether or not a given input | |
223 * row is "dense". A row is "dense" if the number of entries in the row | |
224 * exceeds Control [AMD_DENSE] times sqrt (n), except that rows with 16 or | |
225 * fewer entries are never considered "dense". To turn off the detection | |
226 * of dense rows, set Control [AMD_DENSE] to a negative number, or to a | |
227 * number larger than sqrt (n). The default value of Control [AMD_DENSE] | |
228 * is AMD_DEFAULT_DENSE, which is defined in amd.h as 10. | |
229 * | |
230 * Control [AMD_AGGRESSIVE] is used to determine whether or not aggressive | |
231 * absorption is to be performed. If nonzero, then aggressive absorption | |
232 * is performed (this is the default). | |
233 | |
234 * ---------------------------------------------------------------------------- | |
235 * INPUT/OUPUT ARGUMENTS: | |
236 * ---------------------------------------------------------------------------- | |
237 * | |
238 * Pe: An integer array of size n. On input, Pe [i] is the index in Iw of | |
239 * the start of row i. Pe [i] is ignored if row i has no off-diagonal | |
240 * entries. Thus Pe [i] must be in the range 0 to pfree-1 for non-empty | |
241 * rows. | |
242 * | |
243 * During execution, it is used for both supervariables and elements: | |
244 * | |
245 * Principal supervariable i: index into Iw of the description of | |
246 * supervariable i. A supervariable represents one or more rows of | |
247 * the matrix with identical nonzero pattern. In this case, | |
248 * Pe [i] >= 0. | |
249 * | |
250 * Non-principal supervariable i: if i has been absorbed into another | |
251 * supervariable j, then Pe [i] = FLIP (j), where FLIP (j) is defined | |
252 * as (-(j)-2). Row j has the same pattern as row i. Note that j | |
253 * might later be absorbed into another supervariable j2, in which | |
254 * case Pe [i] is still FLIP (j), and Pe [j] = FLIP (j2) which is | |
255 * < EMPTY, where EMPTY is defined as (-1) in amd_internal.h. | |
256 * | |
257 * Unabsorbed element e: the index into Iw of the description of element | |
258 * e, if e has not yet been absorbed by a subsequent element. Element | |
259 * e is created when the supervariable of the same name is selected as | |
260 * the pivot. In this case, Pe [i] >= 0. | |
261 * | |
262 * Absorbed element e: if element e is absorbed into element e2, then | |
263 * Pe [e] = FLIP (e2). This occurs when the pattern of e (which we | |
264 * refer to as Le) is found to be a subset of the pattern of e2 (that | |
265 * is, Le2). In this case, Pe [i] < EMPTY. If element e is "null" | |
266 * (it has no nonzeros outside its pivot block), then Pe [e] = EMPTY, | |
267 * and e is the root of an assembly subtree (or the whole tree if | |
268 * there is just one such root). | |
269 * | |
270 * Dense variable i: if i is "dense", then Pe [i] = EMPTY. | |
271 * | |
272 * On output, Pe holds the assembly tree/forest, which implicitly | |
273 * represents a pivot order with identical fill-in as the actual order | |
274 * (via a depth-first search of the tree), as follows. If Nv [i] > 0, | |
275 * then i represents a node in the assembly tree, and the parent of i is | |
276 * Pe [i], or EMPTY if i is a root. If Nv [i] = 0, then (i, Pe [i]) | |
277 * represents an edge in a subtree, the root of which is a node in the | |
278 * assembly tree. Note that i refers to a row/column in the original | |
279 * matrix, not the permuted matrix. | |
280 * | |
281 * Info: A double array of size AMD_INFO. If present, (that is, not NULL), | |
282 * then statistics about the ordering are returned in the Info array. | |
283 * See amd.h for a description. | |
284 | |
285 * ---------------------------------------------------------------------------- | |
286 * INPUT/MODIFIED (undefined on output): | |
287 * ---------------------------------------------------------------------------- | |
288 * | |
289 * Len: An integer array of size n. On input, Len [i] holds the number of | |
290 * entries in row i of the matrix, excluding the diagonal. The contents | |
291 * of Len are undefined on output. | |
292 * | |
293 * Iw: An integer array of size iwlen. On input, Iw [0..pfree-1] holds the | |
294 * description of each row i in the matrix. The matrix must be symmetric, | |
295 * and both upper and lower triangular parts must be present. The | |
296 * diagonal must not be present. Row i is held as follows: | |
297 * | |
298 * Len [i]: the length of the row i data structure in the Iw array. | |
299 * Iw [Pe [i] ... Pe [i] + Len [i] - 1]: | |
300 * the list of column indices for nonzeros in row i (simple | |
301 * supervariables), excluding the diagonal. All supervariables | |
302 * start with one row/column each (supervariable i is just row i). | |
303 * If Len [i] is zero on input, then Pe [i] is ignored on input. | |
304 * | |
305 * Note that the rows need not be in any particular order, and there | |
306 * may be empty space between the rows. | |
307 * | |
308 * During execution, the supervariable i experiences fill-in. This is | |
309 * represented by placing in i a list of the elements that cause fill-in | |
310 * in supervariable i: | |
311 * | |
312 * Len [i]: the length of supervariable i in the Iw array. | |
313 * Iw [Pe [i] ... Pe [i] + Elen [i] - 1]: | |
314 * the list of elements that contain i. This list is kept short | |
315 * by removing absorbed elements. | |
316 * Iw [Pe [i] + Elen [i] ... Pe [i] + Len [i] - 1]: | |
317 * the list of supervariables in i. This list is kept short by | |
318 * removing nonprincipal variables, and any entry j that is also | |
319 * contained in at least one of the elements (j in Le) in the list | |
320 * for i (e in row i). | |
321 * | |
322 * When supervariable i is selected as pivot, we create an element e of | |
323 * the same name (e=i): | |
324 * | |
325 * Len [e]: the length of element e in the Iw array. | |
326 * Iw [Pe [e] ... Pe [e] + Len [e] - 1]: | |
327 * the list of supervariables in element e. | |
328 * | |
329 * An element represents the fill-in that occurs when supervariable i is | |
330 * selected as pivot (which represents the selection of row i and all | |
331 * non-principal variables whose principal variable is i). We use the | |
332 * term Le to denote the set of all supervariables in element e. Absorbed | |
333 * supervariables and elements are pruned from these lists when | |
334 * computationally convenient. | |
335 * | |
336 * CAUTION: THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION. | |
337 * The contents of Iw are undefined on output. | |
338 | |
339 * ---------------------------------------------------------------------------- | |
340 * OUTPUT (need not be set on input): | |
341 * ---------------------------------------------------------------------------- | |
342 * | |
343 * Nv: An integer array of size n. During execution, ABS (Nv [i]) is equal to | |
344 * the number of rows that are represented by the principal supervariable | |
345 * i. If i is a nonprincipal or dense variable, then Nv [i] = 0. | |
346 * Initially, Nv [i] = 1 for all i. Nv [i] < 0 signifies that i is a | |
347 * principal variable in the pattern Lme of the current pivot element me. | |
348 * After element me is constructed, Nv [i] is set back to a positive | |
349 * value. | |
350 * | |
351 * On output, Nv [i] holds the number of pivots represented by super | |
352 * row/column i of the original matrix, or Nv [i] = 0 for non-principal | |
353 * rows/columns. Note that i refers to a row/column in the original | |
354 * matrix, not the permuted matrix. | |
355 * | |
356 * Elen: An integer array of size n. See the description of Iw above. At the | |
357 * start of execution, Elen [i] is set to zero for all rows i. During | |
358 * execution, Elen [i] is the number of elements in the list for | |
359 * supervariable i. When e becomes an element, Elen [e] = FLIP (esize) is | |
360 * set, where esize is the size of the element (the number of pivots, plus | |
361 * the number of nonpivotal entries). Thus Elen [e] < EMPTY. | |
362 * Elen (i) = EMPTY set when variable i becomes nonprincipal. | |
363 * | |
364 * For variables, Elen (i) >= EMPTY holds until just before the | |
365 * postordering and permutation vectors are computed. For elements, | |
366 * Elen [e] < EMPTY holds. | |
367 * | |
368 * On output, Elen [i] is the degree of the row/column in the Cholesky | |
369 * factorization of the permuted matrix, corresponding to the original row | |
370 * i, if i is a super row/column. It is equal to EMPTY if i is | |
371 * non-principal. Note that i refers to a row/column in the original | |
372 * matrix, not the permuted matrix. | |
373 * | |
374 * Note that the contents of Elen on output differ from the Fortran | |
375 * version (Elen holds the inverse permutation in the Fortran version, | |
376 * which is instead returned in the Next array in this C version, | |
377 * described below). | |
378 * | |
379 * Last: In a degree list, Last [i] is the supervariable preceding i, or EMPTY | |
380 * if i is the head of the list. In a hash bucket, Last [i] is the hash | |
381 * key for i. | |
382 * | |
383 * Last [Head [hash]] is also used as the head of a hash bucket if | |
384 * Head [hash] contains a degree list (see the description of Head, | |
385 * below). | |
386 * | |
387 * On output, Last [0..n-1] holds the permutation. That is, if | |
388 * i = Last [k], then row i is the kth pivot row (where k ranges from 0 to | |
389 * n-1). Row Last [k] of A is the kth row in the permuted matrix, PAP'. | |
390 * | |
391 * Next: Next [i] is the supervariable following i in a link list, or EMPTY if | |
392 * i is the last in the list. Used for two kinds of lists: degree lists | |
393 * and hash buckets (a supervariable can be in only one kind of list at a | |
394 * time). | |
395 * | |
396 * On output Next [0..n-1] holds the inverse permutation. That is, if | |
397 * k = Next [i], then row i is the kth pivot row. Row i of A appears as | |
398 * the (Next[i])-th row in the permuted matrix, PAP'. | |
399 * | |
400 * Note that the contents of Next on output differ from the Fortran | |
401 * version (Next is undefined on output in the Fortran version). | |
402 | |
403 * ---------------------------------------------------------------------------- | |
404 * LOCAL WORKSPACE (not input or output - used only during execution): | |
405 * ---------------------------------------------------------------------------- | |
406 * | |
407 * Degree: An integer array of size n. If i is a supervariable, then | |
408 * Degree [i] holds the current approximation of the external degree of | |
409 * row i (an upper bound). The external degree is the number of nonzeros | |
410 * in row i, minus ABS (Nv [i]), the diagonal part. The bound is equal to | |
411 * the exact external degree if Elen [i] is less than or equal to two. | |
412 * | |
413 * We also use the term "external degree" for elements e to refer to | |
414 * |Le \ Lme|. If e is an element, then Degree [e] is |Le|, which is the | |
415 * degree of the off-diagonal part of the element e (not including the | |
416 * diagonal part). | |
417 * | |
418 * Head: An integer array of size n. Head is used for degree lists. | |
419 * Head [deg] is the first supervariable in a degree list. All | |
420 * supervariables i in a degree list Head [deg] have the same approximate | |
421 * degree, namely, deg = Degree [i]. If the list Head [deg] is empty then | |
422 * Head [deg] = EMPTY. | |
423 * | |
424 * During supervariable detection Head [hash] also serves as a pointer to | |
425 * a hash bucket. If Head [hash] >= 0, there is a degree list of degree | |
426 * hash. The hash bucket head pointer is Last [Head [hash]]. If | |
427 * Head [hash] = EMPTY, then the degree list and hash bucket are both | |
428 * empty. If Head [hash] < EMPTY, then the degree list is empty, and | |
429 * FLIP (Head [hash]) is the head of the hash bucket. After supervariable | |
430 * detection is complete, all hash buckets are empty, and the | |
431 * (Last [Head [hash]] = EMPTY) condition is restored for the non-empty | |
432 * degree lists. | |
433 * | |
434 * W: An integer array of size n. The flag array W determines the status of | |
435 * elements and variables, and the external degree of elements. | |
436 * | |
437 * for elements: | |
438 * if W [e] = 0, then the element e is absorbed. | |
439 * if W [e] >= wflg, then W [e] - wflg is the size of the set | |
440 * |Le \ Lme|, in terms of nonzeros (the sum of ABS (Nv [i]) for | |
441 * each principal variable i that is both in the pattern of | |
442 * element e and NOT in the pattern of the current pivot element, | |
443 * me). | |
444 * if wflg > W [e] > 0, then e is not absorbed and has not yet been | |
445 * seen in the scan of the element lists in the computation of | |
446 * |Le\Lme| in Scan 1 below. | |
447 * | |
448 * for variables: | |
449 * during supervariable detection, if W [j] != wflg then j is | |
450 * not in the pattern of variable i. | |
451 * | |
452 * The W array is initialized by setting W [i] = 1 for all i, and by | |
453 * setting wflg = 2. It is reinitialized if wflg becomes too large (to | |
454 * ensure that wflg+n does not cause integer overflow). | |
455 | |
456 * ---------------------------------------------------------------------------- | |
457 * LOCAL INTEGERS: | |
458 * ---------------------------------------------------------------------------- | |
459 */ | |
460 | |
461 Int deg, degme, dext, lemax, e, elenme, eln, i, ilast, inext, j, | |
462 jlast, jnext, k, knt1, knt2, knt3, lenj, ln, me, mindeg, nel, nleft, | |
463 nvi, nvj, nvpiv, slenme, wbig, we, wflg, wnvi, ok, ndense, ncmpa, | |
464 dense, aggressive ; | |
465 | |
466 unsigned Int hash ; /* unsigned, so that hash % n is well defined.*/ | |
467 | |
468 /* | |
469 * deg: the degree of a variable or element | |
470 * degme: size, |Lme|, of the current element, me (= Degree [me]) | |
471 * dext: external degree, |Le \ Lme|, of some element e | |
472 * lemax: largest |Le| seen so far (called dmax in Fortran version) | |
473 * e: an element | |
474 * elenme: the length, Elen [me], of element list of pivotal variable | |
475 * eln: the length, Elen [...], of an element list | |
476 * hash: the computed value of the hash function | |
477 * i: a supervariable | |
478 * ilast: the entry in a link list preceding i | |
479 * inext: the entry in a link list following i | |
480 * j: a supervariable | |
481 * jlast: the entry in a link list preceding j | |
482 * jnext: the entry in a link list, or path, following j | |
483 * k: the pivot order of an element or variable | |
484 * knt1: loop counter used during element construction | |
485 * knt2: loop counter used during element construction | |
486 * knt3: loop counter used during compression | |
487 * lenj: Len [j] | |
488 * ln: length of a supervariable list | |
489 * me: current supervariable being eliminated, and the current | |
490 * element created by eliminating that supervariable | |
491 * mindeg: current minimum degree | |
492 * nel: number of pivots selected so far | |
493 * nleft: n - nel, the number of nonpivotal rows/columns remaining | |
494 * nvi: the number of variables in a supervariable i (= Nv [i]) | |
495 * nvj: the number of variables in a supervariable j (= Nv [j]) | |
496 * nvpiv: number of pivots in current element | |
497 * slenme: number of variables in variable list of pivotal variable | |
498 * wbig: = INT_MAX - n for the int version, UF_long_max - n for the | |
499 * UF_long version. wflg is not allowed to be >= wbig. | |
500 * we: W [e] | |
501 * wflg: used for flagging the W array. See description of Iw. | |
502 * wnvi: wflg - Nv [i] | |
503 * x: either a supervariable or an element | |
504 * | |
505 * ok: true if supervariable j can be absorbed into i | |
506 * ndense: number of "dense" rows/columns | |
507 * dense: rows/columns with initial degree > dense are considered "dense" | |
508 * aggressive: true if aggressive absorption is being performed | |
509 * ncmpa: number of garbage collections | |
510 | |
511 * ---------------------------------------------------------------------------- | |
512 * LOCAL DOUBLES, used for statistical output only (except for alpha): | |
513 * ---------------------------------------------------------------------------- | |
514 */ | |
515 | |
516 double f, r, ndiv, s, nms_lu, nms_ldl, dmax, alpha, lnz, lnzme ; | |
517 | |
518 /* | |
519 * f: nvpiv | |
520 * r: degme + nvpiv | |
521 * ndiv: number of divisions for LU or LDL' factorizations | |
522 * s: number of multiply-subtract pairs for LU factorization, for the | |
523 * current element me | |
524 * nms_lu number of multiply-subtract pairs for LU factorization | |
525 * nms_ldl number of multiply-subtract pairs for LDL' factorization | |
526 * dmax: the largest number of entries in any column of L, including the | |
527 * diagonal | |
528 * alpha: "dense" degree ratio | |
529 * lnz: the number of nonzeros in L (excluding the diagonal) | |
530 * lnzme: the number of nonzeros in L (excl. the diagonal) for the | |
531 * current element me | |
532 | |
533 * ---------------------------------------------------------------------------- | |
534 * LOCAL "POINTERS" (indices into the Iw array) | |
535 * ---------------------------------------------------------------------------- | |
536 */ | |
537 | |
538 Int p, p1, p2, p3, p4, pdst, pend, pj, pme, pme1, pme2, pn, psrc ; | |
539 | |
540 /* | |
541 * Any parameter (Pe [...] or pfree) or local variable starting with "p" (for | |
542 * Pointer) is an index into Iw, and all indices into Iw use variables starting | |
543 * with "p." The only exception to this rule is the iwlen input argument. | |
544 * | |
545 * p: pointer into lots of things | |
546 * p1: Pe [i] for some variable i (start of element list) | |
547 * p2: Pe [i] + Elen [i] - 1 for some variable i | |
548 * p3: index of first supervariable in clean list | |
549 * p4: | |
550 * pdst: destination pointer, for compression | |
551 * pend: end of memory to compress | |
552 * pj: pointer into an element or variable | |
553 * pme: pointer into the current element (pme1...pme2) | |
554 * pme1: the current element, me, is stored in Iw [pme1...pme2] | |
555 * pme2: the end of the current element | |
556 * pn: pointer into a "clean" variable, also used to compress | |
557 * psrc: source pointer, for compression | |
558 */ | |
559 | |
560 /* ========================================================================= */ | |
561 /* INITIALIZATIONS */ | |
562 /* ========================================================================= */ | |
563 | |
564 /* Note that this restriction on iwlen is slightly more restrictive than | |
565 * what is actually required in AMD_2. AMD_2 can operate with no elbow | |
566 * room at all, but it will be slow. For better performance, at least | |
567 * size-n elbow room is enforced. */ | |
568 ASSERT (iwlen >= pfree + n) ; | |
569 ASSERT (n > 0) ; | |
570 | |
571 /* initialize output statistics */ | |
572 lnz = 0 ; | |
573 ndiv = 0 ; | |
574 nms_lu = 0 ; | |
575 nms_ldl = 0 ; | |
576 dmax = 1 ; | |
577 me = EMPTY ; | |
578 | |
579 mindeg = 0 ; | |
580 ncmpa = 0 ; | |
581 nel = 0 ; | |
582 lemax = 0 ; | |
583 | |
584 /* get control parameters */ | |
585 if (Control != (double *) NULL) | |
586 { | |
587 alpha = Control [AMD_DENSE] ; | |
588 aggressive = (Control [AMD_AGGRESSIVE] != 0) ; | |
589 } | |
590 else | |
591 { | |
592 alpha = AMD_DEFAULT_DENSE ; | |
593 aggressive = AMD_DEFAULT_AGGRESSIVE ; | |
594 } | |
595 /* Note: if alpha is NaN, this is undefined: */ | |
596 if (alpha < 0) | |
597 { | |
598 /* only remove completely dense rows/columns */ | |
599 dense = n-2 ; | |
600 } | |
601 else | |
602 { | |
603 dense = alpha * sqrt ((double) n) ; | |
604 } | |
605 dense = MAX (16, dense) ; | |
606 dense = MIN (n, dense) ; | |
607 AMD_DEBUG1 (("\n\nAMD (debug), alpha %g, aggr. "ID"\n", | |
608 alpha, aggressive)) ; | |
609 | |
610 for (i = 0 ; i < n ; i++) | |
611 { | |
612 Last [i] = EMPTY ; | |
613 Head [i] = EMPTY ; | |
614 Next [i] = EMPTY ; | |
615 /* if separate Hhead array is used for hash buckets: * | |
616 Hhead [i] = EMPTY ; | |
617 */ | |
618 Nv [i] = 1 ; | |
619 W [i] = 1 ; | |
620 Elen [i] = 0 ; | |
621 Degree [i] = Len [i] ; | |
622 } | |
623 | |
624 #ifndef NDEBUG | |
625 AMD_DEBUG1 (("\n======Nel "ID" initial\n", nel)) ; | |
626 AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next, Last, | |
627 Head, Elen, Degree, W, -1) ; | |
628 #endif | |
629 | |
630 /* initialize wflg */ | |
631 wbig = Int_MAX - n ; | |
632 wflg = clear_flag (0, wbig, W, n) ; | |
633 | |
634 /* --------------------------------------------------------------------- */ | |
635 /* initialize degree lists and eliminate dense and empty rows */ | |
636 /* --------------------------------------------------------------------- */ | |
637 | |
638 ndense = 0 ; | |
639 | |
640 for (i = 0 ; i < n ; i++) | |
641 { | |
642 deg = Degree [i] ; | |
643 ASSERT (deg >= 0 && deg < n) ; | |
644 if (deg == 0) | |
645 { | |
646 | |
647 /* ------------------------------------------------------------- | |
648 * we have a variable that can be eliminated at once because | |
649 * there is no off-diagonal non-zero in its row. Note that | |
650 * Nv [i] = 1 for an empty variable i. It is treated just | |
651 * the same as an eliminated element i. | |
652 * ------------------------------------------------------------- */ | |
653 | |
654 Elen [i] = FLIP (1) ; | |
655 nel++ ; | |
656 Pe [i] = EMPTY ; | |
657 W [i] = 0 ; | |
658 | |
659 } | |
660 else if (deg > dense) | |
661 { | |
662 | |
663 /* ------------------------------------------------------------- | |
664 * Dense variables are not treated as elements, but as unordered, | |
665 * non-principal variables that have no parent. They do not take | |
666 * part in the postorder, since Nv [i] = 0. Note that the Fortran | |
667 * version does not have this option. | |
668 * ------------------------------------------------------------- */ | |
669 | |
670 AMD_DEBUG1 (("Dense node "ID" degree "ID"\n", i, deg)) ; | |
671 ndense++ ; | |
672 Nv [i] = 0 ; /* do not postorder this node */ | |
673 Elen [i] = EMPTY ; | |
674 nel++ ; | |
675 Pe [i] = EMPTY ; | |
676 | |
677 } | |
678 else | |
679 { | |
680 | |
681 /* ------------------------------------------------------------- | |
682 * place i in the degree list corresponding to its degree | |
683 * ------------------------------------------------------------- */ | |
684 | |
685 inext = Head [deg] ; | |
686 ASSERT (inext >= EMPTY && inext < n) ; | |
687 if (inext != EMPTY) Last [inext] = i ; | |
688 Next [i] = inext ; | |
689 Head [deg] = i ; | |
690 | |
691 } | |
692 } | |
693 | |
694 /* ========================================================================= */ | |
695 /* WHILE (selecting pivots) DO */ | |
696 /* ========================================================================= */ | |
697 | |
698 while (nel < n) | |
699 { | |
700 | |
701 #ifndef NDEBUG | |
702 AMD_DEBUG1 (("\n======Nel "ID"\n", nel)) ; | |
703 if (AMD_debug >= 2) | |
704 { | |
705 AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next, | |
706 Last, Head, Elen, Degree, W, nel) ; | |
707 } | |
708 #endif | |
709 | |
710 /* ========================================================================= */ | |
711 /* GET PIVOT OF MINIMUM DEGREE */ | |
712 /* ========================================================================= */ | |
713 | |
714 /* ----------------------------------------------------------------- */ | |
715 /* find next supervariable for elimination */ | |
716 /* ----------------------------------------------------------------- */ | |
717 | |
718 ASSERT (mindeg >= 0 && mindeg < n) ; | |
719 for (deg = mindeg ; deg < n ; deg++) | |
720 { | |
721 me = Head [deg] ; | |
722 if (me != EMPTY) break ; | |
723 } | |
724 mindeg = deg ; | |
725 ASSERT (me >= 0 && me < n) ; | |
726 AMD_DEBUG1 (("=================me: "ID"\n", me)) ; | |
727 | |
728 /* ----------------------------------------------------------------- */ | |
729 /* remove chosen variable from link list */ | |
730 /* ----------------------------------------------------------------- */ | |
731 | |
732 inext = Next [me] ; | |
733 ASSERT (inext >= EMPTY && inext < n) ; | |
734 if (inext != EMPTY) Last [inext] = EMPTY ; | |
735 Head [deg] = inext ; | |
736 | |
737 /* ----------------------------------------------------------------- */ | |
738 /* me represents the elimination of pivots nel to nel+Nv[me]-1. */ | |
739 /* place me itself as the first in this set. */ | |
740 /* ----------------------------------------------------------------- */ | |
741 | |
742 elenme = Elen [me] ; | |
743 nvpiv = Nv [me] ; | |
744 ASSERT (nvpiv > 0) ; | |
745 nel += nvpiv ; | |
746 | |
747 /* ========================================================================= */ | |
748 /* CONSTRUCT NEW ELEMENT */ | |
749 /* ========================================================================= */ | |
750 | |
751 /* ----------------------------------------------------------------- | |
752 * At this point, me is the pivotal supervariable. It will be | |
753 * converted into the current element. Scan list of the pivotal | |
754 * supervariable, me, setting tree pointers and constructing new list | |
755 * of supervariables for the new element, me. p is a pointer to the | |
756 * current position in the old list. | |
757 * ----------------------------------------------------------------- */ | |
758 | |
759 /* flag the variable "me" as being in Lme by negating Nv [me] */ | |
760 Nv [me] = -nvpiv ; | |
761 degme = 0 ; | |
762 ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ; | |
763 | |
764 if (elenme == 0) | |
765 { | |
766 | |
767 /* ------------------------------------------------------------- */ | |
768 /* construct the new element in place */ | |
769 /* ------------------------------------------------------------- */ | |
770 | |
771 pme1 = Pe [me] ; | |
772 pme2 = pme1 - 1 ; | |
773 | |
774 for (p = pme1 ; p <= pme1 + Len [me] - 1 ; p++) | |
775 { | |
776 i = Iw [p] ; | |
777 ASSERT (i >= 0 && i < n && Nv [i] >= 0) ; | |
778 nvi = Nv [i] ; | |
779 if (nvi > 0) | |
780 { | |
781 | |
782 /* ----------------------------------------------------- */ | |
783 /* i is a principal variable not yet placed in Lme. */ | |
784 /* store i in new list */ | |
785 /* ----------------------------------------------------- */ | |
786 | |
787 /* flag i as being in Lme by negating Nv [i] */ | |
788 degme += nvi ; | |
789 Nv [i] = -nvi ; | |
790 Iw [++pme2] = i ; | |
791 | |
792 /* ----------------------------------------------------- */ | |
793 /* remove variable i from degree list. */ | |
794 /* ----------------------------------------------------- */ | |
795 | |
796 ilast = Last [i] ; | |
797 inext = Next [i] ; | |
798 ASSERT (ilast >= EMPTY && ilast < n) ; | |
799 ASSERT (inext >= EMPTY && inext < n) ; | |
800 if (inext != EMPTY) Last [inext] = ilast ; | |
801 if (ilast != EMPTY) | |
802 { | |
803 Next [ilast] = inext ; | |
804 } | |
805 else | |
806 { | |
807 /* i is at the head of the degree list */ | |
808 ASSERT (Degree [i] >= 0 && Degree [i] < n) ; | |
809 Head [Degree [i]] = inext ; | |
810 } | |
811 } | |
812 } | |
813 } | |
814 else | |
815 { | |
816 | |
817 /* ------------------------------------------------------------- */ | |
818 /* construct the new element in empty space, Iw [pfree ...] */ | |
819 /* ------------------------------------------------------------- */ | |
820 | |
821 p = Pe [me] ; | |
822 pme1 = pfree ; | |
823 slenme = Len [me] - elenme ; | |
824 | |
825 for (knt1 = 1 ; knt1 <= elenme + 1 ; knt1++) | |
826 { | |
827 | |
828 if (knt1 > elenme) | |
829 { | |
830 /* search the supervariables in me. */ | |
831 e = me ; | |
832 pj = p ; | |
833 ln = slenme ; | |
834 AMD_DEBUG2 (("Search sv: "ID" "ID" "ID"\n", me,pj,ln)) ; | |
835 } | |
836 else | |
837 { | |
838 /* search the elements in me. */ | |
839 e = Iw [p++] ; | |
840 ASSERT (e >= 0 && e < n) ; | |
841 pj = Pe [e] ; | |
842 ln = Len [e] ; | |
843 AMD_DEBUG2 (("Search element e "ID" in me "ID"\n", e,me)) ; | |
844 ASSERT (Elen [e] < EMPTY && W [e] > 0 && pj >= 0) ; | |
845 } | |
846 ASSERT (ln >= 0 && (ln == 0 || (pj >= 0 && pj < iwlen))) ; | |
847 | |
848 /* --------------------------------------------------------- | |
849 * search for different supervariables and add them to the | |
850 * new list, compressing when necessary. this loop is | |
851 * executed once for each element in the list and once for | |
852 * all the supervariables in the list. | |
853 * --------------------------------------------------------- */ | |
854 | |
855 for (knt2 = 1 ; knt2 <= ln ; knt2++) | |
856 { | |
857 i = Iw [pj++] ; | |
858 ASSERT (i >= 0 && i < n && (i == me || Elen [i] >= EMPTY)); | |
859 nvi = Nv [i] ; | |
860 AMD_DEBUG2 ((": "ID" "ID" "ID" "ID"\n", | |
861 i, Elen [i], Nv [i], wflg)) ; | |
862 | |
863 if (nvi > 0) | |
864 { | |
865 | |
866 /* ------------------------------------------------- */ | |
867 /* compress Iw, if necessary */ | |
868 /* ------------------------------------------------- */ | |
869 | |
870 if (pfree >= iwlen) | |
871 { | |
872 | |
873 AMD_DEBUG1 (("GARBAGE COLLECTION\n")) ; | |
874 | |
875 /* prepare for compressing Iw by adjusting pointers | |
876 * and lengths so that the lists being searched in | |
877 * the inner and outer loops contain only the | |
878 * remaining entries. */ | |
879 | |
880 Pe [me] = p ; | |
881 Len [me] -= knt1 ; | |
882 /* check if nothing left of supervariable me */ | |
883 if (Len [me] == 0) Pe [me] = EMPTY ; | |
884 Pe [e] = pj ; | |
885 Len [e] = ln - knt2 ; | |
886 /* nothing left of element e */ | |
887 if (Len [e] == 0) Pe [e] = EMPTY ; | |
888 | |
889 ncmpa++ ; /* one more garbage collection */ | |
890 | |
891 /* store first entry of each object in Pe */ | |
892 /* FLIP the first entry in each object */ | |
893 for (j = 0 ; j < n ; j++) | |
894 { | |
895 pn = Pe [j] ; | |
896 if (pn >= 0) | |
897 { | |
898 ASSERT (pn >= 0 && pn < iwlen) ; | |
899 Pe [j] = Iw [pn] ; | |
900 Iw [pn] = FLIP (j) ; | |
901 } | |
902 } | |
903 | |
904 /* psrc/pdst point to source/destination */ | |
905 psrc = 0 ; | |
906 pdst = 0 ; | |
907 pend = pme1 - 1 ; | |
908 | |
909 while (psrc <= pend) | |
910 { | |
911 /* search for next FLIP'd entry */ | |
912 j = FLIP (Iw [psrc++]) ; | |
913 if (j >= 0) | |
914 { | |
915 AMD_DEBUG2 (("Got object j: "ID"\n", j)) ; | |
916 Iw [pdst] = Pe [j] ; | |
917 Pe [j] = pdst++ ; | |
918 lenj = Len [j] ; | |
919 /* copy from source to destination */ | |
920 for (knt3 = 0 ; knt3 <= lenj - 2 ; knt3++) | |
921 { | |
922 Iw [pdst++] = Iw [psrc++] ; | |
923 } | |
924 } | |
925 } | |
926 | |
927 /* move the new partially-constructed element */ | |
928 p1 = pdst ; | |
929 for (psrc = pme1 ; psrc <= pfree-1 ; psrc++) | |
930 { | |
931 Iw [pdst++] = Iw [psrc] ; | |
932 } | |
933 pme1 = p1 ; | |
934 pfree = pdst ; | |
935 pj = Pe [e] ; | |
936 p = Pe [me] ; | |
937 | |
938 } | |
939 | |
940 /* ------------------------------------------------- */ | |
941 /* i is a principal variable not yet placed in Lme */ | |
942 /* store i in new list */ | |
943 /* ------------------------------------------------- */ | |
944 | |
945 /* flag i as being in Lme by negating Nv [i] */ | |
946 degme += nvi ; | |
947 Nv [i] = -nvi ; | |
948 Iw [pfree++] = i ; | |
949 AMD_DEBUG2 ((" s: "ID" nv "ID"\n", i, Nv [i])); | |
950 | |
951 /* ------------------------------------------------- */ | |
952 /* remove variable i from degree link list */ | |
953 /* ------------------------------------------------- */ | |
954 | |
955 ilast = Last [i] ; | |
956 inext = Next [i] ; | |
957 ASSERT (ilast >= EMPTY && ilast < n) ; | |
958 ASSERT (inext >= EMPTY && inext < n) ; | |
959 if (inext != EMPTY) Last [inext] = ilast ; | |
960 if (ilast != EMPTY) | |
961 { | |
962 Next [ilast] = inext ; | |
963 } | |
964 else | |
965 { | |
966 /* i is at the head of the degree list */ | |
967 ASSERT (Degree [i] >= 0 && Degree [i] < n) ; | |
968 Head [Degree [i]] = inext ; | |
969 } | |
970 } | |
971 } | |
972 | |
973 if (e != me) | |
974 { | |
975 /* set tree pointer and flag to indicate element e is | |
976 * absorbed into new element me (the parent of e is me) */ | |
977 AMD_DEBUG1 ((" Element "ID" => "ID"\n", e, me)) ; | |
978 Pe [e] = FLIP (me) ; | |
979 W [e] = 0 ; | |
980 } | |
981 } | |
982 | |
983 pme2 = pfree - 1 ; | |
984 } | |
985 | |
986 /* ----------------------------------------------------------------- */ | |
987 /* me has now been converted into an element in Iw [pme1..pme2] */ | |
988 /* ----------------------------------------------------------------- */ | |
989 | |
990 /* degme holds the external degree of new element */ | |
991 Degree [me] = degme ; | |
992 Pe [me] = pme1 ; | |
993 Len [me] = pme2 - pme1 + 1 ; | |
994 ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ; | |
995 | |
996 Elen [me] = FLIP (nvpiv + degme) ; | |
997 /* FLIP (Elen (me)) is now the degree of pivot (including | |
998 * diagonal part). */ | |
999 | |
1000 #ifndef NDEBUG | |
1001 AMD_DEBUG2 (("New element structure: length= "ID"\n", pme2-pme1+1)) ; | |
1002 for (pme = pme1 ; pme <= pme2 ; pme++) AMD_DEBUG3 ((" "ID"", Iw[pme])); | |
1003 AMD_DEBUG3 (("\n")) ; | |
1004 #endif | |
1005 | |
1006 /* ----------------------------------------------------------------- */ | |
1007 /* make sure that wflg is not too large. */ | |
1008 /* ----------------------------------------------------------------- */ | |
1009 | |
1010 /* With the current value of wflg, wflg+n must not cause integer | |
1011 * overflow */ | |
1012 | |
1013 wflg = clear_flag (wflg, wbig, W, n) ; | |
1014 | |
1015 /* ========================================================================= */ | |
1016 /* COMPUTE (W [e] - wflg) = |Le\Lme| FOR ALL ELEMENTS */ | |
1017 /* ========================================================================= */ | |
1018 | |
1019 /* ----------------------------------------------------------------- | |
1020 * Scan 1: compute the external degrees of previous elements with | |
1021 * respect to the current element. That is: | |
1022 * (W [e] - wflg) = |Le \ Lme| | |
1023 * for each element e that appears in any supervariable in Lme. The | |
1024 * notation Le refers to the pattern (list of supervariables) of a | |
1025 * previous element e, where e is not yet absorbed, stored in | |
1026 * Iw [Pe [e] + 1 ... Pe [e] + Len [e]]. The notation Lme | |
1027 * refers to the pattern of the current element (stored in | |
1028 * Iw [pme1..pme2]). If aggressive absorption is enabled, and | |
1029 * (W [e] - wflg) becomes zero, then the element e will be absorbed | |
1030 * in Scan 2. | |
1031 * ----------------------------------------------------------------- */ | |
1032 | |
1033 AMD_DEBUG2 (("me: ")) ; | |
1034 for (pme = pme1 ; pme <= pme2 ; pme++) | |
1035 { | |
1036 i = Iw [pme] ; | |
1037 ASSERT (i >= 0 && i < n) ; | |
1038 eln = Elen [i] ; | |
1039 AMD_DEBUG3 ((""ID" Elen "ID": \n", i, eln)) ; | |
1040 if (eln > 0) | |
1041 { | |
1042 /* note that Nv [i] has been negated to denote i in Lme: */ | |
1043 nvi = -Nv [i] ; | |
1044 ASSERT (nvi > 0 && Pe [i] >= 0 && Pe [i] < iwlen) ; | |
1045 wnvi = wflg - nvi ; | |
1046 for (p = Pe [i] ; p <= Pe [i] + eln - 1 ; p++) | |
1047 { | |
1048 e = Iw [p] ; | |
1049 ASSERT (e >= 0 && e < n) ; | |
1050 we = W [e] ; | |
1051 AMD_DEBUG4 ((" e "ID" we "ID" ", e, we)) ; | |
1052 if (we >= wflg) | |
1053 { | |
1054 /* unabsorbed element e has been seen in this loop */ | |
1055 AMD_DEBUG4 ((" unabsorbed, first time seen")) ; | |
1056 we -= nvi ; | |
1057 } | |
1058 else if (we != 0) | |
1059 { | |
1060 /* e is an unabsorbed element */ | |
1061 /* this is the first we have seen e in all of Scan 1 */ | |
1062 AMD_DEBUG4 ((" unabsorbed")) ; | |
1063 we = Degree [e] + wnvi ; | |
1064 } | |
1065 AMD_DEBUG4 (("\n")) ; | |
1066 W [e] = we ; | |
1067 } | |
1068 } | |
1069 } | |
1070 AMD_DEBUG2 (("\n")) ; | |
1071 | |
1072 /* ========================================================================= */ | |
1073 /* DEGREE UPDATE AND ELEMENT ABSORPTION */ | |
1074 /* ========================================================================= */ | |
1075 | |
1076 /* ----------------------------------------------------------------- | |
1077 * Scan 2: for each i in Lme, sum up the degree of Lme (which is | |
1078 * degme), plus the sum of the external degrees of each Le for the | |
1079 * elements e appearing within i, plus the supervariables in i. | |
1080 * Place i in hash list. | |
1081 * ----------------------------------------------------------------- */ | |
1082 | |
1083 for (pme = pme1 ; pme <= pme2 ; pme++) | |
1084 { | |
1085 i = Iw [pme] ; | |
1086 ASSERT (i >= 0 && i < n && Nv [i] < 0 && Elen [i] >= 0) ; | |
1087 AMD_DEBUG2 (("Updating: i "ID" "ID" "ID"\n", i, Elen[i], Len [i])); | |
1088 p1 = Pe [i] ; | |
1089 p2 = p1 + Elen [i] - 1 ; | |
1090 pn = p1 ; | |
1091 hash = 0 ; | |
1092 deg = 0 ; | |
1093 ASSERT (p1 >= 0 && p1 < iwlen && p2 >= -1 && p2 < iwlen) ; | |
1094 | |
1095 /* ------------------------------------------------------------- */ | |
1096 /* scan the element list associated with supervariable i */ | |
1097 /* ------------------------------------------------------------- */ | |
1098 | |
1099 /* UMFPACK/MA38-style approximate degree: */ | |
1100 if (aggressive) | |
1101 { | |
1102 for (p = p1 ; p <= p2 ; p++) | |
1103 { | |
1104 e = Iw [p] ; | |
1105 ASSERT (e >= 0 && e < n) ; | |
1106 we = W [e] ; | |
1107 if (we != 0) | |
1108 { | |
1109 /* e is an unabsorbed element */ | |
1110 /* dext = | Le \ Lme | */ | |
1111 dext = we - wflg ; | |
1112 if (dext > 0) | |
1113 { | |
1114 deg += dext ; | |
1115 Iw [pn++] = e ; | |
1116 hash += e ; | |
1117 AMD_DEBUG4 ((" e: "ID" hash = "ID"\n",e,hash)) ; | |
1118 } | |
1119 else | |
1120 { | |
1121 /* external degree of e is zero, absorb e into me*/ | |
1122 AMD_DEBUG1 ((" Element "ID" =>"ID" (aggressive)\n", | |
1123 e, me)) ; | |
1124 ASSERT (dext == 0) ; | |
1125 Pe [e] = FLIP (me) ; | |
1126 W [e] = 0 ; | |
1127 } | |
1128 } | |
1129 } | |
1130 } | |
1131 else | |
1132 { | |
1133 for (p = p1 ; p <= p2 ; p++) | |
1134 { | |
1135 e = Iw [p] ; | |
1136 ASSERT (e >= 0 && e < n) ; | |
1137 we = W [e] ; | |
1138 if (we != 0) | |
1139 { | |
1140 /* e is an unabsorbed element */ | |
1141 dext = we - wflg ; | |
1142 ASSERT (dext >= 0) ; | |
1143 deg += dext ; | |
1144 Iw [pn++] = e ; | |
1145 hash += e ; | |
1146 AMD_DEBUG4 ((" e: "ID" hash = "ID"\n",e,hash)) ; | |
1147 } | |
1148 } | |
1149 } | |
1150 | |
1151 /* count the number of elements in i (including me): */ | |
1152 Elen [i] = pn - p1 + 1 ; | |
1153 | |
1154 /* ------------------------------------------------------------- */ | |
1155 /* scan the supervariables in the list associated with i */ | |
1156 /* ------------------------------------------------------------- */ | |
1157 | |
1158 /* The bulk of the AMD run time is typically spent in this loop, | |
1159 * particularly if the matrix has many dense rows that are not | |
1160 * removed prior to ordering. */ | |
1161 p3 = pn ; | |
1162 p4 = p1 + Len [i] ; | |
1163 for (p = p2 + 1 ; p < p4 ; p++) | |
1164 { | |
1165 j = Iw [p] ; | |
1166 ASSERT (j >= 0 && j < n) ; | |
1167 nvj = Nv [j] ; | |
1168 if (nvj > 0) | |
1169 { | |
1170 /* j is unabsorbed, and not in Lme. */ | |
1171 /* add to degree and add to new list */ | |
1172 deg += nvj ; | |
1173 Iw [pn++] = j ; | |
1174 hash += j ; | |
1175 AMD_DEBUG4 ((" s: "ID" hash "ID" Nv[j]= "ID"\n", | |
1176 j, hash, nvj)) ; | |
1177 } | |
1178 } | |
1179 | |
1180 /* ------------------------------------------------------------- */ | |
1181 /* update the degree and check for mass elimination */ | |
1182 /* ------------------------------------------------------------- */ | |
1183 | |
1184 /* with aggressive absorption, deg==0 is identical to the | |
1185 * Elen [i] == 1 && p3 == pn test, below. */ | |
1186 ASSERT (IMPLIES (aggressive, (deg==0) == (Elen[i]==1 && p3==pn))) ; | |
1187 | |
1188 if (Elen [i] == 1 && p3 == pn) | |
1189 { | |
1190 | |
1191 /* --------------------------------------------------------- */ | |
1192 /* mass elimination */ | |
1193 /* --------------------------------------------------------- */ | |
1194 | |
1195 /* There is nothing left of this node except for an edge to | |
1196 * the current pivot element. Elen [i] is 1, and there are | |
1197 * no variables adjacent to node i. Absorb i into the | |
1198 * current pivot element, me. Note that if there are two or | |
1199 * more mass eliminations, fillin due to mass elimination is | |
1200 * possible within the nvpiv-by-nvpiv pivot block. It is this | |
1201 * step that causes AMD's analysis to be an upper bound. | |
1202 * | |
1203 * The reason is that the selected pivot has a lower | |
1204 * approximate degree than the true degree of the two mass | |
1205 * eliminated nodes. There is no edge between the two mass | |
1206 * eliminated nodes. They are merged with the current pivot | |
1207 * anyway. | |
1208 * | |
1209 * No fillin occurs in the Schur complement, in any case, | |
1210 * and this effect does not decrease the quality of the | |
1211 * ordering itself, just the quality of the nonzero and | |
1212 * flop count analysis. It also means that the post-ordering | |
1213 * is not an exact elimination tree post-ordering. */ | |
1214 | |
1215 AMD_DEBUG1 ((" MASS i "ID" => parent e "ID"\n", i, me)) ; | |
1216 Pe [i] = FLIP (me) ; | |
1217 nvi = -Nv [i] ; | |
1218 degme -= nvi ; | |
1219 nvpiv += nvi ; | |
1220 nel += nvi ; | |
1221 Nv [i] = 0 ; | |
1222 Elen [i] = EMPTY ; | |
1223 | |
1224 } | |
1225 else | |
1226 { | |
1227 | |
1228 /* --------------------------------------------------------- */ | |
1229 /* update the upper-bound degree of i */ | |
1230 /* --------------------------------------------------------- */ | |
1231 | |
1232 /* the following degree does not yet include the size | |
1233 * of the current element, which is added later: */ | |
1234 | |
1235 Degree [i] = MIN (Degree [i], deg) ; | |
1236 | |
1237 /* --------------------------------------------------------- */ | |
1238 /* add me to the list for i */ | |
1239 /* --------------------------------------------------------- */ | |
1240 | |
1241 /* move first supervariable to end of list */ | |
1242 Iw [pn] = Iw [p3] ; | |
1243 /* move first element to end of element part of list */ | |
1244 Iw [p3] = Iw [p1] ; | |
1245 /* add new element, me, to front of list. */ | |
1246 Iw [p1] = me ; | |
1247 /* store the new length of the list in Len [i] */ | |
1248 Len [i] = pn - p1 + 1 ; | |
1249 | |
1250 /* --------------------------------------------------------- */ | |
1251 /* place in hash bucket. Save hash key of i in Last [i]. */ | |
1252 /* --------------------------------------------------------- */ | |
1253 | |
1254 /* NOTE: this can fail if hash is negative, because the ANSI C | |
1255 * standard does not define a % b when a and/or b are negative. | |
1256 * That's why hash is defined as an unsigned Int, to avoid this | |
1257 * problem. */ | |
1258 hash = hash % n ; | |
1259 ASSERT (((Int) hash) >= 0 && ((Int) hash) < n) ; | |
1260 | |
1261 /* if the Hhead array is not used: */ | |
1262 j = Head [hash] ; | |
1263 if (j <= EMPTY) | |
1264 { | |
1265 /* degree list is empty, hash head is FLIP (j) */ | |
1266 Next [i] = FLIP (j) ; | |
1267 Head [hash] = FLIP (i) ; | |
1268 } | |
1269 else | |
1270 { | |
1271 /* degree list is not empty, use Last [Head [hash]] as | |
1272 * hash head. */ | |
1273 Next [i] = Last [j] ; | |
1274 Last [j] = i ; | |
1275 } | |
1276 | |
1277 /* if a separate Hhead array is used: * | |
1278 Next [i] = Hhead [hash] ; | |
1279 Hhead [hash] = i ; | |
1280 */ | |
1281 | |
1282 Last [i] = hash ; | |
1283 } | |
1284 } | |
1285 | |
1286 Degree [me] = degme ; | |
1287 | |
1288 /* ----------------------------------------------------------------- */ | |
1289 /* Clear the counter array, W [...], by incrementing wflg. */ | |
1290 /* ----------------------------------------------------------------- */ | |
1291 | |
1292 /* make sure that wflg+n does not cause integer overflow */ | |
1293 lemax = MAX (lemax, degme) ; | |
1294 wflg += lemax ; | |
1295 wflg = clear_flag (wflg, wbig, W, n) ; | |
1296 /* at this point, W [0..n-1] < wflg holds */ | |
1297 | |
1298 /* ========================================================================= */ | |
1299 /* SUPERVARIABLE DETECTION */ | |
1300 /* ========================================================================= */ | |
1301 | |
1302 AMD_DEBUG1 (("Detecting supervariables:\n")) ; | |
1303 for (pme = pme1 ; pme <= pme2 ; pme++) | |
1304 { | |
1305 i = Iw [pme] ; | |
1306 ASSERT (i >= 0 && i < n) ; | |
1307 AMD_DEBUG2 (("Consider i "ID" nv "ID"\n", i, Nv [i])) ; | |
1308 if (Nv [i] < 0) | |
1309 { | |
1310 /* i is a principal variable in Lme */ | |
1311 | |
1312 /* --------------------------------------------------------- | |
1313 * examine all hash buckets with 2 or more variables. We do | |
1314 * this by examing all unique hash keys for supervariables in | |
1315 * the pattern Lme of the current element, me | |
1316 * --------------------------------------------------------- */ | |
1317 | |
1318 /* let i = head of hash bucket, and empty the hash bucket */ | |
1319 ASSERT (Last [i] >= 0 && Last [i] < n) ; | |
1320 hash = Last [i] ; | |
1321 | |
1322 /* if Hhead array is not used: */ | |
1323 j = Head [hash] ; | |
1324 if (j == EMPTY) | |
1325 { | |
1326 /* hash bucket and degree list are both empty */ | |
1327 i = EMPTY ; | |
1328 } | |
1329 else if (j < EMPTY) | |
1330 { | |
1331 /* degree list is empty */ | |
1332 i = FLIP (j) ; | |
1333 Head [hash] = EMPTY ; | |
1334 } | |
1335 else | |
1336 { | |
1337 /* degree list is not empty, restore Last [j] of head j */ | |
1338 i = Last [j] ; | |
1339 Last [j] = EMPTY ; | |
1340 } | |
1341 | |
1342 /* if separate Hhead array is used: * | |
1343 i = Hhead [hash] ; | |
1344 Hhead [hash] = EMPTY ; | |
1345 */ | |
1346 | |
1347 ASSERT (i >= EMPTY && i < n) ; | |
1348 AMD_DEBUG2 (("----i "ID" hash "ID"\n", i, hash)) ; | |
1349 | |
1350 while (i != EMPTY && Next [i] != EMPTY) | |
1351 { | |
1352 | |
1353 /* ----------------------------------------------------- | |
1354 * this bucket has one or more variables following i. | |
1355 * scan all of them to see if i can absorb any entries | |
1356 * that follow i in hash bucket. Scatter i into w. | |
1357 * ----------------------------------------------------- */ | |
1358 | |
1359 ln = Len [i] ; | |
1360 eln = Elen [i] ; | |
1361 ASSERT (ln >= 0 && eln >= 0) ; | |
1362 ASSERT (Pe [i] >= 0 && Pe [i] < iwlen) ; | |
1363 /* do not flag the first element in the list (me) */ | |
1364 for (p = Pe [i] + 1 ; p <= Pe [i] + ln - 1 ; p++) | |
1365 { | |
1366 ASSERT (Iw [p] >= 0 && Iw [p] < n) ; | |
1367 W [Iw [p]] = wflg ; | |
1368 } | |
1369 | |
1370 /* ----------------------------------------------------- */ | |
1371 /* scan every other entry j following i in bucket */ | |
1372 /* ----------------------------------------------------- */ | |
1373 | |
1374 jlast = i ; | |
1375 j = Next [i] ; | |
1376 ASSERT (j >= EMPTY && j < n) ; | |
1377 | |
1378 while (j != EMPTY) | |
1379 { | |
1380 /* ------------------------------------------------- */ | |
1381 /* check if j and i have identical nonzero pattern */ | |
1382 /* ------------------------------------------------- */ | |
1383 | |
1384 AMD_DEBUG3 (("compare i "ID" and j "ID"\n", i,j)) ; | |
1385 | |
1386 /* check if i and j have the same Len and Elen */ | |
1387 ASSERT (Len [j] >= 0 && Elen [j] >= 0) ; | |
1388 ASSERT (Pe [j] >= 0 && Pe [j] < iwlen) ; | |
1389 ok = (Len [j] == ln) && (Elen [j] == eln) ; | |
1390 /* skip the first element in the list (me) */ | |
1391 for (p = Pe [j] + 1 ; ok && p <= Pe [j] + ln - 1 ; p++) | |
1392 { | |
1393 ASSERT (Iw [p] >= 0 && Iw [p] < n) ; | |
1394 if (W [Iw [p]] != wflg) ok = 0 ; | |
1395 } | |
1396 if (ok) | |
1397 { | |
1398 /* --------------------------------------------- */ | |
1399 /* found it! j can be absorbed into i */ | |
1400 /* --------------------------------------------- */ | |
1401 | |
1402 AMD_DEBUG1 (("found it! j "ID" => i "ID"\n", j,i)); | |
1403 Pe [j] = FLIP (i) ; | |
1404 /* both Nv [i] and Nv [j] are negated since they */ | |
1405 /* are in Lme, and the absolute values of each */ | |
1406 /* are the number of variables in i and j: */ | |
1407 Nv [i] += Nv [j] ; | |
1408 Nv [j] = 0 ; | |
1409 Elen [j] = EMPTY ; | |
1410 /* delete j from hash bucket */ | |
1411 ASSERT (j != Next [j]) ; | |
1412 j = Next [j] ; | |
1413 Next [jlast] = j ; | |
1414 | |
1415 } | |
1416 else | |
1417 { | |
1418 /* j cannot be absorbed into i */ | |
1419 jlast = j ; | |
1420 ASSERT (j != Next [j]) ; | |
1421 j = Next [j] ; | |
1422 } | |
1423 ASSERT (j >= EMPTY && j < n) ; | |
1424 } | |
1425 | |
1426 /* ----------------------------------------------------- | |
1427 * no more variables can be absorbed into i | |
1428 * go to next i in bucket and clear flag array | |
1429 * ----------------------------------------------------- */ | |
1430 | |
1431 wflg++ ; | |
1432 i = Next [i] ; | |
1433 ASSERT (i >= EMPTY && i < n) ; | |
1434 | |
1435 } | |
1436 } | |
1437 } | |
1438 AMD_DEBUG2 (("detect done\n")) ; | |
1439 | |
1440 /* ========================================================================= */ | |
1441 /* RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVARIABLES FROM ELEMENT */ | |
1442 /* ========================================================================= */ | |
1443 | |
1444 p = pme1 ; | |
1445 nleft = n - nel ; | |
1446 for (pme = pme1 ; pme <= pme2 ; pme++) | |
1447 { | |
1448 i = Iw [pme] ; | |
1449 ASSERT (i >= 0 && i < n) ; | |
1450 nvi = -Nv [i] ; | |
1451 AMD_DEBUG3 (("Restore i "ID" "ID"\n", i, nvi)) ; | |
1452 if (nvi > 0) | |
1453 { | |
1454 /* i is a principal variable in Lme */ | |
1455 /* restore Nv [i] to signify that i is principal */ | |
1456 Nv [i] = nvi ; | |
1457 | |
1458 /* --------------------------------------------------------- */ | |
1459 /* compute the external degree (add size of current element) */ | |
1460 /* --------------------------------------------------------- */ | |
1461 | |
1462 deg = Degree [i] + degme - nvi ; | |
1463 deg = MIN (deg, nleft - nvi) ; | |
1464 ASSERT (IMPLIES (aggressive, deg > 0) && deg >= 0 && deg < n) ; | |
1465 | |
1466 /* --------------------------------------------------------- */ | |
1467 /* place the supervariable at the head of the degree list */ | |
1468 /* --------------------------------------------------------- */ | |
1469 | |
1470 inext = Head [deg] ; | |
1471 ASSERT (inext >= EMPTY && inext < n) ; | |
1472 if (inext != EMPTY) Last [inext] = i ; | |
1473 Next [i] = inext ; | |
1474 Last [i] = EMPTY ; | |
1475 Head [deg] = i ; | |
1476 | |
1477 /* --------------------------------------------------------- */ | |
1478 /* save the new degree, and find the minimum degree */ | |
1479 /* --------------------------------------------------------- */ | |
1480 | |
1481 mindeg = MIN (mindeg, deg) ; | |
1482 Degree [i] = deg ; | |
1483 | |
1484 /* --------------------------------------------------------- */ | |
1485 /* place the supervariable in the element pattern */ | |
1486 /* --------------------------------------------------------- */ | |
1487 | |
1488 Iw [p++] = i ; | |
1489 | |
1490 } | |
1491 } | |
1492 AMD_DEBUG2 (("restore done\n")) ; | |
1493 | |
1494 /* ========================================================================= */ | |
1495 /* FINALIZE THE NEW ELEMENT */ | |
1496 /* ========================================================================= */ | |
1497 | |
1498 AMD_DEBUG2 (("ME = "ID" DONE\n", me)) ; | |
1499 Nv [me] = nvpiv ; | |
1500 /* save the length of the list for the new element me */ | |
1501 Len [me] = p - pme1 ; | |
1502 if (Len [me] == 0) | |
1503 { | |
1504 /* there is nothing left of the current pivot element */ | |
1505 /* it is a root of the assembly tree */ | |
1506 Pe [me] = EMPTY ; | |
1507 W [me] = 0 ; | |
1508 } | |
1509 if (elenme != 0) | |
1510 { | |
1511 /* element was not constructed in place: deallocate part of */ | |
1512 /* it since newly nonprincipal variables may have been removed */ | |
1513 pfree = p ; | |
1514 } | |
1515 | |
1516 /* The new element has nvpiv pivots and the size of the contribution | |
1517 * block for a multifrontal method is degme-by-degme, not including | |
1518 * the "dense" rows/columns. If the "dense" rows/columns are included, | |
1519 * the frontal matrix is no larger than | |
1520 * (degme+ndense)-by-(degme+ndense). | |
1521 */ | |
1522 | |
1523 if (Info != (double *) NULL) | |
1524 { | |
1525 f = nvpiv ; | |
1526 r = degme + ndense ; | |
1527 dmax = MAX (dmax, f + r) ; | |
1528 | |
1529 /* number of nonzeros in L (excluding the diagonal) */ | |
1530 lnzme = f*r + (f-1)*f/2 ; | |
1531 lnz += lnzme ; | |
1532 | |
1533 /* number of divide operations for LDL' and for LU */ | |
1534 ndiv += lnzme ; | |
1535 | |
1536 /* number of multiply-subtract pairs for LU */ | |
1537 s = f*r*r + r*(f-1)*f + (f-1)*f*(2*f-1)/6 ; | |
1538 nms_lu += s ; | |
1539 | |
1540 /* number of multiply-subtract pairs for LDL' */ | |
1541 nms_ldl += (s + lnzme)/2 ; | |
1542 } | |
1543 | |
1544 #ifndef NDEBUG | |
1545 AMD_DEBUG2 (("finalize done nel "ID" n "ID"\n ::::\n", nel, n)) ; | |
1546 for (pme = Pe [me] ; pme <= Pe [me] + Len [me] - 1 ; pme++) | |
1547 { | |
1548 AMD_DEBUG3 ((" "ID"", Iw [pme])) ; | |
1549 } | |
1550 AMD_DEBUG3 (("\n")) ; | |
1551 #endif | |
1552 | |
1553 } | |
1554 | |
1555 /* ========================================================================= */ | |
1556 /* DONE SELECTING PIVOTS */ | |
1557 /* ========================================================================= */ | |
1558 | |
1559 if (Info != (double *) NULL) | |
1560 { | |
1561 | |
1562 /* count the work to factorize the ndense-by-ndense submatrix */ | |
1563 f = ndense ; | |
1564 dmax = MAX (dmax, (double) ndense) ; | |
1565 | |
1566 /* number of nonzeros in L (excluding the diagonal) */ | |
1567 lnzme = (f-1)*f/2 ; | |
1568 lnz += lnzme ; | |
1569 | |
1570 /* number of divide operations for LDL' and for LU */ | |
1571 ndiv += lnzme ; | |
1572 | |
1573 /* number of multiply-subtract pairs for LU */ | |
1574 s = (f-1)*f*(2*f-1)/6 ; | |
1575 nms_lu += s ; | |
1576 | |
1577 /* number of multiply-subtract pairs for LDL' */ | |
1578 nms_ldl += (s + lnzme)/2 ; | |
1579 | |
1580 /* number of nz's in L (excl. diagonal) */ | |
1581 Info [AMD_LNZ] = lnz ; | |
1582 | |
1583 /* number of divide ops for LU and LDL' */ | |
1584 Info [AMD_NDIV] = ndiv ; | |
1585 | |
1586 /* number of multiply-subtract pairs for LDL' */ | |
1587 Info [AMD_NMULTSUBS_LDL] = nms_ldl ; | |
1588 | |
1589 /* number of multiply-subtract pairs for LU */ | |
1590 Info [AMD_NMULTSUBS_LU] = nms_lu ; | |
1591 | |
1592 /* number of "dense" rows/columns */ | |
1593 Info [AMD_NDENSE] = ndense ; | |
1594 | |
1595 /* largest front is dmax-by-dmax */ | |
1596 Info [AMD_DMAX] = dmax ; | |
1597 | |
1598 /* number of garbage collections in AMD */ | |
1599 Info [AMD_NCMPA] = ncmpa ; | |
1600 | |
1601 /* successful ordering */ | |
1602 Info [AMD_STATUS] = AMD_OK ; | |
1603 } | |
1604 | |
1605 /* ========================================================================= */ | |
1606 /* POST-ORDERING */ | |
1607 /* ========================================================================= */ | |
1608 | |
1609 /* ------------------------------------------------------------------------- | |
1610 * Variables at this point: | |
1611 * | |
1612 * Pe: holds the elimination tree. The parent of j is FLIP (Pe [j]), | |
1613 * or EMPTY if j is a root. The tree holds both elements and | |
1614 * non-principal (unordered) variables absorbed into them. | |
1615 * Dense variables are non-principal and unordered. | |
1616 * | |
1617 * Elen: holds the size of each element, including the diagonal part. | |
1618 * FLIP (Elen [e]) > 0 if e is an element. For unordered | |
1619 * variables i, Elen [i] is EMPTY. | |
1620 * | |
1621 * Nv: Nv [e] > 0 is the number of pivots represented by the element e. | |
1622 * For unordered variables i, Nv [i] is zero. | |
1623 * | |
1624 * Contents no longer needed: | |
1625 * W, Iw, Len, Degree, Head, Next, Last. | |
1626 * | |
1627 * The matrix itself has been destroyed. | |
1628 * | |
1629 * n: the size of the matrix. | |
1630 * No other scalars needed (pfree, iwlen, etc.) | |
1631 * ------------------------------------------------------------------------- */ | |
1632 | |
1633 /* restore Pe */ | |
1634 for (i = 0 ; i < n ; i++) | |
1635 { | |
1636 Pe [i] = FLIP (Pe [i]) ; | |
1637 } | |
1638 | |
1639 /* restore Elen, for output information, and for postordering */ | |
1640 for (i = 0 ; i < n ; i++) | |
1641 { | |
1642 Elen [i] = FLIP (Elen [i]) ; | |
1643 } | |
1644 | |
1645 /* Now the parent of j is Pe [j], or EMPTY if j is a root. Elen [e] > 0 | |
1646 * is the size of element e. Elen [i] is EMPTY for unordered variable i. */ | |
1647 | |
1648 #ifndef NDEBUG | |
1649 AMD_DEBUG2 (("\nTree:\n")) ; | |
1650 for (i = 0 ; i < n ; i++) | |
1651 { | |
1652 AMD_DEBUG2 ((" "ID" parent: "ID" ", i, Pe [i])) ; | |
1653 ASSERT (Pe [i] >= EMPTY && Pe [i] < n) ; | |
1654 if (Nv [i] > 0) | |
1655 { | |
1656 /* this is an element */ | |
1657 e = i ; | |
1658 AMD_DEBUG2 ((" element, size is "ID"\n", Elen [i])) ; | |
1659 ASSERT (Elen [e] > 0) ; | |
1660 } | |
1661 AMD_DEBUG2 (("\n")) ; | |
1662 } | |
1663 AMD_DEBUG2 (("\nelements:\n")) ; | |
1664 for (e = 0 ; e < n ; e++) | |
1665 { | |
1666 if (Nv [e] > 0) | |
1667 { | |
1668 AMD_DEBUG3 (("Element e= "ID" size "ID" nv "ID" \n", e, | |
1669 Elen [e], Nv [e])) ; | |
1670 } | |
1671 } | |
1672 AMD_DEBUG2 (("\nvariables:\n")) ; | |
1673 for (i = 0 ; i < n ; i++) | |
1674 { | |
1675 Int cnt ; | |
1676 if (Nv [i] == 0) | |
1677 { | |
1678 AMD_DEBUG3 (("i unordered: "ID"\n", i)) ; | |
1679 j = Pe [i] ; | |
1680 cnt = 0 ; | |
1681 AMD_DEBUG3 ((" j: "ID"\n", j)) ; | |
1682 if (j == EMPTY) | |
1683 { | |
1684 AMD_DEBUG3 ((" i is a dense variable\n")) ; | |
1685 } | |
1686 else | |
1687 { | |
1688 ASSERT (j >= 0 && j < n) ; | |
1689 while (Nv [j] == 0) | |
1690 { | |
1691 AMD_DEBUG3 ((" j : "ID"\n", j)) ; | |
1692 j = Pe [j] ; | |
1693 AMD_DEBUG3 ((" j:: "ID"\n", j)) ; | |
1694 cnt++ ; | |
1695 if (cnt > n) break ; | |
1696 } | |
1697 e = j ; | |
1698 AMD_DEBUG3 ((" got to e: "ID"\n", e)) ; | |
1699 } | |
1700 } | |
1701 } | |
1702 #endif | |
1703 | |
1704 /* ========================================================================= */ | |
1705 /* compress the paths of the variables */ | |
1706 /* ========================================================================= */ | |
1707 | |
1708 for (i = 0 ; i < n ; i++) | |
1709 { | |
1710 if (Nv [i] == 0) | |
1711 { | |
1712 | |
1713 /* ------------------------------------------------------------- | |
1714 * i is an un-ordered row. Traverse the tree from i until | |
1715 * reaching an element, e. The element, e, was the principal | |
1716 * supervariable of i and all nodes in the path from i to when e | |
1717 * was selected as pivot. | |
1718 * ------------------------------------------------------------- */ | |
1719 | |
1720 AMD_DEBUG1 (("Path compression, i unordered: "ID"\n", i)) ; | |
1721 j = Pe [i] ; | |
1722 ASSERT (j >= EMPTY && j < n) ; | |
1723 AMD_DEBUG3 ((" j: "ID"\n", j)) ; | |
1724 if (j == EMPTY) | |
1725 { | |
1726 /* Skip a dense variable. It has no parent. */ | |
1727 AMD_DEBUG3 ((" i is a dense variable\n")) ; | |
1728 continue ; | |
1729 } | |
1730 | |
1731 /* while (j is a variable) */ | |
1732 while (Nv [j] == 0) | |
1733 { | |
1734 AMD_DEBUG3 ((" j : "ID"\n", j)) ; | |
1735 j = Pe [j] ; | |
1736 AMD_DEBUG3 ((" j:: "ID"\n", j)) ; | |
1737 ASSERT (j >= 0 && j < n) ; | |
1738 } | |
1739 /* got to an element e */ | |
1740 e = j ; | |
1741 AMD_DEBUG3 (("got to e: "ID"\n", e)) ; | |
1742 | |
1743 /* ------------------------------------------------------------- | |
1744 * traverse the path again from i to e, and compress the path | |
1745 * (all nodes point to e). Path compression allows this code to | |
1746 * compute in O(n) time. | |
1747 * ------------------------------------------------------------- */ | |
1748 | |
1749 j = i ; | |
1750 /* while (j is a variable) */ | |
1751 while (Nv [j] == 0) | |
1752 { | |
1753 jnext = Pe [j] ; | |
1754 AMD_DEBUG3 (("j "ID" jnext "ID"\n", j, jnext)) ; | |
1755 Pe [j] = e ; | |
1756 j = jnext ; | |
1757 ASSERT (j >= 0 && j < n) ; | |
1758 } | |
1759 } | |
1760 } | |
1761 | |
1762 /* ========================================================================= */ | |
1763 /* postorder the assembly tree */ | |
1764 /* ========================================================================= */ | |
1765 | |
1766 AMD_postorder (n, Pe, Nv, Elen, | |
1767 W, /* output order */ | |
1768 Head, Next, Last) ; /* workspace */ | |
1769 | |
1770 /* ========================================================================= */ | |
1771 /* compute output permutation and inverse permutation */ | |
1772 /* ========================================================================= */ | |
1773 | |
1774 /* W [e] = k means that element e is the kth element in the new | |
1775 * order. e is in the range 0 to n-1, and k is in the range 0 to | |
1776 * the number of elements. Use Head for inverse order. */ | |
1777 | |
1778 for (k = 0 ; k < n ; k++) | |
1779 { | |
1780 Head [k] = EMPTY ; | |
1781 Next [k] = EMPTY ; | |
1782 } | |
1783 for (e = 0 ; e < n ; e++) | |
1784 { | |
1785 k = W [e] ; | |
1786 ASSERT ((k == EMPTY) == (Nv [e] == 0)) ; | |
1787 if (k != EMPTY) | |
1788 { | |
1789 ASSERT (k >= 0 && k < n) ; | |
1790 Head [k] = e ; | |
1791 } | |
1792 } | |
1793 | |
1794 /* construct output inverse permutation in Next, | |
1795 * and permutation in Last */ | |
1796 nel = 0 ; | |
1797 for (k = 0 ; k < n ; k++) | |
1798 { | |
1799 e = Head [k] ; | |
1800 if (e == EMPTY) break ; | |
1801 ASSERT (e >= 0 && e < n && Nv [e] > 0) ; | |
1802 Next [e] = nel ; | |
1803 nel += Nv [e] ; | |
1804 } | |
1805 ASSERT (nel == n - ndense) ; | |
1806 | |
1807 /* order non-principal variables (dense, & those merged into supervar's) */ | |
1808 for (i = 0 ; i < n ; i++) | |
1809 { | |
1810 if (Nv [i] == 0) | |
1811 { | |
1812 e = Pe [i] ; | |
1813 ASSERT (e >= EMPTY && e < n) ; | |
1814 if (e != EMPTY) | |
1815 { | |
1816 /* This is an unordered variable that was merged | |
1817 * into element e via supernode detection or mass | |
1818 * elimination of i when e became the pivot element. | |
1819 * Place i in order just before e. */ | |
1820 ASSERT (Next [i] == EMPTY && Nv [e] > 0) ; | |
1821 Next [i] = Next [e] ; | |
1822 Next [e]++ ; | |
1823 } | |
1824 else | |
1825 { | |
1826 /* This is a dense unordered variable, with no parent. | |
1827 * Place it last in the output order. */ | |
1828 Next [i] = nel++ ; | |
1829 } | |
1830 } | |
1831 } | |
1832 ASSERT (nel == n) ; | |
1833 | |
1834 AMD_DEBUG2 (("\n\nPerm:\n")) ; | |
1835 for (i = 0 ; i < n ; i++) | |
1836 { | |
1837 k = Next [i] ; | |
1838 ASSERT (k >= 0 && k < n) ; | |
1839 Last [k] = i ; | |
1840 AMD_DEBUG2 ((" perm ["ID"] = "ID"\n", k, i)) ; | |
1841 } | |
1842 } |