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