      SUBROUTINE REDUCE(NDSTK, NR, IOLD, RENUM, NDEG, LVL, LVLS1,       
     * LVLS2, CCSTOR, IBW2, IPF2)
C  SUBROUTINE REDUCE DETERMINES A ROW AND COLUMN PERMUTATION WHICH,
C  WHEN APPLIED TO A GIVEN SPARSE MATRIX, PRODUCES A PERMUTED
C  MATRIX WITH A SMALLER BANDWIDTH AND PROFILE.
C  THE INPUT ARRAY IS A CONNECTION TABLE WHICH REPRESENTS THE
C  INDICES OF THE NONZERO ELEMENTS OF THE MATRIX, A.  THE ALGO-
C  RITHM IS DESCRIBED IN TERMS OF THE ADJACENCY GRAPH WHICH
C  HAS THE CHARACTERISTIC THAT THERE IS AN EDGE (CONNECTION)
C  BETWEEN NODES I AND J IF A(I,J) .NE. 0 AND I .NE. J.
C  DIMENSIONING INFORMATION--THE FOLLOWING INTEGER ARRAYS MUST BE
C  DIMENSIONED IN THE CALLING ROUTINE.
C    NDSTK(NR,D1)        D1 IS .GE. MAXIMUM DEGREE OF ALL NODES.
C    IOLD(D2)            D2 AND NR ARE .GE. THE TOTAL NUMBER OF
C    RENUM(D2+1)         NODES IN THE GRAPH.
C    NDEG(D2)            STORAGE REQUIREMENTS CAN BE SIGNIFICANTLY
C    LVL(D2)             DECREASED FOR IBM 360 AND 370 COMPUTERS
C    LVLS1(D2)           BY REPLACING INTEGER NDSTK BY
C    LVLS2(D2)           INTEGER*2 NDSTK IN SUBROUTINES REDUCE,
C    CCSTOR(D2)          DGREE, FNDIAM, TREE AND NUMBER.
C  COMMON INFORMATION--THE FOLLOWING COMMON BLOCK MUST BE IN THE
C  CALLING ROUTINE.
C    COMMON/GRA/N,IDPTH,IDEG
C  EXPLANATION OF INPUT VARIABLES--
C    NDSTK-     CONNECTION TABLE REPRESENTING GRAPH.
C               NDSTK(I,J)=NODE NUMBER OF JTH CONNECTION TO NODE
C               NUMBER I.  A CONNECTION OF A NODE TO ITSELF IS NOT
C               LISTED.  EXTRA POSITIONS MUST HAVE ZERO FILL.
C    NR-        ROW DIMENSION ASSIGNED NDSTK IN CALLING PROGRAM.
C    IOLD(I)-   NUMBERING OF ITH NODE UPON INPUT.
C               IF NO NUMBERING EXISTS THEN IOLD(I)=I.
C    N-         NUMBER OF NODES IN GRAPH (EQUAL TO ORDER OF MATRIX).
C    IDEG-      MAXIMUM DEGREE OF ANY NODE IN THE GRAPH.
C  EXPLANATION OF OUTPUT VARIABLES--
C    RENUM(I)-  THE NEW NUMBER FOR THE ITH NODE.
C    NDEG(I)-   THE DEGREE OF THE ITH NODE.
C    IBW2-      THE BANDWIDTH AFTER RENUMBERING.
C    IPF2-      THE PROFILE AFTER RENUMBERING.
C    IDPTH-     NUMBER OF LEVELS IN REDUCE LEVEL STRUCTURE.
C  THE FOLLOWING ONLY HAVE MEANING IF THE GRAPH WAS CONNECTED--
C    LVL(I)-    INDEX INTO LVLS1 TO THE FIRST NODE IN LEVEL I.
C               LVL(I+1)-LVL(I)= NUMBER OF NODES IN ITH LEVEL
C    LVLS1-     NODE NUMBERS LISTED BY LEVEL.
C    LVLS2(I)-  THE LEVEL ASSIGNED TO NODE I BY REDUCE.
C  WORKING STORAGE VARIABLE--
C    CCSTOR
C  LOCAL STORAGE--
C    COMMON/CC/-SUBROUTINES REDUCE, SORT2 AND PIKLVL ASSUME THAT
C               THE GRAPH HAS AT MOST 50 CONNECTED COMPONENTS.
C               SUBROUTINE FNDIAM ASSUMES THAT THERE ARE AT MOST
C               100 NODES IN THE LAST LEVEL.
C    COMMON/LVLW/-SUBROUTINES SETUP AND PIKLVL ASSUME THAT THERE
C               ARE AT MOST 100 LEVELS.
C USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
      INTEGER NDSTK
      INTEGER STNODE, RVNODE, RENUM, XC, SORT2, STNUM, CCSTOR,
     * SIZE, STPT, SBNUM
      COMMON /GRA/ N, IDPTH, IDEG
