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