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