C IT IS ASSUMED THAT THE GRAPH HAS AT MOST 50 CONNECTED COMPONENTS.
      COMMON /CC/ XC, SIZE(5000), STPT(5000)
      COMMON /LVLW/ NHIGH(10000), NLOW(10000), NACUM(10000)
      save gra, cc, lvlw
      DIMENSION CCSTOR(1), IOLD(1)
      DIMENSION NDSTK(NR,1), LVL(1), LVLS1(1), LVLS2(1), RENUM(1),
     * NDEG(1)
      IBW2 = 0
      IPF2 = 0
C SET RENUM(I)=0 FOR ALL I TO INDICATE NODE I IS UNNUMBERED
      DO 10 I=1,N
        RENUM(I) = 0
   10 CONTINUE
C COMPUTE DEGREE OF EACH NODE AND ORIGINAL BANDWIDTH AND PROFILE
      CALL DGREE(NDSTK, NR, NDEG, IOLD, IBW1, IPF1)
C SBNUM= LOW END OF AVAILABLE NUMBERS FOR RENUMBERING
C STNUM= HIGH END OF AVAILABLE NUMBERS FOR RENUMBERING
      SBNUM = 1
      STNUM = N
C NUMBER THE NODES OF DEGREE ZERO
      DO 20 I=1,N
        IF (NDEG(I).GT.0) GO TO 20
        RENUM(I) = STNUM
        STNUM = STNUM - 1
   20 CONTINUE
C FIND AN UNNUMBERED NODE OF MIN DEGREE TO START ON
   30 LOWDG = IDEG + 1
      NFLG = 1
      ISDIR = 1
      DO 40 I=1,N
        IF (NDEG(I).GE.LOWDG) GO TO 40
        IF (RENUM(I).GT.0) GO TO 40
        LOWDG = NDEG(I)
        STNODE = I
   40 CONTINUE
C FIND PSEUDO-DIAMETER AND ASSOCIATED LEVEL STRUCTURES.
C STNODE AND RVNODE ARE THE ENDS OF THE DIAM AND LVLS1 AND LVLS2
C ARE THE RESPECTIVE LEVEL STRUCTURES.
      CALL FNDIAM(STNODE, RVNODE, NDSTK, NR, NDEG, LVL, LVLS1,
     * LVLS2, CCSTOR, IDFLT)
      IF (NDEG(STNODE).LE.NDEG(RVNODE)) GO TO 50
C NFLG INDICATES THE END TO BEGIN NUMBERING ON
      NFLG = -1
      STNODE = RVNODE
   50 CALL GPS_SETUP(LVL, LVLS1, LVLS2)
C FIND ALL THE CONNECTED COMPONENTS  (XC COUNTS THEM)
      XC = 0
      LROOT = 1
      LVLN = 1
      DO 60 I=1,N
        IF (LVL(I).NE.0) GO TO 60
        XC = XC + 1
        STPT(XC) = LROOT
        CALL TREE(I, NDSTK, NR, LVL, CCSTOR, NDEG, LVLWTH, LVLBOT,
     *   LVLN, MAXLW, N)
        SIZE(XC) = LVLBOT + LVLWTH - LROOT
        LROOT = LVLBOT + LVLWTH
        LVLN = LROOT
   60 CONTINUE
      IF (SORT2(DMY).EQ.0) GO TO 70
      CALL PIKLVL(LVLS1, LVLS2, CCSTOR, IDFLT, ISDIR)
