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 | } |
---|