C-----------------------------------------------------------------------
C AMD: approximate minimum degree, with aggressive absorption
C-----------------------------------------------------------------------
SUBROUTINE AMD
$ (N, PE, IW, LEN, IWLEN, PFREE, NV, NEXT,
$ LAST, HEAD, ELEN, DEGREE, NCMPA, W)
INTEGER N, IWLEN, PFREE, NCMPA, IW (IWLEN), PE (N),
$ DEGREE (N), NV (N), NEXT (N), LAST (N), HEAD (N),
$ ELEN (N), W (N), LEN (N)
C Given a representation of the nonzero pattern of a symmetric matrix,
C A, (excluding the diagonal) perform an approximate minimum
C (UMFPACK/MA38-style) degree ordering to compute a pivot order
C such that the introduction of nonzeros (fill-in) in the Cholesky
C factors A = LL^T are kept low. At each step, the pivot
C selected is the one with the minimum UMFPACK/MA38-style
C upper-bound on the external degree.
C
C Aggresive absorption is used to tighten the bound on the degree.
C **********************************************************************
C ***** CAUTION: ARGUMENTS ARE NOT CHECKED FOR ERRORS ON INPUT. ******
C **********************************************************************
C References:
C
C [1] Timothy A. Davis and Iain Duff, "An unsymmetric-pattern
C multifrontal method for sparse LU factorization", SIAM J.
C Matrix Analysis and Applications, vol. 18, no. 1, pp.
C 140-158. Discusses UMFPACK / MA38, which first introduced
C the approximate minimum degree used by this routine.
C
C [2] Patrick Amestoy, Timothy A. Davis, and Iain S. Duff, "An
C approximate degree ordering algorithm," SIAM J. Matrix
C Analysis and Applications, vol. 17, no. 4, pp. 886-905,
C 1996. Discusses AMD, AMDBAR, and MC47B.
C
C [3] Alan George and Joseph Liu, "The evolution of the minimum
C degree ordering algorithm," SIAM Review, vol. 31, no. 1,
C pp. 1-19, 1989. We list below the features mentioned in
C that paper that this code includes:
C
C mass elimination:
C Yes. MA27 relied on supervariable detection for mass
C elimination.
C indistinguishable nodes:
C Yes (we call these "supervariables"). This was also in
C the MA27 code - although we modified the method of
C detecting them (the previous hash was the true degree,
C which we no longer keep track of). A supervariable is
C a set of rows with identical nonzero pattern. All
C variables in a supervariable are eliminated together.
C Each supervariable has as its numerical name that of
C one of its variables (its principal variable).
C quotient graph representation:
C Yes. We use the term "element" for the cliques formed
C during elimination. This was also in the MA27 code.
C The algorithm can operate in place, but it will work
C more efficiently if given some "elbow room."
C element absorption:
C Yes. This was also in the MA27 code.
C external degree:
C Yes. The MA27 code was based on the true degree.
C incomplete degree update and multiple elimination:
C No. This was not in MA27, either. Our method of
C degree update within MC47B/BD is element-based, not
C variable-based. It is thus not well-suited for use
C with incomplete degree update or multiple elimination.
C-----------------------------------------------------------------------
C Authors, and Copyright (C) 1995 by:
C Timothy A. Davis, Patrick Amestoy, Iain S. Duff, & John K. Reid.
C
C Acknowledgements:
C This work (and the UMFPACK package) was supported by the
C National Science Foundation (ASC-9111263 and DMS-9223088).
C The UMFPACK/MA38 approximate degree update algorithm, the
C unsymmetric analog which forms the basis of MC47B/BD, was
C developed while Tim Davis was supported by CERFACS (Toulouse,
C France) in a post-doctoral position.
C
C Date: September, 1995
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C INPUT ARGUMENTS (unaltered):
C-----------------------------------------------------------------------
C n: The matrix order.
C
C Restriction: 1 .le. n .lt. (iovflo/2)-2, where iovflo is
C the largest positive integer that your computer can represent.
C iwlen: The length of iw (1..iwlen). On input, the matrix is
C stored in iw (1..pfree-1). However, iw (1..iwlen) should be
C slightly larger than what is required to hold the matrix, at
C least iwlen .ge. pfree + n is recommended. Otherwise,
C excessive compressions will take place.
C *** We do not recommend running this algorithm with ***
C *** iwlen .lt. pfree + n. ***
C *** Better performance will be obtained if ***
C *** iwlen .ge. pfree + n ***
C *** or better yet ***
C *** iwlen .gt. 1.2 * pfree ***
C *** (where pfree is its value on input). ***
C The algorithm will not run at all if iwlen .lt. pfree-1.
C
C Restriction: iwlen .ge. pfree-1
C-----------------------------------------------------------------------
C INPUT/OUPUT ARGUMENTS:
C-----------------------------------------------------------------------
C pe: On input, pe (i) is the index in iw of the start of row i, or
C zero if row i has no off-diagonal non-zeros.
C
C During execution, it is used for both supervariables and
C elements:
C
C * Principal supervariable i: index into iw of the
C description of supervariable i. A supervariable
C represents one or more rows of the matrix
C with identical nonzero pattern.
C * Non-principal supervariable i: if i has been absorbed
C into another supervariable j, then pe (i) = -j.
C That is, j has the same pattern as i.
C Note that j might later be absorbed into another
C supervariable j2, in which case pe (i) is still -j,
C and pe (j) = -j2.
C * Unabsorbed element e: the index into iw of the description
C of element e, if e has not yet been absorbed by a
C subsequent element. Element e is created when
C the supervariable of the same name is selected as
C the pivot.
C * Absorbed element e: if element e is absorbed into element
C e2, then pe (e) = -e2. This occurs when the pattern of
C e (that is, Le) is found to be a subset of the pattern
C of e2 (that is, Le2). If element e is "null" (it has
C no nonzeros outside its pivot block), then pe (e) = 0.
C
C On output, pe holds the assembly tree/forest, which implicitly
C represents a pivot order with identical fill-in as the actual
C order (via a depth-first search of the tree).
C
C On output:
C If nv (i) .gt. 0, then i represents a node in the assembly tree,
C and the parent of i is -pe (i), or zero if i is a root.
C If nv (i) = 0, then (i,-pe (i)) represents an edge in a
C subtree, the root of which is a node in the assembly tree.
C pfree: On input the tail end of the array, iw (pfree..iwlen),
C is empty, and the matrix is stored in iw (1..pfree-1).
C During execution, additional data is placed in iw, and pfree
C is modified so that iw (pfree..iwlen) is always the unused part
C of iw. On output, pfree is set equal to the size of iw that
C would have been needed for no compressions to occur. If
C ncmpa is zero, then pfree (on output) is less than or equal to
C iwlen, and the space iw (pfree+1 ... iwlen) was not used.
C Otherwise, pfree (on output) is greater than iwlen, and all the
C memory in iw was used.
C-----------------------------------------------------------------------
C INPUT/MODIFIED (undefined on output):
C-----------------------------------------------------------------------
C len: On input, len (i) holds the number of entries in row i of the
C matrix, excluding the diagonal. The contents of len (1..n)
C are undefined on output.
C iw: On input, iw (1..pfree-1) holds the description of each row i
C in the matrix. The matrix must be symmetric, and both upper
C and lower triangular parts must be present. The diagonal must
C not be present. Row i is held as follows:
C
C len (i): the length of the row i data structure
C iw (pe (i) ... pe (i) + len (i) - 1):
C the list of column indices for nonzeros
C in row i (simple supervariables), excluding
C the diagonal. All supervariables start with
C one row/column each (supervariable i is just
C row i).
C if len (i) is zero on input, then pe (i) is ignored
C on input.
C
C Note that the rows need not be in any particular order,
C and there may be empty space between the rows.
C
C During execution, the supervariable i experiences fill-in.
C This is represented by placing in i a list of the elements
C that cause fill-in in supervariable i:
C
C len (i): the length of supervariable i
C iw (pe (i) ... pe (i) + elen (i) - 1):
C the list of elements that contain i. This list
C is kept short by removing absorbed elements.
C iw (pe (i) + elen (i) ... pe (i) + len (i) - 1):
C the list of supervariables in i. This list
C is kept short by removing nonprincipal
C variables, and any entry j that is also
C contained in at least one of the elements
C (j in Le) in the list for i (e in row i).
C
C When supervariable i is selected as pivot, we create an
C element e of the same name (e=i):
C
C len (e): the length of element e
C iw (pe (e) ... pe (e) + len (e) - 1):
C the list of supervariables in element e.
C
C An element represents the fill-in that occurs when supervariable
C i is selected as pivot (which represents the selection of row i
C and all non-principal variables whose principal variable is i).
C We use the term Le to denote the set of all supervariables
C in element e. Absorbed supervariables and elements are pruned
C from these lists when computationally convenient.
C
C CAUTION: THE INPUT MATRIX IS OVERWRITTEN DURING COMPUTATION.
C The contents of iw are undefined on output.
C-----------------------------------------------------------------------
C OUTPUT (need not be set on input):
C-----------------------------------------------------------------------
C nv: During execution, abs (nv (i)) is equal to the number of rows
C that are represented by the principal supervariable i. If i is
C a nonprincipal variable, then nv (i) = 0. Initially,
C nv (i) = 1 for all i. nv (i) .lt. 0 signifies that i is a
C principal variable in the pattern Lme of the current pivot
C element me. On output, nv (e) holds the true degree of element
C e at the time it was created (including the diagonal part).
C ncmpa: The number of times iw was compressed. If this is
C excessive, then the execution took longer than what could have
C been. To reduce ncmpa, try increasing iwlen to be 10% or 20%
C larger than the value of pfree on input (or at least
C iwlen .ge. pfree + n). The fastest performance will be
C obtained when ncmpa is returned as zero. If iwlen is set to
C the value returned by pfree on *output*, then no compressions
C will occur.
C elen: See the description of iw above. At the start of execution,
C elen (i) is set to zero. During execution, elen (i) is the
C number of elements in the list for supervariable i. When e
C becomes an element, elen (e) = -nel is set, where nel is the
C current step of factorization. elen (i) = 0 is done when i
C becomes nonprincipal.
C
C For variables, elen (i) .ge. 0 holds until just before the
C permutation vectors are computed. For elements,
C elen (e) .lt. 0 holds.
C
C On output elen (1..n) holds the inverse permutation (the same
C as the 'INVP' argument in Sparspak). That is, if k = elen (i),
C then row i is the kth pivot row. Row i of A appears as the
C (elen(i))-th row in the permuted matrix, PAP^T.
C last: In a degree list, last (i) is the supervariable preceding i,
C or zero if i is the head of the list. In a hash bucket,
C last (i) is the hash key for i. last (head (hash)) is also
C used as the head of a hash bucket if head (hash) contains a
C degree list (see head, below).
C
C On output, last (1..n) holds the permutation (the same as the
C 'PERM' argument in Sparspak). That is, if i = last (k), then
C row i is the kth pivot row. Row last (k) of A is the k-th row
C in the permuted matrix, PAP^T.
C-----------------------------------------------------------------------
C LOCAL (not input or output - used only during execution):
C-----------------------------------------------------------------------
C degree: If i is a supervariable, then degree (i) holds the
C current approximation of the external degree of row i (an upper
C bound). The external degree is the number of nonzeros in row i,
C minus abs (nv (i)) (the diagonal part). The bound is equal to
C the external degree if elen (i) is less than or equal to two.
C
C We also use the term "external degree" for elements e to refer
C to |Le \ Lme|. If e is an element, then degree (e) holds |Le|,
C which is the degree of the off-diagonal part of the element e
C (not including the diagonal part).
C head: head is used for degree lists. head (deg) is the first
C supervariable in a degree list (all supervariables i in a
C degree list deg have the same approximate degree, namely,
C deg = degree (i)). If the list deg is empty then
C head (deg) = 0.
C
C During supervariable detection head (hash) also serves as a
C pointer to a hash bucket.
C If head (hash) .gt. 0, there is a degree list of degree hash.
C The hash bucket head pointer is last (head (hash)).
C If head (hash) = 0, then the degree list and hash bucket are
C both empty.
C If head (hash) .lt. 0, then the degree list is empty, and
C -head (hash) is the head of the hash bucket.
C After supervariable detection is complete, all hash buckets
C are empty, and the (last (head (hash)) = 0) condition is
C restored for the non-empty degree lists.
C next: next (i) is the supervariable following i in a link list, or
C zero if i is the last in the list. Used for two kinds of
C lists: degree lists and hash buckets (a supervariable can be
C in only one kind of list at a time).
C w: The flag array w determines the status of elements and
C variables, and the external degree of elements.
C
C for elements:
C if w (e) = 0, then the element e is absorbed
C if w (e) .ge. wflg, then w (e) - wflg is the size of
C the set |Le \ Lme|, in terms of nonzeros (the
C sum of abs (nv (i)) for each principal variable i that
C is both in the pattern of element e and NOT in the
C pattern of the current pivot element, me).
C if wflg .gt. w (e) .gt. 0, then e is not absorbed and has
C not yet been seen in the scan of the element lists in
C the computation of |Le\Lme| in loop 150 below.
C
C for variables:
C during supervariable detection, if w (j) .ne. wflg then j is
C not in the pattern of variable i
C
C The w array is initialized by setting w (i) = 1 for all i,
C and by setting wflg = 2. It is reinitialized if wflg becomes
C too large (to ensure that wflg+n does not cause integer
C overflow).
C-----------------------------------------------------------------------
C LOCAL INTEGERS:
C-----------------------------------------------------------------------
INTEGER DEG, DEGME, DEXT, DMAX, E, ELENME, ELN, HASH, HMOD, I,
$ ILAST, INEXT, J, JLAST, JNEXT, K, KNT1, KNT2, KNT3,
$ LENJ, LN, MAXMEM, ME, MEM, MINDEG, NEL, NEWMEM,
$ NLEFT, NVI, NVJ, NVPIV, SLENME, WE, WFLG, WNVI, X
C deg: the degree of a variable or element
C degme: size, |Lme|, of the current element, me (= degree (me))
C dext: external degree, |Le \ Lme|, of some element e
C dmax: largest |Le| seen so far
C e: an element
C elenme: the length, elen (me), of element list of pivotal var.
C eln: the length, elen (...), of an element list
C hash: the computed value of the hash function
C hmod: the hash function is computed modulo hmod = max (1,n-1)
C i: a supervariable
C ilast: the entry in a link list preceding i
C inext: the entry in a link list following i
C j: a supervariable
C jlast: the entry in a link list preceding j
C jnext: the entry in a link list, or path, following j
C k: the pivot order of an element or variable
C knt1: loop counter used during element construction
C knt2: loop counter used during element construction
C knt3: loop counter used during compression
C lenj: len (j)
C ln: length of a supervariable list
C maxmem: amount of memory needed for no compressions
C me: current supervariable being eliminated, and the
C current element created by eliminating that
C supervariable
C mem: memory in use assuming no compressions have occurred
C mindeg: current minimum degree
C nel: number of pivots selected so far
C newmem: amount of new memory needed for current pivot element
C nleft: n - nel, the number of nonpivotal rows/columns remaining
C nvi: the number of variables in a supervariable i (= nv (i))
C nvj: the number of variables in a supervariable j (= nv (j))
C nvpiv: number of pivots in current element
C slenme: number of variables in variable list of pivotal variable
C we: w (e)
C wflg: used for flagging the w array. See description of iw.
C wnvi: wflg - nv (i)
C x: either a supervariable or an element
C-----------------------------------------------------------------------
C LOCAL POINTERS:
C-----------------------------------------------------------------------
INTEGER P, P1, P2, P3, PDST, PEND, PJ, PME, PME1, PME2, PN, PSRC
C Any parameter (pe (...) or pfree) or local variable
C starting with "p" (for Pointer) is an index into iw,
C and all indices into iw use variables starting with
C "p." The only exception to this rule is the iwlen
C input argument.
C p: pointer into lots of things
C p1: pe (i) for some variable i (start of element list)
C p2: pe (i) + elen (i) - 1 for some var. i (end of el. list)
C p3: index of first supervariable in clean list
C pdst: destination pointer, for compression
C pend: end of memory to compress
C pj: pointer into an element or variable
C pme: pointer into the current element (pme1...pme2)
C pme1: the current element, me, is stored in iw (pme1...pme2)
C pme2: the end of the current element
C pn: pointer into a "clean" variable, also used to compress
C psrc: source pointer, for compression
C-----------------------------------------------------------------------
C FUNCTIONS CALLED:
C-----------------------------------------------------------------------
INTRINSIC MAX, MIN, MOD
C=======================================================================
C INITIALIZATIONS
C=======================================================================
WFLG = 2
MINDEG = 1
NCMPA = 0
NEL = 0
HMOD = MAX (1, N-1)
DMAX = 0
MEM = PFREE - 1
MAXMEM = MEM
ME = 0
DO 10 I = 1, N
LAST (I) = 0
HEAD (I) = 0
NV (I) = 1
W (I) = 1
ELEN (I) = 0
DEGREE (I) = LEN (I)
10 CONTINUE
C ----------------------------------------------------------------
C initialize degree lists and eliminate rows with no off-diag. nz.
C ----------------------------------------------------------------
DO 20 I = 1, N
DEG = DEGREE (I)
IF (DEG .GT. 0) THEN
C ----------------------------------------------------------
C place i in the degree list corresponding to its degree
C ----------------------------------------------------------
INEXT = HEAD (DEG)
IF (INEXT .NE. 0) LAST (INEXT) = I
NEXT (I) = INEXT
HEAD (DEG) = I
ELSE
C ----------------------------------------------------------
C we have a variable that can be eliminated at once because
C there is no off-diagonal non-zero in its row.
C ----------------------------------------------------------
NEL = NEL + 1
ELEN (I) = -NEL
PE (I) = 0
W (I) = 0
ENDIF
20 CONTINUE
C=======================================================================
C WHILE (selecting pivots) DO
C=======================================================================
30 CONTINUE
IF (NEL .LT. N) THEN
C=======================================================================
C GET PIVOT OF MINIMUM DEGREE
C=======================================================================
C -------------------------------------------------------------
C find next supervariable for elimination
C -------------------------------------------------------------
DO 40 DEG = MINDEG, N
ME = HEAD (DEG)
IF (ME .GT. 0) GOTO 50
40 CONTINUE
50 CONTINUE
MINDEG = DEG
C -------------------------------------------------------------
C remove chosen variable from link list
C -------------------------------------------------------------
INEXT = NEXT (ME)
IF (INEXT .NE. 0) LAST (INEXT) = 0
HEAD (DEG) = INEXT
C -------------------------------------------------------------
C me represents the elimination of pivots nel+1 to nel+nv(me).
C place me itself as the first in this set. It will be moved
C to the nel+nv(me) position when the permutation vectors are
C computed.
C -------------------------------------------------------------
ELENME = ELEN (ME)
ELEN (ME) = - (NEL + 1)
NVPIV = NV (ME)
NEL = NEL + NVPIV
C=======================================================================
C CONSTRUCT NEW ELEMENT
C=======================================================================
C -------------------------------------------------------------
C At this point, me is the pivotal supervariable. It will be
C converted into the current element. Scan list of the
C pivotal supervariable, me, setting tree pointers and
C constructing new list of supervariables for the new element,
C me. p is a pointer to the current position in the old list.
C -------------------------------------------------------------
C flag the variable "me" as being in Lme by negating nv (me)
NV (ME) = -NVPIV
DEGME = 0
IF (ELENME .EQ. 0) THEN
C ----------------------------------------------------------
C construct the new element in place
C ----------------------------------------------------------
PME1 = PE (ME)
PME2 = PME1 - 1
DO 60 P = PME1, PME1 + LEN (ME) - 1
I = IW (P)
NVI = NV (I)
IF (NVI .GT. 0) THEN
C ----------------------------------------------------
C i is a principal variable not yet placed in Lme.
C store i in new list
C ----------------------------------------------------
DEGME = DEGME + NVI
C flag i as being in Lme by negating nv (i)
NV (I) = -NVI
PME2 = PME2 + 1
IW (PME2) = I
C ----------------------------------------------------
C remove variable i from degree list.
C ----------------------------------------------------
ILAST = LAST (I)
INEXT = NEXT (I)
IF (INEXT .NE. 0) LAST (INEXT) = ILAST
IF (ILAST .NE. 0) THEN
NEXT (ILAST) = INEXT
ELSE
C i is at the head of the degree list
HEAD (DEGREE (I)) = INEXT
ENDIF
ENDIF
60 CONTINUE
C this element takes no new memory in iw:
NEWMEM = 0
ELSE
C ----------------------------------------------------------
C construct the new element in empty space, iw (pfree ...)
C ----------------------------------------------------------
P = PE (ME)
PME1 = PFREE
SLENME = LEN (ME) - ELENME
DO 120 KNT1 = 1, ELENME + 1
IF (KNT1 .GT. ELENME) THEN
C search the supervariables in me.
E = ME
PJ = P
LN = SLENME
ELSE
C search the elements in me.
E = IW (P)
P = P + 1
PJ = PE (E)
LN = LEN (E)
ENDIF
C -------------------------------------------------------
C search for different supervariables and add them to the
C new list, compressing when necessary. this loop is
C executed once for each element in the list and once for
C all the supervariables in the list.
C -------------------------------------------------------
DO 110 KNT2 = 1, LN
I = IW (PJ)
PJ = PJ + 1
NVI = NV (I)
IF (NVI .GT. 0) THEN
C -------------------------------------------------
C compress iw, if necessary
C -------------------------------------------------
IF (PFREE .GT. IWLEN) THEN
C prepare for compressing iw by adjusting
C pointers and lengths so that the lists being
C searched in the inner and outer loops contain
C only the remaining entries.
PE (ME) = P
LEN (ME) = LEN (ME) - KNT1
IF (LEN (ME) .EQ. 0) THEN
C nothing left of supervariable me
PE (ME) = 0
ENDIF
PE (E) = PJ
LEN (E) = LN - KNT2
IF (LEN (E) .EQ. 0) THEN
C nothing left of element e
PE (E) = 0
ENDIF
NCMPA = NCMPA + 1
C store first item in pe
C set first entry to -item
DO 70 J = 1, N
PN = PE (J)
IF (PN .GT. 0) THEN
PE (J) = IW (PN)
IW (PN) = -J
ENDIF
70 CONTINUE
C psrc/pdst point to source/destination
PDST = 1
PSRC = 1
PEND = PME1 - 1
C while loop:
80 CONTINUE
IF (PSRC .LE. PEND) THEN
C search for next negative entry
J = -IW (PSRC)
PSRC = PSRC + 1
IF (J .GT. 0) THEN
IW (PDST) = PE (J)
PE (J) = PDST
PDST = PDST + 1
C copy from source to destination
LENJ = LEN (J)
DO 90 KNT3 = 0, LENJ - 2
IW (PDST + KNT3) = IW (PSRC + KNT3)
90 CONTINUE
PDST = PDST + LENJ - 1
PSRC = PSRC + LENJ - 1
ENDIF
GOTO 80
ENDIF
C move the new partially-constructed element
P1 = PDST
DO 100 PSRC = PME1, PFREE - 1
IW (PDST) = IW (PSRC)
PDST = PDST + 1
100 CONTINUE
PME1 = P1
PFREE = PDST
PJ = PE (E)
P = PE (ME)
ENDIF
C -------------------------------------------------
C i is a principal variable not yet placed in Lme
C store i in new list
C -------------------------------------------------
DEGME = DEGME + NVI
C flag i as being in Lme by negating nv (i)
NV (I) = -NVI
IW (PFREE) = I
PFREE = PFREE + 1
C -------------------------------------------------
C remove variable i from degree link list
C -------------------------------------------------
ILAST = LAST (I)
INEXT = NEXT (I)
IF (INEXT .NE. 0) LAST (INEXT) = ILAST
IF (ILAST .NE. 0) THEN
NEXT (ILAST) = INEXT
ELSE
C i is at the head of the degree list
HEAD (DEGREE (I)) = INEXT
ENDIF
ENDIF
110 CONTINUE
IF (E .NE. ME) THEN
C set tree pointer and flag to indicate element e is
C absorbed into new element me (the parent of e is me)
PE (E) = -ME
W (E) = 0
ENDIF
120 CONTINUE
PME2 = PFREE - 1
C this element takes newmem new memory in iw (possibly zero)
NEWMEM = PFREE - PME1
MEM = MEM + NEWMEM
MAXMEM = MAX (MAXMEM, MEM)
ENDIF
C -------------------------------------------------------------
C me has now been converted into an element in iw (pme1..pme2)
C -------------------------------------------------------------
C degme holds the external degree of new element
DEGREE (ME) = DEGME
PE (ME) = PME1
LEN (ME) = PME2 - PME1 + 1
C -------------------------------------------------------------
C make sure that wflg is not too large. With the current
C value of wflg, wflg+n must not cause integer overflow
C -------------------------------------------------------------
IF (WFLG + N .LE. WFLG) THEN
DO 130 X = 1, N
IF (W (X) .NE. 0) W (X) = 1
130 CONTINUE
WFLG = 2
ENDIF
C=======================================================================
C COMPUTE (w (e) - wflg) = |Le\Lme| FOR ALL ELEMENTS
C=======================================================================
C -------------------------------------------------------------
C Scan 1: compute the external degrees of previous elements
C with respect to the current element. That is:
C (w (e) - wflg) = |Le \ Lme|
C for each element e that appears in any supervariable in Lme.
C The notation Le refers to the pattern (list of
C supervariables) of a previous element e, where e is not yet
C absorbed, stored in iw (pe (e) + 1 ... pe (e) + iw (pe (e))).
C The notation Lme refers to the pattern of the current element
C (stored in iw (pme1..pme2)). If (w (e) - wflg) becomes
C zero, then the element e will be absorbed in scan 2.
C -------------------------------------------------------------
DO 150 PME = PME1, PME2
I = IW (PME)
ELN = ELEN (I)
IF (ELN .GT. 0) THEN
C note that nv (i) has been negated to denote i in Lme:
NVI = -NV (I)
WNVI = WFLG - NVI
DO 140 P = PE (I), PE (I) + ELN - 1
E = IW (P)
WE = W (E)
IF (WE .GE. WFLG) THEN
C unabsorbed element e has been seen in this loop
WE = WE - NVI
ELSE IF (WE .NE. 0) THEN
C e is an unabsorbed element
C this is the first we have seen e in all of Scan 1
WE = DEGREE (E) + WNVI
ENDIF
W (E) = WE
140 CONTINUE
ENDIF
150 CONTINUE
C=======================================================================
C DEGREE UPDATE AND ELEMENT ABSORPTION
C=======================================================================
C -------------------------------------------------------------
C Scan 2: for each i in Lme, sum up the degree of Lme (which
C is degme), plus the sum of the external degrees of each Le
C for the elements e appearing within i, plus the
C supervariables in i. Place i in hash list.
C -------------------------------------------------------------
DO 180 PME = PME1, PME2
I = IW (PME)
P1 = PE (I)
P2 = P1 + ELEN (I) - 1
PN = P1
HASH = 0
DEG = 0
C ----------------------------------------------------------
C scan the element list associated with supervariable i
C ----------------------------------------------------------
DO 160 P = P1, P2
E = IW (P)
C dext = | Le \ Lme |
DEXT = W (E) - WFLG
IF (DEXT .GT. 0) THEN
DEG = DEG + DEXT
IW (PN) = E
PN = PN + 1
HASH = HASH + E
ELSE IF (DEXT .EQ. 0) THEN
C aggressive absorption: e is not adjacent to me, but
C the |Le \ Lme| is 0, so absorb it into me
PE (E) = -ME
W (E) = 0
ELSE
C element e has already been absorbed, due to
C regular absorption, in do loop 120 above. Ignore it.
CONTINUE
ENDIF
160 CONTINUE
C count the number of elements in i (including me):
ELEN (I) = PN - P1 + 1
C ----------------------------------------------------------
C scan the supervariables in the list associated with i
C ----------------------------------------------------------
P3 = PN
DO 170 P = P2 + 1, P1 + LEN (I) - 1
J = IW (P)
NVJ = NV (J)
IF (NVJ .GT. 0) THEN
C j is unabsorbed, and not in Lme.
C add to degree and add to new list
DEG = DEG + NVJ
IW (PN) = J
PN = PN + 1
HASH = HASH + J
ENDIF
170 CONTINUE
C ----------------------------------------------------------
C update the degree and check for mass elimination
C ----------------------------------------------------------
IF (DEG .EQ. 0) THEN
C -------------------------------------------------------
C mass elimination
C -------------------------------------------------------
C There is nothing left of this node except for an
C edge to the current pivot element. elen (i) is 1,
C and there are no variables adjacent to node i.
C Absorb i into the current pivot element, me.
PE (I) = -ME
NVI = -NV (I)
DEGME = DEGME - NVI
NVPIV = NVPIV + NVI
NEL = NEL + NVI
NV (I) = 0
ELEN (I) = 0
ELSE
C -------------------------------------------------------
C update the upper-bound degree of i
C -------------------------------------------------------
C the following degree does not yet include the size
C of the current element, which is added later:
DEGREE (I) = MIN (DEGREE (I), DEG)
C -------------------------------------------------------
C add me to the list for i
C -------------------------------------------------------
C move first supervariable to end of list
IW (PN) = IW (P3)
C move first element to end of element part of list
IW (P3) = IW (P1)
C add new element to front of list.
IW (P1) = ME
C store the new length of the list in len (i)
LEN (I) = PN - P1 + 1
C -------------------------------------------------------
C place in hash bucket. Save hash key of i in last (i).
C -------------------------------------------------------
HASH = MOD (HASH, HMOD) + 1
J = HEAD (HASH)
IF (J .LE. 0) THEN
C the degree list is empty, hash head is -j
NEXT (I) = -J
HEAD (HASH) = -I
ELSE
C degree list is not empty
C use last (head (hash)) as hash head
NEXT (I) = LAST (J)
LAST (J) = I
ENDIF
LAST (I) = HASH
ENDIF
180 CONTINUE
DEGREE (ME) = DEGME
C -------------------------------------------------------------
C Clear the counter array, w (...), by incrementing wflg.
C -------------------------------------------------------------
DMAX = MAX (DMAX, DEGME)
WFLG = WFLG + DMAX
C make sure that wflg+n does not cause integer overflow
IF (WFLG + N .LE. WFLG) THEN
DO 190 X = 1, N
IF (W (X) .NE. 0) W (X) = 1
190 CONTINUE
WFLG = 2
ENDIF
C at this point, w (1..n) .lt. wflg holds
C=======================================================================
C SUPERVARIABLE DETECTION
C=======================================================================
DO 250 PME = PME1, PME2
I = IW (PME)
IF (NV (I) .LT. 0) THEN
C i is a principal variable in Lme
C -------------------------------------------------------
C examine all hash buckets with 2 or more variables. We
C do this by examing all unique hash keys for super-
C variables in the pattern Lme of the current element, me
C -------------------------------------------------------
HASH = LAST (I)
C let i = head of hash bucket, and empty the hash bucket
J = HEAD (HASH)
IF (J .EQ. 0) GOTO 250
IF (J .LT. 0) THEN
C degree list is empty
I = -J
HEAD (HASH) = 0
ELSE
C degree list is not empty, restore last () of head
I = LAST (J)
LAST (J) = 0
ENDIF
IF (I .EQ. 0) GOTO 250
C while loop:
200 CONTINUE
IF (NEXT (I) .NE. 0) THEN
C ----------------------------------------------------
C this bucket has one or more variables following i.
C scan all of them to see if i can absorb any entries
C that follow i in hash bucket. Scatter i into w.
C ----------------------------------------------------
LN = LEN (I)
ELN = ELEN (I)
C do not flag the first element in the list (me)
DO 210 P = PE (I) + 1, PE (I) + LN - 1
W (IW (P)) = WFLG
210 CONTINUE
C ----------------------------------------------------
C scan every other entry j following i in bucket
C ----------------------------------------------------
JLAST = I
J = NEXT (I)
C while loop:
220 CONTINUE
IF (J .NE. 0) THEN
C -------------------------------------------------
C check if j and i have identical nonzero pattern
C -------------------------------------------------
IF (LEN (J) .NE. LN) THEN
C i and j do not have same size data structure
GOTO 240
ENDIF
IF (ELEN (J) .NE. ELN) THEN
C i and j do not have same number of adjacent el
GOTO 240
ENDIF
C do not flag the first element in the list (me)
DO 230 P = PE (J) + 1, PE (J) + LN - 1
IF (W (IW (P)) .NE. WFLG) THEN
C an entry (iw(p)) is in j but not in i
GOTO 240
ENDIF
230 CONTINUE
C -------------------------------------------------
C found it
C -------------------------------------------------
PE (J) = -I
C both nv (i) and nv (j) are negated since they
C are in Lme, and the absolute values of each
C are the number of variables in i and j:
NV (I) = NV (I) + NV (J)
NV (J) = 0
ELEN (J) = 0
C delete j from hash bucket
J = NEXT (J)
NEXT (JLAST) = J
GOTO 220
C -------------------------------------------------
240 CONTINUE
C j cannot be absorbed into i
C -------------------------------------------------
JLAST = J
J = NEXT (J)
GOTO 220
ENDIF
C ----------------------------------------------------
C no more variables can be absorbed into i
C go to next i in bucket and clear flag array
C ----------------------------------------------------
WFLG = WFLG + 1
I = NEXT (I)
IF (I .NE. 0) GOTO 200
ENDIF
ENDIF
250 CONTINUE
C=======================================================================
C RESTORE DEGREE LISTS AND REMOVE NONPRINCIPAL SUPERVAR. FROM ELEMENT
C=======================================================================
P = PME1
NLEFT = N - NEL
DO 260 PME = PME1, PME2
I = IW (PME)
NVI = -NV (I)
IF (NVI .GT. 0) THEN
C i is a principal variable in Lme
C restore nv (i) to signify that i is principal
NV (I) = NVI
C -------------------------------------------------------
C compute the external degree (add size of current elem)
C -------------------------------------------------------
DEG = MIN (DEGREE (I) + DEGME - NVI, NLEFT - NVI)
C -------------------------------------------------------
C place the supervariable at the head of the degree list
C -------------------------------------------------------
INEXT = HEAD (DEG)
IF (INEXT .NE. 0) LAST (INEXT) = I
NEXT (I) = INEXT
LAST (I) = 0
HEAD (DEG) = I
C -------------------------------------------------------
C save the new degree, and find the minimum degree
C -------------------------------------------------------
MINDEG = MIN (MINDEG, DEG)
DEGREE (I) = DEG
C -------------------------------------------------------
C place the supervariable in the element pattern
C -------------------------------------------------------
IW (P) = I
P = P + 1
ENDIF
260 CONTINUE
C=======================================================================
C FINALIZE THE NEW ELEMENT
C=======================================================================
NV (ME) = NVPIV + DEGME
C nv (me) is now the degree of pivot (including diagonal part)
C save the length of the list for the new element me
LEN (ME) = P - PME1
IF (LEN (ME) .EQ. 0) THEN
C there is nothing left of the current pivot element
PE (ME) = 0
W (ME) = 0
ENDIF
IF (NEWMEM .NE. 0) THEN
C element was not constructed in place: deallocate part
C of it (final size is less than or equal to newmem,
C since newly nonprincipal variables have been removed).
PFREE = P
MEM = MEM - NEWMEM + LEN (ME)
ENDIF
C=======================================================================
C END WHILE (selecting pivots)
GOTO 30
ENDIF
C=======================================================================
C=======================================================================
C COMPUTE THE PERMUTATION VECTORS
C=======================================================================
C ----------------------------------------------------------------
C The time taken by the following code is O(n). At this
C point, elen (e) = -k has been done for all elements e,
C and elen (i) = 0 has been done for all nonprincipal
C variables i. At this point, there are no principal
C supervariables left, and all elements are absorbed.
C ----------------------------------------------------------------
C ----------------------------------------------------------------
C compute the ordering of unordered nonprincipal variables
C ----------------------------------------------------------------
DO 290 I = 1, N
IF (ELEN (I) .EQ. 0) THEN
C ----------------------------------------------------------
C i is an un-ordered row. Traverse the tree from i until
C reaching an element, e. The element, e, was the
C principal supervariable of i and all nodes in the path
C from i to when e was selected as pivot.
C ----------------------------------------------------------
J = -PE (I)
C while (j is a variable) do:
270 CONTINUE
IF (ELEN (J) .GE. 0) THEN
J = -PE (J)
GOTO 270
ENDIF
E = J
C ----------------------------------------------------------
C get the current pivot ordering of e
C ----------------------------------------------------------
K = -ELEN (E)
C ----------------------------------------------------------
C traverse the path again from i to e, and compress the
C path (all nodes point to e). Path compression allows
C this code to compute in O(n) time. Order the unordered
C nodes in the path, and place the element e at the end.
C ----------------------------------------------------------
J = I
C while (j is a variable) do:
280 CONTINUE
IF (ELEN (J) .GE. 0) THEN
JNEXT = -PE (J)
PE (J) = -E
IF (ELEN (J) .EQ. 0) THEN
C j is an unordered row
ELEN (J) = K
K = K + 1
ENDIF
J = JNEXT
GOTO 280
ENDIF
C leave elen (e) negative, so we know it is an element
ELEN (E) = -K
ENDIF
290 CONTINUE
C ----------------------------------------------------------------
C reset the inverse permutation (elen (1..n)) to be positive,
C and compute the permutation (last (1..n)).
C ----------------------------------------------------------------
DO 300 I = 1, N
K = ABS (ELEN (I))
LAST (K) = I
ELEN (I) = K
300 CONTINUE
C=======================================================================
C RETURN THE MEMORY USAGE IN IW
C=======================================================================
C If maxmem is less than or equal to iwlen, then no compressions
C occurred, and iw (maxmem+1 ... iwlen) was unused. Otherwise
C compressions did occur, and iwlen would have had to have been
C greater than or equal to maxmem for no compressions to occur.
C Return the value of maxmem in the pfree argument.
PFREE = MAXMEM
RETURN
END
