lemon-project-template-glpk: deps/glpk/src/amd/amd

lemon-project-template-glpk

annotate deps/glpk/src/amd/amd_2.c @ 9:33de93886c88

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

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