ma27ad.f

      NVROOT = 0
C
C Link together variables having same degree
      DO 30 IS = 1,N
        K = IPE(IS)
        IF (K.GT.0) THEN
          ID = IW(K) + 1
          NS = IPD(ID)
          IF (NS.GT.0) LST(NS) = IS
          NXT(IS) = NS
          IPD(ID) = IS
          LST(IS) = -ID
        ELSE
C We have a variable that can be eliminated at once because there is
C     no off-diagonal nonzero in its row.
          NEL = NEL + 1
          FLAG(IS) = -1
          NXT(IS) = 0
          LST(IS) = 0
        ENDIF
   30 CONTINUE

C
C Start of main loop
C
      DO 340 ML = 1,N
C Leave loop if all variables have been eliminated.
        IF (NEL+NVROOT+1.GE.N) GO TO 350
C
C Find next supervariable for elimination.
        DO 40 ID = MD,N
          MS = IPD(ID)
          IF (MS.GT.0) GO TO 50
   40   CONTINUE
   50   MD = ID
C Nvpiv holds the number of variables in the pivot.
        NVPIV = NV(MS)
C
C Remove chosen variable from linked list
        NS = NXT(MS)
        NXT(MS) = 0
        LST(MS) = 0
        IF (NS.GT.0) LST(NS) = -ID
        IPD(ID) = NS
        ME = MS
        NEL = NEL + NVPIV
C IDN holds the degree of the new element.
        IDN = 0
C
C Run through the list of the pivotal supervariable, setting tree
C     pointers and constructing new list of supervariables.
C KP is a pointer to the current position in the old list.
        KP = IPE(ME)
        FLAG(MS) = -1
C IP points to the start of the new list.
        IP = IWFR
C LEN holds the length of the list associated with the pivot.
        LEN = IW(KP)
        DO 140 KP1 = 1,LEN
          KP = KP + 1
          KE = IW(KP)
C Jump if KE is an element that has not been merged into another.
          IF (FLAG(KE).LE.-2) GO TO 60
C Jump if KE is an element that has been merged into another or is
C     a supervariable that has been eliminated.
          IF (FLAG(KE).LE.0) THEN
             IF (IPE(KE).NE.-ROOT) GO TO 140
C KE has been merged into the root
             KE = ROOT
             IF (FLAG(KE).LE.0) GO TO 140
          END IF
C We have a supervariable. Prepare to search rest of list.
          JP = KP - 1
          LN = LEN - KP1 + 1
          IE = MS
          GO TO 70
C Search variable list of element KE, using JP as a pointer to it.
   60     IE = KE
          JP = IPE(IE)
          LN = IW(JP)
C
C Search for different supervariables and add them to the new list,
C     compressing when necessary. This loop is executed once for
C     each element in the list and once for all the supervariables
C     in the list.
   70     DO 130 JP1 = 1,LN
            JP = JP + 1
            IS = IW(JP)
C Jump if IS is not a principal variable or has already been counted.
            IF (FLAG(IS).LE.0) THEN
               IF (IPE(IS).EQ.-ROOT) THEN
C IS has been merged into the root
                  IS = ROOT
                  IW(JP) = ROOT
                  IF (FLAG(IS).LE.0) GO TO 130
               ELSE
                  GO TO 130
               END IF
            END IF
            FLAG(IS) = 0

C To fix Nick bug need to add one here to store (eventually) length
C     of new row
            IF (IWFR .GE. LW-1) THEN
C Logic was previously as below
CCC         IF (IWFR.LT.LW) GO TO 100
C Prepare for compressing IW by adjusting pointers and
C     lengths so that the lists being searched in the inner and outer
C     loops contain only the remaining entries.
              IPE(MS) = KP
              IW(KP) = LEN - KP1
              IPE(IE) = JP
              IW(JP) = LN - JP1
C Compress IW
              CALL MA27UD(N,IPE,IW,IP-1,LWFR,NCMPA)
C Copy new list forward
              JP2 = IWFR - 1
              IWFR = LWFR
              IF (IP.GT.JP2) GO TO 90
              DO 80 JP = IP,JP2
                IW(IWFR) = IW(JP)
                IWFR = IWFR + 1
   80         CONTINUE
C Adjust pointers for the new list and the lists being searched.
   90         IP = LWFR
              JP = IPE(IE)
              KP = IPE(ME)
            ENDIF

C Store IS in new list.
            IW(IWFR) = IS
            IDN = IDN + NV(IS)
            IWFR = IWFR + 1
C Remove IS from degree linked list
            LS = LST(IS)
            LST(IS) = 0
            NS = NXT(IS)
            NXT(IS) = 0
            IF (NS.GT.0) LST(NS) = LS
            IF (LS.LT.0) THEN
              LS = -LS
              IPD(LS) = NS
            ELSE IF (LS.GT.0) THEN
              NXT(LS) = NS
            END IF
  130     CONTINUE
C Jump if we have just been searching the variables at the end of
C     the list of the pivot.
          IF (IE.EQ.MS) GO TO 150
C Set tree pointer and flag to indicate element IE is absorbed into
C     new element ME.
          IPE(IE) = -ME
          FLAG(IE) = -1
  140   CONTINUE

C Store the degree of the pivot.
  150   NV(MS) = IDN + NVPIV

C Jump if new element is null.
        IF (IWFR.EQ.IP) THEN
          IPE(ME) = 0
          GO TO 340
        ENDIF

        K1 = IP
        K2 = IWFR - 1
C
C Run through new list of supervariables revising each associated list,
C     recalculating degrees and removing duplicates.
        LIMIT = NINT(FRATIO*(N-NEL))
        DO 310 K = K1,K2
          IS = IW(K)
          IF (IS.EQ.ROOT) GO TO 310
          IF (NFLG.GT.2) GO TO 170