C ON RETURN FROM PIKLVL, ISDIR INDICATES THE DIRECTION THE LARGEST
C COMPONENT FELL.  ISDIR IS MODIFIED NOW TO INDICATE THE NUMBERING
C DIRECTION.  NUM IS SET TO THE PROPER VALUE FOR THIS DIRECTION.
   70 ISDIR = ISDIR*NFLG
      NUM = SBNUM
      IF (ISDIR.LT.0) NUM = STNUM
      CALL NUMBER(STNODE, NUM, NDSTK, LVLS2, NDEG, RENUM, LVLS1,
     * LVL, NR, NFLG, IBW2, IPF2, CCSTOR, ISDIR)
C UPDATE STNUM OR SBNUM AFTER NUMBERING
      IF (ISDIR.LT.0) STNUM = NUM
      IF (ISDIR.GT.0) SBNUM = NUM
      IF (SBNUM.LE.STNUM) GO TO 30
      IF (IBW2.LE.IBW1) RETURN
C IF ORIGINAL NUMBERING IS BETTER THAN NEW ONE, SET UP TO RETURN IT
      DO 80 I=1,N
        RENUM(I) = IOLD(I)
   80 CONTINUE
      IBW2 = IBW1
      IPF2 = IPF1
      RETURN
      END
      SUBROUTINE DGREE(NDSTK, NR, NDEG, IOLD, IBW1, IPF1)               
C  DGREE COMPUTES THE DEGREE OF EACH NODE IN NDSTK AND STORES
C  IT IN THE ARRAY NDEG.  THE BANDWIDTH AND PROFILE FOR THE ORIGINAL
C  OR INPUT RENUMBERING OF THE GRAPH IS COMPUTED ALSO.
C USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
      INTEGER NDSTK
      COMMON /GRA/ N, IDPTH, IDEG
      DIMENSION NDSTK(NR,1), NDEG(1), IOLD(1)
      IBW1 = 0
      IPF1 = 0
      DO 40 I=1,N
        NDEG(I) = 0
        IRW = 0
        DO 20 J=1,IDEG
          ITST = NDSTK(I,J)
          IF (ITST) 30, 30, 10
   10     NDEG(I) = NDEG(I) + 1
          IDIF = IOLD(I) - IOLD(ITST)
          IF (IRW.LT.IDIF) IRW = IDIF
   20   CONTINUE
   30   IPF1 = IPF1 + IRW
        IF (IRW.GT.IBW1) IBW1 = IRW
   40 CONTINUE
      RETURN
      END
      SUBROUTINE FNDIAM(SND1, SND2, NDSTK, NR, NDEG, LVL, LVLS1,        
     * LVLS2, IWK, IDFLT)
C  FNDIAM IS THE CONTROL PROCEDURE FOR FINDING THE PSEUDO-DIAMETER OF
C  NDSTK AS WELL AS THE LEVEL STRUCTURE FROM EACH END
C  SND1-        ON INPUT THIS IS THE NODE NUMBER OF THE FIRST
C               ATTEMPT AT FINDING A DIAMETER.  ON OUTPUT IT
C               CONTAINS THE ACTUAL NUMBER USED.
C  SND2-        ON OUTPUT CONTAINS OTHER END OF DIAMETER
C  LVLS1-       ARRAY CONTAINING LEVEL STRUCTURE WITH SND1 AS ROOT
C  LVLS2-       ARRAY CONTAINING LEVEL STRUCTURE WITH SND2 AS ROOT
C  IDFLT-       FLAG USED IN PICKING FINAL LEVEL STRUCTURE, SET
C               =1 IF WIDTH OF LVLS1 .LE. WIDTH OF LVLS2, OTHERWISE =2
C  LVL,IWK-     WORKING STORAGE
C USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
      INTEGER NDSTK
      INTEGER FLAG, SND, SND1, SND2
      COMMON /GRA/ N, IDPTH, IDEG
C IT IS ASSUMED THAT THE LAST LEVEL HAS AT MOST 100 NODES.
      COMMON /CC/ NDLST(10001)
      DIMENSION NDSTK(NR,1), NDEG(1), LVL(1), LVLS1(1), LVLS2(1),
     * IWK(1)
      FLAG = 0
      MTW2 = N
      SND = SND1
