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