C Reset FLAG values to +/-IOVFLO.
          DO 160 I = 1,N
            IF (FLAG(I).GT.0) FLAG(I) = IOVFLO
            IF (FLAG(I).LE.-2) FLAG(I) = -IOVFLO
  160     CONTINUE
          NFLG = IOVFLO
C Reduce NFLG by one to cater for this supervariable.
  170     NFLG = NFLG - 1
C Begin with the degree of the new element. Its variables must always
C     be counted during the degree calculation and they are already
C     flagged with the value 0.
          ID = IDN
C Run through the list associated with supervariable IS
          KP1 = IPE(IS) + 1
C NP points to the next entry in the revised list.
          NP = KP1
          KP2 = IW(KP1-1) + KP1 - 1
          DO 220 KP = KP1,KP2
            KE = IW(KP)
C Test whether KE is an element, a redundant entry or a supervariable.
            IF (FLAG(KE).EQ.-1) THEN
              IF (IPE(KE).NE.-ROOT) GO TO 220
C KE has been merged into the root
              KE = ROOT
              IW(KP) = ROOT
              IF (FLAG(KE).EQ.-1) GO TO 220
            END IF
            IF (FLAG(KE).GE.0) GO TO 230
C Search list of element KE, revising the degree when new variables
C     found.
            JP1 = IPE(KE) + 1
            JP2 = IW(JP1-1) + JP1 - 1
            IDL = ID
            DO 190 JP = JP1,JP2
              JS = IW(JP)
C Jump if JS has already been counted.
              IF (FLAG(JS).LE.NFLG) GO TO 190
              ID = ID + NV(JS)
              FLAG(JS) = NFLG
  190       CONTINUE
C Jump if one or more new supervariables were found.
            IF (ID.GT.IDL) GO TO 210
C Check whether every variable of element KE is in new element ME.
            DO 200 JP = JP1,JP2
              JS = IW(JP)
              IF (FLAG(JS).NE.0) GO TO 210
  200       CONTINUE
C Set tree pointer and FLAG to indicate that element KE is absorbed
C     into new element ME.
            IPE(KE) = -ME
            FLAG(KE) = -1
            GO TO 220
C Store element KE in the revised list for supervariable IS and flag it.
  210       IW(NP) = KE
            FLAG(KE) = -NFLG
            NP = NP + 1
  220     CONTINUE
          NP0 = NP
          GO TO 250
C Treat the rest of the list associated with supervariable IS. It
C     consists entirely of supervariables.
  230     KP0 = KP
          NP0 = NP
          DO 240 KP = KP0,KP2
            KS = IW(KP)
            IF (FLAG(KS).LE.NFLG) THEN
               IF (IPE(KS).EQ.-ROOT) THEN
                  KS = ROOT
                  IW(KP) = ROOT
                  IF (FLAG(KS).LE.NFLG) GO TO 240
               ELSE
                  GO TO 240
               END IF
            END IF
C Add to degree, flag supervariable KS and add it to new list.
            ID = ID + NV(KS)
            FLAG(KS) = NFLG
            IW(NP) = KS
            NP = NP + 1
  240     CONTINUE
C Move first supervariable to end of list, move first element to end
C     of element part of list and add new element to front of list.
  250     IF (ID.GE.LIMIT) GO TO 295
          IW(NP) = IW(NP0)
          IW(NP0) = IW(KP1)
          IW(KP1) = ME
C Store the new length of the list.
          IW(KP1-1) = NP - KP1 + 1
C
C Check whether row is is identical to another by looking in linked
C     list of supervariables with degree ID at those whose lists have
C     first entry ME. Note that those containing ME come first so the
C     search can be terminated when a list not starting with ME is
C     found.
          JS = IPD(ID)
          DO 280 L = 1,N
            IF (JS.LE.0) GO TO 300
            KP1 = IPE(JS) + 1
            IF (IW(KP1).NE.ME) GO TO 300
C JS has same degree and is active. Check if identical to IS.
            KP2 = KP1 - 1 + IW(KP1-1)
            DO 260 KP = KP1,KP2
              IE = IW(KP)
C Jump if IE is a supervariable or an element not in the list of IS.
              IF (ABS(FLAG(IE)+0).GT.NFLG) GO TO 270
  260       CONTINUE
            GO TO 290

  270       JS = NXT(JS)
  280     CONTINUE
C Supervariable amalgamation. Row IS is identical to row JS.
C Regard all variables in the two supervariables as being in IS. Set
C     tree pointer, FLAG and NV entries.
  290     IPE(JS) = -IS
          NV(IS) = NV(IS) + NV(JS)
          NV(JS) = 0
          FLAG(JS) = -1
C Replace JS by IS in linked list.
          NS = NXT(JS)
          LS = LST(JS)
          IF (NS.GT.0) LST(NS) = IS
          IF (LS.GT.0) NXT(LS) = IS
          LST(IS) = LS
          NXT(IS) = NS
          LST(JS) = 0
          NXT(JS) = 0
          IF (IPD(ID).EQ.JS) IPD(ID) = IS
          GO TO 310
C Treat IS as full. Merge it into the root node.
  295     IF (NVROOT.EQ.0) THEN
            ROOT = IS
            IPE(IS) = 0
          ELSE
            IW(K) = ROOT
            IPE(IS) = -ROOT
            NV(ROOT) = NV(ROOT) + NV(IS)
            NV(IS) = 0
            FLAG(IS) = -1
          END IF
          NVROOT = NV(ROOT)
          GO TO 310