C ZERO LVL TO INDICATE ALL NODES ARE AVAILABLE TO TREE
   10 DO 20 I=1,N
        LVL(I) = 0
   20 CONTINUE
      LVLN = 1
C DROP A TREE FROM SND
      CALL TREE(SND, NDSTK, NR, LVL, IWK, NDEG, LVLWTH, LVLBOT,
     * LVLN, MAXLW, MTW2)
      IF (FLAG.GE.1) GO TO 50
      FLAG = 1
   30 IDPTH = LVLN - 1
      MTW1 = MAXLW
C COPY LEVEL STRUCTURE INTO LVLS1
      DO 40 I=1,N
        LVLS1(I) = LVL(I)
   40 CONTINUE
      NDXN = 1
      NDXL = 0
      MTW2 = N
C SORT LAST LEVEL BY DEGREE  AND STORE IN NDLST
      CALL SORTDG(NDLST, IWK(LVLBOT), NDXL, LVLWTH, NDEG)
      SND = NDLST(1)
      GO TO 10
   50 IF (IDPTH.GE.LVLN-1) GO TO 60
C START AGAIN WITH NEW STARTING NODE
      SND1 = SND
      GO TO 30
   60 IF (MAXLW.GE.MTW2) GO TO 80
      MTW2 = MAXLW
      SND2 = SND
C STORE NARROWEST REVERSE LEVEL STRUCTURE IN LVLS2
      DO 70 I=1,N
        LVLS2(I) = LVL(I)
   70 CONTINUE
   80 IF (NDXN.EQ.NDXL) GO TO 90
C TRY NEXT NODE IN NDLST
      NDXN = NDXN + 1
      SND = NDLST(NDXN)
      GO TO 10
   90 IDFLT = 1
      IF (MTW2.LE.MTW1) IDFLT = 2
      RETURN
      END
      SUBROUTINE TREE(IROOT, NDSTK, NR, LVL, IWK, NDEG, LVLWTH,         
     * LVLBOT, LVLN, MAXLW, IBORT)
C  TREE DROPS A TREE IN NDSTK FROM IROOT
C  LVL-         ARRAY INDICATING AVAILABLE NODES IN NDSTK WITH ZERO
C               ENTRIES. TREE ENTERS LEVEL NUMBERS ASSIGNED
C               DURING EXECUTION OF THIS PROCEDURE
C  IWK-         ON OUTPUT CONTAINS NODE NUMBERS USED IN TREE
C               ARRANGED BY LEVELS (IWK(LVLN) CONTAINS IROOT
C               AND IWK(LVLBOT+LVLWTH-1) CONTAINS LAST NODE ENTERED)
C  LVLWTH-      ON OUTPUT CONTAINS WIDTH OF LAST LEVEL
C  LVLBOT-      ON OUTPUT CONTAINS INDEX INTO IWK OF FIRST
C               NODE IN LAST LEVEL
C  MAXLW-       ON OUTPUT CONTAINS THE MAXIMUM LEVEL WIDTH
C  LVLN-        ON INPUT THE FIRST AVAILABLE LOCATION IN IWK
C               USUALLY ONE BUT IF IWK IS USED TO STORE PREVIOUS
C               CONNECTED COMPONENTS, LVLN IS NEXT AVAILABLE LOCATION.
C               ON OUTPUT THE TOTAL NUMBER OF LEVELS + 1
C  IBORT-       INPUT PARAM WHICH TRIGGERS EARLY RETURN IF
C               MAXLW BECOMES .GE. IBORT
C USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
      INTEGER NDSTK
      DIMENSION NDSTK(NR,1), LVL(1), IWK(1), NDEG(1)
      MAXLW = 0
      ITOP = LVLN
      INOW = LVLN
      LVLBOT = LVLN
      LVLTOP = LVLN + 1
      LVLN = 1
      LVL(IROOT) = 1
      IWK(ITOP) = IROOT
   10 LVLN = LVLN + 1
   20 IWKNOW = IWK(INOW)
      NDROW = NDEG(IWKNOW)
      DO 30 J=1,NDROW
        ITEST = NDSTK(IWKNOW,J)
        IF (LVL(ITEST).NE.0) GO TO 30
        LVL(ITEST) = LVLN
        ITOP = ITOP + 1
        IWK(ITOP) = ITEST
   30 CONTINUE
      INOW = INOW + 1
      IF (INOW.LT.LVLTOP) GO TO 20
      LVLWTH = LVLTOP - LVLBOT
      IF (MAXLW.LT.LVLWTH) MAXLW = LVLWTH
      IF (MAXLW.GE.IBORT) RETURN
      IF (ITOP.LT.LVLTOP) RETURN
      LVLBOT = INOW
      LVLTOP = ITOP + 1
      GO TO 10
      END
      SUBROUTINE SORTDG(STK1, STK2, X1, X2, NDEG)                       
