Mercurial > octave-nkf
diff liboctave/UMFPACK/AMD/Source/amdbar.f @ 5164:57077d0ddc8e
[project @ 2005-02-25 19:55:24 by jwe]
| author | jwe |
|---|---|
| date | Fri, 25 Feb 2005 19:55:28 +0000 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/UMFPACK/AMD/Source/amdbar.f Fri Feb 25 19:55:28 2005 +0000 @@ -0,0 +1,1206 @@ +C----------------------------------------------------------------------- +C AMDBAR: approximate minimum degree, without aggressive absorption +C----------------------------------------------------------------------- + + SUBROUTINE AMDBAR + $ (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 This routine does not do aggresive absorption (as done by AMD). + +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, 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 ---------------------------------------------------------- + +C UMFPACK/MA38-style approximate degree: + DO 160 P = P1, P2 + E = IW (P) + WE = W (E) + IF (WE .NE. 0) THEN +C e is an unabsorbed element + DEG = DEG + WE - WFLG + IW (PN) = E + PN = PN + 1 + HASH = HASH + E + 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 (ELEN (I) .EQ. 1 .AND. P3 .EQ. PN) 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! j can be absorbed into i +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 = MAX (1, 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 +