C Insert IS into linked list of supervariables of same degree.
  300     NS = IPD(ID)
          IF (NS.GT.0) LST(NS) = IS
          NXT(IS) = NS
          IPD(ID) = IS
          LST(IS) = -ID
          MD = MIN(MD,ID)
  310   CONTINUE

C
C Reset flags for supervariables in newly created element and
C     remove those absorbed into others.
        DO 320 K = K1,K2
          IS = IW(K)
          IF (NV(IS).EQ.0) GO TO 320
          FLAG(IS) = NFLG
          IW(IP) = IS
          IP = IP + 1
  320   CONTINUE

        FLAG(ME) = -NFLG
C Move first entry to end to make room for length.
        IW(IP) = IW(K1)
        IW(K1) = IP - K1
C Set pointer for new element and reset IWFR.
        IPE(ME) = K1
        IWFR = IP + 1

C  End of main loop
  340 CONTINUE
C

C Absorb any remaining variables into the root
  350 DO 360 IS = 1,N
        IF(NXT(IS).NE.0 .OR. LST(IS).NE.0) THEN
          IF (NVROOT.EQ.0) THEN
            ROOT = IS
            IPE(IS) = 0
          ELSE
            IPE(IS) = -ROOT
          END IF
          NVROOT = NVROOT + NV(IS)
          NV(IS) = 0
         END IF
  360 CONTINUE
C Link any remaining elements to the root
      DO 370 IE = 1,N
        IF (IPE(IE).GT.0) IPE(IE) = -ROOT
  370 CONTINUE
      IF(NVROOT.GT.0)NV(ROOT)=NVROOT
      END
      SUBROUTINE MA27UD(N,IPE,IW,LW,IWFR,NCMPA)
C COMPRESS LISTS HELD BY MA27H/HD AND MA27K/KD IN IW AND ADJUST POINTERS
C     IN IPE TO CORRESPOND.
C N IS THE MATRIX ORDER. IT IS NOT ALTERED.
C IPE(I) POINTS TO THE POSITION IN IW OF THE START OF LIST I OR IS
C     ZERO IF THERE IS NO LIST I. ON EXIT IT POINTS TO THE NEW POSITION.
C IW HOLDS THE LISTS, EACH HEADED BY ITS LENGTH. ON OUTPUT THE SAME
C     LISTS ARE HELD, BUT THEY ARE NOW COMPRESSED TOGETHER.
C LW HOLDS THE LENGTH OF IW. IT IS NOT ALTERED.
C IWFR NEED NOT BE SET ON ENTRY. ON EXIT IT POINTS TO THE FIRST FREE
C     LOCATION IN IW.
C     ON RETURN IT IS SET TO THE FIRST FREE LOCATION IN IW.
C NCMPA see INFO(11) in MA27A/AD.
C
C     .. Scalar Arguments ..
      INTEGER IWFR,LW,N,NCMPA
C     ..
C     .. Array Arguments ..
      INTEGER IPE(N),IW(LW)
C     ..
C     .. Local Scalars ..
      INTEGER I,IR,K,K1,K2,LWFR
C     ..
C     .. Executable Statements ..
      NCMPA = NCMPA + 1
C PREPARE FOR COMPRESSING BY STORING THE LENGTHS OF THE
C     LISTS IN IPE AND SETTING THE FIRST ENTRY OF EACH LIST TO
C     -(LIST NUMBER).
      DO 10 I = 1,N
        K1 = IPE(I)
        IF (K1.LE.0) GO TO 10
        IPE(I) = IW(K1)
        IW(K1) = -I
   10 CONTINUE
C
C COMPRESS
C IWFR POINTS JUST BEYOND THE END OF THE COMPRESSED FILE.
C LWFR POINTS JUST BEYOND THE END OF THE UNCOMPRESSED FILE.
      IWFR = 1
      LWFR = IWFR
      DO 60 IR = 1,N
        IF (LWFR.GT.LW) GO TO 70
C SEARCH FOR THE NEXT NEGATIVE ENTRY.
        DO 20 K = LWFR,LW
          IF (IW(K).LT.0) GO TO 30
   20   CONTINUE
        GO TO 70
C PICK UP ENTRY NUMBER, STORE LENGTH IN NEW POSITION, SET NEW POINTER
C     AND PREPARE TO COPY LIST.
   30   I = -IW(K)
        IW(IWFR) = IPE(I)
        IPE(I) = IWFR
        K1 = K + 1
        K2 = K + IW(IWFR)
        IWFR = IWFR + 1
        IF (K1.GT.K2) GO TO 50
C COPY LIST TO NEW POSITION.
        DO 40 K = K1,K2
          IW(IWFR) = IW(K)
          IWFR = IWFR + 1
   40   CONTINUE
   50   LWFR = K2 + 1
   60 CONTINUE
   70 RETURN

      END
      SUBROUTINE MA27JD(N,NZ,IRN,ICN,PERM,IW,LW,IPE,IQ,FLAG,IWFR,
     +                 ICNTL,INFO)