C SORTDG SORTS STK2 BY DEGREE OF THE NODE AND ADDS IT TO THE END
C OF STK1 IN ORDER OF LOWEST TO HIGHEST DEGREE.  X1 AND X2 ARE THE
C NUMBER OF NODES IN STK1 AND STK2 RESPECTIVELY.
      INTEGER X1, X2, STK1, STK2, TEMP
      COMMON /GRA/ N, IDPTH, IDEG
      DIMENSION NDEG(1), STK1(1), STK2(1)
      IND = X2
   10 ITEST = 0
      IND = IND - 1
      IF (IND.LT.1) GO TO 30
      DO 20 I=1,IND
        J = I + 1
        ISTK2 = STK2(I)
        JSTK2 = STK2(J)
        IF (NDEG(ISTK2).LE.NDEG(JSTK2)) GO TO 20
        ITEST = 1
        TEMP = STK2(I)
        STK2(I) = STK2(J)
        STK2(J) = TEMP
   20 CONTINUE
      IF (ITEST.EQ.1) GO TO 10
   30 DO 40 I=1,X2
        X1 = X1 + 1
        STK1(X1) = STK2(I)
   40 CONTINUE
      RETURN
      END
      SUBROUTINE GPS_SETUP(LVL, LVLS1, LVLS2)                               
C SETUP COMPUTES THE REVERSE LEVELING INFO FROM LVLS2 AND STORES
C IT INTO LVLS2.  NACUM(I) IS INITIALIZED TO NODES/ITH LEVEL FOR NODES
C ON THE PSEUDO-DIAMETER OF THE GRAPH.  LVL IS INITIALIZED TO NON-
C ZERO FOR NODES ON THE PSEUDO-DIAM AND NODES IN A DIFFERENT
C COMPONENT OF THE GRAPH.
      COMMON /GRA/ N, IDPTH, IDEG
C IT IS ASSUMED THAT THERE ARE AT MOST 100 LEVELS.
      COMMON /LVLW/ NHIGH(10000), NLOW(10000), NACUM(10000)
      DIMENSION LVL(1), LVLS1(1), LVLS2(1)
      DO 10 I=1,IDPTH
        NACUM(I) = 0
   10 CONTINUE
      DO 30 I=1,N
        LVL(I) = 1
        LVLS2(I) = IDPTH + 1 - LVLS2(I)
        ITEMP = LVLS2(I)
        IF (ITEMP.GT.IDPTH) GO TO 30
        IF (ITEMP.NE.LVLS1(I)) GO TO 20
        NACUM(ITEMP) = NACUM(ITEMP) + 1
        GO TO 30
   20   LVL(I) = 0
   30 CONTINUE
      RETURN
      END
      INTEGER FUNCTION SORT2(DMY)                                       
C SORT2 SORTS SIZE AND STPT INTO DESCENDING ORDER ACCORDING TO
C VALUES OF SIZE. XC=NUMBER OF ENTRIES IN EACH ARRAY
      INTEGER TEMP, CCSTOR, SIZE, STPT, XC
