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