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