C IT IS ASSUMED THAT THE GRAPH HAS AT MOST 50 CONNECTED COMPONENTS.
      COMMON /CC/ XC, SIZE(5000), STPT(5000)
      SORT2 = 0
      IF (XC.EQ.0) RETURN
      SORT2 = 1
      IND = XC
   10 ITEST = 0
      IND = IND - 1
      IF (IND.LT.1) RETURN
      DO 20 I=1,IND
        J = I + 1
        IF (SIZE(I).GE.SIZE(J)) GO TO 20
        ITEST = 1
        TEMP = SIZE(I)
        SIZE(I) = SIZE(J)
        SIZE(J) = TEMP
        TEMP = STPT(I)
        STPT(I) = STPT(J)
        STPT(J) = TEMP
   20 CONTINUE
      IF (ITEST.EQ.1) GO TO 10
      RETURN
      END
      SUBROUTINE PIKLVL(LVLS1, LVLS2, CCSTOR, IDFLT, ISDIR)             
C PIKLVL CHOOSES THE LEVEL STRUCTURE  USED IN NUMBERING GRAPH
C LVLS1-       ON INPUT CONTAINS FORWARD LEVELING INFO
C LVLS2-       ON INPUT CONTAINS REVERSE LEVELING INFO
C              ON OUTPUT THE FINAL LEVEL STRUCTURE CHOSEN
C CCSTOR-      ON INPUT CONTAINS CONNECTED COMPONENT INFO
C IDFLT-       ON INPUT =1 IF WDTH LVLS1.LE.WDTH LVLS2, =2 OTHERWISE
C NHIGH        KEEPS TRACK OF LEVEL WIDTHS FOR HIGH NUMBERING
C NLOW-        KEEPS TRACK OF LEVEL WIDTHS FOR LOW NUMBERING
C NACUM-       KEEPS TRACK OF LEVEL WIDTHS FOR CHOSEN LEVEL STRUCTURE
C XC-          NUMBER OF CONNECTED COMPONENTS
C SIZE(I)-     SIZE OF ITH CONNECTED COMPONENT
C STPT(I)-     INDEX INTO CCSTORE OF 1ST NODE IN ITH CON COMPT
C ISDIR-       FLAG WHICH INDICATES WHICH WAY THE LARGEST CONNECTED
C              COMPONENT FELL.  =+1 IF LOW AND -1 IF HIGH
      INTEGER CCSTOR, SIZE, STPT, XC, END
      COMMON /GRA/ N, IDPTH, IDEG
C IT IS ASSUMED THAT THE GRAPH HAS AT MOST 50 COMPONENTS AND
C THAT THERE ARE AT MOST 100 LEVELS.
      COMMON /LVLW/ NHIGH(10000), NLOW(10000), NACUM(10000)
      COMMON /CC/ XC, SIZE(5000), STPT(5000)
      DIMENSION LVLS1(1), LVLS2(1), CCSTOR(1)
C FOR EACH CONNECTED COMPONENT DO
      DO 80 I=1,XC
        J = STPT(I)
        END = SIZE(I) + J - 1
C SET NHIGH AND NLOW EQUAL TO NACUM
        DO 10 K=1,IDPTH
          NHIGH(K) = NACUM(K)
          NLOW(K) = NACUM(K)
   10   CONTINUE
C UPDATE NHIGH AND NLOW FOR EACH NODE IN CONNECTED COMPONENT
        DO 20 K=J,END
          INODE = CCSTOR(K)
          LVLNH = LVLS1(INODE)
          NHIGH(LVLNH) = NHIGH(LVLNH) + 1
          LVLNL = LVLS2(INODE)
          NLOW(LVLNL) = NLOW(LVLNL) + 1
   20   CONTINUE
        MAX1 = 0
        MAX2 = 0
C SET MAX1=LARGEST NEW NUMBER IN NHIGH
C SET MAX2=LARGEST NEW NUMBER IN NLOW
        DO 30 K=1,IDPTH
          IF (2*NACUM(K).EQ.NLOW(K)+NHIGH(K)) GO TO 30
          IF (NHIGH(K).GT.MAX1) MAX1 = NHIGH(K)
          IF (NLOW(K).GT.MAX2) MAX2 = NLOW(K)
   30   CONTINUE
C SET IT= NUMBER OF LEVEL STRUCTURE TO BE USED
        IT = 1
        IF (MAX1.GT.MAX2) IT = 2
        IF (MAX1.EQ.MAX2) IT = IDFLT
        IF (IT.EQ.2) GO TO 60
        IF (I.EQ.1) ISDIR = -1
