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 = f*r + (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 = f*r*r + 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 }