C
C SORT PRIOR TO CALLING ANALYSIS ROUTINE MA27K/KD.
C
C GIVEN THE POSITIONS OF THE OFF-DIAGONAL NON-ZEROS OF A SYMMETRIC
C     MATRIX AND A PERMUTATION, CONSTRUCT THE SPARSITY PATTERN
C     OF THE STRICTLY UPPER TRIANGULAR PART OF THE PERMUTED MATRIX.
C     EITHER ONE OF A PAIR (I,J),(J,I) MAY BE USED TO REPRESENT
C     THE PAIR. DIAGONAL ELEMENTS ARE IGNORED. NO CHECK IS MADE
C     FOR DUPLICATE ELEMENTS UNLESS ANY ROW HAS MORE THAN ICNTL(4)
C     NON-ZEROS, IN WHICH CASE DUPLICATES ARE REMOVED.
C
C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
C NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT. IT IS NOT
C     ALTERED.
C IRN(I),I=1,2,...,NZ MUST BE SET TO THE ROW INDICES OF THE
C     NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS EQUIVALENCED WITH IW.
C     IRN(1) MAY BE EQUIVALENCED WITH IW(1).
C ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN INDICES OF THE
C     NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS EQUIVALENCED
C     WITH IW.ICN(1) MAY BE EQUIVELENCED WITH IW(K),K.GT.NZ.
C PERM(I) MUST BE SET TO HOLD THE POSITION OF VARIABLE I IN THE
C     PERMUTED ORDER. IT IS NOT ALTERED.
C IW NEED NOT BE SET ON INPUT. ON OUTPUT IT CONTAINS LISTS OF
C     COLUMN NUMBERS, EACH LIST BEING HEADED BY ITS LENGTH.
C LW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST
C     MAX(NZ,N+(NO. OF OFF-DIAGONAL NON-ZEROS)). IT IS NOT ALTERED.
C IPE NEED NOT BE SET ON INPUT. ON OUTPUT IPE(I) POINTS TO THE START OF
C     THE ENTRY IN IW FOR ROW I, OR IS ZERO IF THERE IS NO ENTRY.
C IQ NEED NOT BE SET. ON OUTPUT IQ(I) CONTAINS THE NUMBER OF
C     OFF-DIAGONAL NON-ZEROS IN ROW I, INCLUDING DUPLICATES.
C FLAG IS USED FOR WORKSPACE TO HOLD FLAGS TO PERMIT DUPLICATE
C     ENTRIES TO BE IDENTIFIED QUICKLY.
C IWFR NEED NOT BE SET ON INPUT. ON OUTPUT IT POINTS TO THE FIRST
C     UNUSED LOCATION IN IW.
C ICNTL is an INTEGER array of length 30, see MA27A/AD.
C INFO is an INTEGER array of length 20, see MA27A/AD.
C
C     .. Scalar Arguments ..
      INTEGER IWFR,LW,N,NZ
C     ..
C     .. Array Arguments ..
      INTEGER FLAG(N),ICN(*),IPE(N),IQ(N),IRN(*),IW(LW),PERM(N)
      INTEGER ICNTL(30),INFO(20)
C     ..
C     .. Local Scalars ..
      INTEGER I,ID,IN,J,JDUMMY,K,K1,K2,L,LBIG,LEN
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC MAX
C     ..
C     .. Executable Statements ..
C
C INITIALIZE INFO(1), INFO(2) AND IQ
      INFO(1) = 0
      INFO(2) = 0
      DO 10 I = 1,N
        IQ(I) = 0
   10 CONTINUE
C
C COUNT THE NUMBERS OF NON-ZEROS IN THE ROWS, PRINT WARNINGS ABOUT
C     OUT-OF-RANGE INDICES AND TRANSFER GENUINE ROW NUMBERS
C     (NEGATED) INTO IW.
      IF (NZ.EQ.0) GO TO 110
      DO 100 K = 1,NZ
        I = IRN(K)
        J = ICN(K)
        IW(K) = -I
        IF(I.LT.J) THEN
          IF (I.GE.1 .AND. J.LE.N) GO TO 80
        ELSE IF(I.GT.J) THEN
          IF (J.GE.1 .AND. I.LE.N) GO TO 80
        ELSE
          IW(K) = 0
          IF (I.GE.1 .AND. I.LE.N) GO TO 100
        END IF
        INFO(2) = INFO(2) + 1
        INFO(1) = 1
        IW(K) = 0
        IF (INFO(2).LE.1 .AND. ICNTL(2).GT.0) THEN
          WRITE (ICNTL(2),FMT=60) INFO(1)
        END IF

   60   FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MA27AD',
     +          '  *** INFO(1) =',I2)

        IF (INFO(2).LE.10 .AND. ICNTL(2).GT.0) THEN
          WRITE (ICNTL(2),FMT=70) K,I,J
        END IF

   70   FORMAT (I6,'TH NON-ZERO (IN ROW',I6,' AND COLUMN',I6,
     +         ') IGNORED')

        GO TO 100

   80   IF (PERM(J).GT.PERM(I)) GO TO 90
        IQ(J) = IQ(J) + 1
        GO TO 100

   90   IQ(I) = IQ(I) + 1
  100 CONTINUE
C
C ACCUMULATE ROW COUNTS TO GET POINTERS TO ROW ENDS
C     IN IPE.
  110 IWFR = 1
      LBIG = 0
      DO 120 I = 1,N
        L = IQ(I)
        LBIG = MAX(L,LBIG)
        IWFR = IWFR + L
        IPE(I) = IWFR - 1
  120 CONTINUE
C
C PERFORM IN-PLACE SORT
      IF (NZ.EQ.0) GO TO 250
      DO 160 K = 1,NZ
        I = -IW(K)
        IF (I.LE.0) GO TO 160
        L = K
        IW(K) = 0
        DO 150 ID = 1,NZ
          J = ICN(L)
          IF (PERM(I).LT.PERM(J)) GO TO 130
          L = IPE(J)
          IPE(J) = L - 1
          IN = IW(L)
          IW(L) = I
          GO TO 140

  130     L = IPE(I)
          IPE(I) = L - 1
          IN = IW(L)
          IW(L) = J
  140     I = -IN
          IF (I.LE.0) GO TO 160
  150   CONTINUE
  160 CONTINUE