C COPY LVLS1 INTO LVLS2 FOR EACH NODE IN CONNECTED COMPONENT
        DO 40 K=J,END
          INODE = CCSTOR(K)
          LVLS2(INODE) = LVLS1(INODE)
   40   CONTINUE
C UPDATE NACUM TO BE THE SAME AS NHIGH
        DO 50 K=1,IDPTH
          NACUM(K) = NHIGH(K)
   50   CONTINUE
        GO TO 80
C UPDATE NACUM TO BE THE SAME AS NLOW
   60   DO 70 K=1,IDPTH
          NACUM(K) = NLOW(K)
   70   CONTINUE
   80 CONTINUE
      RETURN
      END
      SUBROUTINE NUMBER(SND, NUM, NDSTK, LVLS2, NDEG, RENUM, LVLST,     
     * LSTPT, NR, NFLG, IBW2, IPF2, IPFA, ISDIR)
C  NUMBER PRODUCES THE NUMBERING OF THE GRAPH FOR MIN BANDWIDTH
C  SND-         ON INPUT THE NODE TO BEGIN NUMBERING ON
C  NUM-         ON INPUT AND OUTPUT, THE NEXT AVAILABLE NUMBER
C  LVLS2-       THE LEVEL STRUCTURE TO BE USED IN NUMBERING
C  RENUM-       THE ARRAY USED TO STORE THE NEW NUMBERING
C  LVLST-       ON OUTPUT CONTAINS LEVEL STRUCTURE
C  LSTPT(I)-    ON OUTPUT, INDEX INTO LVLST TO FIRST NODE IN ITH LVL
C               LSTPT(I+1) - LSTPT(I) = NUMBER OF NODES IN ITH LVL
C  NFLG-        =+1 IF SND IS FORWARD END OF PSEUDO-DIAM
C               =-1 IF SND IS REVERSE END OF PSEUDO-DIAM
C  IBW2-        BANDWIDTH OF NEW NUMBERING COMPUTED BY NUMBER
C  IPF2-        PROFILE OF NEW NUMBERING COMPUTED BY NUMBER
C  IPFA-        WORKING STORAGE USED TO COMPUTE PROFILE AND BANDWIDTH
C  ISDIR-       INDICATES STEP DIRECTION USED IN NUMBERING(+1 OR -1)
C USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
      INTEGER NDSTK
      INTEGER SND, STKA, STKB, STKC, STKD, XA, XB, XC, XD, CX, END,
     * RENUM, TEST
      COMMON /GRA/ N, IDPTH, IDEG
C THE STORAGE IN COMMON BLOCKS CC AND LVLW IS NOW FREE AND CAN
C BE USED FOR STACKS.
      COMMON /LVLW/ STKA(10000), STKB(10000), STKC(10000)
      COMMON /CC/ STKD(10001)
      DIMENSION IPFA(1)
      DIMENSION NDSTK(NR,1), LVLS2(1), NDEG(1), RENUM(1), LVLST(1),
     * LSTPT(1)
C SET UP LVLST AND LSTPT FROM LVLS2
      DO 10 I=1,N
        IPFA(I) = 0
   10 CONTINUE
      write(0,*) 'NUMBER: initialization on NSTPT'
      NSTPT = 1
      DO 30 I=1,IDPTH
        LSTPT(I) = NSTPT
        DO 20 J=1,N
          IF (LVLS2(J).NE.I) GO TO 20
          LVLST(NSTPT) = J
          NSTPT = NSTPT + 1
   20   CONTINUE
   30 CONTINUE
      LSTPT(IDPTH+1) = NSTPT
      write(0,*) 'NUMBER: initialization completed', idpth,nstpt
C STKA, STKB, STKC AND STKD ARE STACKS WITH POINTERS
C XA,XB,XC, AND XD.  CX IS A SPECIAL POINTER INTO STKC WHICH
C INDICATES THE PARTICULAR NODE BEING PROCESSED.
C LVLN KEEPS TRACK OF THE LEVEL WE ARE WORKING AT.
C INITIALLY STKC CONTAINS ONLY THE INITIAL NODE, SND.
      LVLN = 0
      IF (NFLG.LT.0) LVLN = IDPTH + 1
      XC = 1
      STKC(XC) = SND
   40 CX = 1
      XD = 0
      LVLN = LVLN + NFLG
      LST = LSTPT(LVLN)
      LND = LSTPT(LVLN+1) - 1
