alpar@1: /* ========================================================================= */ alpar@1: /* === AMD_2 =============================================================== */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: /* ------------------------------------------------------------------------- */ alpar@1: /* AMD, Copyright (c) Timothy A. Davis, */ alpar@1: /* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */ alpar@1: /* email: davis at cise.ufl.edu CISE Department, Univ. of Florida. */ alpar@1: /* web: http://www.cise.ufl.edu/research/sparse/amd */ alpar@1: /* ------------------------------------------------------------------------- */ alpar@1: alpar@1: /* AMD_2: performs the AMD ordering on a symmetric sparse matrix A, followed alpar@1: * by a postordering (via depth-first search) of the assembly tree using the alpar@1: * AMD_postorder routine. alpar@1: */ alpar@1: alpar@1: #include "amd_internal.h" alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* === clear_flag ========================================================== */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: static Int clear_flag (Int wflg, Int wbig, Int W [ ], Int n) alpar@1: { alpar@1: Int x ; alpar@1: if (wflg < 2 || wflg >= wbig) alpar@1: { alpar@1: for (x = 0 ; x < n ; x++) alpar@1: { alpar@1: if (W [x] != 0) W [x] = 1 ; alpar@1: } alpar@1: wflg = 2 ; alpar@1: } alpar@1: /* at this point, W [0..n-1] < wflg holds */ alpar@1: return (wflg) ; alpar@1: } alpar@1: alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* === AMD_2 =============================================================== */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: GLOBAL void AMD_2 alpar@1: ( alpar@1: Int n, /* A is n-by-n, where n > 0 */ alpar@1: Int Pe [ ], /* Pe [0..n-1]: index in Iw of row i on input */ alpar@1: Int Iw [ ], /* workspace of size iwlen. Iw [0..pfree-1] alpar@1: * holds the matrix on input */ alpar@1: Int Len [ ], /* Len [0..n-1]: length for row/column i on input */ alpar@1: Int iwlen, /* length of Iw. iwlen >= pfree + n */ alpar@1: Int pfree, /* Iw [pfree ... iwlen-1] is empty on input */ alpar@1: alpar@1: /* 7 size-n workspaces, not defined on input: */ alpar@1: Int Nv [ ], /* the size of each supernode on output */ alpar@1: Int Next [ ], /* the output inverse permutation */ alpar@1: Int Last [ ], /* the output permutation */ alpar@1: Int Head [ ], alpar@1: Int Elen [ ], /* the size columns of L for each supernode */ alpar@1: Int Degree [ ], alpar@1: Int W [ ], alpar@1: alpar@1: /* control parameters and output statistics */ alpar@1: double Control [ ], /* array of size AMD_CONTROL */ alpar@1: double Info [ ] /* array of size AMD_INFO */ alpar@1: ) alpar@1: { alpar@1: alpar@1: /* alpar@1: * Given a representation of the nonzero pattern of a symmetric matrix, A, alpar@1: * (excluding the diagonal) perform an approximate minimum (UMFPACK/MA38-style) alpar@1: * degree ordering to compute a pivot order such that the introduction of alpar@1: * nonzeros (fill-in) in the Cholesky factors A = LL' is kept low. At each alpar@1: * step, the pivot selected is the one with the minimum UMFAPACK/MA38-style alpar@1: * upper-bound on the external degree. This routine can optionally perform alpar@1: * aggresive absorption (as done by MC47B in the Harwell Subroutine alpar@1: * Library). alpar@1: * alpar@1: * The approximate degree algorithm implemented here is the symmetric analog of alpar@1: * the degree update algorithm in MA38 and UMFPACK (the Unsymmetric-pattern alpar@1: * MultiFrontal PACKage, both by Davis and Duff). The routine is based on the alpar@1: * MA27 minimum degree ordering algorithm by Iain Duff and John Reid. alpar@1: * alpar@1: * This routine is a translation of the original AMDBAR and MC47B routines, alpar@1: * in Fortran, with the following modifications: alpar@1: * alpar@1: * (1) dense rows/columns are removed prior to ordering the matrix, and placed alpar@1: * last in the output order. The presence of a dense row/column can alpar@1: * increase the ordering time by up to O(n^2), unless they are removed alpar@1: * prior to ordering. alpar@1: * alpar@1: * (2) the minimum degree ordering is followed by a postordering (depth-first alpar@1: * search) of the assembly tree. Note that mass elimination (discussed alpar@1: * below) combined with the approximate degree update can lead to the mass alpar@1: * elimination of nodes with lower exact degree than the current pivot alpar@1: * element. No additional fill-in is caused in the representation of the alpar@1: * Schur complement. The mass-eliminated nodes merge with the current alpar@1: * pivot element. They are ordered prior to the current pivot element. alpar@1: * Because they can have lower exact degree than the current element, the alpar@1: * merger of two or more of these nodes in the current pivot element can alpar@1: * lead to a single element that is not a "fundamental supernode". The alpar@1: * diagonal block can have zeros in it. Thus, the assembly tree used here alpar@1: * is not guaranteed to be the precise supernodal elemination tree (with alpar@1: * "funadmental" supernodes), and the postordering performed by this alpar@1: * routine is not guaranteed to be a precise postordering of the alpar@1: * elimination tree. alpar@1: * alpar@1: * (3) input parameters are added, to control aggressive absorption and the alpar@1: * detection of "dense" rows/columns of A. alpar@1: * alpar@1: * (4) additional statistical information is returned, such as the number of alpar@1: * nonzeros in L, and the flop counts for subsequent LDL' and LU alpar@1: * factorizations. These are slight upper bounds, because of the mass alpar@1: * elimination issue discussed above. alpar@1: * alpar@1: * (5) additional routines are added to interface this routine to MATLAB alpar@1: * to provide a simple C-callable user-interface, to check inputs for alpar@1: * errors, compute the symmetry of the pattern of A and the number of alpar@1: * nonzeros in each row/column of A+A', to compute the pattern of A+A', alpar@1: * to perform the assembly tree postordering, and to provide debugging alpar@1: * ouput. Many of these functions are also provided by the Fortran alpar@1: * Harwell Subroutine Library routine MC47A. alpar@1: * alpar@1: * (6) both int and UF_long versions are provided. In the descriptions below alpar@1: * and integer is and int or UF_long depending on which version is alpar@1: * being used. alpar@1: alpar@1: ********************************************************************** alpar@1: ***** CAUTION: ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT. ****** alpar@1: ********************************************************************** alpar@1: ** If you want error checking, a more versatile input format, and a ** alpar@1: ** simpler user interface, use amd_order or amd_l_order instead. ** alpar@1: ** This routine is not meant to be user-callable. ** alpar@1: ********************************************************************** alpar@1: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * References: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * alpar@1: * [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern multifrontal alpar@1: * method for sparse LU factorization", SIAM J. Matrix Analysis and alpar@1: * Applications, vol. 18, no. 1, pp. 140-158. Discusses UMFPACK / MA38, alpar@1: * which first introduced the approximate minimum degree used by this alpar@1: * routine. alpar@1: * alpar@1: * [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, "An approximate alpar@1: * minimum degree ordering algorithm," SIAM J. Matrix Analysis and alpar@1: * Applications, vol. 17, no. 4, pp. 886-905, 1996. Discusses AMDBAR and alpar@1: * MC47B, which are the Fortran versions of this routine. alpar@1: * alpar@1: * [3] Alan George and Joseph Liu, "The evolution of the minimum degree alpar@1: * ordering algorithm," SIAM Review, vol. 31, no. 1, pp. 1-19, 1989. alpar@1: * We list below the features mentioned in that paper that this code alpar@1: * includes: alpar@1: * alpar@1: * mass elimination: alpar@1: * Yes. MA27 relied on supervariable detection for mass elimination. alpar@1: * alpar@1: * indistinguishable nodes: alpar@1: * Yes (we call these "supervariables"). This was also in the MA27 alpar@1: * code - although we modified the method of detecting them (the alpar@1: * previous hash was the true degree, which we no longer keep track alpar@1: * of). A supervariable is a set of rows with identical nonzero alpar@1: * pattern. All variables in a supervariable are eliminated together. alpar@1: * Each supervariable has as its numerical name that of one of its alpar@1: * variables (its principal variable). alpar@1: * alpar@1: * quotient graph representation: alpar@1: * Yes. We use the term "element" for the cliques formed during alpar@1: * elimination. This was also in the MA27 code. The algorithm can alpar@1: * operate in place, but it will work more efficiently if given some alpar@1: * "elbow room." alpar@1: * alpar@1: * element absorption: alpar@1: * Yes. This was also in the MA27 code. alpar@1: * alpar@1: * external degree: alpar@1: * Yes. The MA27 code was based on the true degree. alpar@1: * alpar@1: * incomplete degree update and multiple elimination: alpar@1: * No. This was not in MA27, either. Our method of degree update alpar@1: * within MC47B is element-based, not variable-based. It is thus alpar@1: * not well-suited for use with incomplete degree update or multiple alpar@1: * elimination. alpar@1: * alpar@1: * Authors, and Copyright (C) 2004 by: alpar@1: * Timothy A. Davis, Patrick Amestoy, Iain S. Duff, John K. Reid. alpar@1: * alpar@1: * Acknowledgements: This work (and the UMFPACK package) was supported by the alpar@1: * National Science Foundation (ASC-9111263, DMS-9223088, and CCR-0203270). alpar@1: * The UMFPACK/MA38 approximate degree update algorithm, the unsymmetric analog alpar@1: * which forms the basis of AMD, was developed while Tim Davis was supported by alpar@1: * CERFACS (Toulouse, France) in a post-doctoral position. This C version, and alpar@1: * the etree postorder, were written while Tim Davis was on sabbatical at alpar@1: * Stanford University and Lawrence Berkeley National Laboratory. alpar@1: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * INPUT ARGUMENTS (unaltered): alpar@1: * ---------------------------------------------------------------------------- alpar@1: alpar@1: * n: The matrix order. Restriction: n >= 1. alpar@1: * alpar@1: * iwlen: The size of the Iw array. On input, the matrix is stored in alpar@1: * Iw [0..pfree-1]. However, Iw [0..iwlen-1] should be slightly larger alpar@1: * than what is required to hold the matrix, at least iwlen >= pfree + n. alpar@1: * Otherwise, excessive compressions will take place. The recommended alpar@1: * value of iwlen is 1.2 * pfree + n, which is the value used in the alpar@1: * user-callable interface to this routine (amd_order.c). The algorithm alpar@1: * will not run at all if iwlen < pfree. Restriction: iwlen >= pfree + n. alpar@1: * Note that this is slightly more restrictive than the actual minimum alpar@1: * (iwlen >= pfree), but AMD_2 will be very slow with no elbow room. alpar@1: * Thus, this routine enforces a bare minimum elbow room of size n. alpar@1: * alpar@1: * pfree: On input the tail end of the array, Iw [pfree..iwlen-1], is empty, alpar@1: * and the matrix is stored in Iw [0..pfree-1]. During execution, alpar@1: * additional data is placed in Iw, and pfree is modified so that alpar@1: * Iw [pfree..iwlen-1] is always the unused part of Iw. alpar@1: * alpar@1: * Control: A double array of size AMD_CONTROL containing input parameters alpar@1: * that affect how the ordering is computed. If NULL, then default alpar@1: * settings are used. alpar@1: * alpar@1: * Control [AMD_DENSE] is used to determine whether or not a given input alpar@1: * row is "dense". A row is "dense" if the number of entries in the row alpar@1: * exceeds Control [AMD_DENSE] times sqrt (n), except that rows with 16 or alpar@1: * fewer entries are never considered "dense". To turn off the detection alpar@1: * of dense rows, set Control [AMD_DENSE] to a negative number, or to a alpar@1: * number larger than sqrt (n). The default value of Control [AMD_DENSE] alpar@1: * is AMD_DEFAULT_DENSE, which is defined in amd.h as 10. alpar@1: * alpar@1: * Control [AMD_AGGRESSIVE] is used to determine whether or not aggressive alpar@1: * absorption is to be performed. If nonzero, then aggressive absorption alpar@1: * is performed (this is the default). alpar@1: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * INPUT/OUPUT ARGUMENTS: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * alpar@1: * Pe: An integer array of size n. On input, Pe [i] is the index in Iw of alpar@1: * the start of row i. Pe [i] is ignored if row i has no off-diagonal alpar@1: * entries. Thus Pe [i] must be in the range 0 to pfree-1 for non-empty alpar@1: * rows. alpar@1: * alpar@1: * During execution, it is used for both supervariables and elements: alpar@1: * alpar@1: * Principal supervariable i: index into Iw of the description of alpar@1: * supervariable i. A supervariable represents one or more rows of alpar@1: * the matrix with identical nonzero pattern. In this case, alpar@1: * Pe [i] >= 0. alpar@1: * alpar@1: * Non-principal supervariable i: if i has been absorbed into another alpar@1: * supervariable j, then Pe [i] = FLIP (j), where FLIP (j) is defined alpar@1: * as (-(j)-2). Row j has the same pattern as row i. Note that j alpar@1: * might later be absorbed into another supervariable j2, in which alpar@1: * case Pe [i] is still FLIP (j), and Pe [j] = FLIP (j2) which is alpar@1: * < EMPTY, where EMPTY is defined as (-1) in amd_internal.h. alpar@1: * alpar@1: * Unabsorbed element e: the index into Iw of the description of element alpar@1: * e, if e has not yet been absorbed by a subsequent element. Element alpar@1: * e is created when the supervariable of the same name is selected as alpar@1: * the pivot. In this case, Pe [i] >= 0. alpar@1: * alpar@1: * Absorbed element e: if element e is absorbed into element e2, then alpar@1: * Pe [e] = FLIP (e2). This occurs when the pattern of e (which we alpar@1: * refer to as Le) is found to be a subset of the pattern of e2 (that alpar@1: * is, Le2). In this case, Pe [i] < EMPTY. If element e is "null" alpar@1: * (it has no nonzeros outside its pivot block), then Pe [e] = EMPTY, alpar@1: * and e is the root of an assembly subtree (or the whole tree if alpar@1: * there is just one such root). alpar@1: * alpar@1: * Dense variable i: if i is "dense", then Pe [i] = EMPTY. alpar@1: * alpar@1: * On output, Pe holds the assembly tree/forest, which implicitly alpar@1: * represents a pivot order with identical fill-in as the actual order alpar@1: * (via a depth-first search of the tree), as follows. If Nv [i] > 0, alpar@1: * then i represents a node in the assembly tree, and the parent of i is alpar@1: * Pe [i], or EMPTY if i is a root. If Nv [i] = 0, then (i, Pe [i]) alpar@1: * represents an edge in a subtree, the root of which is a node in the alpar@1: * assembly tree. Note that i refers to a row/column in the original alpar@1: * matrix, not the permuted matrix. alpar@1: * alpar@1: * Info: A double array of size AMD_INFO. If present, (that is, not NULL), alpar@1: * then statistics about the ordering are returned in the Info array. alpar@1: * See amd.h for a description. alpar@1: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * INPUT/MODIFIED (undefined on output): alpar@1: * ---------------------------------------------------------------------------- alpar@1: * alpar@1: * Len: An integer array of size n. On input, Len [i] holds the number of alpar@1: * entries in row i of the matrix, excluding the diagonal. The contents alpar@1: * of Len are undefined on output. alpar@1: * alpar@1: * Iw: An integer array of size iwlen. On input, Iw [0..pfree-1] holds the alpar@1: * description of each row i in the matrix. The matrix must be symmetric, alpar@1: * and both upper and lower triangular parts must be present. The alpar@1: * diagonal must not be present. Row i is held as follows: alpar@1: * alpar@1: * Len [i]: the length of the row i data structure in the Iw array. alpar@1: * Iw [Pe [i] ... Pe [i] + Len [i] - 1]: alpar@1: * the list of column indices for nonzeros in row i (simple alpar@1: * supervariables), excluding the diagonal. All supervariables alpar@1: * start with one row/column each (supervariable i is just row i). alpar@1: * If Len [i] is zero on input, then Pe [i] is ignored on input. alpar@1: * alpar@1: * Note that the rows need not be in any particular order, and there alpar@1: * may be empty space between the rows. alpar@1: * alpar@1: * During execution, the supervariable i experiences fill-in. This is alpar@1: * represented by placing in i a list of the elements that cause fill-in alpar@1: * in supervariable i: alpar@1: * alpar@1: * Len [i]: the length of supervariable i in the Iw array. alpar@1: * Iw [Pe [i] ... Pe [i] + Elen [i] - 1]: alpar@1: * the list of elements that contain i. This list is kept short alpar@1: * by removing absorbed elements. alpar@1: * Iw [Pe [i] + Elen [i] ... Pe [i] + Len [i] - 1]: alpar@1: * the list of supervariables in i. This list is kept short by alpar@1: * removing nonprincipal variables, and any entry j that is also alpar@1: * contained in at least one of the elements (j in Le) in the list alpar@1: * for i (e in row i). alpar@1: * alpar@1: * When supervariable i is selected as pivot, we create an element e of alpar@1: * the same name (e=i): alpar@1: * alpar@1: * Len [e]: the length of element e in the Iw array. alpar@1: * Iw [Pe [e] ... Pe [e] + Len [e] - 1]: alpar@1: * the list of supervariables in element e. alpar@1: * alpar@1: * An element represents the fill-in that occurs when supervariable i is alpar@1: * selected as pivot (which represents the selection of row i and all alpar@1: * non-principal variables whose principal variable is i). We use the alpar@1: * term Le to denote the set of all supervariables in element e. Absorbed alpar@1: * supervariables and elements are pruned from these lists when alpar@1: * computationally convenient. alpar@1: * alpar@1: * CAUTION: THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION. alpar@1: * The contents of Iw are undefined on output. alpar@1: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * OUTPUT (need not be set on input): alpar@1: * ---------------------------------------------------------------------------- alpar@1: * alpar@1: * Nv: An integer array of size n. During execution, ABS (Nv [i]) is equal to alpar@1: * the number of rows that are represented by the principal supervariable alpar@1: * i. If i is a nonprincipal or dense variable, then Nv [i] = 0. alpar@1: * Initially, Nv [i] = 1 for all i. Nv [i] < 0 signifies that i is a alpar@1: * principal variable in the pattern Lme of the current pivot element me. alpar@1: * After element me is constructed, Nv [i] is set back to a positive alpar@1: * value. alpar@1: * alpar@1: * On output, Nv [i] holds the number of pivots represented by super alpar@1: * row/column i of the original matrix, or Nv [i] = 0 for non-principal alpar@1: * rows/columns. Note that i refers to a row/column in the original alpar@1: * matrix, not the permuted matrix. alpar@1: * alpar@1: * Elen: An integer array of size n. See the description of Iw above. At the alpar@1: * start of execution, Elen [i] is set to zero for all rows i. During alpar@1: * execution, Elen [i] is the number of elements in the list for alpar@1: * supervariable i. When e becomes an element, Elen [e] = FLIP (esize) is alpar@1: * set, where esize is the size of the element (the number of pivots, plus alpar@1: * the number of nonpivotal entries). Thus Elen [e] < EMPTY. alpar@1: * Elen (i) = EMPTY set when variable i becomes nonprincipal. alpar@1: * alpar@1: * For variables, Elen (i) >= EMPTY holds until just before the alpar@1: * postordering and permutation vectors are computed. For elements, alpar@1: * Elen [e] < EMPTY holds. alpar@1: * alpar@1: * On output, Elen [i] is the degree of the row/column in the Cholesky alpar@1: * factorization of the permuted matrix, corresponding to the original row alpar@1: * i, if i is a super row/column. It is equal to EMPTY if i is alpar@1: * non-principal. Note that i refers to a row/column in the original alpar@1: * matrix, not the permuted matrix. alpar@1: * alpar@1: * Note that the contents of Elen on output differ from the Fortran alpar@1: * version (Elen holds the inverse permutation in the Fortran version, alpar@1: * which is instead returned in the Next array in this C version, alpar@1: * described below). alpar@1: * alpar@1: * Last: In a degree list, Last [i] is the supervariable preceding i, or EMPTY alpar@1: * if i is the head of the list. In a hash bucket, Last [i] is the hash alpar@1: * key for i. alpar@1: * alpar@1: * Last [Head [hash]] is also used as the head of a hash bucket if alpar@1: * Head [hash] contains a degree list (see the description of Head, alpar@1: * below). alpar@1: * alpar@1: * On output, Last [0..n-1] holds the permutation. That is, if alpar@1: * i = Last [k], then row i is the kth pivot row (where k ranges from 0 to alpar@1: * n-1). Row Last [k] of A is the kth row in the permuted matrix, PAP'. alpar@1: * alpar@1: * Next: Next [i] is the supervariable following i in a link list, or EMPTY if alpar@1: * i is the last in the list. Used for two kinds of lists: degree lists alpar@1: * and hash buckets (a supervariable can be in only one kind of list at a alpar@1: * time). alpar@1: * alpar@1: * On output Next [0..n-1] holds the inverse permutation. That is, if alpar@1: * k = Next [i], then row i is the kth pivot row. Row i of A appears as alpar@1: * the (Next[i])-th row in the permuted matrix, PAP'. alpar@1: * alpar@1: * Note that the contents of Next on output differ from the Fortran alpar@1: * version (Next is undefined on output in the Fortran version). alpar@1: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * LOCAL WORKSPACE (not input or output - used only during execution): alpar@1: * ---------------------------------------------------------------------------- alpar@1: * alpar@1: * Degree: An integer array of size n. If i is a supervariable, then alpar@1: * Degree [i] holds the current approximation of the external degree of alpar@1: * row i (an upper bound). The external degree is the number of nonzeros alpar@1: * in row i, minus ABS (Nv [i]), the diagonal part. The bound is equal to alpar@1: * the exact external degree if Elen [i] is less than or equal to two. alpar@1: * alpar@1: * We also use the term "external degree" for elements e to refer to alpar@1: * |Le \ Lme|. If e is an element, then Degree [e] is |Le|, which is the alpar@1: * degree of the off-diagonal part of the element e (not including the alpar@1: * diagonal part). alpar@1: * alpar@1: * Head: An integer array of size n. Head is used for degree lists. alpar@1: * Head [deg] is the first supervariable in a degree list. All alpar@1: * supervariables i in a degree list Head [deg] have the same approximate alpar@1: * degree, namely, deg = Degree [i]. If the list Head [deg] is empty then alpar@1: * Head [deg] = EMPTY. alpar@1: * alpar@1: * During supervariable detection Head [hash] also serves as a pointer to alpar@1: * a hash bucket. If Head [hash] >= 0, there is a degree list of degree alpar@1: * hash. The hash bucket head pointer is Last [Head [hash]]. If alpar@1: * Head [hash] = EMPTY, then the degree list and hash bucket are both alpar@1: * empty. If Head [hash] < EMPTY, then the degree list is empty, and alpar@1: * FLIP (Head [hash]) is the head of the hash bucket. After supervariable alpar@1: * detection is complete, all hash buckets are empty, and the alpar@1: * (Last [Head [hash]] = EMPTY) condition is restored for the non-empty alpar@1: * degree lists. alpar@1: * alpar@1: * W: An integer array of size n. The flag array W determines the status of alpar@1: * elements and variables, and the external degree of elements. alpar@1: * alpar@1: * for elements: alpar@1: * if W [e] = 0, then the element e is absorbed. alpar@1: * if W [e] >= wflg, then W [e] - wflg is the size of the set alpar@1: * |Le \ Lme|, in terms of nonzeros (the sum of ABS (Nv [i]) for alpar@1: * each principal variable i that is both in the pattern of alpar@1: * element e and NOT in the pattern of the current pivot element, alpar@1: * me). alpar@1: * if wflg > W [e] > 0, then e is not absorbed and has not yet been alpar@1: * seen in the scan of the element lists in the computation of alpar@1: * |Le\Lme| in Scan 1 below. alpar@1: * alpar@1: * for variables: alpar@1: * during supervariable detection, if W [j] != wflg then j is alpar@1: * not in the pattern of variable i. alpar@1: * alpar@1: * The W array is initialized by setting W [i] = 1 for all i, and by alpar@1: * setting wflg = 2. It is reinitialized if wflg becomes too large (to alpar@1: * ensure that wflg+n does not cause integer overflow). alpar@1: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * LOCAL INTEGERS: alpar@1: * ---------------------------------------------------------------------------- alpar@1: */ alpar@1: alpar@1: Int deg, degme, dext, lemax, e, elenme, eln, i, ilast, inext, j, alpar@1: jlast, jnext, k, knt1, knt2, knt3, lenj, ln, me, mindeg, nel, nleft, alpar@1: nvi, nvj, nvpiv, slenme, wbig, we, wflg, wnvi, ok, ndense, ncmpa, alpar@1: dense, aggressive ; alpar@1: alpar@1: unsigned Int hash ; /* unsigned, so that hash % n is well defined.*/ alpar@1: alpar@1: /* alpar@1: * deg: the degree of a variable or element alpar@1: * degme: size, |Lme|, of the current element, me (= Degree [me]) alpar@1: * dext: external degree, |Le \ Lme|, of some element e alpar@1: * lemax: largest |Le| seen so far (called dmax in Fortran version) alpar@1: * e: an element alpar@1: * elenme: the length, Elen [me], of element list of pivotal variable alpar@1: * eln: the length, Elen [...], of an element list alpar@1: * hash: the computed value of the hash function alpar@1: * i: a supervariable alpar@1: * ilast: the entry in a link list preceding i alpar@1: * inext: the entry in a link list following i alpar@1: * j: a supervariable alpar@1: * jlast: the entry in a link list preceding j alpar@1: * jnext: the entry in a link list, or path, following j alpar@1: * k: the pivot order of an element or variable alpar@1: * knt1: loop counter used during element construction alpar@1: * knt2: loop counter used during element construction alpar@1: * knt3: loop counter used during compression alpar@1: * lenj: Len [j] alpar@1: * ln: length of a supervariable list alpar@1: * me: current supervariable being eliminated, and the current alpar@1: * element created by eliminating that supervariable alpar@1: * mindeg: current minimum degree alpar@1: * nel: number of pivots selected so far alpar@1: * nleft: n - nel, the number of nonpivotal rows/columns remaining alpar@1: * nvi: the number of variables in a supervariable i (= Nv [i]) alpar@1: * nvj: the number of variables in a supervariable j (= Nv [j]) alpar@1: * nvpiv: number of pivots in current element alpar@1: * slenme: number of variables in variable list of pivotal variable alpar@1: * wbig: = INT_MAX - n for the int version, UF_long_max - n for the alpar@1: * UF_long version. wflg is not allowed to be >= wbig. alpar@1: * we: W [e] alpar@1: * wflg: used for flagging the W array. See description of Iw. alpar@1: * wnvi: wflg - Nv [i] alpar@1: * x: either a supervariable or an element alpar@1: * alpar@1: * ok: true if supervariable j can be absorbed into i alpar@1: * ndense: number of "dense" rows/columns alpar@1: * dense: rows/columns with initial degree > dense are considered "dense" alpar@1: * aggressive: true if aggressive absorption is being performed alpar@1: * ncmpa: number of garbage collections alpar@1: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * LOCAL DOUBLES, used for statistical output only (except for alpha): alpar@1: * ---------------------------------------------------------------------------- alpar@1: */ alpar@1: alpar@1: double f, r, ndiv, s, nms_lu, nms_ldl, dmax, alpha, lnz, lnzme ; alpar@1: alpar@1: /* alpar@1: * f: nvpiv alpar@1: * r: degme + nvpiv alpar@1: * ndiv: number of divisions for LU or LDL' factorizations alpar@1: * s: number of multiply-subtract pairs for LU factorization, for the alpar@1: * current element me alpar@1: * nms_lu number of multiply-subtract pairs for LU factorization alpar@1: * nms_ldl number of multiply-subtract pairs for LDL' factorization alpar@1: * dmax: the largest number of entries in any column of L, including the alpar@1: * diagonal alpar@1: * alpha: "dense" degree ratio alpar@1: * lnz: the number of nonzeros in L (excluding the diagonal) alpar@1: * lnzme: the number of nonzeros in L (excl. the diagonal) for the alpar@1: * current element me alpar@1: alpar@1: * ---------------------------------------------------------------------------- alpar@1: * LOCAL "POINTERS" (indices into the Iw array) alpar@1: * ---------------------------------------------------------------------------- alpar@1: */ alpar@1: alpar@1: Int p, p1, p2, p3, p4, pdst, pend, pj, pme, pme1, pme2, pn, psrc ; alpar@1: alpar@1: /* alpar@1: * Any parameter (Pe [...] or pfree) or local variable starting with "p" (for alpar@1: * Pointer) is an index into Iw, and all indices into Iw use variables starting alpar@1: * with "p." The only exception to this rule is the iwlen input argument. alpar@1: * alpar@1: * p: pointer into lots of things alpar@1: * p1: Pe [i] for some variable i (start of element list) alpar@1: * p2: Pe [i] + Elen [i] - 1 for some variable i alpar@1: * p3: index of first supervariable in clean list alpar@1: * p4: alpar@1: * pdst: destination pointer, for compression alpar@1: * pend: end of memory to compress alpar@1: * pj: pointer into an element or variable alpar@1: * pme: pointer into the current element (pme1...pme2) alpar@1: * pme1: the current element, me, is stored in Iw [pme1...pme2] alpar@1: * pme2: the end of the current element alpar@1: * pn: pointer into a "clean" variable, also used to compress alpar@1: * psrc: source pointer, for compression alpar@1: */ alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* INITIALIZATIONS */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: /* Note that this restriction on iwlen is slightly more restrictive than alpar@1: * what is actually required in AMD_2. AMD_2 can operate with no elbow alpar@1: * room at all, but it will be slow. For better performance, at least alpar@1: * size-n elbow room is enforced. */ alpar@1: ASSERT (iwlen >= pfree + n) ; alpar@1: ASSERT (n > 0) ; alpar@1: alpar@1: /* initialize output statistics */ alpar@1: lnz = 0 ; alpar@1: ndiv = 0 ; alpar@1: nms_lu = 0 ; alpar@1: nms_ldl = 0 ; alpar@1: dmax = 1 ; alpar@1: me = EMPTY ; alpar@1: alpar@1: mindeg = 0 ; alpar@1: ncmpa = 0 ; alpar@1: nel = 0 ; alpar@1: lemax = 0 ; alpar@1: alpar@1: /* get control parameters */ alpar@1: if (Control != (double *) NULL) alpar@1: { alpar@1: alpha = Control [AMD_DENSE] ; alpar@1: aggressive = (Control [AMD_AGGRESSIVE] != 0) ; alpar@1: } alpar@1: else alpar@1: { alpar@1: alpha = AMD_DEFAULT_DENSE ; alpar@1: aggressive = AMD_DEFAULT_AGGRESSIVE ; alpar@1: } alpar@1: /* Note: if alpha is NaN, this is undefined: */ alpar@1: if (alpha < 0) alpar@1: { alpar@1: /* only remove completely dense rows/columns */ alpar@1: dense = n-2 ; alpar@1: } alpar@1: else alpar@1: { alpar@1: dense = alpha * sqrt ((double) n) ; alpar@1: } alpar@1: dense = MAX (16, dense) ; alpar@1: dense = MIN (n, dense) ; alpar@1: AMD_DEBUG1 (("\n\nAMD (debug), alpha %g, aggr. "ID"\n", alpar@1: alpha, aggressive)) ; alpar@1: alpar@1: for (i = 0 ; i < n ; i++) alpar@1: { alpar@1: Last [i] = EMPTY ; alpar@1: Head [i] = EMPTY ; alpar@1: Next [i] = EMPTY ; alpar@1: /* if separate Hhead array is used for hash buckets: * alpar@1: Hhead [i] = EMPTY ; alpar@1: */ alpar@1: Nv [i] = 1 ; alpar@1: W [i] = 1 ; alpar@1: Elen [i] = 0 ; alpar@1: Degree [i] = Len [i] ; alpar@1: } alpar@1: alpar@1: #ifndef NDEBUG alpar@1: AMD_DEBUG1 (("\n======Nel "ID" initial\n", nel)) ; alpar@1: AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next, Last, alpar@1: Head, Elen, Degree, W, -1) ; alpar@1: #endif alpar@1: alpar@1: /* initialize wflg */ alpar@1: wbig = Int_MAX - n ; alpar@1: wflg = clear_flag (0, wbig, W, n) ; alpar@1: alpar@1: /* --------------------------------------------------------------------- */ alpar@1: /* initialize degree lists and eliminate dense and empty rows */ alpar@1: /* --------------------------------------------------------------------- */ alpar@1: alpar@1: ndense = 0 ; alpar@1: alpar@1: for (i = 0 ; i < n ; i++) alpar@1: { alpar@1: deg = Degree [i] ; alpar@1: ASSERT (deg >= 0 && deg < n) ; alpar@1: if (deg == 0) alpar@1: { alpar@1: alpar@1: /* ------------------------------------------------------------- alpar@1: * we have a variable that can be eliminated at once because alpar@1: * there is no off-diagonal non-zero in its row. Note that alpar@1: * Nv [i] = 1 for an empty variable i. It is treated just alpar@1: * the same as an eliminated element i. alpar@1: * ------------------------------------------------------------- */ alpar@1: alpar@1: Elen [i] = FLIP (1) ; alpar@1: nel++ ; alpar@1: Pe [i] = EMPTY ; alpar@1: W [i] = 0 ; alpar@1: alpar@1: } alpar@1: else if (deg > dense) alpar@1: { alpar@1: alpar@1: /* ------------------------------------------------------------- alpar@1: * Dense variables are not treated as elements, but as unordered, alpar@1: * non-principal variables that have no parent. They do not take alpar@1: * part in the postorder, since Nv [i] = 0. Note that the Fortran alpar@1: * version does not have this option. alpar@1: * ------------------------------------------------------------- */ alpar@1: alpar@1: AMD_DEBUG1 (("Dense node "ID" degree "ID"\n", i, deg)) ; alpar@1: ndense++ ; alpar@1: Nv [i] = 0 ; /* do not postorder this node */ alpar@1: Elen [i] = EMPTY ; alpar@1: nel++ ; alpar@1: Pe [i] = EMPTY ; alpar@1: alpar@1: } alpar@1: else alpar@1: { alpar@1: alpar@1: /* ------------------------------------------------------------- alpar@1: * place i in the degree list corresponding to its degree alpar@1: * ------------------------------------------------------------- */ alpar@1: alpar@1: inext = Head [deg] ; alpar@1: ASSERT (inext >= EMPTY && inext < n) ; alpar@1: if (inext != EMPTY) Last [inext] = i ; alpar@1: Next [i] = inext ; alpar@1: Head [deg] = i ; alpar@1: alpar@1: } alpar@1: } alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* WHILE (selecting pivots) DO */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: while (nel < n) alpar@1: { alpar@1: alpar@1: #ifndef NDEBUG alpar@1: AMD_DEBUG1 (("\n======Nel "ID"\n", nel)) ; alpar@1: if (AMD_debug >= 2) alpar@1: { alpar@1: AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next, alpar@1: Last, Head, Elen, Degree, W, nel) ; alpar@1: } alpar@1: #endif alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* GET PIVOT OF MINIMUM DEGREE */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: /* ----------------------------------------------------------------- */ alpar@1: /* find next supervariable for elimination */ alpar@1: /* ----------------------------------------------------------------- */ alpar@1: alpar@1: ASSERT (mindeg >= 0 && mindeg < n) ; alpar@1: for (deg = mindeg ; deg < n ; deg++) alpar@1: { alpar@1: me = Head [deg] ; alpar@1: if (me != EMPTY) break ; alpar@1: } alpar@1: mindeg = deg ; alpar@1: ASSERT (me >= 0 && me < n) ; alpar@1: AMD_DEBUG1 (("=================me: "ID"\n", me)) ; alpar@1: alpar@1: /* ----------------------------------------------------------------- */ alpar@1: /* remove chosen variable from link list */ alpar@1: /* ----------------------------------------------------------------- */ alpar@1: alpar@1: inext = Next [me] ; alpar@1: ASSERT (inext >= EMPTY && inext < n) ; alpar@1: if (inext != EMPTY) Last [inext] = EMPTY ; alpar@1: Head [deg] = inext ; alpar@1: alpar@1: /* ----------------------------------------------------------------- */ alpar@1: /* me represents the elimination of pivots nel to nel+Nv[me]-1. */ alpar@1: /* place me itself as the first in this set. */ alpar@1: /* ----------------------------------------------------------------- */ alpar@1: alpar@1: elenme = Elen [me] ; alpar@1: nvpiv = Nv [me] ; alpar@1: ASSERT (nvpiv > 0) ; alpar@1: nel += nvpiv ; alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* CONSTRUCT NEW ELEMENT */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: /* ----------------------------------------------------------------- alpar@1: * At this point, me is the pivotal supervariable. It will be alpar@1: * converted into the current element. Scan list of the pivotal alpar@1: * supervariable, me, setting tree pointers and constructing new list alpar@1: * of supervariables for the new element, me. p is a pointer to the alpar@1: * current position in the old list. alpar@1: * ----------------------------------------------------------------- */ alpar@1: alpar@1: /* flag the variable "me" as being in Lme by negating Nv [me] */ alpar@1: Nv [me] = -nvpiv ; alpar@1: degme = 0 ; alpar@1: ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ; alpar@1: alpar@1: if (elenme == 0) alpar@1: { alpar@1: alpar@1: /* ------------------------------------------------------------- */ alpar@1: /* construct the new element in place */ alpar@1: /* ------------------------------------------------------------- */ alpar@1: alpar@1: pme1 = Pe [me] ; alpar@1: pme2 = pme1 - 1 ; alpar@1: alpar@1: for (p = pme1 ; p <= pme1 + Len [me] - 1 ; p++) alpar@1: { alpar@1: i = Iw [p] ; alpar@1: ASSERT (i >= 0 && i < n && Nv [i] >= 0) ; alpar@1: nvi = Nv [i] ; alpar@1: if (nvi > 0) alpar@1: { alpar@1: alpar@1: /* ----------------------------------------------------- */ alpar@1: /* i is a principal variable not yet placed in Lme. */ alpar@1: /* store i in new list */ alpar@1: /* ----------------------------------------------------- */ alpar@1: alpar@1: /* flag i as being in Lme by negating Nv [i] */ alpar@1: degme += nvi ; alpar@1: Nv [i] = -nvi ; alpar@1: Iw [++pme2] = i ; alpar@1: alpar@1: /* ----------------------------------------------------- */ alpar@1: /* remove variable i from degree list. */ alpar@1: /* ----------------------------------------------------- */ alpar@1: alpar@1: ilast = Last [i] ; alpar@1: inext = Next [i] ; alpar@1: ASSERT (ilast >= EMPTY && ilast < n) ; alpar@1: ASSERT (inext >= EMPTY && inext < n) ; alpar@1: if (inext != EMPTY) Last [inext] = ilast ; alpar@1: if (ilast != EMPTY) alpar@1: { alpar@1: Next [ilast] = inext ; alpar@1: } alpar@1: else alpar@1: { alpar@1: /* i is at the head of the degree list */ alpar@1: ASSERT (Degree [i] >= 0 && Degree [i] < n) ; alpar@1: Head [Degree [i]] = inext ; alpar@1: } alpar@1: } alpar@1: } alpar@1: } alpar@1: else alpar@1: { alpar@1: alpar@1: /* ------------------------------------------------------------- */ alpar@1: /* construct the new element in empty space, Iw [pfree ...] */ alpar@1: /* ------------------------------------------------------------- */ alpar@1: alpar@1: p = Pe [me] ; alpar@1: pme1 = pfree ; alpar@1: slenme = Len [me] - elenme ; alpar@1: alpar@1: for (knt1 = 1 ; knt1 <= elenme + 1 ; knt1++) alpar@1: { alpar@1: alpar@1: if (knt1 > elenme) alpar@1: { alpar@1: /* search the supervariables in me. */ alpar@1: e = me ; alpar@1: pj = p ; alpar@1: ln = slenme ; alpar@1: AMD_DEBUG2 (("Search sv: "ID" "ID" "ID"\n", me,pj,ln)) ; alpar@1: } alpar@1: else alpar@1: { alpar@1: /* search the elements in me. */ alpar@1: e = Iw [p++] ; alpar@1: ASSERT (e >= 0 && e < n) ; alpar@1: pj = Pe [e] ; alpar@1: ln = Len [e] ; alpar@1: AMD_DEBUG2 (("Search element e "ID" in me "ID"\n", e,me)) ; alpar@1: ASSERT (Elen [e] < EMPTY && W [e] > 0 && pj >= 0) ; alpar@1: } alpar@1: ASSERT (ln >= 0 && (ln == 0 || (pj >= 0 && pj < iwlen))) ; alpar@1: alpar@1: /* --------------------------------------------------------- alpar@1: * search for different supervariables and add them to the alpar@1: * new list, compressing when necessary. this loop is alpar@1: * executed once for each element in the list and once for alpar@1: * all the supervariables in the list. alpar@1: * --------------------------------------------------------- */ alpar@1: alpar@1: for (knt2 = 1 ; knt2 <= ln ; knt2++) alpar@1: { alpar@1: i = Iw [pj++] ; alpar@1: ASSERT (i >= 0 && i < n && (i == me || Elen [i] >= EMPTY)); alpar@1: nvi = Nv [i] ; alpar@1: AMD_DEBUG2 ((": "ID" "ID" "ID" "ID"\n", alpar@1: i, Elen [i], Nv [i], wflg)) ; alpar@1: alpar@1: if (nvi > 0) alpar@1: { alpar@1: alpar@1: /* ------------------------------------------------- */ alpar@1: /* compress Iw, if necessary */ alpar@1: /* ------------------------------------------------- */ alpar@1: alpar@1: if (pfree >= iwlen) alpar@1: { alpar@1: alpar@1: AMD_DEBUG1 (("GARBAGE COLLECTION\n")) ; alpar@1: alpar@1: /* prepare for compressing Iw by adjusting pointers alpar@1: * and lengths so that the lists being searched in alpar@1: * the inner and outer loops contain only the alpar@1: * remaining entries. */ alpar@1: alpar@1: Pe [me] = p ; alpar@1: Len [me] -= knt1 ; alpar@1: /* check if nothing left of supervariable me */ alpar@1: if (Len [me] == 0) Pe [me] = EMPTY ; alpar@1: Pe [e] = pj ; alpar@1: Len [e] = ln - knt2 ; alpar@1: /* nothing left of element e */ alpar@1: if (Len [e] == 0) Pe [e] = EMPTY ; alpar@1: alpar@1: ncmpa++ ; /* one more garbage collection */ alpar@1: alpar@1: /* store first entry of each object in Pe */ alpar@1: /* FLIP the first entry in each object */ alpar@1: for (j = 0 ; j < n ; j++) alpar@1: { alpar@1: pn = Pe [j] ; alpar@1: if (pn >= 0) alpar@1: { alpar@1: ASSERT (pn >= 0 && pn < iwlen) ; alpar@1: Pe [j] = Iw [pn] ; alpar@1: Iw [pn] = FLIP (j) ; alpar@1: } alpar@1: } alpar@1: alpar@1: /* psrc/pdst point to source/destination */ alpar@1: psrc = 0 ; alpar@1: pdst = 0 ; alpar@1: pend = pme1 - 1 ; alpar@1: alpar@1: while (psrc <= pend) alpar@1: { alpar@1: /* search for next FLIP'd entry */ alpar@1: j = FLIP (Iw [psrc++]) ; alpar@1: if (j >= 0) alpar@1: { alpar@1: AMD_DEBUG2 (("Got object j: "ID"\n", j)) ; alpar@1: Iw [pdst] = Pe [j] ; alpar@1: Pe [j] = pdst++ ; alpar@1: lenj = Len [j] ; alpar@1: /* copy from source to destination */ alpar@1: for (knt3 = 0 ; knt3 <= lenj - 2 ; knt3++) alpar@1: { alpar@1: Iw [pdst++] = Iw [psrc++] ; alpar@1: } alpar@1: } alpar@1: } alpar@1: alpar@1: /* move the new partially-constructed element */ alpar@1: p1 = pdst ; alpar@1: for (psrc = pme1 ; psrc <= pfree-1 ; psrc++) alpar@1: { alpar@1: Iw [pdst++] = Iw [psrc] ; alpar@1: } alpar@1: pme1 = p1 ; alpar@1: pfree = pdst ; alpar@1: pj = Pe [e] ; alpar@1: p = Pe [me] ; alpar@1: alpar@1: } alpar@1: alpar@1: /* ------------------------------------------------- */ alpar@1: /* i is a principal variable not yet placed in Lme */ alpar@1: /* store i in new list */ alpar@1: /* ------------------------------------------------- */ alpar@1: alpar@1: /* flag i as being in Lme by negating Nv [i] */ alpar@1: degme += nvi ; alpar@1: Nv [i] = -nvi ; alpar@1: Iw [pfree++] = i ; alpar@1: AMD_DEBUG2 ((" s: "ID" nv "ID"\n", i, Nv [i])); alpar@1: alpar@1: /* ------------------------------------------------- */ alpar@1: /* remove variable i from degree link list */ alpar@1: /* ------------------------------------------------- */ alpar@1: alpar@1: ilast = Last [i] ; alpar@1: inext = Next [i] ; alpar@1: ASSERT (ilast >= EMPTY && ilast < n) ; alpar@1: ASSERT (inext >= EMPTY && inext < n) ; alpar@1: if (inext != EMPTY) Last [inext] = ilast ; alpar@1: if (ilast != EMPTY) alpar@1: { alpar@1: Next [ilast] = inext ; alpar@1: } alpar@1: else alpar@1: { alpar@1: /* i is at the head of the degree list */ alpar@1: ASSERT (Degree [i] >= 0 && Degree [i] < n) ; alpar@1: Head [Degree [i]] = inext ; alpar@1: } alpar@1: } alpar@1: } alpar@1: alpar@1: if (e != me) alpar@1: { alpar@1: /* set tree pointer and flag to indicate element e is alpar@1: * absorbed into new element me (the parent of e is me) */ alpar@1: AMD_DEBUG1 ((" Element "ID" => "ID"\n", e, me)) ; alpar@1: Pe [e] = FLIP (me) ; alpar@1: W [e] = 0 ; alpar@1: } alpar@1: } alpar@1: alpar@1: pme2 = pfree - 1 ; alpar@1: } alpar@1: alpar@1: /* ----------------------------------------------------------------- */ alpar@1: /* me has now been converted into an element in Iw [pme1..pme2] */ alpar@1: /* ----------------------------------------------------------------- */ alpar@1: alpar@1: /* degme holds the external degree of new element */ alpar@1: Degree [me] = degme ; alpar@1: Pe [me] = pme1 ; alpar@1: Len [me] = pme2 - pme1 + 1 ; alpar@1: ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ; alpar@1: alpar@1: Elen [me] = FLIP (nvpiv + degme) ; alpar@1: /* FLIP (Elen (me)) is now the degree of pivot (including alpar@1: * diagonal part). */ alpar@1: alpar@1: #ifndef NDEBUG alpar@1: AMD_DEBUG2 (("New element structure: length= "ID"\n", pme2-pme1+1)) ; alpar@1: for (pme = pme1 ; pme <= pme2 ; pme++) AMD_DEBUG3 ((" "ID"", Iw[pme])); alpar@1: AMD_DEBUG3 (("\n")) ; alpar@1: #endif alpar@1: alpar@1: /* ----------------------------------------------------------------- */ alpar@1: /* make sure that wflg is not too large. */ alpar@1: /* ----------------------------------------------------------------- */ alpar@1: alpar@1: /* With the current value of wflg, wflg+n must not cause integer alpar@1: * overflow */ alpar@1: alpar@1: wflg = clear_flag (wflg, wbig, W, n) ; alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* COMPUTE (W [e] - wflg) = |Le\Lme| FOR ALL ELEMENTS */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: /* ----------------------------------------------------------------- alpar@1: * Scan 1: compute the external degrees of previous elements with alpar@1: * respect to the current element. That is: alpar@1: * (W [e] - wflg) = |Le \ Lme| alpar@1: * for each element e that appears in any supervariable in Lme. The alpar@1: * notation Le refers to the pattern (list of supervariables) of a alpar@1: * previous element e, where e is not yet absorbed, stored in alpar@1: * Iw [Pe [e] + 1 ... Pe [e] + Len [e]]. The notation Lme alpar@1: * refers to the pattern of the current element (stored in alpar@1: * Iw [pme1..pme2]). If aggressive absorption is enabled, and alpar@1: * (W [e] - wflg) becomes zero, then the element e will be absorbed alpar@1: * in Scan 2. alpar@1: * ----------------------------------------------------------------- */ alpar@1: alpar@1: AMD_DEBUG2 (("me: ")) ; alpar@1: for (pme = pme1 ; pme <= pme2 ; pme++) alpar@1: { alpar@1: i = Iw [pme] ; alpar@1: ASSERT (i >= 0 && i < n) ; alpar@1: eln = Elen [i] ; alpar@1: AMD_DEBUG3 ((""ID" Elen "ID": \n", i, eln)) ; alpar@1: if (eln > 0) alpar@1: { alpar@1: /* note that Nv [i] has been negated to denote i in Lme: */ alpar@1: nvi = -Nv [i] ; alpar@1: ASSERT (nvi > 0 && Pe [i] >= 0 && Pe [i] < iwlen) ; alpar@1: wnvi = wflg - nvi ; alpar@1: for (p = Pe [i] ; p <= Pe [i] + eln - 1 ; p++) alpar@1: { alpar@1: e = Iw [p] ; alpar@1: ASSERT (e >= 0 && e < n) ; alpar@1: we = W [e] ; alpar@1: AMD_DEBUG4 ((" e "ID" we "ID" ", e, we)) ; alpar@1: if (we >= wflg) alpar@1: { alpar@1: /* unabsorbed element e has been seen in this loop */ alpar@1: AMD_DEBUG4 ((" unabsorbed, first time seen")) ; alpar@1: we -= nvi ; alpar@1: } alpar@1: else if (we != 0) alpar@1: { alpar@1: /* e is an unabsorbed element */ alpar@1: /* this is the first we have seen e in all of Scan 1 */ alpar@1: AMD_DEBUG4 ((" unabsorbed")) ; alpar@1: we = Degree [e] + wnvi ; alpar@1: } alpar@1: AMD_DEBUG4 (("\n")) ; alpar@1: W [e] = we ; alpar@1: } alpar@1: } alpar@1: } alpar@1: AMD_DEBUG2 (("\n")) ; alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* DEGREE UPDATE AND ELEMENT ABSORPTION */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: /* ----------------------------------------------------------------- alpar@1: * Scan 2: for each i in Lme, sum up the degree of Lme (which is alpar@1: * degme), plus the sum of the external degrees of each Le for the alpar@1: * elements e appearing within i, plus the supervariables in i. alpar@1: * Place i in hash list. alpar@1: * ----------------------------------------------------------------- */ alpar@1: alpar@1: for (pme = pme1 ; pme <= pme2 ; pme++) alpar@1: { alpar@1: i = Iw [pme] ; alpar@1: ASSERT (i >= 0 && i < n && Nv [i] < 0 && Elen [i] >= 0) ; alpar@1: AMD_DEBUG2 (("Updating: i "ID" "ID" "ID"\n", i, Elen[i], Len [i])); alpar@1: p1 = Pe [i] ; alpar@1: p2 = p1 + Elen [i] - 1 ; alpar@1: pn = p1 ; alpar@1: hash = 0 ; alpar@1: deg = 0 ; alpar@1: ASSERT (p1 >= 0 && p1 < iwlen && p2 >= -1 && p2 < iwlen) ; alpar@1: alpar@1: /* ------------------------------------------------------------- */ alpar@1: /* scan the element list associated with supervariable i */ alpar@1: /* ------------------------------------------------------------- */ alpar@1: alpar@1: /* UMFPACK/MA38-style approximate degree: */ alpar@1: if (aggressive) alpar@1: { alpar@1: for (p = p1 ; p <= p2 ; p++) alpar@1: { alpar@1: e = Iw [p] ; alpar@1: ASSERT (e >= 0 && e < n) ; alpar@1: we = W [e] ; alpar@1: if (we != 0) alpar@1: { alpar@1: /* e is an unabsorbed element */ alpar@1: /* dext = | Le \ Lme | */ alpar@1: dext = we - wflg ; alpar@1: if (dext > 0) alpar@1: { alpar@1: deg += dext ; alpar@1: Iw [pn++] = e ; alpar@1: hash += e ; alpar@1: AMD_DEBUG4 ((" e: "ID" hash = "ID"\n",e,hash)) ; alpar@1: } alpar@1: else alpar@1: { alpar@1: /* external degree of e is zero, absorb e into me*/ alpar@1: AMD_DEBUG1 ((" Element "ID" =>"ID" (aggressive)\n", alpar@1: e, me)) ; alpar@1: ASSERT (dext == 0) ; alpar@1: Pe [e] = FLIP (me) ; alpar@1: W [e] = 0 ; alpar@1: } alpar@1: } alpar@1: } alpar@1: } alpar@1: else alpar@1: { alpar@1: for (p = p1 ; p <= p2 ; p++) alpar@1: { alpar@1: e = Iw [p] ; alpar@1: ASSERT (e >= 0 && e < n) ; alpar@1: we = W [e] ; alpar@1: if (we != 0) alpar@1: { alpar@1: /* e is an unabsorbed element */ alpar@1: dext = we - wflg ; alpar@1: ASSERT (dext >= 0) ; alpar@1: deg += dext ; alpar@1: Iw [pn++] = e ; alpar@1: hash += e ; alpar@1: AMD_DEBUG4 ((" e: "ID" hash = "ID"\n",e,hash)) ; alpar@1: } alpar@1: } alpar@1: } alpar@1: alpar@1: /* count the number of elements in i (including me): */ alpar@1: Elen [i] = pn - p1 + 1 ; alpar@1: alpar@1: /* ------------------------------------------------------------- */ alpar@1: /* scan the supervariables in the list associated with i */ alpar@1: /* ------------------------------------------------------------- */ alpar@1: alpar@1: /* The bulk of the AMD run time is typically spent in this loop, alpar@1: * particularly if the matrix has many dense rows that are not alpar@1: * removed prior to ordering. */ alpar@1: p3 = pn ; alpar@1: p4 = p1 + Len [i] ; alpar@1: for (p = p2 + 1 ; p < p4 ; p++) alpar@1: { alpar@1: j = Iw [p] ; alpar@1: ASSERT (j >= 0 && j < n) ; alpar@1: nvj = Nv [j] ; alpar@1: if (nvj > 0) alpar@1: { alpar@1: /* j is unabsorbed, and not in Lme. */ alpar@1: /* add to degree and add to new list */ alpar@1: deg += nvj ; alpar@1: Iw [pn++] = j ; alpar@1: hash += j ; alpar@1: AMD_DEBUG4 ((" s: "ID" hash "ID" Nv[j]= "ID"\n", alpar@1: j, hash, nvj)) ; alpar@1: } alpar@1: } alpar@1: alpar@1: /* ------------------------------------------------------------- */ alpar@1: /* update the degree and check for mass elimination */ alpar@1: /* ------------------------------------------------------------- */ alpar@1: alpar@1: /* with aggressive absorption, deg==0 is identical to the alpar@1: * Elen [i] == 1 && p3 == pn test, below. */ alpar@1: ASSERT (IMPLIES (aggressive, (deg==0) == (Elen[i]==1 && p3==pn))) ; alpar@1: alpar@1: if (Elen [i] == 1 && p3 == pn) alpar@1: { alpar@1: alpar@1: /* --------------------------------------------------------- */ alpar@1: /* mass elimination */ alpar@1: /* --------------------------------------------------------- */ alpar@1: alpar@1: /* There is nothing left of this node except for an edge to alpar@1: * the current pivot element. Elen [i] is 1, and there are alpar@1: * no variables adjacent to node i. Absorb i into the alpar@1: * current pivot element, me. Note that if there are two or alpar@1: * more mass eliminations, fillin due to mass elimination is alpar@1: * possible within the nvpiv-by-nvpiv pivot block. It is this alpar@1: * step that causes AMD's analysis to be an upper bound. alpar@1: * alpar@1: * The reason is that the selected pivot has a lower alpar@1: * approximate degree than the true degree of the two mass alpar@1: * eliminated nodes. There is no edge between the two mass alpar@1: * eliminated nodes. They are merged with the current pivot alpar@1: * anyway. alpar@1: * alpar@1: * No fillin occurs in the Schur complement, in any case, alpar@1: * and this effect does not decrease the quality of the alpar@1: * ordering itself, just the quality of the nonzero and alpar@1: * flop count analysis. It also means that the post-ordering alpar@1: * is not an exact elimination tree post-ordering. */ alpar@1: alpar@1: AMD_DEBUG1 ((" MASS i "ID" => parent e "ID"\n", i, me)) ; alpar@1: Pe [i] = FLIP (me) ; alpar@1: nvi = -Nv [i] ; alpar@1: degme -= nvi ; alpar@1: nvpiv += nvi ; alpar@1: nel += nvi ; alpar@1: Nv [i] = 0 ; alpar@1: Elen [i] = EMPTY ; alpar@1: alpar@1: } alpar@1: else alpar@1: { alpar@1: alpar@1: /* --------------------------------------------------------- */ alpar@1: /* update the upper-bound degree of i */ alpar@1: /* --------------------------------------------------------- */ alpar@1: alpar@1: /* the following degree does not yet include the size alpar@1: * of the current element, which is added later: */ alpar@1: alpar@1: Degree [i] = MIN (Degree [i], deg) ; alpar@1: alpar@1: /* --------------------------------------------------------- */ alpar@1: /* add me to the list for i */ alpar@1: /* --------------------------------------------------------- */ alpar@1: alpar@1: /* move first supervariable to end of list */ alpar@1: Iw [pn] = Iw [p3] ; alpar@1: /* move first element to end of element part of list */ alpar@1: Iw [p3] = Iw [p1] ; alpar@1: /* add new element, me, to front of list. */ alpar@1: Iw [p1] = me ; alpar@1: /* store the new length of the list in Len [i] */ alpar@1: Len [i] = pn - p1 + 1 ; alpar@1: alpar@1: /* --------------------------------------------------------- */ alpar@1: /* place in hash bucket. Save hash key of i in Last [i]. */ alpar@1: /* --------------------------------------------------------- */ alpar@1: alpar@1: /* NOTE: this can fail if hash is negative, because the ANSI C alpar@1: * standard does not define a % b when a and/or b are negative. alpar@1: * That's why hash is defined as an unsigned Int, to avoid this alpar@1: * problem. */ alpar@1: hash = hash % n ; alpar@1: ASSERT (((Int) hash) >= 0 && ((Int) hash) < n) ; alpar@1: alpar@1: /* if the Hhead array is not used: */ alpar@1: j = Head [hash] ; alpar@1: if (j <= EMPTY) alpar@1: { alpar@1: /* degree list is empty, hash head is FLIP (j) */ alpar@1: Next [i] = FLIP (j) ; alpar@1: Head [hash] = FLIP (i) ; alpar@1: } alpar@1: else alpar@1: { alpar@1: /* degree list is not empty, use Last [Head [hash]] as alpar@1: * hash head. */ alpar@1: Next [i] = Last [j] ; alpar@1: Last [j] = i ; alpar@1: } alpar@1: alpar@1: /* if a separate Hhead array is used: * alpar@1: Next [i] = Hhead [hash] ; alpar@1: Hhead [hash] = i ; alpar@1: */ alpar@1: alpar@1: Last [i] = hash ; alpar@1: } alpar@1: } alpar@1: alpar@1: Degree [me] = degme ; alpar@1: alpar@1: /* ----------------------------------------------------------------- */ alpar@1: /* Clear the counter array, W [...], by incrementing wflg. */ alpar@1: /* ----------------------------------------------------------------- */ alpar@1: alpar@1: /* make sure that wflg+n does not cause integer overflow */ alpar@1: lemax = MAX (lemax, degme) ; alpar@1: wflg += lemax ; alpar@1: wflg = clear_flag (wflg, wbig, W, n) ; alpar@1: /* at this point, W [0..n-1] < wflg holds */ alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* SUPERVARIABLE DETECTION */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: AMD_DEBUG1 (("Detecting supervariables:\n")) ; alpar@1: for (pme = pme1 ; pme <= pme2 ; pme++) alpar@1: { alpar@1: i = Iw [pme] ; alpar@1: ASSERT (i >= 0 && i < n) ; alpar@1: AMD_DEBUG2 (("Consider i "ID" nv "ID"\n", i, Nv [i])) ; alpar@1: if (Nv [i] < 0) alpar@1: { alpar@1: /* i is a principal variable in Lme */ alpar@1: alpar@1: /* --------------------------------------------------------- alpar@1: * examine all hash buckets with 2 or more variables. We do alpar@1: * this by examing all unique hash keys for supervariables in alpar@1: * the pattern Lme of the current element, me alpar@1: * --------------------------------------------------------- */ alpar@1: alpar@1: /* let i = head of hash bucket, and empty the hash bucket */ alpar@1: ASSERT (Last [i] >= 0 && Last [i] < n) ; alpar@1: hash = Last [i] ; alpar@1: alpar@1: /* if Hhead array is not used: */ alpar@1: j = Head [hash] ; alpar@1: if (j == EMPTY) alpar@1: { alpar@1: /* hash bucket and degree list are both empty */ alpar@1: i = EMPTY ; alpar@1: } alpar@1: else if (j < EMPTY) alpar@1: { alpar@1: /* degree list is empty */ alpar@1: i = FLIP (j) ; alpar@1: Head [hash] = EMPTY ; alpar@1: } alpar@1: else alpar@1: { alpar@1: /* degree list is not empty, restore Last [j] of head j */ alpar@1: i = Last [j] ; alpar@1: Last [j] = EMPTY ; alpar@1: } alpar@1: alpar@1: /* if separate Hhead array is used: * alpar@1: i = Hhead [hash] ; alpar@1: Hhead [hash] = EMPTY ; alpar@1: */ alpar@1: alpar@1: ASSERT (i >= EMPTY && i < n) ; alpar@1: AMD_DEBUG2 (("----i "ID" hash "ID"\n", i, hash)) ; alpar@1: alpar@1: while (i != EMPTY && Next [i] != EMPTY) alpar@1: { alpar@1: alpar@1: /* ----------------------------------------------------- alpar@1: * this bucket has one or more variables following i. alpar@1: * scan all of them to see if i can absorb any entries alpar@1: * that follow i in hash bucket. Scatter i into w. alpar@1: * ----------------------------------------------------- */ alpar@1: alpar@1: ln = Len [i] ; alpar@1: eln = Elen [i] ; alpar@1: ASSERT (ln >= 0 && eln >= 0) ; alpar@1: ASSERT (Pe [i] >= 0 && Pe [i] < iwlen) ; alpar@1: /* do not flag the first element in the list (me) */ alpar@1: for (p = Pe [i] + 1 ; p <= Pe [i] + ln - 1 ; p++) alpar@1: { alpar@1: ASSERT (Iw [p] >= 0 && Iw [p] < n) ; alpar@1: W [Iw [p]] = wflg ; alpar@1: } alpar@1: alpar@1: /* ----------------------------------------------------- */ alpar@1: /* scan every other entry j following i in bucket */ alpar@1: /* ----------------------------------------------------- */ alpar@1: alpar@1: jlast = i ; alpar@1: j = Next [i] ; alpar@1: ASSERT (j >= EMPTY && j < n) ; alpar@1: alpar@1: while (j != EMPTY) alpar@1: { alpar@1: /* ------------------------------------------------- */ alpar@1: /* check if j and i have identical nonzero pattern */ alpar@1: /* ------------------------------------------------- */ alpar@1: alpar@1: AMD_DEBUG3 (("compare i "ID" and j "ID"\n", i,j)) ; alpar@1: alpar@1: /* check if i and j have the same Len and Elen */ alpar@1: ASSERT (Len [j] >= 0 && Elen [j] >= 0) ; alpar@1: ASSERT (Pe [j] >= 0 && Pe [j] < iwlen) ; alpar@1: ok = (Len [j] == ln) && (Elen [j] == eln) ; alpar@1: /* skip the first element in the list (me) */ alpar@1: for (p = Pe [j] + 1 ; ok && p <= Pe [j] + ln - 1 ; p++) alpar@1: { alpar@1: ASSERT (Iw [p] >= 0 && Iw [p] < n) ; alpar@1: if (W [Iw [p]] != wflg) ok = 0 ; alpar@1: } alpar@1: if (ok) alpar@1: { alpar@1: /* --------------------------------------------- */ alpar@1: /* found it! j can be absorbed into i */ alpar@1: /* --------------------------------------------- */ alpar@1: alpar@1: AMD_DEBUG1 (("found it! j "ID" => i "ID"\n", j,i)); alpar@1: Pe [j] = FLIP (i) ; alpar@1: /* both Nv [i] and Nv [j] are negated since they */ alpar@1: /* are in Lme, and the absolute values of each */ alpar@1: /* are the number of variables in i and j: */ alpar@1: Nv [i] += Nv [j] ; alpar@1: Nv [j] = 0 ; alpar@1: Elen [j] = EMPTY ; alpar@1: /* delete j from hash bucket */ alpar@1: ASSERT (j != Next [j]) ; alpar@1: j = Next [j] ; alpar@1: Next [jlast] = j ; alpar@1: alpar@1: } alpar@1: else alpar@1: { alpar@1: /* j cannot be absorbed into i */ alpar@1: jlast = j ; alpar@1: ASSERT (j != Next [j]) ; alpar@1: j = Next [j] ; alpar@1: } alpar@1: ASSERT (j >= EMPTY && j < n) ; alpar@1: } alpar@1: alpar@1: /* ----------------------------------------------------- alpar@1: * no more variables can be absorbed into i alpar@1: * go to next i in bucket and clear flag array alpar@1: * ----------------------------------------------------- */ alpar@1: alpar@1: wflg++ ; alpar@1: i = Next [i] ; alpar@1: ASSERT (i >= EMPTY && i < n) ; alpar@1: alpar@1: } alpar@1: } alpar@1: } alpar@1: AMD_DEBUG2 (("detect done\n")) ; alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVARIABLES FROM ELEMENT */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: p = pme1 ; alpar@1: nleft = n - nel ; alpar@1: for (pme = pme1 ; pme <= pme2 ; pme++) alpar@1: { alpar@1: i = Iw [pme] ; alpar@1: ASSERT (i >= 0 && i < n) ; alpar@1: nvi = -Nv [i] ; alpar@1: AMD_DEBUG3 (("Restore i "ID" "ID"\n", i, nvi)) ; alpar@1: if (nvi > 0) alpar@1: { alpar@1: /* i is a principal variable in Lme */ alpar@1: /* restore Nv [i] to signify that i is principal */ alpar@1: Nv [i] = nvi ; alpar@1: alpar@1: /* --------------------------------------------------------- */ alpar@1: /* compute the external degree (add size of current element) */ alpar@1: /* --------------------------------------------------------- */ alpar@1: alpar@1: deg = Degree [i] + degme - nvi ; alpar@1: deg = MIN (deg, nleft - nvi) ; alpar@1: ASSERT (IMPLIES (aggressive, deg > 0) && deg >= 0 && deg < n) ; alpar@1: alpar@1: /* --------------------------------------------------------- */ alpar@1: /* place the supervariable at the head of the degree list */ alpar@1: /* --------------------------------------------------------- */ alpar@1: alpar@1: inext = Head [deg] ; alpar@1: ASSERT (inext >= EMPTY && inext < n) ; alpar@1: if (inext != EMPTY) Last [inext] = i ; alpar@1: Next [i] = inext ; alpar@1: Last [i] = EMPTY ; alpar@1: Head [deg] = i ; alpar@1: alpar@1: /* --------------------------------------------------------- */ alpar@1: /* save the new degree, and find the minimum degree */ alpar@1: /* --------------------------------------------------------- */ alpar@1: alpar@1: mindeg = MIN (mindeg, deg) ; alpar@1: Degree [i] = deg ; alpar@1: alpar@1: /* --------------------------------------------------------- */ alpar@1: /* place the supervariable in the element pattern */ alpar@1: /* --------------------------------------------------------- */ alpar@1: alpar@1: Iw [p++] = i ; alpar@1: alpar@1: } alpar@1: } alpar@1: AMD_DEBUG2 (("restore done\n")) ; alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* FINALIZE THE NEW ELEMENT */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: AMD_DEBUG2 (("ME = "ID" DONE\n", me)) ; alpar@1: Nv [me] = nvpiv ; alpar@1: /* save the length of the list for the new element me */ alpar@1: Len [me] = p - pme1 ; alpar@1: if (Len [me] == 0) alpar@1: { alpar@1: /* there is nothing left of the current pivot element */ alpar@1: /* it is a root of the assembly tree */ alpar@1: Pe [me] = EMPTY ; alpar@1: W [me] = 0 ; alpar@1: } alpar@1: if (elenme != 0) alpar@1: { alpar@1: /* element was not constructed in place: deallocate part of */ alpar@1: /* it since newly nonprincipal variables may have been removed */ alpar@1: pfree = p ; alpar@1: } alpar@1: alpar@1: /* The new element has nvpiv pivots and the size of the contribution alpar@1: * block for a multifrontal method is degme-by-degme, not including alpar@1: * the "dense" rows/columns. If the "dense" rows/columns are included, alpar@1: * the frontal matrix is no larger than alpar@1: * (degme+ndense)-by-(degme+ndense). alpar@1: */ alpar@1: alpar@1: if (Info != (double *) NULL) alpar@1: { alpar@1: f = nvpiv ; alpar@1: r = degme + ndense ; alpar@1: dmax = MAX (dmax, f + r) ; alpar@1: alpar@1: /* number of nonzeros in L (excluding the diagonal) */ alpar@1: lnzme = f*r + (f-1)*f/2 ; alpar@1: lnz += lnzme ; alpar@1: alpar@1: /* number of divide operations for LDL' and for LU */ alpar@1: ndiv += lnzme ; alpar@1: alpar@1: /* number of multiply-subtract pairs for LU */ alpar@1: s = f*r*r + r*(f-1)*f + (f-1)*f*(2*f-1)/6 ; alpar@1: nms_lu += s ; alpar@1: alpar@1: /* number of multiply-subtract pairs for LDL' */ alpar@1: nms_ldl += (s + lnzme)/2 ; alpar@1: } alpar@1: alpar@1: #ifndef NDEBUG alpar@1: AMD_DEBUG2 (("finalize done nel "ID" n "ID"\n ::::\n", nel, n)) ; alpar@1: for (pme = Pe [me] ; pme <= Pe [me] + Len [me] - 1 ; pme++) alpar@1: { alpar@1: AMD_DEBUG3 ((" "ID"", Iw [pme])) ; alpar@1: } alpar@1: AMD_DEBUG3 (("\n")) ; alpar@1: #endif alpar@1: alpar@1: } alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* DONE SELECTING PIVOTS */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: if (Info != (double *) NULL) alpar@1: { alpar@1: alpar@1: /* count the work to factorize the ndense-by-ndense submatrix */ alpar@1: f = ndense ; alpar@1: dmax = MAX (dmax, (double) ndense) ; alpar@1: alpar@1: /* number of nonzeros in L (excluding the diagonal) */ alpar@1: lnzme = (f-1)*f/2 ; alpar@1: lnz += lnzme ; alpar@1: alpar@1: /* number of divide operations for LDL' and for LU */ alpar@1: ndiv += lnzme ; alpar@1: alpar@1: /* number of multiply-subtract pairs for LU */ alpar@1: s = (f-1)*f*(2*f-1)/6 ; alpar@1: nms_lu += s ; alpar@1: alpar@1: /* number of multiply-subtract pairs for LDL' */ alpar@1: nms_ldl += (s + lnzme)/2 ; alpar@1: alpar@1: /* number of nz's in L (excl. diagonal) */ alpar@1: Info [AMD_LNZ] = lnz ; alpar@1: alpar@1: /* number of divide ops for LU and LDL' */ alpar@1: Info [AMD_NDIV] = ndiv ; alpar@1: alpar@1: /* number of multiply-subtract pairs for LDL' */ alpar@1: Info [AMD_NMULTSUBS_LDL] = nms_ldl ; alpar@1: alpar@1: /* number of multiply-subtract pairs for LU */ alpar@1: Info [AMD_NMULTSUBS_LU] = nms_lu ; alpar@1: alpar@1: /* number of "dense" rows/columns */ alpar@1: Info [AMD_NDENSE] = ndense ; alpar@1: alpar@1: /* largest front is dmax-by-dmax */ alpar@1: Info [AMD_DMAX] = dmax ; alpar@1: alpar@1: /* number of garbage collections in AMD */ alpar@1: Info [AMD_NCMPA] = ncmpa ; alpar@1: alpar@1: /* successful ordering */ alpar@1: Info [AMD_STATUS] = AMD_OK ; alpar@1: } alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* POST-ORDERING */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: /* ------------------------------------------------------------------------- alpar@1: * Variables at this point: alpar@1: * alpar@1: * Pe: holds the elimination tree. The parent of j is FLIP (Pe [j]), alpar@1: * or EMPTY if j is a root. The tree holds both elements and alpar@1: * non-principal (unordered) variables absorbed into them. alpar@1: * Dense variables are non-principal and unordered. alpar@1: * alpar@1: * Elen: holds the size of each element, including the diagonal part. alpar@1: * FLIP (Elen [e]) > 0 if e is an element. For unordered alpar@1: * variables i, Elen [i] is EMPTY. alpar@1: * alpar@1: * Nv: Nv [e] > 0 is the number of pivots represented by the element e. alpar@1: * For unordered variables i, Nv [i] is zero. alpar@1: * alpar@1: * Contents no longer needed: alpar@1: * W, Iw, Len, Degree, Head, Next, Last. alpar@1: * alpar@1: * The matrix itself has been destroyed. alpar@1: * alpar@1: * n: the size of the matrix. alpar@1: * No other scalars needed (pfree, iwlen, etc.) alpar@1: * ------------------------------------------------------------------------- */ alpar@1: alpar@1: /* restore Pe */ alpar@1: for (i = 0 ; i < n ; i++) alpar@1: { alpar@1: Pe [i] = FLIP (Pe [i]) ; alpar@1: } alpar@1: alpar@1: /* restore Elen, for output information, and for postordering */ alpar@1: for (i = 0 ; i < n ; i++) alpar@1: { alpar@1: Elen [i] = FLIP (Elen [i]) ; alpar@1: } alpar@1: alpar@1: /* Now the parent of j is Pe [j], or EMPTY if j is a root. Elen [e] > 0 alpar@1: * is the size of element e. Elen [i] is EMPTY for unordered variable i. */ alpar@1: alpar@1: #ifndef NDEBUG alpar@1: AMD_DEBUG2 (("\nTree:\n")) ; alpar@1: for (i = 0 ; i < n ; i++) alpar@1: { alpar@1: AMD_DEBUG2 ((" "ID" parent: "ID" ", i, Pe [i])) ; alpar@1: ASSERT (Pe [i] >= EMPTY && Pe [i] < n) ; alpar@1: if (Nv [i] > 0) alpar@1: { alpar@1: /* this is an element */ alpar@1: e = i ; alpar@1: AMD_DEBUG2 ((" element, size is "ID"\n", Elen [i])) ; alpar@1: ASSERT (Elen [e] > 0) ; alpar@1: } alpar@1: AMD_DEBUG2 (("\n")) ; alpar@1: } alpar@1: AMD_DEBUG2 (("\nelements:\n")) ; alpar@1: for (e = 0 ; e < n ; e++) alpar@1: { alpar@1: if (Nv [e] > 0) alpar@1: { alpar@1: AMD_DEBUG3 (("Element e= "ID" size "ID" nv "ID" \n", e, alpar@1: Elen [e], Nv [e])) ; alpar@1: } alpar@1: } alpar@1: AMD_DEBUG2 (("\nvariables:\n")) ; alpar@1: for (i = 0 ; i < n ; i++) alpar@1: { alpar@1: Int cnt ; alpar@1: if (Nv [i] == 0) alpar@1: { alpar@1: AMD_DEBUG3 (("i unordered: "ID"\n", i)) ; alpar@1: j = Pe [i] ; alpar@1: cnt = 0 ; alpar@1: AMD_DEBUG3 ((" j: "ID"\n", j)) ; alpar@1: if (j == EMPTY) alpar@1: { alpar@1: AMD_DEBUG3 ((" i is a dense variable\n")) ; alpar@1: } alpar@1: else alpar@1: { alpar@1: ASSERT (j >= 0 && j < n) ; alpar@1: while (Nv [j] == 0) alpar@1: { alpar@1: AMD_DEBUG3 ((" j : "ID"\n", j)) ; alpar@1: j = Pe [j] ; alpar@1: AMD_DEBUG3 ((" j:: "ID"\n", j)) ; alpar@1: cnt++ ; alpar@1: if (cnt > n) break ; alpar@1: } alpar@1: e = j ; alpar@1: AMD_DEBUG3 ((" got to e: "ID"\n", e)) ; alpar@1: } alpar@1: } alpar@1: } alpar@1: #endif alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* compress the paths of the variables */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: for (i = 0 ; i < n ; i++) alpar@1: { alpar@1: if (Nv [i] == 0) alpar@1: { alpar@1: alpar@1: /* ------------------------------------------------------------- alpar@1: * i is an un-ordered row. Traverse the tree from i until alpar@1: * reaching an element, e. The element, e, was the principal alpar@1: * supervariable of i and all nodes in the path from i to when e alpar@1: * was selected as pivot. alpar@1: * ------------------------------------------------------------- */ alpar@1: alpar@1: AMD_DEBUG1 (("Path compression, i unordered: "ID"\n", i)) ; alpar@1: j = Pe [i] ; alpar@1: ASSERT (j >= EMPTY && j < n) ; alpar@1: AMD_DEBUG3 ((" j: "ID"\n", j)) ; alpar@1: if (j == EMPTY) alpar@1: { alpar@1: /* Skip a dense variable. It has no parent. */ alpar@1: AMD_DEBUG3 ((" i is a dense variable\n")) ; alpar@1: continue ; alpar@1: } alpar@1: alpar@1: /* while (j is a variable) */ alpar@1: while (Nv [j] == 0) alpar@1: { alpar@1: AMD_DEBUG3 ((" j : "ID"\n", j)) ; alpar@1: j = Pe [j] ; alpar@1: AMD_DEBUG3 ((" j:: "ID"\n", j)) ; alpar@1: ASSERT (j >= 0 && j < n) ; alpar@1: } alpar@1: /* got to an element e */ alpar@1: e = j ; alpar@1: AMD_DEBUG3 (("got to e: "ID"\n", e)) ; alpar@1: alpar@1: /* ------------------------------------------------------------- alpar@1: * traverse the path again from i to e, and compress the path alpar@1: * (all nodes point to e). Path compression allows this code to alpar@1: * compute in O(n) time. alpar@1: * ------------------------------------------------------------- */ alpar@1: alpar@1: j = i ; alpar@1: /* while (j is a variable) */ alpar@1: while (Nv [j] == 0) alpar@1: { alpar@1: jnext = Pe [j] ; alpar@1: AMD_DEBUG3 (("j "ID" jnext "ID"\n", j, jnext)) ; alpar@1: Pe [j] = e ; alpar@1: j = jnext ; alpar@1: ASSERT (j >= 0 && j < n) ; alpar@1: } alpar@1: } alpar@1: } alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* postorder the assembly tree */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: AMD_postorder (n, Pe, Nv, Elen, alpar@1: W, /* output order */ alpar@1: Head, Next, Last) ; /* workspace */ alpar@1: alpar@1: /* ========================================================================= */ alpar@1: /* compute output permutation and inverse permutation */ alpar@1: /* ========================================================================= */ alpar@1: alpar@1: /* W [e] = k means that element e is the kth element in the new alpar@1: * order. e is in the range 0 to n-1, and k is in the range 0 to alpar@1: * the number of elements. Use Head for inverse order. */ alpar@1: alpar@1: for (k = 0 ; k < n ; k++) alpar@1: { alpar@1: Head [k] = EMPTY ; alpar@1: Next [k] = EMPTY ; alpar@1: } alpar@1: for (e = 0 ; e < n ; e++) alpar@1: { alpar@1: k = W [e] ; alpar@1: ASSERT ((k == EMPTY) == (Nv [e] == 0)) ; alpar@1: if (k != EMPTY) alpar@1: { alpar@1: ASSERT (k >= 0 && k < n) ; alpar@1: Head [k] = e ; alpar@1: } alpar@1: } alpar@1: alpar@1: /* construct output inverse permutation in Next, alpar@1: * and permutation in Last */ alpar@1: nel = 0 ; alpar@1: for (k = 0 ; k < n ; k++) alpar@1: { alpar@1: e = Head [k] ; alpar@1: if (e == EMPTY) break ; alpar@1: ASSERT (e >= 0 && e < n && Nv [e] > 0) ; alpar@1: Next [e] = nel ; alpar@1: nel += Nv [e] ; alpar@1: } alpar@1: ASSERT (nel == n - ndense) ; alpar@1: alpar@1: /* order non-principal variables (dense, & those merged into supervar's) */ alpar@1: for (i = 0 ; i < n ; i++) alpar@1: { alpar@1: if (Nv [i] == 0) alpar@1: { alpar@1: e = Pe [i] ; alpar@1: ASSERT (e >= EMPTY && e < n) ; alpar@1: if (e != EMPTY) alpar@1: { alpar@1: /* This is an unordered variable that was merged alpar@1: * into element e via supernode detection or mass alpar@1: * elimination of i when e became the pivot element. alpar@1: * Place i in order just before e. */ alpar@1: ASSERT (Next [i] == EMPTY && Nv [e] > 0) ; alpar@1: Next [i] = Next [e] ; alpar@1: Next [e]++ ; alpar@1: } alpar@1: else alpar@1: { alpar@1: /* This is a dense unordered variable, with no parent. alpar@1: * Place it last in the output order. */ alpar@1: Next [i] = nel++ ; alpar@1: } alpar@1: } alpar@1: } alpar@1: ASSERT (nel == n) ; alpar@1: alpar@1: AMD_DEBUG2 (("\n\nPerm:\n")) ; alpar@1: for (i = 0 ; i < n ; i++) alpar@1: { alpar@1: k = Next [i] ; alpar@1: ASSERT (k >= 0 && k < n) ; alpar@1: Last [k] = i ; alpar@1: AMD_DEBUG2 ((" perm ["ID"] = "ID"\n", k, i)) ; alpar@1: } alpar@1: }