C
C MAKE ROOM IN IW FOR ROW LENGTHS AND INITIALIZE FLAG.
      K = IWFR - 1
      L = K + N
      IWFR = L + 1
      DO 190 I = 1,N
        FLAG(I) = 0
        J = N + 1 - I
        LEN = IQ(J)
        IF (LEN.LE.0) GO TO 180
        DO 170 JDUMMY = 1,LEN
          IW(L) = IW(K)
          K = K - 1
          L = L - 1
  170   CONTINUE
  180   IPE(J) = L
        L = L - 1
  190 CONTINUE
      IF (LBIG.GE.ICNTL(4)) GO TO 210
C
C PLACE ROW LENGTHS IN IW
      DO 200 I = 1,N
        K = IPE(I)
        IW(K) = IQ(I)
        IF (IQ(I).EQ.0) IPE(I) = 0
  200 CONTINUE
      GO TO 250
C
C
C REMOVE DUPLICATE ENTRIES
  210 IWFR = 1
      DO 240 I = 1,N
        K1 = IPE(I) + 1
        K2 = IPE(I) + IQ(I)
        IF (K1.LE.K2) GO TO 220
        IPE(I) = 0
        GO TO 240

  220   IPE(I) = IWFR
        IWFR = IWFR + 1
        DO 230 K = K1,K2
          J = IW(K)
          IF (FLAG(J).EQ.I) GO TO 230
          IW(IWFR) = J
          IWFR = IWFR + 1
          FLAG(J) = I
  230   CONTINUE
        K = IPE(I)
        IW(K) = IWFR - K - 1
  240 CONTINUE
  250 RETURN

      END
      SUBROUTINE MA27KD(N,IPE,IW,LW,IWFR,IPS,IPV,NV,FLAG,NCMPA)
C
C USING A GIVEN PIVOTAL SEQUENCE AND A REPRESENTATION OF THE MATRIX THAT
C     INCLUDES ONLY NON-ZEROS OF THE STRICTLY UPPER-TRIANGULAR PART
C     OF THE PERMUTED MATRIX, CONSTRUCT TREE POINTERS.
C
C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
C IPE(I) MUST BE SET TO POINT TO THE POSITION IN IW OF THE
C     START OF ROW I OR HAVE THE VALUE ZERO IF ROW I HAS NO OFF-
C     DIAGONAL NON-ZEROS. DURING EXECUTION IT IS USED AS FOLLOWS.
C     IF VARIABLE I IS ELIMINATED THEN IPE(I) POINTS TO THE LIST
C     OF VARIABLES FOR CREATED ELEMENT I. IF ELEMENT I IS
C     ABSORBED INTO NEWLY CREATED ELEMENT J THEN IPE(I)=-J.
C IW MUST BE SET ON ENTRY TO HOLD LISTS OF VARIABLES BY
C     ROWS, EACH LIST BEING HEADED BY ITS LENGTH. WHEN A VARIABLE
C     IS ELIMINATED ITS LIST IS REPLACED BY A LIST OF VARIABLES
C     IN THE NEW ELEMENT.
C LW MUST BE SET TO THE LENGTH OF IW. IT IS NOT ALTERED.
C IWFR MUST BE SET TO THE POSITION IN IW OF THE FIRST FREE VARIABLE.
C     IT IS REVISED DURING EXECUTION, CONTINUING TO HAVE THIS MEANING.
C IPS(I) MUST BE SET TO THE POSITION OF VARIABLE I IN THE REQUIRED
C     ORDERING. IT IS NOT ALTERED.
C IPV NEED NOT BE SET. IPV(K) IS SET TO HOLD THE K TH VARIABLE
C     IN PIVOT ORDER.
C NV NEED NOT BE SET. IF VARIABLE J HAS NOT BEEN ELIMINATED THEN
C     THE LAST ELEMENT WHOSE LEADING VARIABLE (VARIABLE EARLIEST
C     IN THE PIVOT SEQUENCE) IS J IS ELEMENT NV(J). IF ELEMENT J
C     EXISTS THEN THE LAST ELEMENT HAVING THE SAME LEADING
C     VARIABLE IS NV(J). IN BOTH CASES NV(J)=0 IF THERE IS NO SUCH
C     ELEMENT. IF ELEMENT J HAS BEEN MERGED INTO A LATER ELEMENT
C     THEN NV(J) IS THE DEGREE AT THE TIME OF ELIMINATION.
C FLAG IS USED AS WORKSPACE FOR VARIABLE FLAGS.
C     FLAG(JS)=ME IF JS HAS BEEN INCLUDED IN THE LIST FOR ME.
C NCMPA see INFO(11) in MA27A/AD.
C
C     .. Scalar Arguments ..
      INTEGER IWFR,LW,N,NCMPA
C     ..
C     .. Array Arguments ..
      INTEGER FLAG(N),IPE(N),IPS(N),IPV(N),IW(LW),NV(N)
C     ..
C     .. Local Scalars ..
      INTEGER I,IE,IP,J,JE,JP,JP1,JP2,JS,KDUMMY,LN,LWFR,ME,MINJS,ML,MS
C     ..
C     .. External Subroutines ..
      EXTERNAL MA27UD
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC MIN
C     ..
C     .. Executable Statements ..
C
C INITIALIZATIONS
      DO 10 I = 1,N
        FLAG(I) = 0
        NV(I) = 0
        J = IPS(I)
        IPV(J) = I
   10 CONTINUE
      NCMPA = 0
C
C START OF MAIN LOOP
C
      DO 100 ML = 1,N
C ME=MS IS THE NAME OF THE VARIABLE ELIMINATED AND
C     OF THE ELEMENT CREATED IN THE MAIN LOOP.
        MS = IPV(ML)
        ME = MS
        FLAG(MS) = ME
C
C MERGE ROW MS WITH ALL THE ELEMENTS HAVING MS AS LEADING VARIABLE.
C IP POINTS TO THE START OF THE NEW LIST.
        IP = IWFR
C MINJS IS SET TO THE POSITION IN THE ORDER OF THE LEADING VARIABLE
C     IN THE NEW LIST.
        MINJS = N
        IE = ME
        DO 70 KDUMMY = 1,N
C SEARCH VARIABLE LIST OF ELEMENT IE.
C JP POINTS TO THE CURRENT POSITION IN THE LIST BEING SEARCHED.
          JP = IPE(IE)
C LN IS THE LENGTH OF THE LIST BEING SEARCHED.
          LN = 0
          IF (JP.LE.0) GO TO 60
          LN = IW(JP)
C
C SEARCH FOR DIFFERENT VARIABLES AND ADD THEM TO LIST,
C     COMPRESSING WHEN NECESSARY
          DO 50 JP1 = 1,LN
            JP = JP + 1
C PLACE NEXT VARIABLE IN JS.
            JS = IW(JP)
C JUMP IF VARIABLE HAS ALREADY BEEN INCLUDED.
            IF (FLAG(JS).EQ.ME) GO TO 50
            FLAG(JS) = ME
            IF (IWFR.LT.LW) GO TO 40
C PREPARE FOR COMPRESSING IW BY ADJUSTING POINTER TO AND LENGTH OF
C     THE LIST FOR IE TO REFER TO THE REMAINING ENTRIES.
            IPE(IE) = JP
            IW(JP) = LN - JP1
C COMPRESS IW.
            CALL MA27UD(N,IPE,IW,IP-1,LWFR,NCMPA)
C COPY NEW LIST FORWARD
            JP2 = IWFR - 1
            IWFR = LWFR
            IF (IP.GT.JP2) GO TO 30
            DO 20 JP = IP,JP2
              IW(IWFR) = IW(JP)
              IWFR = IWFR + 1
   20       CONTINUE
   30       IP = LWFR
            JP = IPE(IE)
C ADD VARIABLE JS TO NEW LIST.
   40       IW(IWFR) = JS
            MINJS = MIN(MINJS,IPS(JS)+0)
            IWFR = IWFR + 1
   50     CONTINUE
C RECORD ABSORPTION OF ELEMENT IE INTO NEW ELEMENT.
   60     IPE(IE) = -ME
C PICK UP NEXT ELEMENT WITH LEADING VARIABLE MS.
          JE = NV(IE)
C STORE DEGREE OF IE.
          NV(IE) = LN + 1
          IE = JE
C LEAVE LOOP IF THERE ARE NO MORE ELEMENTS.
          IF (IE.EQ.0) GO TO 80
   70   CONTINUE
   80   IF (IWFR.GT.IP) GO TO 90
C DEAL WITH NULL NEW ELEMENT.
        IPE(ME) = 0
        NV(ME) = 1
        GO TO 100
C LINK NEW ELEMENT WITH OTHERS HAVING SAME LEADING VARIABLE.
   90   MINJS = IPV(MINJS)
        NV(ME) = NV(MINJS)
        NV(MINJS) = ME
C MOVE FIRST ENTRY IN NEW LIST TO END TO ALLOW ROOM FOR LENGTH AT
C     FRONT. SET POINTER TO FRONT.
        IW(IWFR) = IW(IP)
        IW(IP) = IWFR - IP
        IPE(ME) = IP
        IWFR = IWFR + 1
  100 CONTINUE
      RETURN

      END
      SUBROUTINE MA27LD(N,IPE,NV,IPS,NE,NA,ND,NSTEPS,NEMIN)
C
C TREE SEARCH
C
C GIVEN SON TO FATHER TREE POINTERS, PERFORM DEPTH-FIRST
C     SEARCH TO FIND PIVOT ORDER AND NUMBER OF ELIMINATIONS
C     AND ASSEMBLIES AT EACH STAGE.
C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
C IPE(I) MUST BE SET EQUAL TO -(FATHER OF NODE I) OR ZERO IF
C      NODE IS A ROOT. IT IS ALTERED TO POINT TO ITS NEXT
C      YOUNGER BROTHER IF IT HAS ONE, BUT OTHERWISE IS NOT
C      CHANGED.
C NV(I) MUST BE SET TO ZERO IF NO VARIABLES ARE ELIMINATED AT NODE
C      I AND TO THE DEGREE OTHERWISE. ONLY LEAF NODES CAN HAVE
C      ZERO VALUES OF NV(I). NV IS NOT ALTERED.
C IPS(I) NEED NOT BE SET. IT IS USED TEMPORARILY TO HOLD
C      -(ELDEST SON OF NODE I) IF IT HAS ONE AND 0 OTHERWISE. IT IS
C      EVENTUALLY SET TO HOLD THE POSITION OF NODE I IN THE ORDER.
C NE(IS) NEED NOT BE SET. IT IS SET TO THE NUMBER OF VARIABLES
C      ELIMINATED AT STAGE IS OF THE ELIMINATION.
C NA(IS) NEED NOT BE SET. IT IS SET TO THE NUMBER OF ELEMENTS
C      ASSEMBLED AT STAGE IS OF THE ELIMINATION.
C ND(IS) NEED NOT BE SET. IT IS SET TO THE DEGREE AT STAGE IS OF
C     THE ELIMINATION.
C NSTEPS NEED NOT BE SET. IT IS SET TO  THE NUMBER OF ELIMINATION
C      STEPS.
C NEMIN see ICNTL(5) in MA27A/AD.
C
C     .. Scalar Arguments ..
      INTEGER N,NSTEPS,NEMIN
C     ..
C     .. Array Arguments ..
      INTEGER IPE(N),IPS(N),NA(N),ND(N),NE(N),NV(N)
C     ..
C     .. Local Scalars ..
      INTEGER I,IB,IF,IL,IS,ISON,K,L,NR
C     ..
C     .. Executable Statements ..
C INITIALIZE IPS AND NE.
      DO 10 I = 1,N
        IPS(I) = 0
        NE(I) = 0
   10 CONTINUE
C
C SET IPS(I) TO -(ELDEST SON OF NODE I) AND IPE(I) TO NEXT YOUNGER
C     BROTHER OF NODE I IF IT HAS ONE.
C FIRST PASS IS FOR NODES WITHOUT ELIMINATIONS.
      DO 20 I = 1,N
        IF (NV(I).GT.0) GO TO 20
        IF = -IPE(I)
        IS = -IPS(IF)
        IF (IS.GT.0) IPE(I) = IS
        IPS(IF) = -I
   20 CONTINUE
C NR IS DECREMENTED FOR EACH ROOT NODE. THESE ARE STORED IN
C     NE(I),I=NR,N.
      NR = N + 1
C SECOND PASS TO ADD NODES WITH ELIMINATIONS.
      DO 50 I = 1,N
        IF (NV(I).LE.0) GO TO 50
C NODE IF IS THE FATHER OF NODE I.
        IF = -IPE(I)
        IF (IF.EQ.0) GO TO 40
        IS = -IPS(IF)
C JUMP IF NODE IF HAS NO SONS YET.
        IF (IS.LE.0) GO TO 30
C SET POINTER TO NEXT BROTHER
        IPE(I) = IS
C NODE I IS ELDEST SON OF NODE IF.
   30   IPS(IF) = -I
        GO TO 50
C WE HAVE A ROOT
   40   NR = NR - 1
        NE(NR) = I
   50 CONTINUE
C
C DEPTH-FIRST SEARCH.
C IL HOLDS THE CURRENT TREE LEVEL. ROOTS ARE AT LEVEL N, THEIR SONS
C     ARE AT LEVEL N-1, ETC.
C IS HOLDS THE CURRENT ELIMINATION STAGE. WE ACCUMULATE THE NUMBER
C     OF ELIMINATIONS AT STAGE IS DIRECTLY IN NE(IS). THE NUMBER OF
C     ASSEMBLIES IS ACCUMULATED TEMPORARILY IN NA(IL), FOR TREE
C     LEVEL IL, AND IS TRANSFERED TO NA(IS) WHEN WE REACH THE
C     APPROPRIATE STAGE IS.
      IS = 1
C I IS THE CURRENT NODE.
      I = 0
      DO 160 K = 1,N
        IF (I.GT.0) GO TO 60
C PICK UP NEXT ROOT.
        I = NE(NR)
        NE(NR) = 0
        NR = NR + 1
        IL = N
        NA(N) = 0
C GO TO SON FOR AS LONG AS POSSIBLE, CLEARING FATHER-SON POINTERS
C     IN IPS AS EACH IS USED AND SETTING NA(IL)=0 FOR ALL LEVELS
C     REACHED.
   60   DO 70 L = 1,N
          IF (IPS(I).GE.0) GO TO 80
          ISON = -IPS(I)
          IPS(I) = 0
          I = ISON
          IL = IL - 1
          NA(IL) = 0
   70   CONTINUE
C RECORD POSITION OF NODE I IN THE ORDER.
   80   IPS(I) = K
        NE(IS) = NE(IS) + 1
C JUMP IF NODE HAS NO ELIMINATIONS.
        IF (NV(I).LE.0) GO TO 120
        IF (IL.LT.N) NA(IL+1) = NA(IL+1) + 1
        NA(IS) = NA(IL)
        ND(IS) = NV(I)
C CHECK FOR STATIC CONDENSATION
        IF (NA(IS).NE.1) GO TO 90
        IF (ND(IS-1)-NE(IS-1).EQ.ND(IS)) GO TO 100
C CHECK FOR SMALL NUMBERS OF ELIMINATIONS IN BOTH LAST TWO STEPS.
   90   IF (NE(IS).GE.NEMIN) GO TO 110
        IF (NA(IS).EQ.0) GO TO 110
        IF (NE(IS-1).GE.NEMIN) GO TO 110
C COMBINE THE LAST TWO STEPS
  100   NA(IS-1) = NA(IS-1) + NA(IS) - 1
        ND(IS-1) = ND(IS) + NE(IS-1)
        NE(IS-1) = NE(IS) + NE(IS-1)
        NE(IS) = 0
        GO TO 120

  110   IS = IS + 1
  120   IB = IPE(I)
        IF (IB.GE.0) THEN
C NODE I HAS A BROTHER OR IS A ROOT
          IF (IB.GT.0) NA(IL) = 0
          I = IB
        ELSE
C GO TO FATHER OF NODE I
          I = -IB
          IL = IL + 1
        END IF
  160 CONTINUE
      NSTEPS = IS - 1
      RETURN

      END
      SUBROUTINE MA27MD(N,NZ,IRN,ICN,PERM,NA,NE,ND,NSTEPS,LSTKI,LSTKR,
     +                 IW,INFO,OPS)
C
C STORAGE AND OPERATION COUNT EVALUATION.
C
C EVALUATE NUMBER OF OPERATIONS AND SPACE REQUIRED BY FACTORIZATION
C     USING MA27B/BD.  THE VALUES GIVEN ARE EXACT ONLY IF NO NUMERICAL
C     PIVOTING IS PERFORMED AND THEN ONLY IF IRN(1) WAS NOT
C     EQUIVALENCED TO IW(1) BY THE USER BEFORE CALLING MA27A/AD.  IF
C     THE EQUIVALENCE HAS BEEN MADE ONLY AN UPPER BOUND FOR NIRNEC
C     AND NRLNEC CAN BE CALCULATED ALTHOUGH THE OTHER COUNTS WILL
C     STILL BE EXACT.
C
C N MUST BE SET TO THE MATRIX ORDER. IT IS NOT ALTERED.
C NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT.  IT IS NOT ALTERED.
C IRN,ICN.  UNLESS IRN(1) HAS BEEN EQUIVALENCED TO IW(1)
C     IRN,ICN MUST BE SET TO THE ROW AND COLUMN INDICES OF THE
C     NON-ZEROS ON INPUT.  THEY ARE NOT ALTERED BY MA27M/MD.
C PERM MUST BE SET TO THE POSITION IN THE PIVOT ORDER OF EACH ROW.
C     IT IS NOT ALTERED.
C NA,NE,ND MUST BE SET TO HOLD, FOR EACH TREE NODE, THE NUMBER OF STACK
C     ELEMENTS ASSEMBLED, THE NUMBER OF ELIMINATIONS AND THE SIZE OF
C     THE ASSEMBLED FRONT MATRIX RESPECTIVELY.  THEY ARE NOT ALTERED.
C NSTEPS MUST BE SET TO HOLD THE NUMBER OF TREE NODES. IT IS NOT
C     ALTERED.
C LSTKI IS USED AS A WORK ARRAY BY MA27M/MD.
C LSTKR.  IF IRN(1) IS EQUIVALENCED TO IW(1)  THEN LSTKR(I)
C     MUST HOLD THE TOTAL NUMBER OF OFF-DIAGONAL ENTRIES (INCLUDING
C     DUPLICATES) IN ROW I (I=1,..,N) OF THE ORIGINAL MATRIX.  IT
C     IS USED AS WORKSPACE BY MA27M/MD.
C IW IS A WORKSPACE ARRAY USED BY OTHER SUBROUTINES AND PASSED TO THIS
C     SUBROUTINE ONLY SO THAT A TEST FOR EQUIVALENCE WITH IRN CAN BE
C     MADE.
C
C COUNTS FOR OPERATIONS AND STORAGE ARE ACCUMULATED IN VARIABLES
C     OPS,NRLTOT,NIRTOT,NRLNEC,NIRNEC,NRLADU,NRLNEC,NIRNEC.
C OPS NUMBER OF MULTIPLICATIONS AND ADDITIONS DURING FACTORIZATION.
C NRLADU,NIRADU REAL AND INTEGER STORAGE RESPECTIVELY FOR THE
C     MATRIX FACTORS.
C NRLTOT,NIRTOT REAL AND INTEGER STRORAGE RESPECTIVELY REQUIRED
C     FOR THE FACTORIZATION IF NO COMPRESSES ARE ALLOWED.
C NRLNEC,NIRNEC REAL AND INTEGER STORAGE RESPECTIVELY REQUIRED FOR
C     THE FACTORIZATION IF COMPRESSES ARE ALLOWED.
C INFO is an INTEGER array of length 20, see MA27A/AD.
C OPS ACCUMULATES THE NO. OF MULTIPLY/ADD PAIRS NEEDED TO CREATE THE
C     TRIANGULAR FACTORIZATION, IN THE DEFINITE CASE.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION OPS
      INTEGER N,NSTEPS,NZ
C     ..
C     .. Array Arguments ..
      INTEGER ICN(*),IRN(*),IW(*),LSTKI(N),LSTKR(N),NA(NSTEPS),
     +        ND(NSTEPS),NE(NSTEPS),PERM(N),INFO(20)
C     ..
C     .. Local Scalars ..
      INTEGER I,INEW,IOLD,IORG,IROW,ISTKI,ISTKR,ITOP,ITREE,JOLD,JORG,K,
     +        LSTK,NASSR,NELIM,NFR,NSTK,NUMORG,NZ1,NZ2
      DOUBLE PRECISION DELIM
      INTEGER NRLADU,NIRADU,NIRTOT,NRLTOT,NIRNEC,NRLNEC
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC MAX,MIN
C     ..
C     .. Executable Statements ..
C
      IF (NZ.EQ.0) GO TO 20
C JUMP IF IW AND IRN HAVE NOT BEEN EQUIVALENCED.
      IF (IRN(1).NE.IW(1)) GO TO 20
C RESET IRN(1).
      IRN(1) = IW(1) - 1
C THE TOTAL NUMBER OF OFF-DIAGONAL ENTRIES IS ACCUMULATED IN NZ2.
C LSTKI IS SET TO HOLD THE TOTAL NUMBER OF ENTRIES (INCUDING
C     THE DIAGONAL) IN EACH ROW IN PERMUTED ORDER.
      NZ2 = 0
      DO 10 IOLD = 1,N
        INEW = PERM(IOLD)
        LSTKI(INEW) = LSTKR(IOLD) + 1
        NZ2 = NZ2 + LSTKR(IOLD)
   10 CONTINUE
C NZ1 IS THE NUMBER OF ENTRIES IN ONE TRIANGLE INCLUDING THE DIAGONAL.
C NZ2 IS THE TOTAL NUMBER OF ENTRIES INCLUDING THE DIAGONAL.
      NZ1 = NZ2/2 + N
      NZ2 = NZ2 + N
      GO TO 60
C COUNT (IN LSTKI) NON-ZEROS IN ORIGINAL MATRIX BY PERMUTED ROW (UPPER
C     TRIANGLE ONLY). INITIALIZE COUNTS.