C BEGIN PROCESSING NODE STKC(CX)
   50 IPRO = STKC(CX)
      RENUM(IPRO) = NUM
      NUM = NUM + ISDIR
      END = NDEG(IPRO)
      XA = 0
      XB = 0
C CHECK ALL ADJACENT NODES
      DO 80 I=1,END
c$$$        write(0,*) 'NUMBER: loop 80 ',i,end, lvln
        TEST = NDSTK(IPRO,I)
        INX = RENUM(TEST)
C ONLY NODES NOT NUMBERED OR ALREADY ON A STACK ARE ADDED
        IF (INX.EQ.0) GO TO 60
        IF (INX.LT.0) GO TO 80
C DO PRELIMINARY BANDWIDTH AND PROFILE CALCULATIONS
        NBW = (RENUM(IPRO)-INX)*ISDIR
        IF (ISDIR.GT.0) INX = RENUM(IPRO)
        IF (IPFA(INX).LT.NBW) IPFA(INX) = NBW
        GO TO 80
   60   RENUM(TEST) = -1
C PUT NODES ON SAME LEVEL ON STKA, ALL OTHERS ON STKB
        IF (LVLS2(TEST).EQ.LVLS2(IPRO)) GO TO 70
        XB = XB + 1
        if (xb>10000) write(0,*) 'XB>10000 in NUMBER'
        STKB(XB) = TEST
        GO TO 80
   70   XA = XA + 1
        if (xa>10000) write(0,*) 'XA>10000 in NUMBER'
        STKA(XA) = TEST
   80 CONTINUE
C SORT STKA AND STKB INTO INCREASING DEGREE AND ADD STKA TO STKC
C AND STKB TO STKD
      IF (XA.EQ.0) GO TO 100
      IF (XA.EQ.1) GO TO 90
      CALL SORTDG(STKC, STKA, XC, XA, NDEG)
      GO TO 100
   90 XC = XC + 1
      if (xc>10000) write(0,*) 'XC>10000 in NUMBER'
      STKC(XC) = STKA(XA)

  100 IF (XB.EQ.0) GO TO 120
      IF (XB.EQ.1) GO TO 110
      CALL SORTDG(STKD, STKB, XD, XB, NDEG)
      GO TO 120
  110 XD = XD + 1
      if (xd>10000) write(0,*) 'XD>10000 in NUMBER'
      STKD(XD) = STKB(XB)
C BE SURE TO PROCESS ALL NODES IN STKC
  120 CX = CX + 1
      if (cx>10000) write(0,*) 'CX>10000 in NUMBER'
      IF (XC.GE.CX) GO TO 50
C WHEN STKC IS EXHAUSTED LOOK FOR MIN DEGREE NODE IN SAME LEVEL
C WHICH HAS NOT BEEN PROCESSED
      MAX = IDEG + 1
      SND = N + 1
      DO 130 I=LST,LND
        TEST = LVLST(I)
        IF (RENUM(TEST).NE.0) GO TO 130
        IF (NDEG(TEST).GE.MAX) GO TO 130
        RENUM(SND) = 0
        RENUM(TEST) = -1
        MAX = NDEG(TEST)
        SND = TEST
  130 CONTINUE
      IF (SND.EQ.N+1) GO TO 140
      XC = XC + 1
      if (xc>10000) write(0,*) 'XC>10000 ...2... in NUMBER'
      STKC(XC) = SND
      GO TO 50
C IF STKD IS EMPTY WE ARE DONE, OTHERWISE COPY STKD ONTO STKC
C AND BEGIN PROCESSING NEW STKC
  140 IF (XD.EQ.0) GO TO 160
      DO 150 I=1,XD
        STKC(I) = STKD(I)
  150 CONTINUE
      XC = XD
      GO TO 40
C DO FINAL BANDWIDTH AND PROFILE CALCULATIONS
  160 DO 170 I=1,N
        IF (IPFA(I).GT.IBW2) IBW2 = IPFA(I)
        IPF2 = IPF2 + IPFA(I)
  170 CONTINUE
      RETURN
      END
