Amesos Package Browser (Single Doxygen Collection) Development
amesos_amd_2.c
Go to the documentation of this file.
00001 /* ========================================================================= */
00002 /* === AMD_2 =============================================================== */
00003 /* ========================================================================= */
00004 
00005 /* ------------------------------------------------------------------------- */
00006 /* AMD, Copyright (c) Timothy A. Davis,              */
00007 /* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
00008 /* email: davis at cise.ufl.edu    CISE Department, Univ. of Florida.        */
00009 /* web: http://www.cise.ufl.edu/research/sparse/amd                          */
00010 /* ------------------------------------------------------------------------- */
00011 
00012 /* AMD_2:  performs the AMD ordering on a symmetric sparse matrix A, followed
00013  * by a postordering (via depth-first search) of the assembly tree using the
00014  * AMD_postorder routine.
00015  */
00016 
00017 #include "amesos_amd_internal.h"
00018 
00019 /* ========================================================================= */
00020 /* === clear_flag ========================================================== */
00021 /* ========================================================================= */
00022 
00023 static Int amesos_clear_flag (Int wflg, Int wbig, Int W [ ], Int n)
00024 {
00025     Int x ;
00026     if (wflg < 2 || wflg >= wbig)
00027     {
00028   for (x = 0 ; x < n ; x++)
00029   {
00030       if (W [x] != 0) W [x] = 1 ;
00031   }
00032   wflg = 2 ;
00033     }
00034     /*  at this point, W [0..n-1] < wflg holds */
00035     return (wflg) ;
00036 }
00037 
00038 
00039 /* ========================================================================= */
00040 /* === AMD_2 =============================================================== */
00041 /* ========================================================================= */
00042 
00043 GLOBAL void AMD_2
00044 (
00045     Int n,    /* A is n-by-n, where n > 0 */
00046     Int Pe [ ],   /* Pe [0..n-1]: index in Iw of row i on input */
00047     Int Iw [ ],   /* workspace of size iwlen. Iw [0..pfree-1]
00048        * holds the matrix on input */
00049     Int Len [ ],  /* Len [0..n-1]: length for row/column i on input */
00050     Int iwlen,    /* length of Iw. iwlen >= pfree + n */
00051     Int pfree,    /* Iw [pfree ... iwlen-1] is empty on input */
00052 
00053     /* 7 size-n workspaces, not defined on input: */
00054     Int Nv [ ],   /* the size of each supernode on output */
00055     Int Next [ ], /* the output inverse permutation */
00056     Int Last [ ], /* the output permutation */
00057     Int Head [ ],
00058     Int Elen [ ], /* the size columns of L for each supernode */
00059     Int Degree [ ],
00060     Int W [ ],
00061 
00062     /* control parameters and output statistics */
00063     double Control [ ], /* array of size AMD_CONTROL */
00064     double Info [ ] /* array of size AMD_INFO */
00065 )
00066 {
00067 
00068 /*
00069  * Given a representation of the nonzero pattern of a symmetric matrix, A,
00070  * (excluding the diagonal) perform an approximate minimum (UMFPACK/MA38-style)
00071  * degree ordering to compute a pivot order such that the introduction of
00072  * nonzeros (fill-in) in the Cholesky factors A = LL' is kept low.  At each
00073  * step, the pivot selected is the one with the minimum UMFAPACK/MA38-style
00074  * upper-bound on the external degree.  This routine can optionally perform
00075  * aggresive absorption (as done by MC47B in the Harwell Subroutine
00076  * Library).
00077  *
00078  * The approximate degree algorithm implemented here is the symmetric analog of
00079  * the degree update algorithm in MA38 and UMFPACK (the Unsymmetric-pattern
00080  * MultiFrontal PACKage, both by Davis and Duff).  The routine is based on the
00081  * MA27 minimum degree ordering algorithm by Iain Duff and John Reid.
00082  *
00083  * This routine is a translation of the original AMDBAR and MC47B routines,
00084  * in Fortran, with the following modifications:
00085  *
00086  * (1) dense rows/columns are removed prior to ordering the matrix, and placed
00087  *  last in the output order.  The presence of a dense row/column can
00088  *  increase the ordering time by up to O(n^2), unless they are removed
00089  *  prior to ordering.
00090  *
00091  * (2) the minimum degree ordering is followed by a postordering (depth-first
00092  *  search) of the assembly tree.  Note that mass elimination (discussed
00093  *  below) combined with the approximate degree update can lead to the mass
00094  *  elimination of nodes with lower exact degree than the current pivot
00095  *  element.  No additional fill-in is caused in the representation of the
00096  *  Schur complement.  The mass-eliminated nodes merge with the current
00097  *  pivot element.  They are ordered prior to the current pivot element.
00098  *  Because they can have lower exact degree than the current element, the
00099  *  merger of two or more of these nodes in the current pivot element can
00100  *  lead to a single element that is not a "fundamental supernode".  The
00101  *  diagonal block can have zeros in it.  Thus, the assembly tree used here
00102  *  is not guaranteed to be the precise supernodal elemination tree (with
00103  *  "funadmental" supernodes), and the postordering performed by this
00104  *  routine is not guaranteed to be a precise postordering of the
00105  *  elimination tree.
00106  *
00107  * (3) input parameters are added, to control aggressive absorption and the
00108  *  detection of "dense" rows/columns of A.
00109  *
00110  * (4) additional statistical information is returned, such as the number of
00111  *  nonzeros in L, and the flop counts for subsequent LDL' and LU
00112  *  factorizations.  These are slight upper bounds, because of the mass
00113  *  elimination issue discussed above.
00114  *
00115  * (5) additional routines are added to interface this routine to MATLAB
00116  *  to provide a simple C-callable user-interface, to check inputs for
00117  *  errors, compute the symmetry of the pattern of A and the number of
00118  *  nonzeros in each row/column of A+A', to compute the pattern of A+A',
00119  *  to perform the assembly tree postordering, and to provide debugging
00120  *  ouput.  Many of these functions are also provided by the Fortran
00121  *  Harwell Subroutine Library routine MC47A.
00122  *
00123  * (6) both int and UF_long versions are provided.  In the descriptions below
00124  *  and integer is and int or UF_long depending on which version is
00125  *  being used.
00126 
00127  **********************************************************************
00128  ***** CAUTION:  ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT.  ******
00129  **********************************************************************
00130  ** If you want error checking, a more versatile input format, and a **
00131  ** simpler user interface, use amd_order or amd_l_order instead.    **
00132  ** This routine is not meant to be user-callable.                   **
00133  **********************************************************************
00134 
00135  * ----------------------------------------------------------------------------
00136  * References:
00137  * ----------------------------------------------------------------------------
00138  *
00139  *  [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern multifrontal
00140  *  method for sparse LU factorization", SIAM J. Matrix Analysis and
00141  *  Applications, vol. 18, no. 1, pp. 140-158.  Discusses UMFPACK / MA38,
00142  *  which first introduced the approximate minimum degree used by this
00143  *  routine.
00144  *
00145  *  [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, "An approximate
00146  *  minimum degree ordering algorithm," SIAM J. Matrix Analysis and
00147  *  Applications, vol. 17, no. 4, pp. 886-905, 1996.  Discusses AMDBAR and
00148  *  MC47B, which are the Fortran versions of this routine.
00149  *
00150  *  [3] Alan George and Joseph Liu, "The evolution of the minimum degree
00151  *  ordering algorithm," SIAM Review, vol. 31, no. 1, pp. 1-19, 1989.
00152  *  We list below the features mentioned in that paper that this code
00153  *  includes:
00154  *
00155  *  mass elimination:
00156  *      Yes.  MA27 relied on supervariable detection for mass elimination.
00157  *
00158  *  indistinguishable nodes:
00159  *      Yes (we call these "supervariables").  This was also in the MA27
00160  *      code - although we modified the method of detecting them (the
00161  *      previous hash was the true degree, which we no longer keep track
00162  *      of).  A supervariable is a set of rows with identical nonzero
00163  *      pattern.  All variables in a supervariable are eliminated together.
00164  *      Each supervariable has as its numerical name that of one of its
00165  *      variables (its principal variable).
00166  *
00167  *  quotient graph representation:
00168  *      Yes.  We use the term "element" for the cliques formed during
00169  *      elimination.  This was also in the MA27 code.  The algorithm can
00170  *      operate in place, but it will work more efficiently if given some
00171  *      "elbow room."
00172  *
00173  *  element absorption:
00174  *      Yes.  This was also in the MA27 code.
00175  *
00176  *  external degree:
00177  *      Yes.  The MA27 code was based on the true degree.
00178  *
00179  *  incomplete degree update and multiple elimination:
00180  *      No.  This was not in MA27, either.  Our method of degree update
00181  *      within MC47B is element-based, not variable-based.  It is thus
00182  *      not well-suited for use with incomplete degree update or multiple
00183  *      elimination.
00184  *
00185  * Authors, and Copyright (C) 2004 by:
00186  * Timothy A. Davis, Patrick Amestoy, Iain S. Duff, John K. Reid.
00187  *
00188  * Acknowledgements: This work (and the UMFPACK package) was supported by the
00189  * National Science Foundation (ASC-9111263, DMS-9223088, and CCR-0203270).
00190  * The UMFPACK/MA38 approximate degree update algorithm, the unsymmetric analog
00191  * which forms the basis of AMD, was developed while Tim Davis was supported by
00192  * CERFACS (Toulouse, France) in a post-doctoral position.  This C version, and
00193  * the etree postorder, were written while Tim Davis was on sabbatical at
00194  * Stanford University and Lawrence Berkeley National Laboratory.
00195 
00196  * ----------------------------------------------------------------------------
00197  * INPUT ARGUMENTS (unaltered):
00198  * ----------------------------------------------------------------------------
00199 
00200  * n:  The matrix order.  Restriction:  n >= 1.
00201  *
00202  * iwlen:  The size of the Iw array.  On input, the matrix is stored in
00203  *  Iw [0..pfree-1].  However, Iw [0..iwlen-1] should be slightly larger
00204  *  than what is required to hold the matrix, at least iwlen >= pfree + n.
00205  *  Otherwise, excessive compressions will take place.  The recommended
00206  *  value of iwlen is 1.2 * pfree + n, which is the value used in the
00207  *  user-callable interface to this routine (amd_order.c).  The algorithm
00208  *  will not run at all if iwlen < pfree.  Restriction: iwlen >= pfree + n.
00209  *  Note that this is slightly more restrictive than the actual minimum
00210  *  (iwlen >= pfree), but AMD_2 will be very slow with no elbow room.
00211  *  Thus, this routine enforces a bare minimum elbow room of size n.
00212  *
00213  * pfree: On input the tail end of the array, Iw [pfree..iwlen-1], is empty,
00214  *  and the matrix is stored in Iw [0..pfree-1].  During execution,
00215  *  additional data is placed in Iw, and pfree is modified so that
00216  *  Iw [pfree..iwlen-1] is always the unused part of Iw.
00217  *
00218  * Control:  A double array of size AMD_CONTROL containing input parameters
00219  *  that affect how the ordering is computed.  If NULL, then default
00220  *  settings are used.
00221  *
00222  *  Control [AMD_DENSE] is used to determine whether or not a given input
00223  *  row is "dense".  A row is "dense" if the number of entries in the row
00224  *  exceeds Control [AMD_DENSE] times sqrt (n), except that rows with 16 or
00225  *  fewer entries are never considered "dense".  To turn off the detection
00226  *  of dense rows, set Control [AMD_DENSE] to a negative number, or to a
00227  *  number larger than sqrt (n).  The default value of Control [AMD_DENSE]
00228  *  is AMD_DEFAULT_DENSE, which is defined in amd.h as 10.
00229  *
00230  *  Control [AMD_AGGRESSIVE] is used to determine whether or not aggressive
00231  *  absorption is to be performed.  If nonzero, then aggressive absorption
00232  *  is performed (this is the default).
00233 
00234  * ----------------------------------------------------------------------------
00235  * INPUT/OUPUT ARGUMENTS:
00236  * ----------------------------------------------------------------------------
00237  *
00238  * Pe:  An integer array of size n.  On input, Pe [i] is the index in Iw of
00239  *  the start of row i.  Pe [i] is ignored if row i has no off-diagonal
00240  *  entries.  Thus Pe [i] must be in the range 0 to pfree-1 for non-empty
00241  *  rows.
00242  *
00243  *  During execution, it is used for both supervariables and elements:
00244  *
00245  *  Principal supervariable i:  index into Iw of the description of
00246  *      supervariable i.  A supervariable represents one or more rows of
00247  *      the matrix with identical nonzero pattern.  In this case,
00248  *      Pe [i] >= 0.
00249  *
00250  *  Non-principal supervariable i:  if i has been absorbed into another
00251  *      supervariable j, then Pe [i] = FLIP (j), where FLIP (j) is defined
00252  *      as (-(j)-2).  Row j has the same pattern as row i.  Note that j
00253  *      might later be absorbed into another supervariable j2, in which
00254  *      case Pe [i] is still FLIP (j), and Pe [j] = FLIP (j2) which is
00255  *      < EMPTY, where EMPTY is defined as (-1) in amd_internal.h.
00256  *
00257  *  Unabsorbed element e:  the index into Iw of the description of element
00258  *      e, if e has not yet been absorbed by a subsequent element.  Element
00259  *      e is created when the supervariable of the same name is selected as
00260  *      the pivot.  In this case, Pe [i] >= 0.
00261  *
00262  *  Absorbed element e:  if element e is absorbed into element e2, then
00263  *      Pe [e] = FLIP (e2).  This occurs when the pattern of e (which we
00264  *      refer to as Le) is found to be a subset of the pattern of e2 (that
00265  *      is, Le2).  In this case, Pe [i] < EMPTY.  If element e is "null"
00266  *      (it has no nonzeros outside its pivot block), then Pe [e] = EMPTY,
00267  *      and e is the root of an assembly subtree (or the whole tree if
00268  *      there is just one such root).
00269  *
00270  *  Dense variable i:  if i is "dense", then Pe [i] = EMPTY.
00271  *
00272  *  On output, Pe holds the assembly tree/forest, which implicitly
00273  *  represents a pivot order with identical fill-in as the actual order
00274  *  (via a depth-first search of the tree), as follows.  If Nv [i] > 0,
00275  *  then i represents a node in the assembly tree, and the parent of i is
00276  *  Pe [i], or EMPTY if i is a root.  If Nv [i] = 0, then (i, Pe [i])
00277  *  represents an edge in a subtree, the root of which is a node in the
00278  *  assembly tree.  Note that i refers to a row/column in the original
00279  *  matrix, not the permuted matrix.
00280  *
00281  * Info:  A double array of size AMD_INFO.  If present, (that is, not NULL),
00282  *  then statistics about the ordering are returned in the Info array.
00283  *  See amd.h for a description.
00284 
00285  * ----------------------------------------------------------------------------
00286  * INPUT/MODIFIED (undefined on output):
00287  * ----------------------------------------------------------------------------
00288  *
00289  * Len:  An integer array of size n.  On input, Len [i] holds the number of
00290  *  entries in row i of the matrix, excluding the diagonal.  The contents
00291  *  of Len are undefined on output.
00292  *
00293  * Iw:  An integer array of size iwlen.  On input, Iw [0..pfree-1] holds the
00294  *  description of each row i in the matrix.  The matrix must be symmetric,
00295  *  and both upper and lower triangular parts must be present.  The
00296  *  diagonal must not be present.  Row i is held as follows:
00297  *
00298  *      Len [i]:  the length of the row i data structure in the Iw array.
00299  *      Iw [Pe [i] ... Pe [i] + Len [i] - 1]:
00300  *    the list of column indices for nonzeros in row i (simple
00301  *    supervariables), excluding the diagonal.  All supervariables
00302  *    start with one row/column each (supervariable i is just row i).
00303  *    If Len [i] is zero on input, then Pe [i] is ignored on input.
00304  *
00305  *      Note that the rows need not be in any particular order, and there
00306  *      may be empty space between the rows.
00307  *
00308  *  During execution, the supervariable i experiences fill-in.  This is
00309  *  represented by placing in i a list of the elements that cause fill-in
00310  *  in supervariable i:
00311  *
00312  *      Len [i]:  the length of supervariable i in the Iw array.
00313  *      Iw [Pe [i] ... Pe [i] + Elen [i] - 1]:
00314  *    the list of elements that contain i.  This list is kept short
00315  *    by removing absorbed elements.
00316  *      Iw [Pe [i] + Elen [i] ... Pe [i] + Len [i] - 1]:
00317  *    the list of supervariables in i.  This list is kept short by
00318  *    removing nonprincipal variables, and any entry j that is also
00319  *    contained in at least one of the elements (j in Le) in the list
00320  *    for i (e in row i).
00321  *
00322  *  When supervariable i is selected as pivot, we create an element e of
00323  *  the same name (e=i):
00324  *
00325  *      Len [e]:  the length of element e in the Iw array.
00326  *      Iw [Pe [e] ... Pe [e] + Len [e] - 1]:
00327  *    the list of supervariables in element e.
00328  *
00329  *  An element represents the fill-in that occurs when supervariable i is
00330  *  selected as pivot (which represents the selection of row i and all
00331  *  non-principal variables whose principal variable is i).  We use the
00332  *  term Le to denote the set of all supervariables in element e.  Absorbed
00333  *  supervariables and elements are pruned from these lists when
00334  *  computationally convenient.
00335  *
00336  *  CAUTION:  THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION.
00337  *  The contents of Iw are undefined on output.
00338 
00339  * ----------------------------------------------------------------------------
00340  * OUTPUT (need not be set on input):
00341  * ----------------------------------------------------------------------------
00342  *
00343  * Nv:  An integer array of size n.  During execution, ABS (Nv [i]) is equal to
00344  *  the number of rows that are represented by the principal supervariable
00345  *  i.  If i is a nonprincipal or dense variable, then Nv [i] = 0.
00346  *  Initially, Nv [i] = 1 for all i.  Nv [i] < 0 signifies that i is a
00347  *  principal variable in the pattern Lme of the current pivot element me.
00348  *  After element me is constructed, Nv [i] is set back to a positive
00349  *  value.
00350  *
00351  *  On output, Nv [i] holds the number of pivots represented by super
00352  *  row/column i of the original matrix, or Nv [i] = 0 for non-principal
00353  *  rows/columns.  Note that i refers to a row/column in the original
00354  *  matrix, not the permuted matrix.
00355  *
00356  * Elen:  An integer array of size n.  See the description of Iw above.  At the
00357  *  start of execution, Elen [i] is set to zero for all rows i.  During
00358  *  execution, Elen [i] is the number of elements in the list for
00359  *  supervariable i.  When e becomes an element, Elen [e] = FLIP (esize) is
00360  *  set, where esize is the size of the element (the number of pivots, plus
00361  *  the number of nonpivotal entries).  Thus Elen [e] < EMPTY.
00362  *  Elen (i) = EMPTY set when variable i becomes nonprincipal.
00363  *
00364  *  For variables, Elen (i) >= EMPTY holds until just before the
00365  *  postordering and permutation vectors are computed.  For elements,
00366  *  Elen [e] < EMPTY holds.
00367  *
00368  *  On output, Elen [i] is the degree of the row/column in the Cholesky
00369  *  factorization of the permuted matrix, corresponding to the original row
00370  *  i, if i is a super row/column.  It is equal to EMPTY if i is
00371  *  non-principal.  Note that i refers to a row/column in the original
00372  *  matrix, not the permuted matrix.
00373  *
00374  *  Note that the contents of Elen on output differ from the Fortran
00375  *  version (Elen holds the inverse permutation in the Fortran version,
00376  *  which is instead returned in the Next array in this C version,
00377  *  described below).
00378  *
00379  * Last: In a degree list, Last [i] is the supervariable preceding i, or EMPTY
00380  *  if i is the head of the list.  In a hash bucket, Last [i] is the hash
00381  *  key for i.
00382  *
00383  *  Last [Head [hash]] is also used as the head of a hash bucket if
00384  *  Head [hash] contains a degree list (see the description of Head,
00385  *  below).
00386  *
00387  *  On output, Last [0..n-1] holds the permutation.  That is, if
00388  *  i = Last [k], then row i is the kth pivot row (where k ranges from 0 to
00389  *  n-1).  Row Last [k] of A is the kth row in the permuted matrix, PAP'.
00390  *
00391  * Next: Next [i] is the supervariable following i in a link list, or EMPTY if
00392  *  i is the last in the list.  Used for two kinds of lists:  degree lists
00393  *  and hash buckets (a supervariable can be in only one kind of list at a
00394  *  time).
00395  *
00396  *  On output Next [0..n-1] holds the inverse permutation.  That is, if
00397  *  k = Next [i], then row i is the kth pivot row. Row i of A appears as
00398  *  the (Next[i])-th row in the permuted matrix, PAP'.
00399  *
00400  *  Note that the contents of Next on output differ from the Fortran
00401  *  version (Next is undefined on output in the Fortran version).
00402 
00403  * ----------------------------------------------------------------------------
00404  * LOCAL WORKSPACE (not input or output - used only during execution):
00405  * ----------------------------------------------------------------------------
00406  *
00407  * Degree:  An integer array of size n.  If i is a supervariable, then
00408  *  Degree [i] holds the current approximation of the external degree of
00409  *  row i (an upper bound).  The external degree is the number of nonzeros
00410  *  in row i, minus ABS (Nv [i]), the diagonal part.  The bound is equal to
00411  *  the exact external degree if Elen [i] is less than or equal to two.
00412  *
00413  *  We also use the term "external degree" for elements e to refer to
00414  *  |Le \ Lme|.  If e is an element, then Degree [e] is |Le|, which is the
00415  *  degree of the off-diagonal part of the element e (not including the
00416  *  diagonal part).
00417  *
00418  * Head:   An integer array of size n.  Head is used for degree lists.
00419  *  Head [deg] is the first supervariable in a degree list.  All
00420  *  supervariables i in a degree list Head [deg] have the same approximate
00421  *  degree, namely, deg = Degree [i].  If the list Head [deg] is empty then
00422  *  Head [deg] = EMPTY.
00423  *
00424  *  During supervariable detection Head [hash] also serves as a pointer to
00425  *  a hash bucket.  If Head [hash] >= 0, there is a degree list of degree
00426  *  hash.  The hash bucket head pointer is Last [Head [hash]].  If
00427  *  Head [hash] = EMPTY, then the degree list and hash bucket are both
00428  *  empty.  If Head [hash] < EMPTY, then the degree list is empty, and
00429  *  FLIP (Head [hash]) is the head of the hash bucket.  After supervariable
00430  *  detection is complete, all hash buckets are empty, and the
00431  *  (Last [Head [hash]] = EMPTY) condition is restored for the non-empty
00432  *  degree lists.
00433  *
00434  * W:  An integer array of size n.  The flag array W determines the status of
00435  *  elements and variables, and the external degree of elements.
00436  *
00437  *  for elements:
00438  *      if W [e] = 0, then the element e is absorbed.
00439  *      if W [e] >= wflg, then W [e] - wflg is the size of the set
00440  *    |Le \ Lme|, in terms of nonzeros (the sum of ABS (Nv [i]) for
00441  *    each principal variable i that is both in the pattern of
00442  *    element e and NOT in the pattern of the current pivot element,
00443  *    me).
00444  *      if wflg > W [e] > 0, then e is not absorbed and has not yet been
00445  *    seen in the scan of the element lists in the computation of
00446  *    |Le\Lme| in Scan 1 below.
00447  *
00448  *  for variables:
00449  *      during supervariable detection, if W [j] != wflg then j is
00450  *      not in the pattern of variable i.
00451  *
00452  *  The W array is initialized by setting W [i] = 1 for all i, and by
00453  *  setting wflg = 2.  It is reinitialized if wflg becomes too large (to
00454  *  ensure that wflg+n does not cause integer overflow).
00455 
00456  * ----------------------------------------------------------------------------
00457  * LOCAL INTEGERS:
00458  * ----------------------------------------------------------------------------
00459  */
00460 
00461     Int deg, degme, dext, lemax, e, elenme, eln, i, ilast, inext, j,
00462   jlast, jnext, k, knt1, knt2, knt3, lenj, ln, me, mindeg, nel, nleft,
00463   nvi, nvj, nvpiv, slenme, wbig, we, wflg, wnvi, ok, ndense, ncmpa,
00464   dense, aggressive ;
00465 
00466     unsigned Int hash ;     /* unsigned, so that hash % n is well defined.*/
00467 
00468 /*
00469  * deg:   the degree of a variable or element
00470  * degme: size, |Lme|, of the current element, me (= Degree [me])
00471  * dext:  external degree, |Le \ Lme|, of some element e
00472  * lemax: largest |Le| seen so far (called dmax in Fortran version)
00473  * e:   an element
00474  * elenme:  the length, Elen [me], of element list of pivotal variable
00475  * eln:   the length, Elen [...], of an element list
00476  * hash:  the computed value of the hash function
00477  * i:   a supervariable
00478  * ilast: the entry in a link list preceding i
00479  * inext: the entry in a link list following i
00480  * j:   a supervariable
00481  * jlast: the entry in a link list preceding j
00482  * jnext: the entry in a link list, or path, following j
00483  * k:   the pivot order of an element or variable
00484  * knt1:  loop counter used during element construction
00485  * knt2:  loop counter used during element construction
00486  * knt3:  loop counter used during compression
00487  * lenj:  Len [j]
00488  * ln:    length of a supervariable list
00489  * me:    current supervariable being eliminated, and the current
00490  *        element created by eliminating that supervariable
00491  * mindeg:  current minimum degree
00492  * nel:   number of pivots selected so far
00493  * nleft: n - nel, the number of nonpivotal rows/columns remaining
00494  * nvi:   the number of variables in a supervariable i (= Nv [i])
00495  * nvj:   the number of variables in a supervariable j (= Nv [j])
00496  * nvpiv: number of pivots in current element
00497  * slenme:  number of variables in variable list of pivotal variable
00498  * wbig:  = INT_MAX - n for the int version, UF_long_max - n for the
00499  *        UF_long version.  wflg is not allowed to be >= wbig.
00500  * we:    W [e]
00501  * wflg:  used for flagging the W array.  See description of Iw.
00502  * wnvi:  wflg - Nv [i]
00503  * x:   either a supervariable or an element
00504  *
00505  * ok:    true if supervariable j can be absorbed into i
00506  * ndense:  number of "dense" rows/columns
00507  * dense: rows/columns with initial degree > dense are considered "dense"
00508  * aggressive:  true if aggressive absorption is being performed
00509  * ncmpa: number of garbage collections
00510 
00511  * ----------------------------------------------------------------------------
00512  * LOCAL DOUBLES, used for statistical output only (except for alpha):
00513  * ----------------------------------------------------------------------------
00514  */
00515 
00516     double f, r, ndiv, s, nms_lu, nms_ldl, dmax, alpha, lnz, lnzme ;
00517 
00518 /*
00519  * f:   nvpiv
00520  * r:   degme + nvpiv
00521  * ndiv:  number of divisions for LU or LDL' factorizations
00522  * s:   number of multiply-subtract pairs for LU factorization, for the
00523  *        current element me
00524  * nms_lu number of multiply-subtract pairs for LU factorization
00525  * nms_ldl  number of multiply-subtract pairs for LDL' factorization
00526  * dmax:  the largest number of entries in any column of L, including the
00527  *        diagonal
00528  * alpha: "dense" degree ratio
00529  * lnz:   the number of nonzeros in L (excluding the diagonal)
00530  * lnzme: the number of nonzeros in L (excl. the diagonal) for the
00531  *        current element me
00532 
00533  * ----------------------------------------------------------------------------
00534  * LOCAL "POINTERS" (indices into the Iw array)
00535  * ----------------------------------------------------------------------------
00536 */
00537 
00538     Int p, p1, p2, p3, p4, pdst, pend, pj, pme, pme1, pme2, pn, psrc ;
00539 
00540 /*
00541  * Any parameter (Pe [...] or pfree) or local variable starting with "p" (for
00542  * Pointer) is an index into Iw, and all indices into Iw use variables starting
00543  * with "p."  The only exception to this rule is the iwlen input argument.
00544  *
00545  * p:           pointer into lots of things
00546  * p1:          Pe [i] for some variable i (start of element list)
00547  * p2:          Pe [i] + Elen [i] -  1 for some variable i
00548  * p3:          index of first supervariable in clean list
00549  * p4:    
00550  * pdst:        destination pointer, for compression
00551  * pend:        end of memory to compress
00552  * pj:          pointer into an element or variable
00553  * pme:         pointer into the current element (pme1...pme2)
00554  * pme1:        the current element, me, is stored in Iw [pme1...pme2]
00555  * pme2:        the end of the current element
00556  * pn:          pointer into a "clean" variable, also used to compress
00557  * psrc:        source pointer, for compression
00558 */
00559 
00560 /* ========================================================================= */
00561 /*  INITIALIZATIONS */
00562 /* ========================================================================= */
00563 
00564     /* Note that this restriction on iwlen is slightly more restrictive than
00565      * what is actually required in AMD_2.  AMD_2 can operate with no elbow
00566      * room at all, but it will be slow.  For better performance, at least
00567      * size-n elbow room is enforced. */
00568     ASSERT (iwlen >= pfree + n) ;
00569     ASSERT (n > 0) ;
00570 
00571     /* initialize output statistics */
00572     lnz = 0 ;
00573     ndiv = 0 ;
00574     nms_lu = 0 ;
00575     nms_ldl = 0 ;
00576     dmax = 1 ;
00577     me = EMPTY ;
00578 
00579     mindeg = 0 ;
00580     ncmpa = 0 ;
00581     nel = 0 ;
00582     lemax = 0 ;
00583 
00584     /* get control parameters */
00585     if (Control != (double *) NULL)
00586     {
00587   alpha = Control [AMD_DENSE] ;
00588   aggressive = (Control [AMD_AGGRESSIVE] != 0) ;
00589     }
00590     else
00591     {
00592   alpha = AMD_DEFAULT_DENSE ;
00593   aggressive = AMD_DEFAULT_AGGRESSIVE ;
00594     }
00595     /* Note: if alpha is NaN, this is undefined: */
00596     if (alpha < 0)
00597     {
00598   /* only remove completely dense rows/columns */
00599   dense = n-2 ;
00600     }
00601     else
00602     {
00603   dense = (int) ( alpha * sqrt ((double) n) ) ;
00604     }
00605     dense = MAX (16, dense) ;
00606     dense = MIN (n,  dense) ;
00607     AMD_DEBUG1 (("\n\nAMD (debug), alpha %g, aggr. "ID"\n",
00608   alpha, aggressive)) ;
00609 
00610     for (i = 0 ; i < n ; i++)
00611     {
00612   Last [i] = EMPTY ;
00613   Head [i] = EMPTY ;
00614   Next [i] = EMPTY ;
00615   /* if separate Hhead array is used for hash buckets: *
00616   Hhead [i] = EMPTY ;
00617   */
00618   Nv [i] = 1 ;
00619   W [i] = 1 ;
00620   Elen [i] = 0 ;
00621   Degree [i] = Len [i] ;
00622     }
00623 
00624 #ifndef NDEBUG
00625     AMD_DEBUG1 (("\n======Nel "ID" initial\n", nel)) ;
00626     AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next, Last,
00627     Head, Elen, Degree, W, -1) ;
00628 #endif
00629 
00630     /* initialize wflg */
00631     wbig = Int_MAX - n ;
00632     wflg = amesos_clear_flag (0, wbig, W, n) ;
00633 
00634     /* --------------------------------------------------------------------- */
00635     /* initialize degree lists and eliminate dense and empty rows */
00636     /* --------------------------------------------------------------------- */
00637 
00638     ndense = 0 ;
00639 
00640     for (i = 0 ; i < n ; i++)
00641     {
00642   deg = Degree [i] ;
00643   ASSERT (deg >= 0 && deg < n) ;
00644   if (deg == 0)
00645   {
00646 
00647       /* -------------------------------------------------------------
00648        * we have a variable that can be eliminated at once because
00649        * there is no off-diagonal non-zero in its row.  Note that
00650        * Nv [i] = 1 for an empty variable i.  It is treated just
00651        * the same as an eliminated element i.
00652        * ------------------------------------------------------------- */
00653 
00654       Elen [i] = FLIP (1) ;
00655       nel++ ;
00656       Pe [i] = EMPTY ;
00657       W [i] = 0 ;
00658 
00659   }
00660   else if (deg > dense)
00661   {
00662 
00663       /* -------------------------------------------------------------
00664        * Dense variables are not treated as elements, but as unordered,
00665        * non-principal variables that have no parent.  They do not take
00666        * part in the postorder, since Nv [i] = 0.  Note that the Fortran
00667        * version does not have this option.
00668        * ------------------------------------------------------------- */
00669 
00670       AMD_DEBUG1 (("Dense node "ID" degree "ID"\n", i, deg)) ;
00671       ndense++ ;
00672       Nv [i] = 0 ;    /* do not postorder this node */
00673       Elen [i] = EMPTY ;
00674       nel++ ;
00675       Pe [i] = EMPTY ;
00676 
00677   }
00678   else
00679   {
00680 
00681       /* -------------------------------------------------------------
00682        * place i in the degree list corresponding to its degree
00683        * ------------------------------------------------------------- */
00684 
00685       inext = Head [deg] ;
00686       ASSERT (inext >= EMPTY && inext < n) ;
00687       if (inext != EMPTY) Last [inext] = i ;
00688       Next [i] = inext ;
00689       Head [deg] = i ;
00690 
00691   }
00692     }
00693 
00694 /* ========================================================================= */
00695 /* WHILE (selecting pivots) DO */
00696 /* ========================================================================= */
00697 
00698     while (nel < n)
00699     {
00700 
00701 #ifndef NDEBUG
00702   AMD_DEBUG1 (("\n======Nel "ID"\n", nel)) ;
00703   if (AMD_debug >= 2)
00704   {
00705       AMD_dump (n, Pe, Iw, Len, iwlen, pfree, Nv, Next,
00706         Last, Head, Elen, Degree, W, nel) ;
00707   }
00708 #endif
00709 
00710 /* ========================================================================= */
00711 /* GET PIVOT OF MINIMUM DEGREE */
00712 /* ========================================================================= */
00713 
00714   /* ----------------------------------------------------------------- */
00715   /* find next supervariable for elimination */
00716   /* ----------------------------------------------------------------- */
00717 
00718   ASSERT (mindeg >= 0 && mindeg < n) ;
00719   for (deg = mindeg ; deg < n ; deg++)
00720   {
00721       me = Head [deg] ;
00722       if (me != EMPTY) break ;
00723   }
00724   mindeg = deg ;
00725   ASSERT (me >= 0 && me < n) ;
00726   AMD_DEBUG1 (("=================me: "ID"\n", me)) ;
00727 
00728   /* ----------------------------------------------------------------- */
00729   /* remove chosen variable from link list */
00730   /* ----------------------------------------------------------------- */
00731 
00732   inext = Next [me] ;
00733   ASSERT (inext >= EMPTY && inext < n) ;
00734   if (inext != EMPTY) Last [inext] = EMPTY ;
00735   Head [deg] = inext ;
00736 
00737   /* ----------------------------------------------------------------- */
00738   /* me represents the elimination of pivots nel to nel+Nv[me]-1. */
00739   /* place me itself as the first in this set. */
00740   /* ----------------------------------------------------------------- */
00741 
00742   elenme = Elen [me] ;
00743   nvpiv = Nv [me] ;
00744   ASSERT (nvpiv > 0) ;
00745   nel += nvpiv ;
00746 
00747 /* ========================================================================= */
00748 /* CONSTRUCT NEW ELEMENT */
00749 /* ========================================================================= */
00750 
00751   /* -----------------------------------------------------------------
00752    * At this point, me is the pivotal supervariable.  It will be
00753    * converted into the current element.  Scan list of the pivotal
00754    * supervariable, me, setting tree pointers and constructing new list
00755    * of supervariables for the new element, me.  p is a pointer to the
00756    * current position in the old list.
00757    * ----------------------------------------------------------------- */
00758 
00759   /* flag the variable "me" as being in Lme by negating Nv [me] */
00760   Nv [me] = -nvpiv ;
00761   degme = 0 ;
00762   ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ;
00763 
00764   if (elenme == 0)
00765   {
00766 
00767       /* ------------------------------------------------------------- */
00768       /* construct the new element in place */
00769       /* ------------------------------------------------------------- */
00770 
00771       pme1 = Pe [me] ;
00772       pme2 = pme1 - 1 ;
00773 
00774       for (p = pme1 ; p <= pme1 + Len [me] - 1 ; p++)
00775       {
00776     i = Iw [p] ;
00777     ASSERT (i >= 0 && i < n && Nv [i] >= 0) ;
00778     nvi = Nv [i] ;
00779     if (nvi > 0)
00780     {
00781 
00782         /* ----------------------------------------------------- */
00783         /* i is a principal variable not yet placed in Lme. */
00784         /* store i in new list */
00785         /* ----------------------------------------------------- */
00786 
00787         /* flag i as being in Lme by negating Nv [i] */
00788         degme += nvi ;
00789         Nv [i] = -nvi ;
00790         Iw [++pme2] = i ;
00791 
00792         /* ----------------------------------------------------- */
00793         /* remove variable i from degree list. */
00794         /* ----------------------------------------------------- */
00795 
00796         ilast = Last [i] ;
00797         inext = Next [i] ;
00798         ASSERT (ilast >= EMPTY && ilast < n) ;
00799         ASSERT (inext >= EMPTY && inext < n) ;
00800         if (inext != EMPTY) Last [inext] = ilast ;
00801         if (ilast != EMPTY)
00802         {
00803       Next [ilast] = inext ;
00804         }
00805         else
00806         {
00807       /* i is at the head of the degree list */
00808       ASSERT (Degree [i] >= 0 && Degree [i] < n) ;
00809       Head [Degree [i]] = inext ;
00810         }
00811     }
00812       }
00813   }
00814   else
00815   {
00816 
00817       /* ------------------------------------------------------------- */
00818       /* construct the new element in empty space, Iw [pfree ...] */
00819       /* ------------------------------------------------------------- */
00820 
00821       p = Pe [me] ;
00822       pme1 = pfree ;
00823       slenme = Len [me] - elenme ;
00824 
00825       for (knt1 = 1 ; knt1 <= elenme + 1 ; knt1++)
00826       {
00827 
00828     if (knt1 > elenme)
00829     {
00830         /* search the supervariables in me. */
00831         e = me ;
00832         pj = p ;
00833         ln = slenme ;
00834         AMD_DEBUG2 (("Search sv: "ID" "ID" "ID"\n", me,pj,ln)) ;
00835     }
00836     else
00837     {
00838         /* search the elements in me. */
00839         e = Iw [p++] ;
00840         ASSERT (e >= 0 && e < n) ;
00841         pj = Pe [e] ;
00842         ln = Len [e] ;
00843         AMD_DEBUG2 (("Search element e "ID" in me "ID"\n", e,me)) ;
00844         ASSERT (Elen [e] < EMPTY && W [e] > 0 && pj >= 0) ;
00845     }
00846     ASSERT (ln >= 0 && (ln == 0 || (pj >= 0 && pj < iwlen))) ;
00847 
00848     /* ---------------------------------------------------------
00849      * search for different supervariables and add them to the
00850      * new list, compressing when necessary. this loop is
00851      * executed once for each element in the list and once for
00852      * all the supervariables in the list.
00853      * --------------------------------------------------------- */
00854 
00855     for (knt2 = 1 ; knt2 <= ln ; knt2++)
00856     {
00857         i = Iw [pj++] ;
00858         ASSERT (i >= 0 && i < n && (i == me || Elen [i] >= EMPTY));
00859         nvi = Nv [i] ;
00860         AMD_DEBUG2 ((": "ID" "ID" "ID" "ID"\n",
00861         i, Elen [i], Nv [i], wflg)) ;
00862 
00863         if (nvi > 0)
00864         {
00865 
00866       /* ------------------------------------------------- */
00867       /* compress Iw, if necessary */
00868       /* ------------------------------------------------- */
00869 
00870       if (pfree >= iwlen)
00871       {
00872 
00873           AMD_DEBUG1 (("GARBAGE COLLECTION\n")) ;
00874 
00875           /* prepare for compressing Iw by adjusting pointers
00876            * and lengths so that the lists being searched in
00877            * the inner and outer loops contain only the
00878            * remaining entries. */
00879 
00880           Pe [me] = p ;
00881           Len [me] -= knt1 ;
00882           /* check if nothing left of supervariable me */
00883           if (Len [me] == 0) Pe [me] = EMPTY ;
00884           Pe [e] = pj ;
00885           Len [e] = ln - knt2 ;
00886           /* nothing left of element e */
00887           if (Len [e] == 0) Pe [e] = EMPTY ;
00888 
00889           ncmpa++ ; /* one more garbage collection */
00890 
00891           /* store first entry of each object in Pe */
00892           /* FLIP the first entry in each object */
00893           for (j = 0 ; j < n ; j++)
00894           {
00895         pn = Pe [j] ;
00896         if (pn >= 0)
00897         {
00898             ASSERT (pn >= 0 && pn < iwlen) ;
00899             Pe [j] = Iw [pn] ;
00900             Iw [pn] = FLIP (j) ;
00901         }
00902           }
00903 
00904           /* psrc/pdst point to source/destination */
00905           psrc = 0 ;
00906           pdst = 0 ;
00907           pend = pme1 - 1 ;
00908 
00909           while (psrc <= pend)
00910           {
00911         /* search for next FLIP'd entry */
00912         j = FLIP (Iw [psrc++]) ;
00913         if (j >= 0)
00914         {
00915             AMD_DEBUG2 (("Got object j: "ID"\n", j)) ;
00916             Iw [pdst] = Pe [j] ;
00917             Pe [j] = pdst++ ;
00918             lenj = Len [j] ;
00919             /* copy from source to destination */
00920             for (knt3 = 0 ; knt3 <= lenj - 2 ; knt3++)
00921             {
00922           Iw [pdst++] = Iw [psrc++] ;
00923             }
00924         }
00925           }
00926 
00927           /* move the new partially-constructed element */
00928           p1 = pdst ;
00929           for (psrc = pme1 ; psrc <= pfree-1 ; psrc++)
00930           {
00931         Iw [pdst++] = Iw [psrc] ;
00932           }
00933           pme1 = p1 ;
00934           pfree = pdst ;
00935           pj = Pe [e] ;
00936           p = Pe [me] ;
00937 
00938       }
00939 
00940       /* ------------------------------------------------- */
00941       /* i is a principal variable not yet placed in Lme */
00942       /* store i in new list */
00943       /* ------------------------------------------------- */
00944 
00945       /* flag i as being in Lme by negating Nv [i] */
00946       degme += nvi ;
00947       Nv [i] = -nvi ;
00948       Iw [pfree++] = i ;
00949       AMD_DEBUG2 (("     s: "ID"     nv "ID"\n", i, Nv [i]));
00950 
00951       /* ------------------------------------------------- */
00952       /* remove variable i from degree link list */
00953       /* ------------------------------------------------- */
00954 
00955       ilast = Last [i] ;
00956       inext = Next [i] ;
00957       ASSERT (ilast >= EMPTY && ilast < n) ;
00958       ASSERT (inext >= EMPTY && inext < n) ;
00959       if (inext != EMPTY) Last [inext] = ilast ;
00960       if (ilast != EMPTY)
00961       {
00962           Next [ilast] = inext ;
00963       }
00964       else
00965       {
00966           /* i is at the head of the degree list */
00967           ASSERT (Degree [i] >= 0 && Degree [i] < n) ;
00968           Head [Degree [i]] = inext ;
00969       }
00970         }
00971     }
00972 
00973     if (e != me)
00974     {
00975         /* set tree pointer and flag to indicate element e is
00976          * absorbed into new element me (the parent of e is me) */
00977         AMD_DEBUG1 ((" Element "ID" => "ID"\n", e, me)) ;
00978         Pe [e] = FLIP (me) ;
00979         W [e] = 0 ;
00980     }
00981       }
00982 
00983       pme2 = pfree - 1 ;
00984   }
00985 
00986   /* ----------------------------------------------------------------- */
00987   /* me has now been converted into an element in Iw [pme1..pme2] */
00988   /* ----------------------------------------------------------------- */
00989 
00990   /* degme holds the external degree of new element */
00991   Degree [me] = degme ;
00992   Pe [me] = pme1 ;
00993   Len [me] = pme2 - pme1 + 1 ;
00994   ASSERT (Pe [me] >= 0 && Pe [me] < iwlen) ;
00995 
00996   Elen [me] = FLIP (nvpiv + degme) ;
00997   /* FLIP (Elen (me)) is now the degree of pivot (including
00998    * diagonal part). */
00999 
01000 #ifndef NDEBUG
01001   AMD_DEBUG2 (("New element structure: length= "ID"\n", pme2-pme1+1)) ;
01002   for (pme = pme1 ; pme <= pme2 ; pme++) AMD_DEBUG3 ((" "ID"", Iw[pme]));
01003   AMD_DEBUG3 (("\n")) ;
01004 #endif
01005 
01006   /* ----------------------------------------------------------------- */
01007   /* make sure that wflg is not too large. */
01008   /* ----------------------------------------------------------------- */
01009 
01010   /* With the current value of wflg, wflg+n must not cause integer
01011    * overflow */
01012 
01013   wflg = amesos_clear_flag (wflg, wbig, W, n) ;
01014 
01015 /* ========================================================================= */
01016 /* COMPUTE (W [e] - wflg) = |Le\Lme| FOR ALL ELEMENTS */
01017 /* ========================================================================= */
01018 
01019   /* -----------------------------------------------------------------
01020    * Scan 1:  compute the external degrees of previous elements with
01021    * respect to the current element.  That is:
01022    *       (W [e] - wflg) = |Le \ Lme|
01023    * for each element e that appears in any supervariable in Lme.  The
01024    * notation Le refers to the pattern (list of supervariables) of a
01025    * previous element e, where e is not yet absorbed, stored in
01026    * Iw [Pe [e] + 1 ... Pe [e] + Len [e]].  The notation Lme
01027    * refers to the pattern of the current element (stored in
01028    * Iw [pme1..pme2]).   If aggressive absorption is enabled, and
01029    * (W [e] - wflg) becomes zero, then the element e will be absorbed
01030    * in Scan 2.
01031    * ----------------------------------------------------------------- */
01032 
01033   AMD_DEBUG2 (("me: ")) ;
01034   for (pme = pme1 ; pme <= pme2 ; pme++)
01035   {
01036       i = Iw [pme] ;
01037       ASSERT (i >= 0 && i < n) ;
01038       eln = Elen [i] ;
01039       AMD_DEBUG3 ((""ID" Elen "ID": \n", i, eln)) ;
01040       if (eln > 0)
01041       {
01042     /* note that Nv [i] has been negated to denote i in Lme: */
01043     nvi = -Nv [i] ;
01044     ASSERT (nvi > 0 && Pe [i] >= 0 && Pe [i] < iwlen) ;
01045     wnvi = wflg - nvi ;
01046     for (p = Pe [i] ; p <= Pe [i] + eln - 1 ; p++)
01047     {
01048         e = Iw [p] ;
01049         ASSERT (e >= 0 && e < n) ;
01050         we = W [e] ;
01051         AMD_DEBUG4 (("    e "ID" we "ID" ", e, we)) ;
01052         if (we >= wflg)
01053         {
01054       /* unabsorbed element e has been seen in this loop */
01055       AMD_DEBUG4 (("    unabsorbed, first time seen")) ;
01056       we -= nvi ;
01057         }
01058         else if (we != 0)
01059         {
01060       /* e is an unabsorbed element */
01061       /* this is the first we have seen e in all of Scan 1 */
01062       AMD_DEBUG4 (("    unabsorbed")) ;
01063       we = Degree [e] + wnvi ;
01064         }
01065         AMD_DEBUG4 (("\n")) ;
01066         W [e] = we ;
01067     }
01068       }
01069   }
01070   AMD_DEBUG2 (("\n")) ;
01071 
01072 /* ========================================================================= */
01073 /* DEGREE UPDATE AND ELEMENT ABSORPTION */
01074 /* ========================================================================= */
01075 
01076   /* -----------------------------------------------------------------
01077    * Scan 2:  for each i in Lme, sum up the degree of Lme (which is
01078    * degme), plus the sum of the external degrees of each Le for the
01079    * elements e appearing within i, plus the supervariables in i.
01080    * Place i in hash list.
01081    * ----------------------------------------------------------------- */
01082 
01083   for (pme = pme1 ; pme <= pme2 ; pme++)
01084   {
01085       i = Iw [pme] ;
01086       ASSERT (i >= 0 && i < n && Nv [i] < 0 && Elen [i] >= 0) ;
01087       AMD_DEBUG2 (("Updating: i "ID" "ID" "ID"\n", i, Elen[i], Len [i]));
01088       p1 = Pe [i] ;
01089       p2 = p1 + Elen [i] - 1 ;
01090       pn = p1 ;
01091       hash = 0 ;
01092       deg = 0 ;
01093       ASSERT (p1 >= 0 && p1 < iwlen && p2 >= -1 && p2 < iwlen) ;
01094 
01095       /* ------------------------------------------------------------- */
01096       /* scan the element list associated with supervariable i */
01097       /* ------------------------------------------------------------- */
01098 
01099       /* UMFPACK/MA38-style approximate degree: */
01100       if (aggressive)
01101       {
01102     for (p = p1 ; p <= p2 ; p++)
01103     {
01104         e = Iw [p] ;
01105         ASSERT (e >= 0 && e < n) ;
01106         we = W [e] ;
01107         if (we != 0)
01108         {
01109       /* e is an unabsorbed element */
01110       /* dext = | Le \ Lme | */
01111       dext = we - wflg ;
01112       if (dext > 0)
01113       {
01114           deg += dext ;
01115           Iw [pn++] = e ;
01116           hash += e ;
01117           AMD_DEBUG4 ((" e: "ID" hash = "ID"\n",e,hash)) ;
01118       }
01119       else
01120       {
01121           /* external degree of e is zero, absorb e into me*/
01122           AMD_DEBUG1 ((" Element "ID" =>"ID" (aggressive)\n",
01123         e, me)) ;
01124           ASSERT (dext == 0) ;
01125           Pe [e] = FLIP (me) ;
01126           W [e] = 0 ;
01127       }
01128         }
01129     }
01130       }
01131       else
01132       {
01133     for (p = p1 ; p <= p2 ; p++)
01134     {
01135         e = Iw [p] ;
01136         ASSERT (e >= 0 && e < n) ;
01137         we = W [e] ;
01138         if (we != 0)
01139         {
01140       /* e is an unabsorbed element */
01141       dext = we - wflg ;
01142       ASSERT (dext >= 0) ;
01143       deg += dext ;
01144       Iw [pn++] = e ;
01145       hash += e ;
01146       AMD_DEBUG4 (("  e: "ID" hash = "ID"\n",e,hash)) ;
01147         }
01148     }
01149       }
01150 
01151       /* count the number of elements in i (including me): */
01152       Elen [i] = pn - p1 + 1 ;
01153 
01154       /* ------------------------------------------------------------- */
01155       /* scan the supervariables in the list associated with i */
01156       /* ------------------------------------------------------------- */
01157 
01158       /* The bulk of the AMD run time is typically spent in this loop,
01159        * particularly if the matrix has many dense rows that are not
01160        * removed prior to ordering. */
01161       p3 = pn ;
01162       p4 = p1 + Len [i] ;
01163       for (p = p2 + 1 ; p < p4 ; p++)
01164       {
01165     j = Iw [p] ;
01166     ASSERT (j >= 0 && j < n) ;
01167     nvj = Nv [j] ;
01168     if (nvj > 0)
01169     {
01170         /* j is unabsorbed, and not in Lme. */
01171         /* add to degree and add to new list */
01172         deg += nvj ;
01173         Iw [pn++] = j ;
01174         hash += j ;
01175         AMD_DEBUG4 (("  s: "ID" hash "ID" Nv[j]= "ID"\n",
01176         j, hash, nvj)) ;
01177     }
01178       }
01179 
01180       /* ------------------------------------------------------------- */
01181       /* update the degree and check for mass elimination */
01182       /* ------------------------------------------------------------- */
01183 
01184       /* with aggressive absorption, deg==0 is identical to the
01185        * Elen [i] == 1 && p3 == pn test, below. */
01186       ASSERT (IMPLIES (aggressive, (deg==0) == (Elen[i]==1 && p3==pn))) ;
01187 
01188       if (Elen [i] == 1 && p3 == pn)
01189       {
01190 
01191     /* --------------------------------------------------------- */
01192     /* mass elimination */
01193     /* --------------------------------------------------------- */
01194 
01195     /* There is nothing left of this node except for an edge to
01196      * the current pivot element.  Elen [i] is 1, and there are
01197      * no variables adjacent to node i.  Absorb i into the
01198      * current pivot element, me.  Note that if there are two or
01199      * more mass eliminations, fillin due to mass elimination is
01200      * possible within the nvpiv-by-nvpiv pivot block.  It is this
01201      * step that causes AMD's analysis to be an upper bound.
01202      *
01203      * The reason is that the selected pivot has a lower
01204      * approximate degree than the true degree of the two mass
01205      * eliminated nodes.  There is no edge between the two mass
01206      * eliminated nodes.  They are merged with the current pivot
01207      * anyway.
01208      *
01209      * No fillin occurs in the Schur complement, in any case,
01210      * and this effect does not decrease the quality of the
01211      * ordering itself, just the quality of the nonzero and
01212      * flop count analysis.  It also means that the post-ordering
01213      * is not an exact elimination tree post-ordering. */
01214 
01215     AMD_DEBUG1 (("  MASS i "ID" => parent e "ID"\n", i, me)) ;
01216     Pe [i] = FLIP (me) ;
01217     nvi = -Nv [i] ;
01218     degme -= nvi ;
01219     nvpiv += nvi ;
01220     nel += nvi ;
01221     Nv [i] = 0 ;
01222     Elen [i] = EMPTY ;
01223 
01224       }
01225       else
01226       {
01227 
01228     /* --------------------------------------------------------- */
01229     /* update the upper-bound degree of i */
01230     /* --------------------------------------------------------- */
01231 
01232     /* the following degree does not yet include the size
01233      * of the current element, which is added later: */
01234 
01235     Degree [i] = MIN (Degree [i], deg) ;
01236 
01237     /* --------------------------------------------------------- */
01238     /* add me to the list for i */
01239     /* --------------------------------------------------------- */
01240 
01241     /* move first supervariable to end of list */
01242     Iw [pn] = Iw [p3] ;
01243     /* move first element to end of element part of list */
01244     Iw [p3] = Iw [p1] ;
01245     /* add new element, me, to front of list. */
01246     Iw [p1] = me ;
01247     /* store the new length of the list in Len [i] */
01248     Len [i] = pn - p1 + 1 ;
01249 
01250     /* --------------------------------------------------------- */
01251     /* place in hash bucket.  Save hash key of i in Last [i]. */
01252     /* --------------------------------------------------------- */
01253 
01254     /* NOTE: this can fail if hash is negative, because the ANSI C
01255      * standard does not define a % b when a and/or b are negative.
01256      * That's why hash is defined as an unsigned Int, to avoid this
01257      * problem. */
01258     hash = hash % n ;
01259     ASSERT (((Int) hash) >= 0 && ((Int) hash) < n) ;
01260 
01261     /* if the Hhead array is not used: */
01262     j = Head [hash] ;
01263     if (j <= EMPTY)
01264     {
01265         /* degree list is empty, hash head is FLIP (j) */
01266         Next [i] = FLIP (j) ;
01267         Head [hash] = FLIP (i) ;
01268     }
01269     else
01270     {
01271         /* degree list is not empty, use Last [Head [hash]] as
01272          * hash head. */
01273         Next [i] = Last [j] ;
01274         Last [j] = i ;
01275     }
01276 
01277     /* if a separate Hhead array is used: *
01278     Next [i] = Hhead [hash] ;
01279     Hhead [hash] = i ;
01280     */
01281 
01282     Last [i] = hash ;
01283       }
01284   }
01285 
01286   Degree [me] = degme ;
01287 
01288   /* ----------------------------------------------------------------- */
01289   /* Clear the counter array, W [...], by incrementing wflg. */
01290   /* ----------------------------------------------------------------- */
01291 
01292   /* make sure that wflg+n does not cause integer overflow */
01293   lemax =  MAX (lemax, degme) ;
01294   wflg += lemax ;
01295   wflg = amesos_clear_flag (wflg, wbig, W, n) ;
01296   /*  at this point, W [0..n-1] < wflg holds */
01297 
01298 /* ========================================================================= */
01299 /* SUPERVARIABLE DETECTION */
01300 /* ========================================================================= */
01301 
01302   AMD_DEBUG1 (("Detecting supervariables:\n")) ;
01303   for (pme = pme1 ; pme <= pme2 ; pme++)
01304   {
01305       i = Iw [pme] ;
01306       ASSERT (i >= 0 && i < n) ;
01307       AMD_DEBUG2 (("Consider i "ID" nv "ID"\n", i, Nv [i])) ;
01308       if (Nv [i] < 0)
01309       {
01310     /* i is a principal variable in Lme */
01311 
01312     /* ---------------------------------------------------------
01313      * examine all hash buckets with 2 or more variables.  We do
01314      * this by examing all unique hash keys for supervariables in
01315      * the pattern Lme of the current element, me
01316      * --------------------------------------------------------- */
01317 
01318     /* let i = head of hash bucket, and empty the hash bucket */
01319     ASSERT (Last [i] >= 0 && Last [i] < n) ;
01320     hash = Last [i] ;
01321 
01322     /* if Hhead array is not used: */
01323     j = Head [hash] ;
01324     if (j == EMPTY)
01325     {
01326         /* hash bucket and degree list are both empty */
01327         i = EMPTY ;
01328     }
01329     else if (j < EMPTY)
01330     {
01331         /* degree list is empty */
01332         i = FLIP (j) ;
01333         Head [hash] = EMPTY ;
01334     }
01335     else
01336     {
01337         /* degree list is not empty, restore Last [j] of head j */
01338         i = Last [j] ;
01339         Last [j] = EMPTY ;
01340     }
01341 
01342     /* if separate Hhead array is used: *
01343     i = Hhead [hash] ;
01344     Hhead [hash] = EMPTY ;
01345     */
01346 
01347     ASSERT (i >= EMPTY && i < n) ;
01348     AMD_DEBUG2 (("----i "ID" hash "ID"\n", i, hash)) ;
01349 
01350     while (i != EMPTY && Next [i] != EMPTY)
01351     {
01352 
01353         /* -----------------------------------------------------
01354          * this bucket has one or more variables following i.
01355          * scan all of them to see if i can absorb any entries
01356          * that follow i in hash bucket.  Scatter i into w.
01357          * ----------------------------------------------------- */
01358 
01359         ln = Len [i] ;
01360         eln = Elen [i] ;
01361         ASSERT (ln >= 0 && eln >= 0) ;
01362         ASSERT (Pe [i] >= 0 && Pe [i] < iwlen) ;
01363         /* do not flag the first element in the list (me) */
01364         for (p = Pe [i] + 1 ; p <= Pe [i] + ln - 1 ; p++)
01365         {
01366       ASSERT (Iw [p] >= 0 && Iw [p] < n) ;
01367       W [Iw [p]] = wflg ;
01368         }
01369 
01370         /* ----------------------------------------------------- */
01371         /* scan every other entry j following i in bucket */
01372         /* ----------------------------------------------------- */
01373 
01374         jlast = i ;
01375         j = Next [i] ;
01376         ASSERT (j >= EMPTY && j < n) ;
01377 
01378         while (j != EMPTY)
01379         {
01380       /* ------------------------------------------------- */
01381       /* check if j and i have identical nonzero pattern */
01382       /* ------------------------------------------------- */
01383 
01384       AMD_DEBUG3 (("compare i "ID" and j "ID"\n", i,j)) ;
01385 
01386       /* check if i and j have the same Len and Elen */
01387       ASSERT (Len [j] >= 0 && Elen [j] >= 0) ;
01388       ASSERT (Pe [j] >= 0 && Pe [j] < iwlen) ;
01389       ok = (Len [j] == ln) && (Elen [j] == eln) ;
01390       /* skip the first element in the list (me) */
01391       for (p = Pe [j] + 1 ; ok && p <= Pe [j] + ln - 1 ; p++)
01392       {
01393           ASSERT (Iw [p] >= 0 && Iw [p] < n) ;
01394           if (W [Iw [p]] != wflg) ok = 0 ;
01395       }
01396       if (ok)
01397       {
01398           /* --------------------------------------------- */
01399           /* found it!  j can be absorbed into i */
01400           /* --------------------------------------------- */
01401 
01402           AMD_DEBUG1 (("found it! j "ID" => i "ID"\n", j,i));
01403           Pe [j] = FLIP (i) ;
01404           /* both Nv [i] and Nv [j] are negated since they */
01405           /* are in Lme, and the absolute values of each */
01406           /* are the number of variables in i and j: */
01407           Nv [i] += Nv [j] ;
01408           Nv [j] = 0 ;
01409           Elen [j] = EMPTY ;
01410           /* delete j from hash bucket */
01411           ASSERT (j != Next [j]) ;
01412           j = Next [j] ;
01413           Next [jlast] = j ;
01414 
01415       }
01416       else
01417       {
01418           /* j cannot be absorbed into i */
01419           jlast = j ;
01420           ASSERT (j != Next [j]) ;
01421           j = Next [j] ;
01422       }
01423       ASSERT (j >= EMPTY && j < n) ;
01424         }
01425 
01426         /* -----------------------------------------------------
01427          * no more variables can be absorbed into i
01428          * go to next i in bucket and clear flag array
01429          * ----------------------------------------------------- */
01430 
01431         wflg++ ;
01432         i = Next [i] ;
01433         ASSERT (i >= EMPTY && i < n) ;
01434 
01435     }
01436       }
01437   }
01438   AMD_DEBUG2 (("detect done\n")) ;
01439 
01440 /* ========================================================================= */
01441 /* RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVARIABLES FROM ELEMENT */
01442 /* ========================================================================= */
01443 
01444   p = pme1 ;
01445   nleft = n - nel ;
01446   for (pme = pme1 ; pme <= pme2 ; pme++)
01447   {
01448       i = Iw [pme] ;
01449       ASSERT (i >= 0 && i < n) ;
01450       nvi = -Nv [i] ;
01451       AMD_DEBUG3 (("Restore i "ID" "ID"\n", i, nvi)) ;
01452       if (nvi > 0)
01453       {
01454     /* i is a principal variable in Lme */
01455     /* restore Nv [i] to signify that i is principal */
01456     Nv [i] = nvi ;
01457 
01458     /* --------------------------------------------------------- */
01459     /* compute the external degree (add size of current element) */
01460     /* --------------------------------------------------------- */
01461 
01462     deg = Degree [i] + degme - nvi ;
01463     deg = MIN (deg, nleft - nvi) ;
01464     ASSERT (IMPLIES (aggressive, deg > 0) && deg >= 0 && deg < n) ;
01465 
01466     /* --------------------------------------------------------- */
01467     /* place the supervariable at the head of the degree list */
01468     /* --------------------------------------------------------- */
01469 
01470     inext = Head [deg] ;
01471     ASSERT (inext >= EMPTY && inext < n) ;
01472     if (inext != EMPTY) Last [inext] = i ;
01473     Next [i] = inext ;
01474     Last [i] = EMPTY ;
01475     Head [deg] = i ;
01476 
01477     /* --------------------------------------------------------- */
01478     /* save the new degree, and find the minimum degree */
01479     /* --------------------------------------------------------- */
01480 
01481     mindeg = MIN (mindeg, deg) ;
01482     Degree [i] = deg ;
01483 
01484     /* --------------------------------------------------------- */
01485     /* place the supervariable in the element pattern */
01486     /* --------------------------------------------------------- */
01487 
01488     Iw [p++] = i ;
01489 
01490       }
01491   }
01492   AMD_DEBUG2 (("restore done\n")) ;
01493 
01494 /* ========================================================================= */
01495 /* FINALIZE THE NEW ELEMENT */
01496 /* ========================================================================= */
01497 
01498   AMD_DEBUG2 (("ME = "ID" DONE\n", me)) ;
01499   Nv [me] = nvpiv ;
01500   /* save the length of the list for the new element me */
01501   Len [me] = p - pme1 ;
01502   if (Len [me] == 0)
01503   {
01504       /* there is nothing left of the current pivot element */
01505       /* it is a root of the assembly tree */
01506       Pe [me] = EMPTY ;
01507       W [me] = 0 ;
01508   }
01509   if (elenme != 0)
01510   {
01511       /* element was not constructed in place: deallocate part of */
01512       /* it since newly nonprincipal variables may have been removed */
01513       pfree = p ;
01514   }
01515 
01516   /* The new element has nvpiv pivots and the size of the contribution
01517    * block for a multifrontal method is degme-by-degme, not including
01518    * the "dense" rows/columns.  If the "dense" rows/columns are included,
01519    * the frontal matrix is no larger than
01520    * (degme+ndense)-by-(degme+ndense).
01521    */
01522 
01523   if (Info != (double *) NULL)
01524   {
01525       f = nvpiv ;
01526       r = degme + ndense ;
01527       dmax = MAX (dmax, f + r) ;
01528 
01529       /* number of nonzeros in L (excluding the diagonal) */
01530       lnzme = f*r + (f-1)*f/2 ;
01531       lnz += lnzme ;
01532 
01533       /* number of divide operations for LDL' and for LU */
01534       ndiv += lnzme ;
01535 
01536       /* number of multiply-subtract pairs for LU */
01537       s = f*r*r + r*(f-1)*f + (f-1)*f*(2*f-1)/6 ;
01538       nms_lu += s ;
01539 
01540       /* number of multiply-subtract pairs for LDL' */
01541       nms_ldl += (s + lnzme)/2 ;
01542   }
01543 
01544 #ifndef NDEBUG
01545   AMD_DEBUG2 (("finalize done nel "ID" n "ID"\n   ::::\n", nel, n)) ;
01546   for (pme = Pe [me] ; pme <= Pe [me] + Len [me] - 1 ; pme++)
01547   {
01548         AMD_DEBUG3 ((" "ID"", Iw [pme])) ;
01549   }
01550   AMD_DEBUG3 (("\n")) ;
01551 #endif
01552 
01553     }
01554 
01555 /* ========================================================================= */
01556 /* DONE SELECTING PIVOTS */
01557 /* ========================================================================= */
01558 
01559     if (Info != (double *) NULL)
01560     {
01561 
01562   /* count the work to factorize the ndense-by-ndense submatrix */
01563   f = ndense ;
01564   dmax = MAX (dmax, (double) ndense) ;
01565 
01566   /* number of nonzeros in L (excluding the diagonal) */
01567   lnzme = (f-1)*f/2 ;
01568   lnz += lnzme ;
01569 
01570   /* number of divide operations for LDL' and for LU */
01571   ndiv += lnzme ;
01572 
01573   /* number of multiply-subtract pairs for LU */
01574   s = (f-1)*f*(2*f-1)/6 ;
01575   nms_lu += s ;
01576 
01577   /* number of multiply-subtract pairs for LDL' */
01578   nms_ldl += (s + lnzme)/2 ;
01579 
01580   /* number of nz's in L (excl. diagonal) */
01581   Info [AMD_LNZ] = lnz ;
01582 
01583   /* number of divide ops for LU and LDL' */
01584   Info [AMD_NDIV] = ndiv ;
01585 
01586   /* number of multiply-subtract pairs for LDL' */
01587   Info [AMD_NMULTSUBS_LDL] = nms_ldl ;
01588 
01589   /* number of multiply-subtract pairs for LU */
01590   Info [AMD_NMULTSUBS_LU] = nms_lu ;
01591 
01592   /* number of "dense" rows/columns */
01593   Info [AMD_NDENSE] = ndense ;
01594 
01595   /* largest front is dmax-by-dmax */
01596   Info [AMD_DMAX] = dmax ;
01597 
01598   /* number of garbage collections in AMD */
01599   Info [AMD_NCMPA] = ncmpa ;
01600 
01601   /* successful ordering */
01602   Info [AMD_STATUS] = AMD_OK ;
01603     }
01604 
01605 /* ========================================================================= */
01606 /* POST-ORDERING */
01607 /* ========================================================================= */
01608 
01609 /* -------------------------------------------------------------------------
01610  * Variables at this point:
01611  *
01612  * Pe: holds the elimination tree.  The parent of j is FLIP (Pe [j]),
01613  *  or EMPTY if j is a root.  The tree holds both elements and
01614  *  non-principal (unordered) variables absorbed into them.
01615  *  Dense variables are non-principal and unordered.
01616  *
01617  * Elen: holds the size of each element, including the diagonal part.
01618  *  FLIP (Elen [e]) > 0 if e is an element.  For unordered
01619  *  variables i, Elen [i] is EMPTY.
01620  *
01621  * Nv: Nv [e] > 0 is the number of pivots represented by the element e.
01622  *  For unordered variables i, Nv [i] is zero.
01623  *
01624  * Contents no longer needed:
01625  *  W, Iw, Len, Degree, Head, Next, Last.
01626  *
01627  * The matrix itself has been destroyed.
01628  *
01629  * n: the size of the matrix.
01630  * No other scalars needed (pfree, iwlen, etc.)
01631  * ------------------------------------------------------------------------- */
01632 
01633     /* restore Pe */
01634     for (i = 0 ; i < n ; i++)
01635     {
01636   Pe [i] = FLIP (Pe [i]) ;
01637     }
01638 
01639     /* restore Elen, for output information, and for postordering */
01640     for (i = 0 ; i < n ; i++)
01641     {
01642   Elen [i] = FLIP (Elen [i]) ;
01643     }
01644 
01645 /* Now the parent of j is Pe [j], or EMPTY if j is a root.  Elen [e] > 0
01646  * is the size of element e.  Elen [i] is EMPTY for unordered variable i. */
01647 
01648 #ifndef NDEBUG
01649     AMD_DEBUG2 (("\nTree:\n")) ;
01650     for (i = 0 ; i < n ; i++)
01651     {
01652   AMD_DEBUG2 ((" "ID" parent: "ID"   ", i, Pe [i])) ;
01653   ASSERT (Pe [i] >= EMPTY && Pe [i] < n) ;
01654   if (Nv [i] > 0)
01655   {
01656       /* this is an element */
01657       e = i ;
01658       AMD_DEBUG2 ((" element, size is "ID"\n", Elen [i])) ;
01659       ASSERT (Elen [e] > 0) ;
01660   }
01661   AMD_DEBUG2 (("\n")) ;
01662     }
01663     AMD_DEBUG2 (("\nelements:\n")) ;
01664     for (e = 0 ; e < n ; e++)
01665     {
01666   if (Nv [e] > 0)
01667   {
01668       AMD_DEBUG3 (("Element e= "ID" size "ID" nv "ID" \n", e,
01669     Elen [e], Nv [e])) ;
01670   }
01671     }
01672     AMD_DEBUG2 (("\nvariables:\n")) ;
01673     for (i = 0 ; i < n ; i++)
01674     {
01675   Int cnt ;
01676   if (Nv [i] == 0)
01677   {
01678       AMD_DEBUG3 (("i unordered: "ID"\n", i)) ;
01679       j = Pe [i] ;
01680       cnt = 0 ;
01681       AMD_DEBUG3 (("  j: "ID"\n", j)) ;
01682       if (j == EMPTY)
01683       {
01684     AMD_DEBUG3 (("  i is a dense variable\n")) ;
01685       }
01686       else
01687       {
01688     ASSERT (j >= 0 && j < n) ;
01689     while (Nv [j] == 0)
01690     {
01691         AMD_DEBUG3 (("  j : "ID"\n", j)) ;
01692         j = Pe [j] ;
01693         AMD_DEBUG3 (("  j:: "ID"\n", j)) ;
01694         cnt++ ;
01695         if (cnt > n) break ;
01696     }
01697     e = j ;
01698     AMD_DEBUG3 (("  got to e: "ID"\n", e)) ;
01699       }
01700   }
01701     }
01702 #endif
01703 
01704 /* ========================================================================= */
01705 /* compress the paths of the variables */
01706 /* ========================================================================= */
01707 
01708     for (i = 0 ; i < n ; i++)
01709     {
01710   if (Nv [i] == 0)
01711   {
01712 
01713       /* -------------------------------------------------------------
01714        * i is an un-ordered row.  Traverse the tree from i until
01715        * reaching an element, e.  The element, e, was the principal
01716        * supervariable of i and all nodes in the path from i to when e
01717        * was selected as pivot.
01718        * ------------------------------------------------------------- */
01719 
01720       AMD_DEBUG1 (("Path compression, i unordered: "ID"\n", i)) ;
01721       j = Pe [i] ;
01722       ASSERT (j >= EMPTY && j < n) ;
01723       AMD_DEBUG3 (("  j: "ID"\n", j)) ;
01724       if (j == EMPTY)
01725       {
01726     /* Skip a dense variable.  It has no parent. */
01727     AMD_DEBUG3 (("      i is a dense variable\n")) ;
01728     continue ;
01729       }
01730 
01731       /* while (j is a variable) */
01732       while (Nv [j] == 0)
01733       {
01734     AMD_DEBUG3 (("    j : "ID"\n", j)) ;
01735     j = Pe [j] ;
01736     AMD_DEBUG3 (("    j:: "ID"\n", j)) ;
01737     ASSERT (j >= 0 && j < n) ;
01738       }
01739       /* got to an element e */
01740       e = j ;
01741       AMD_DEBUG3 (("got to e: "ID"\n", e)) ;
01742 
01743       /* -------------------------------------------------------------
01744        * traverse the path again from i to e, and compress the path
01745        * (all nodes point to e).  Path compression allows this code to
01746        * compute in O(n) time.
01747        * ------------------------------------------------------------- */
01748 
01749       j = i ;
01750       /* while (j is a variable) */
01751       while (Nv [j] == 0)
01752       {
01753     jnext = Pe [j] ;
01754     AMD_DEBUG3 (("j "ID" jnext "ID"\n", j, jnext)) ;
01755     Pe [j] = e ;
01756     j = jnext ;
01757     ASSERT (j >= 0 && j < n) ;
01758       }
01759   }
01760     }
01761 
01762 /* ========================================================================= */
01763 /* postorder the assembly tree */
01764 /* ========================================================================= */
01765 
01766     AMD_postorder (n, Pe, Nv, Elen,
01767   W,      /* output order */
01768   Head, Next, Last) ; /* workspace */
01769 
01770 /* ========================================================================= */
01771 /* compute output permutation and inverse permutation */
01772 /* ========================================================================= */
01773 
01774     /* W [e] = k means that element e is the kth element in the new
01775      * order.  e is in the range 0 to n-1, and k is in the range 0 to
01776      * the number of elements.  Use Head for inverse order. */
01777 
01778     for (k = 0 ; k < n ; k++)
01779     {
01780   Head [k] = EMPTY ;
01781   Next [k] = EMPTY ;
01782     }
01783     for (e = 0 ; e < n ; e++)
01784     {
01785   k = W [e] ;
01786   ASSERT ((k == EMPTY) == (Nv [e] == 0)) ;
01787   if (k != EMPTY)
01788   {
01789       ASSERT (k >= 0 && k < n) ;
01790       Head [k] = e ;
01791   }
01792     }
01793 
01794     /* construct output inverse permutation in Next,
01795      * and permutation in Last */
01796     nel = 0 ;
01797     for (k = 0 ; k < n ; k++)
01798     {
01799   e = Head [k] ;
01800   if (e == EMPTY) break ;
01801   ASSERT (e >= 0 && e < n && Nv [e] > 0) ;
01802   Next [e] = nel ;
01803   nel += Nv [e] ;
01804     }
01805     ASSERT (nel == n - ndense) ;
01806 
01807     /* order non-principal variables (dense, & those merged into supervar's) */
01808     for (i = 0 ; i < n ; i++)
01809     {
01810   if (Nv [i] == 0)
01811   {
01812       e = Pe [i] ;
01813       ASSERT (e >= EMPTY && e < n) ;
01814       if (e != EMPTY)
01815       {
01816     /* This is an unordered variable that was merged
01817      * into element e via supernode detection or mass
01818      * elimination of i when e became the pivot element.
01819      * Place i in order just before e. */
01820     ASSERT (Next [i] == EMPTY && Nv [e] > 0) ;
01821     Next [i] = Next [e] ;
01822     Next [e]++ ;
01823       }
01824       else
01825       {
01826     /* This is a dense unordered variable, with no parent.
01827      * Place it last in the output order. */
01828     Next [i] = nel++ ;
01829       }
01830   }
01831     }
01832     ASSERT (nel == n) ;
01833 
01834     AMD_DEBUG2 (("\n\nPerm:\n")) ;
01835     for (i = 0 ; i < n ; i++)
01836     {
01837   k = Next [i] ;
01838   ASSERT (k >= 0 && k < n) ;
01839   Last [k] = i ;
01840   AMD_DEBUG2 (("   perm ["ID"] = "ID"\n", k, i)) ;
01841     }
01842 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines