Mercurial > octave-nkf
view 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 source
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