ma27ad.f

* COPYRIGHT (c) 1982 AEA Technology
*######DATE 20 September 2001
C  September 2001: threadsafe version of MA27
C  19/3/03. Array ICNTL in MA27GD made assumed size. 
C  17/9/09. Bug corrected in MA27HD.  Also some tidying of code in
C     MA27HD and change of most comments to lower case.
C  16/6/10. Statement function in MA27OD replaced by in-line code.

      SUBROUTINE MA27ID(ICNTL,CNTL)
      INTEGER ICNTL(30)
      DOUBLE PRECISION CNTL(5)

      INTEGER IFRLVL
      PARAMETER ( IFRLVL=5 )

C Stream number for error messages
      ICNTL(1) = 6
C Stream number for diagnostic messages
      ICNTL(2) = 6
C Control the level of diagnostic printing.
C   0 no printing
C   1 printing of scalar parameters and first parts of arrays.
C   2 printing of scalar parameters and whole of arrays.
      ICNTL(3) = 0
C The largest integer such that all integers I in the range
C -ICNTL(4).LE.I.LE.ICNTL(4) can be handled by the shortest integer
C type in use.
      ICNTL(4) = 2139062143
C Minimum number of eliminations in a step that is automatically
C accepted. if two adjacent steps can be combined and each has less
C eliminations then they are combined.
      ICNTL(5) = 1
C Control whether direct or indirect access is used by MA27C/CD.
C Indirect access is employed in forward and back substitution 
C respectively if the size of a block is less than
C ICNTL(IFRLVL+MIN(10,NPIV)) and ICNTL(IFRLVL+10+MIN(10,NPIV))
C respectively, where NPIV is the number of pivots in the block.
      ICNTL(IFRLVL+1)  = 32639
      ICNTL(IFRLVL+2)  = 32639
      ICNTL(IFRLVL+3)  = 32639
      ICNTL(IFRLVL+4)  = 32639
      ICNTL(IFRLVL+5)  = 14
      ICNTL(IFRLVL+6)  = 9
      ICNTL(IFRLVL+7)  = 8
      ICNTL(IFRLVL+8)  = 8
      ICNTL(IFRLVL+9)  = 9
      ICNTL(IFRLVL+10) = 10
      ICNTL(IFRLVL+11) = 32639
      ICNTL(IFRLVL+12) = 32639
      ICNTL(IFRLVL+13) = 32639
      ICNTL(IFRLVL+14) = 32689
      ICNTL(IFRLVL+15) = 24
      ICNTL(IFRLVL+16) = 11
      ICNTL(IFRLVL+17) = 9
      ICNTL(IFRLVL+18) = 8
      ICNTL(IFRLVL+19) = 9
      ICNTL(IFRLVL+20) = 10
C Not used
      ICNTL(26) = 0
      ICNTL(27) = 0
      ICNTL(28) = 0
      ICNTL(29) = 0
      ICNTL(30) = 0

C Control threshold pivoting.
      CNTL(1) = 0.1D0
C If a column of the reduced matrix has relative density greater than
C CNTL(2), it is forced into the root. All such columns are taken to
C have sparsity pattern equal to their merged patterns, so the fill
C and operation counts may be overestimated.
      CNTL(2) = 1.0D0
C An entry with absolute value less than CNTL(3) is never accepted as
C a 1x1 pivot or as the off-diagonal of a 2x2 pivot.
      CNTL(3) = 0.0D0
C Not used
      CNTL(4) = 0.0
      CNTL(5) = 0.0

      RETURN
      END

      SUBROUTINE MA27AD(N,NZ,IRN,ICN,IW,LIW,IKEEP,IW1,NSTEPS,IFLAG,
     +                 ICNTL,CNTL,INFO,OPS)
C THIS SUBROUTINE COMPUTES A MINIMUM DEGREE ORDERING OR ACCEPTS A GIVEN
C     ORDERING. IT COMPUTES ASSOCIATED ASSEMBLY AND ELIMINATION
C     INFORMATION FOR MA27B/BD.
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 NUMBERS OF THE
C     NON-ZEROS. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED
C     TO IW (SEE DESCRIPTION OF IW).
C ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN NUMBERS OF THE
C     NON-ZEROS. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED
C     TO IW (SEE DESCRIPTION OF IW).
C IW NEED NOT BE SET ON INPUT. IT IS USED FOR WORKSPACE.
C     IRN(1) MAY BE EQUIVALENCED TO IW(1) AND ICN(1) MAY BE
C     EQUIVALENCED TO IW(K), WHERE K.GT.NZ.
C LIW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST 2*NZ+3*N
C     FOR THE IFLAG=0 ENTRY AND AT LEAST NZ+3*N FOR THE IFLAG=1
C     ENTRY. IT IS NOT ALTERED.
C IKEEP NEED NOT BE SET UNLESS AN ORDERING IS GIVEN, IN WHICH CASE
C     IKEEP(I,1) MUST BE SET TO THE POSITION OF VARIABLE I IN THE
C     ORDER. ON OUTPUT IKEEP CONTAINS INFORMATION NEEDED BY MA27B/BD.
C     IKEEP(I,1) HOLDS THE POSITION OF VARIABLE I IN THE PIVOT ORDER.
C     IKEEP(I,2), IKEEP(I,3) HOLD THE NUMBER OF ELIMINATIONS, ASSEMBLIES
C     AT MAJOR STEP I, I=1,2,...,NSTEPS. NOTE THAT WHEN AN ORDER IS
C     GIVEN IT MAY BE REPLACED BY ANOTHER ORDER THAT GIVES IDENTICAL
C     NUMERICAL RESULTS.
C IW1 IS USED FOR WORKSPACE.
C NSTEPS NEED NOT BE SET. ON OUTPUT IT CONTAINS THE NUMBER OF MAJOR
C     STEPS NEEDED FOR A DEFINITE MATRIX AND MUST BE PASSED UNCHANGED
C     TO MA27B/BD.
C IFLAG MUST SET TO ZERO IF THE USER WANTS THE PIVOT ORDER CHOSEN
C     AUTOMATICALLY AND TO ONE IF HE WANTS TO SPECIFY IT IN IKEEP.
C ICNTL is an INTEGER array of length 30 containing user options
C     with integer values (defaults set in MA27I/ID)
C   ICNTL(1) (LP) MUST BE SET TO THE STREAM NUMBER FOR ERROR MESSAGES.
C     ERROR MESSAGES ARE SUPPRESSED IF ICNTL(1) IS NOT POSITIVE.
C     IT IS NOT ALTERED.
C   ICNTL(2) (MP) MUST BE SET TO THE STREAM NUMBER FOR DIAGNOSTIC
C     MESSAGES.  DIAGNOSTIC MESSAGES ARE SUPPRESSED IF ICNTL(2) IS NOT
C     POSITIVE.  IT IS NOT ALTERED.
C   ICNTL(3) (LDIAG) CONTROLS THE LEVEL OF DIAGNOSTIC PRINTING.
C     0 NO PRINTING
C     1 PRINTING OF SCALAR PARAMETERS AND FIRST PARTS OF ARRAYS.
C     2 PRINTING OF SCALAR PARAMETERS AND WHOLE OF ARRAYS.
C   ICNTL(4) (IOVFLO) IS THE LARGEST INTEGER SUCH THAT ALL INTEGERS
C     I IN THE RANGE -IOVFLO.LE.I.LE.IOVFLO CAN BE HANDLED BY THE
C     SHORTEST INTEGER TYPE IN USE.
C   ICNT(5) (NEMIN) MUST BE SET TO THE MINIMUM NUMBER OF ELIMINATIONS
C     IN A STEP THAT IS AUTOMATICALLY ACCEPTED. IF TWO ADJACENT STEPS
C     CAN BE COMBINED AND EACH HAS LESS ELIMINATIONS THEN THEY ARE
C     COMBINED.
C   ICNTL(IFRLVL+I) I=1,20, (IFRLVL) MUST BE SET TO CONTROL WHETHER
C     DIRECT OR INDIRECT ACCESS IS USED BY MA27C/CD.  INDIRECT ACCESS
C     IS EMPLOYED IN FORWARD AND BACK SUBSTITUTION RESPECTIVELY IF THE
C     SIZE OF A BLOCK IS LESS THAN ICNTL(IFRLVL+(MIN(10,NPIV)) AND
C     ICNTL(IFRLVL+10+MIN(10,NPIV)) RESPECTIVELY, WHERE NPIV IS THE
C     NUMBER OF PIVOTS IN THE BLOCK.
C   ICNTL(I) I=26,30 are not used.
C CNTL is an DOUBLE PRECISION array of length 5 containing user options
C     with real values (defaults set in MA27I/ID)
C   CNTL(1) (U) IS USED TO CONTROL THRESHOLD PIVOTING. IT IS NOT
C     ALTERED.
C   CNTL(2) (FRATIO) has default value 1.0.  If a column of the
C      reduced matrix has relative density greater than CNTL(2), it
C      is forced into the root. All such columns are taken to have
C      sparsity pattern equal to their merged patterns, so the fill
C      and operation counts may be overestimated.
C   CNTL(3) (PIVTOL) has default value 0.0. An entry with absolute
C      value less than CNTL(3) is never accepted as a 1x1 pivot or
C      as the off-diagonal of a 2x2 pivot.
C   CNTL(I) I=4,5 are not used.
C INFO is an INTEGER array of length 20 which is used to return
C     information to the user.
C   INFO(1) (IFLAG) is an error return code, zero for success, greater
C      than zero for a warning and less than zero for errors, see
C      INFO(2).
C   INFO(2) (IERROR) HOLDS ADDITIONAL INFORMATION IN THE EVENT OF ERRORS.
C     IF INFO(1)=-3 INFO(2) HOLDS A LENGTH THAT MAY SUFFICE FOR IW.
C     IF INFO(1)=-4 INFO(2) HOLDS A LENGTH THAT MAY SUFFICE FOR A.
C     IF INFO(1)=-5 INFO(2) IS SET TO THE PIVOT STEP AT WHICH SINGULARITY
C                 WAS DETECTED.
C     IF INFO(1)=-6 INFO(2) IS SET TO THE PIVOT STEP AT WHICH A CHANGE OF
C                 PIVOT SIGN WAS FOUND.
C     IF INFO(1)= 1 INFO(2) HOLDS THE NUMBER OF FAULTY ENTRIES.
C     IF INFO(1)= 2 INFO(2) IS SET TO THE NUMBER OF SIGNS CHANGES IN
C                 THE PIVOTS.
C     IF INFO(1)=3 INFO(2) IS SET TO THE RANK OF THE MATRIX.
C   INFO(3) and INFO(4) (NRLTOT and NIRTOT) REAL AND INTEGER STRORAGE
C     RESPECTIVELY REQUIRED FOR THE FACTORIZATION IF NO COMPRESSES ARE
C     ALLOWED.
C   INFO(5) and INFO(6) (NRLNEC and NIRNEC) REAL AND INTEGER STORAGE
C     RESPECTIVELY REQUIRED FOR THE FACTORIZATION IF COMPRESSES ARE
C     ALLOWED AND THE MATRIX IS DEFINITE.
C   INFO(7) and INFO(8) (NRLADU and NIRADU) REAL AND INTEGER STORAGE
C     RESPECTIVELY FOR THE MATRIX FACTORS AS CALCULATED BY MA27A/AD
C     FOR THE DEFINITE CASE.  
C   INFO(9) and INFO(10) (NRLBDU and NIRBDU) REAL AND INTEGER STORAGE
C     RESPECTIVELY FOR THE MATRIX FACTORS AS FOUND  BY MA27B/BD.
C   INFO(11) (NCMPA) ACCUMULATES THE NUMBER OF TIMES THE ARRAY IW IS
C     COMPRESSED BY MA27A/AD.
C   INFO(12) and INFO(13) (NCMPBR and NCMPBI) ACCUMULATE THE NUMBER
C     OF COMPRESSES OF THE REALS AND INTEGERS PERFORMED BY MA27B/BD.
C   INFO(14) (NTWO) IS USED BY MA27B/BD TO RECORD THE NUMBER OF 2*2
C     PIVOTS USED.
C   INFO(15) (NEIG) IS USED BY ME27B/BD TO RECORD THE NUMBER OF
C     NEGATIVE EIGENVALUES OF A.
C   INFO(16) to INFO(20) are not used. 
C OPS ACCUMULATES THE NO. OF MULTIPLY/ADD PAIRS NEEDED TO CREATE THE
C     TRIANGULAR FACTORIZATION, IN THE DEFINITE CASE.
C
C     .. Scalar Arguments ..
      INTEGER IFLAG,LIW,N,NSTEPS,NZ
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION CNTL(5),OPS
      INTEGER ICNTL(30),INFO(20)
      INTEGER ICN(*),IKEEP(N,3),IRN(*),IW(LIW),IW1(N,2)
C     ..
C     .. Local Scalars ..
      INTEGER I,IWFR,K,L1,L2,LLIW
C     ..
C     .. External Subroutines ..
      EXTERNAL MA27GD,MA27HD,MA27JD,MA27KD,MA27LD,MA27MD,MA27UD
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC MIN
C     ..
C     .. Executable Statements ..
      DO 5 I = 1,15
        INFO(I) = 0
    5 CONTINUE

      IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0) GO TO 40
C PRINT INPUT VARIABLES.
      WRITE (ICNTL(2),FMT=10) N,NZ,LIW,IFLAG

   10 FORMAT(/,/,' ENTERING MA27AD WITH      N     NZ      LIW  IFLAG',
     +       /,21X,I7,I7,I9,I7)

      NSTEPS = 0
      K = MIN(8,NZ)
      IF (ICNTL(3).GT.1) K = NZ
      IF (K.GT.0) WRITE (ICNTL(2),FMT=20) (IRN(I),ICN(I),I=1,K)

   20 FORMAT (' MATRIX NON-ZEROS',/,4 (I9,I6),/,
     +       (I9,I6,I9,I6,I9,I6,I9,I6))

      K = MIN(10,N)
      IF (ICNTL(3).GT.1) K = N
      IF (IFLAG.EQ.1 .AND. K.GT.0) THEN
        WRITE (ICNTL(2),FMT=30) (IKEEP(I,1),I=1,K)
      END IF

   30 FORMAT (' IKEEP(.,1)=',10I6,/, (12X,10I6))

   40 IF (N.LT.1 .OR. N.GT.ICNTL(4)) GO TO 70
      IF (NZ.LT.0) GO TO 100
      LLIW = LIW - 2*N
      L1 = LLIW + 1
      L2 = L1 + N
      IF (IFLAG.EQ.1) GO TO 50
      IF (LIW.LT.2*NZ+3*N+1) GO TO 130
C SORT
      CALL MA27GD(N,NZ,IRN,ICN,IW,LLIW,IW1,IW1(1,2),IW(L1),IWFR,
     +           ICNTL,INFO)
C ANALYZE USING MINIMUM DEGREE ORDERING
      CALL MA27HD(N,IW1,IW,LLIW,IWFR,IW(L1),IW(L2),IKEEP(1,2),
     +           IKEEP(1,3),IKEEP,ICNTL(4),INFO(11),CNTL(2))
      GO TO 60
C SORT USING GIVEN ORDER
   50 IF (LIW.LT.NZ+3*N+1) GO TO 120
      CALL MA27JD(N,NZ,IRN,ICN,IKEEP,IW,LLIW,IW1,IW1(1,2),IW(L1),IWFR,
     +           ICNTL,INFO)
C ANALYZE USING GIVEN ORDER
      CALL MA27KD(N,IW1,IW,LLIW,IWFR,IKEEP,IKEEP(1,2),IW(L1),IW(L2),
     +           INFO(11))
C PERFORM DEPTH-FIRST SEARCH OF ASSEMBLY TREE
   60 CALL MA27LD(N,IW1,IW(L1),IKEEP,IKEEP(1,2),IKEEP(1,3),IW(L2),
     +           NSTEPS,ICNTL(5))
C EVALUATE STORAGE AND OPERATION COUNT REQUIRED BY MA27B/BD IN THE
C     DEFINITE CASE.
C SET IW(1) SO THAT ARRAYS IW AND IRN CAN BE TESTED FOR EQUIVALENCE
C     IN MA27M/MD.
      IF(NZ.GE.1) IW(1) = IRN(1) + 1
      CALL MA27MD(N,NZ,IRN,ICN,IKEEP,IKEEP(1,3),IKEEP(1,2),IW(L2),
     +           NSTEPS,IW1,IW1(1,2),IW,INFO,OPS)
      GO TO 160

   70 INFO(1) = -1
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)

   80 FORMAT (' **** ERROR RETURN FROM MA27AD **** INFO(1)=',I3)

      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=90) N

   90 FORMAT (' VALUE OF N OUT OF RANGE ... =',I10)

      GO TO 160

  100 INFO(1) = -2
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=110) NZ

  110 FORMAT (' VALUE OF NZ OUT OF RANGE .. =',I10)

      GO TO 160

  120 INFO(2) = NZ + 3*N + 1
      GO TO 140

  130 INFO(2) = 2*NZ + 3*N + 1
  140 INFO(1) = -3
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=150) LIW,INFO(2)

  150 FORMAT (' LIW TOO SMALL, MUST BE INCREASED FROM',I10,
     +       ' TO AT LEAST',I10)

  160 IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0 .OR. INFO(1).LT.0) GO TO 200
C PRINT PARAMETER VALUES ON EXIT.
      WRITE (ICNTL(2),FMT=170) NSTEPS,INFO(1),OPS,INFO(2),INFO(3),
     +  INFO(4),INFO(5),INFO(6),INFO(7),INFO(8),INFO(11)

  170 FORMAT (/,' LEAVING MA27AD WITH NSTEPS  INFO(1)    OPS IERROR',
     +          ' NRLTOT NIRTOT',
     +        /,20X,2I7,F7.0,3I7,
     +        /,20X,' NRLNEC NIRNEC NRLADU NIRADU  NCMPA',
     +        /,20X,6I7)

      K = MIN(9,N)
      IF (ICNTL(3).GT.1) K = N
      IF (K.GT.0) WRITE (ICNTL(2),FMT=30) (IKEEP(I,1),I=1,K)
      K = MIN(K,NSTEPS)
      IF (K.GT.0) WRITE (ICNTL(2),FMT=180) (IKEEP(I,2),I=1,K)

  180 FORMAT (' IKEEP(.,2)=',10I6,/, (12X,10I6))

      IF (K.GT.0) WRITE (ICNTL(2),FMT=190) (IKEEP(I,3),I=1,K)

  190 FORMAT (' IKEEP(.,3)=',10I6,/, (12X,10I6))

  200 CONTINUE

      RETURN
      END

      SUBROUTINE MA27BD(N,NZ,IRN,ICN,A,LA,IW,LIW,IKEEP,NSTEPS,MAXFRT,
     +                 IW1,ICNTL,CNTL,INFO)
C THIS SUBROUTINE COMPUTES THE FACTORISATION OF THE MATRIX INPUT IN
C     A,IRN,ICN USING INFORMATION (IN IKEEP) FROM MA27A/AD.
C N MUST BE SET TO THE ORDER OF THE MATRIX.  IT IS NOT ALTERED.
C NZ MUST BE SET TO THE NUMBER OF NON-ZEROS INPUT.  IT IS NOT
C     ALTERED.
C IRN,ICN,A.  ENTRY K (K=1,NZ) OF IRN,ICN MUST BE SET TO THE ROW
C     AND COLUMN INDEX RESPECTIVELY OF THE NON-ZERO IN A.
C     IRN AND ICN ARE UNALTERED BY MA27B/BD.
C     ON EXIT, ENTRIES 1 TO NRLBDU OF A HOLD REAL INFORMATION
C     ON THE FACTORS AND SHOULD BE PASSED UNCHANGED TO MA27C/CD.
C LA LENGTH OF ARRAY A.  AN INDICATION OF A SUITABLE VALUE,
C     SUFFICIENT FOR FACTORIZATION OF A DEFINITE MATRIX, WILL
C     HAVE BEEN PROVIDED IN NRLNEC AND NRLTOT BY MA27A/AD.
C     IT IS NOT ALTERED BY MA27B/BD.
C IW NEED NOT BE SET ON ENTRY.  USED AS A WORK ARRAY BY MA27B/BD.
C     ON EXIT, ENTRIES 1 TO NIRBDU HOLD INTEGER INFORMATION ON THE
C     FACTORS AND SHOULD BE PASSED UNCHANGED TO MA27C/CD.
C LIW LENGTH OF ARRAY IW.  AN INDICATION OF A SUITABLE VALUE WILL
C     HAVE BEEN PROVIDED IN NIRNEC AND NIRTOT BY MA27A/AD.
C     IT IS NOT ALTERED BY MA27B/BD.
C IKEEP MUST BE UNCHANGED SINCE THE CALL TO MA27A/AD.
C     IT IS NOT ALTERED BY MA27B/BD.
C NSTEPS MUST BE UNCHANGED SINCE THE CALL TO MA27A/AD.
C     IT IS NOT ALTERED BY MA27B/BD.
C MAXFRT NEED NOT BE SET AND MUST BE PASSED UNCHANGED TO MA27C/CD.
C     IT IS THE MAXIMUM SIZE OF THE FRONTAL MATRIX GENERATED DURING
C     THE DECOMPOSITION.
C IW1 USED AS WORKSPACE BY MA27B/BD.
C ICNTL is an INTEGER array of length 30, see MA27A/AD.
C CNTL is a DOUBLE PRECISION array of length 5, see MA27A/AD.
C INFO is an INTEGER array of length 20, see MA27A/AD.
C
C     .. Scalar Arguments ..
      INTEGER LA,LIW,MAXFRT,N,NSTEPS,NZ
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(LA),CNTL(5)
      INTEGER ICN(*),IKEEP(N,3),IRN(*),IW(LIW),IW1(N)
      INTEGER ICNTL(30),INFO(20)
C     ..
C     .. Local Scalars ..
      INTEGER I,IAPOS,IBLK,IPOS,IROWS,J1,J2,JJ,K,KBLK,KZ,LEN,NCOLS,
     +        NROWS,NZ1
C     ..
C     .. External Subroutines ..
      EXTERNAL MA27ND,MA27OD
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,MIN
C     ..
C     .. Executable Statements ..
      INFO(1) = 0

      IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0) GO TO 60
C PRINT INPUT PARAMETERS.
      WRITE (ICNTL(2),FMT=10) N,NZ,LA,LIW,NSTEPS,CNTL(1)

   10 FORMAT (/,/,
     + ' ENTERING MA27BD WITH      N     NZ       LA      LIW',
     +       ' NSTEPS      U',/,21X,I7,I7,I9,I9,I7,1PD10.2)

      KZ = MIN(6,NZ)
      IF (ICNTL(3).GT.1) KZ = NZ
      IF (NZ.GT.0) WRITE (ICNTL(2),FMT=20) (A(K),IRN(K),ICN(K),K=1,KZ)

   20 FORMAT (' MATRIX NON-ZEROS',/,1X,2 (1P,D16.3,2I6),/,
     +       (1X,1P,D16.3,2I6,1P,D16.3,2I6))

      K = MIN(9,N)
      IF (ICNTL(3).GT.1) K = N
      IF (K.GT.0) WRITE (ICNTL(2),FMT=30) (IKEEP(I,1),I=1,K)

   30 FORMAT (' IKEEP(.,1)=',10I6,/, (12X,10I6))

      K = MIN(K,NSTEPS)
      IF (K.GT.0) WRITE (ICNTL(2),FMT=40) (IKEEP(I,2),I=1,K)

   40 FORMAT (' IKEEP(.,2)=',10I6,/, (12X,10I6))

      IF (K.GT.0) WRITE (ICNTL(2),FMT=50) (IKEEP(I,3),I=1,K)

   50 FORMAT (' IKEEP(.,3)=',10I6,/, (12X,10I6))

   60 IF (N.LT.1 .OR. N.GT.ICNTL(4)) GO TO 70
      IF (NZ.LT.0) GO TO 100
      IF (LIW.LT.NZ) GO TO 120
      IF (LA.LT.NZ+N) GO TO 150
      IF (NSTEPS.LT.1 .OR. NSTEPS.GT.N) GO TO 175
C SORT
      CALL MA27ND(N,NZ,NZ1,A,LA,IRN,ICN,IW,LIW,IKEEP,IW1,ICNTL,INFO)
      IF (INFO(1).EQ.-3) GO TO 130
      IF (INFO(1).EQ.-4) GO TO 160
C FACTORIZE
      CALL MA27OD(N,NZ1,A,LA,IW,LIW,IKEEP,IKEEP(1,3),NSTEPS,MAXFRT,
     +           IKEEP(1,2),IW1,ICNTL,CNTL,INFO)
      IF (INFO(1).EQ.-3) GO TO 130
      IF (INFO(1).EQ.-4) GO TO 160
      IF (INFO(1).EQ.-5) GO TO 180
      IF (INFO(1).EQ.-6) GO TO 200
C **** WARNING MESSAGE ****
      IF (INFO(1).EQ.3 .AND. ICNTL(2).GT.0) THEN
        WRITE (ICNTL(2),FMT=65) INFO(1),INFO(2)
      END IF

   65 FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MA27BD',
     +        '  *** INFO(1) =',I2,
     +        /,5X,'MATRIX IS SINGULAR. RANK=',I5)

      GO TO 220
C **** ERROR RETURNS ****
   70 INFO(1) = -1
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)

   80 FORMAT (' **** ERROR RETURN FROM MA27BD **** INFO(1)=',I3)

      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=90) N

   90 FORMAT (' VALUE OF N OUT OF RANGE ... =',I10)

      GO TO 220

  100 INFO(1) = -2
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=110) NZ

  110 FORMAT (' VALUE OF NZ OUT OF RANGE .. =',I10)

      GO TO 220

  120 INFO(1) = -3
      INFO(2) = NZ
  130 IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=140) LIW,INFO(2)

  140 FORMAT (' LIW TOO SMALL, MUST BE INCREASED FROM',I10,' TO',
     +       ' AT LEAST',I10)

      GO TO 220

  150 INFO(1) = -4
      INFO(2) = NZ + N
  160 IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=170) LA,INFO(2)

  170 FORMAT (' LA TOO SMALL, MUST BE INCREASED FROM ',I10,' TO',
     +       ' AT LEAST',I10)

      GO TO 220

  175 INFO(1) = -7
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
      IF (ICNTL(1).GT.0) THEN
        WRITE (ICNTL(1),FMT='(A)') ' NSTEPS is out of range'
      END IF
      GO TO 220

  180 IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=190) INFO(2)

  190 FORMAT (' ZERO PIVOT AT STAGE',I10,
     +        ' WHEN INPUT MATRIX DECLARED DEFINITE')

      GO TO 220

  200 IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=80) INFO(1)
      IF (ICNTL(1).GT.0) WRITE (ICNTL(1),FMT=210)

  210 FORMAT (' CHANGE IN SIGN OF PIVOT ENCOUNTERED',
     +        ' WHEN FACTORING ALLEGEDLY DEFINITE MATRIX')

  220 IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0 .OR. INFO(1).LT.0) GO TO 310
C PRINT OUTPUT PARAMETERS.
      WRITE (ICNTL(2),FMT=230) MAXFRT,INFO(1),INFO(9),INFO(10),INFO(12),
     +  INFO(13),INFO(14),INFO(2)

  230 FORMAT (/,' LEAVING MA27BD WITH',
     +        /,10X,'  MAXFRT  INFO(1) NRLBDU NIRBDU NCMPBR',
     +         ' NCMPBI   NTWO IERROR',
     +        /,11X,8I7)

      IF (INFO(1).LT.0) GO TO 310
C PRINT OUT MATRIX FACTORS FROM MA27B/BD.
      KBLK = ABS(IW(1)+0)
      IF (KBLK.EQ.0) GO TO 310
      IF (ICNTL(3).EQ.1) KBLK = 1
      IPOS = 2
      IAPOS = 1
      DO 300 IBLK = 1,KBLK
        NCOLS = IW(IPOS)
        NROWS = IW(IPOS+1)
        J1 = IPOS + 2
        IF (NCOLS.GT.0) GO TO 240
        NCOLS = -NCOLS
        NROWS = 1
        J1 = J1 - 1
  240   WRITE (ICNTL(2),FMT=250) IBLK,NROWS,NCOLS

  250   FORMAT (' BLOCK PIVOT =',I8,' NROWS =',I8,' NCOLS =',I8)

        J2 = J1 + NCOLS - 1
        IPOS = J2 + 1
        WRITE (ICNTL(2),FMT=260) (IW(JJ),JJ=J1,J2)

  260   FORMAT (' COLUMN INDICES =',10I6,/, (17X,10I6))

        WRITE (ICNTL(2),FMT=270)

  270   FORMAT (' REAL ENTRIES .. EACH ROW STARTS ON A NEW LINE')

        LEN = NCOLS
        DO 290 IROWS = 1,NROWS
          J1 = IAPOS
          J2 = IAPOS + LEN - 1
          WRITE (ICNTL(2),FMT=280) (A(JJ),JJ=J1,J2)

  280     FORMAT (1P,5D13.3)

          LEN = LEN - 1
          IAPOS = J2 + 1
  290   CONTINUE
  300 CONTINUE
  310 RETURN
      END

      SUBROUTINE MA27CD(N,A,LA,IW,LIW,W,MAXFRT,RHS,IW1,NSTEPS,
     + ICNTL,INFO)
C THIS SUBROUTINE USES THE FACTORISATION OF THE MATRIX IN A,IW TO
C     SOLVE A SYSTEM OF EQUATIONS.
C N MUST BE SET TO THE ORDER OF THE MATRIX.  IT IS NOT ALTERED.
C A,IW HOLD INFORMATION ON THE FACTORS AND MUST BE UNCHANGED SINCE
C     THE CALL TO MA27B/BD. THEY ARE NOT ALTERED BY MA27C/CDD.
C LA,LIW MUST BE SET TO THE LENGTHS OF A,IW RESPECTIVELY.  THEY
C     ARE NOT ALTERED.
C W USED AS A WORK ARRAY.
C MAXFRT IS THE LENGTH OF W AND MUST BE PASSED UNCHANGED FROM THE
C     CALL TO MA27B/BD.  IT IS NOT ALTERED.
C RHS MUST BE SET TO THE RIGHT HAND SIDE FOR THE EQUATIONS BEING
C     SOLVED.  ON EXIT, THIS ARRAY WILL HOLD THE SOLUTION.
C IW1 USED AS A WORK ARRAY.
C NSTEPS IS THE LENGTH OF IW1 WHICH MUST BE AT LEAST THE ABSOLUTE
C     VALUE OF IW(1).  IT IS NOT ALTERED.
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 LA,LIW,MAXFRT,N,NSTEPS
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(LA),RHS(N),W(MAXFRT)
      INTEGER IW(LIW),IW1(NSTEPS),ICNTL(30),INFO(20)
C     ..
C     .. Local Scalars ..
      INTEGER I,IAPOS,IBLK,IPOS,IROWS,J1,J2,JJ,K,KBLK,LATOP,LEN,NBLK,
     +        NCOLS,NROWS
C     ..
C     .. External Subroutines ..
      EXTERNAL MA27QD,MA27RD
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,MIN
C     ..
C     .. Executable Statements ..
      INFO(1) = 0

      IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0) GO TO 110
C PRINT INPUT PARAMETERS.
      WRITE (ICNTL(2),FMT=10) N,LA,LIW,MAXFRT,NSTEPS

   10 FORMAT (/,/,' ENTERING MA27CD WITH      N     LA    LIW MAXFRT',
     +       '  NSTEPS',/,21X,5I7)
C PRINT OUT MATRIX FACTORS FROM MA27B/BD.
      KBLK = ABS(IW(1)+0)
      IF (KBLK.EQ.0) GO TO 90
      IF (ICNTL(3).EQ.1) KBLK = 1
      IPOS = 2
      IAPOS = 1
      DO 80 IBLK = 1,KBLK
        NCOLS = IW(IPOS)
        NROWS = IW(IPOS+1)
        J1 = IPOS + 2
        IF (NCOLS.GT.0) GO TO 20
        NCOLS = -NCOLS
        NROWS = 1
        J1 = J1 - 1
   20   WRITE (ICNTL(2),FMT=30) IBLK,NROWS,NCOLS

   30   FORMAT (' BLOCK PIVOT =',I8,' NROWS =',I8,' NCOLS =',I8)

        J2 = J1 + NCOLS - 1
        IPOS = J2 + 1
        WRITE (ICNTL(2),FMT=40) (IW(JJ),JJ=J1,J2)

   40   FORMAT (' COLUMN INDICES =',10I6,/, (17X,10I6))

        WRITE (ICNTL(2),FMT=50)

   50   FORMAT (' REAL ENTRIES .. EACH ROW STARTS ON A NEW LINE')

        LEN = NCOLS
        DO 70 IROWS = 1,NROWS
          J1 = IAPOS
          J2 = IAPOS + LEN - 1
          WRITE (ICNTL(2),FMT=60) (A(JJ),JJ=J1,J2)

   60     FORMAT (1P,5D13.3)

          LEN = LEN - 1
          IAPOS = J2 + 1
   70   CONTINUE
   80 CONTINUE
   90 K = MIN(10,N)
      IF (ICNTL(3).GT.1) K = N
      IF (N.GT.0) WRITE (ICNTL(2),FMT=100) (RHS(I),I=1,K)

  100 FORMAT (' RHS',1P,5D13.3,/, (4X,1P,5D13.3))

  110 IF (IW(1).LT.0) GO TO 130
      NBLK = IW(1)
      IF (NBLK.GT.0) GO TO 140
C WE HAVE A ZERO MATRIX
      DO 120 I = 1,N
        RHS(I) = 0.0D0
  120 CONTINUE
      GO TO 150

  130 NBLK = -IW(1)
C FORWARD SUBSTITUTION
  140 CALL MA27QD(N,A,LA,IW(2),LIW-1,W,MAXFRT,RHS,IW1,NBLK,LATOP,ICNTL)
C BACK SUBSTITUTION.
      CALL MA27RD(N,A,LA,IW(2),LIW-1,W,MAXFRT,RHS,IW1,NBLK,LATOP,ICNTL)
  150 IF (ICNTL(3).LE.0 .OR. ICNTL(2).LE.0) GO TO 170
C PRINT OUTPUT PARAMETERS.
      WRITE (ICNTL(2),FMT=160)

  160 FORMAT (/,/,' LEAVING MA27CD WITH')

      IF (N.GT.0) WRITE (ICNTL(2),FMT=100) (RHS(I),I=1,K)
  170 CONTINUE

      RETURN
      END
      SUBROUTINE MA27GD(N,NZ,IRN,ICN,IW,LW,IPE,IQ,FLAG,IWFR,
     +                 ICNTL,INFO)
C
C SORT PRIOR TO CALLING ANALYSIS ROUTINE MA27H/HD.
C
C GIVEN THE POSITIONS OF THE OFF-DIAGONAL NON-ZEROS OF A SYMMETRIC
C     MATRIX, CONSTRUCT THE SPARSITY PATTERN OF THE OFF-DIAGONAL
C     PART OF THE WHOLE MATRIX (UPPER AND LOWER TRIANGULAR PARTS).
C EITHER ONE OF A PAIR (I,J),(J,I) MAY BE USED TO REPRESENT
C     THE PAIR. DIAGONAL ELEMENTS AND DUPLICATES ARE IGNORED.
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 NUMBERS OF THE
C     NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED
C     TO IW (SEE DESCRIPTION OF IW).
C ICN(I),I=1,2,...,NZ MUST BE SET TO THE COLUMN NUMBERS OF THE
C     NON-ZEROS ON INPUT. IT IS NOT ALTERED UNLESS IT IS EQUIVALENCED
C     TO IW (SEE DESCRIPTION OF IW).
C IW NEED NOT BE SET ON INPUT. ON OUTPUT IT CONTAINS LISTS OF
C     COLUMN INDICES, EACH LIST BEING HEADED BY ITS LENGTH.
C     IRN(1) MAY BE EQUIVALENCED TO IW(1) AND ICN(1) MAY BE
C     EQUIVALENCED TO IW(K), WHERE K.GT.NZ.
C LW MUST BE SET TO THE LENGTH OF IW. IT MUST BE AT LEAST 2*NZ+N.
C     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),I=1,N CONTAINS THE NUMBER OF
C     OFF-DIAGONAL N0N-ZEROS IN ROW I INCLUDING DUPLICATES.
C FLAG IS USED FOR WORKSPACE TO HOLD FLAGS TO PERMIT DUPLICATE ENTRIES
C     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)
      INTEGER ICNTL(*),INFO(20)
C     ..
C     .. Local Scalars ..
      INTEGER I,ID,J,JN,K,K1,K2,L,LAST,LR,N1,NDUP
C     ..
C     .. Executable Statements ..
C
C INITIALIZE INFO(2) AND COUNT IN IPE THE
C     NUMBERS OF NON-ZEROS IN THE ROWS AND MOVE ROW AND COLUMN
C     NUMBERS INTO IW.
      INFO(2) = 0
      DO 10 I = 1,N
        IPE(I) = 0
   10 CONTINUE
      LR = NZ
      IF (NZ.EQ.0) GO TO 120
      DO 110 K = 1,NZ
        I = IRN(K)
        J = ICN(K)
        IF (I.LT.J) THEN
          IF (I.GE.1 .AND. J.LE.N) GO TO 90
        ELSE IF (I.GT.J) THEN
          IF (J.GE.1 .AND. I.LE.N) GO TO 90
        ELSE
          IF (I.GE.1 .AND. I.LE.N) GO TO 80
        END IF
        INFO(2) = INFO(2) + 1
        INFO(1) = 1
        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')

   80   I = 0
        J = 0
        GO TO 100

   90   IPE(I) = IPE(I) + 1
        IPE(J) = IPE(J) + 1
  100   IW(K) = J
        LR = LR + 1
        IW(LR) = I
  110 CONTINUE
C
C ACCUMULATE ROW COUNTS TO GET POINTERS TO ROW STARTS IN BOTH IPE AND IQ
C     AND INITIALIZE FLAG
  120 IQ(1) = 1
      N1 = N - 1
      IF (N1.LE.0) GO TO 140
      DO 130 I = 1,N1
        FLAG(I) = 0
        IF (IPE(I).EQ.0) IPE(I) = -1
        IQ(I+1) = IPE(I) + IQ(I) + 1
        IPE(I) = IQ(I)
  130 CONTINUE
  140 LAST = IPE(N) + IQ(N)
      FLAG(N) = 0
      IF (LR.GE.LAST) GO TO 160
      K1 = LR + 1
      DO 150 K = K1,LAST
        IW(K) = 0
  150 CONTINUE
  160 IPE(N) = IQ(N)
      IWFR = LAST + 1
      IF (NZ.EQ.0) GO TO 230
C
C RUN THROUGH PUTTING THE MATRIX ELEMENTS IN THE RIGHT PLACE
C     BUT WITH SIGNS INVERTED. IQ IS USED FOR HOLDING RUNNING POINTERS
C     AND IS LEFT HOLDING POINTERS TO ROW ENDS.
      DO 220 K = 1,NZ
        J = IW(K)
        IF (J.LE.0) GO TO 220
        L = K
        IW(K) = 0
        DO 210 ID = 1,NZ
          IF (L.GT.NZ) GO TO 170
          L = L + NZ
          GO TO 180

  170     L = L - NZ
  180     I = IW(L)
          IW(L) = 0
          IF (I.LT.J) GO TO 190
          L = IQ(J) + 1
          IQ(J) = L
          JN = IW(L)
          IW(L) = -I
          GO TO 200

  190     L = IQ(I) + 1
          IQ(I) = L
          JN = IW(L)
          IW(L) = -J
  200     J = JN
          IF (J.LE.0) GO TO 220
  210   CONTINUE
  220 CONTINUE
C
C RUN THROUGH RESTORING SIGNS, REMOVING DUPLICATES AND SETTING THE
C     MATE OF EACH NON-ZERO.
C NDUP COUNTS THE NUMBER OF DUPLICATE ELEMENTS.
  230 NDUP = 0
      DO 280 I = 1,N
        K1 = IPE(I) + 1
        K2 = IQ(I)
        IF (K1.LE.K2) GO TO 240
C ROW IS EMPTY. SET POINTER TO ZERO.
        IPE(I) = 0
        IQ(I) = 0
        GO TO 280
C ON ENTRY TO THIS LOOP FLAG(J).LT.I FOR J=1,2,...,N. DURING THE LOOP
C     FLAG(J) IS SET TO I IF A NON-ZERO IN COLUMN J IS FOUND. THIS
C     PERMITS DUPLICATES TO BE RECOGNIZED QUICKLY.
  240   DO 260 K = K1,K2
          J = -IW(K)
          IF (J.LE.0) GO TO 270
          L = IQ(J) + 1
          IQ(J) = L
          IW(L) = I
          IW(K) = J
          IF (FLAG(J).NE.I) GO TO 250
          NDUP = NDUP + 1
          IW(L) = 0
          IW(K) = 0
  250     FLAG(J) = I
  260   CONTINUE
  270   IQ(I) = IQ(I) - IPE(I)
        IF (NDUP.EQ.0) IW(K1-1) = IQ(I)
  280 CONTINUE
      IF (NDUP.EQ.0) GO TO 310
C
C COMPRESS IW TO REMOVE DUMMY ENTRIES CAUSED BY DUPLICATES.
      IWFR = 1
      DO 300 I = 1,N
        K1 = IPE(I) + 1
        IF (K1.EQ.1) GO TO 300
        K2 = IQ(I) + IPE(I)
        L = IWFR
        IPE(I) = IWFR
        IWFR = IWFR + 1
        DO 290 K = K1,K2
          IF (IW(K).EQ.0) GO TO 290
          IW(IWFR) = IW(K)
          IWFR = IWFR + 1
  290   CONTINUE
        IW(L) = IWFR - L - 1
  300 CONTINUE
  310 RETURN

      END

      SUBROUTINE MA27HD(N,IPE,IW,LW,IWFR,NV,NXT,LST,IPD,FLAG,IOVFLO,
     +                 NCMPA,FRATIO)
C Bug was found in the then identical subroutine MA57HD which was then
C     corrected.
C     Changes made in September 2009 because of bug in compress control
C     found by Nick.
C
C ANALYSIS SUBROUTINE
C
C GIVEN REPRESENTATION OF THE WHOLE MATRIX (EXCLUDING DIAGONAL)
C     PERFORM MINIMUM DEGREE ORDERING, CONSTRUCTING TREE POINTERS.
C     IT WORKS WITH SUPERVARIABLES WHICH ARE COLLECTIONS OF ONE OR MORE
C     VARIABLES, STARTING WITH SUPERVARIABLE I CONTAINING VARIABLE I FOR
C     I = 1,2,...,N. ALL VARIABLES IN A SUPERVARIABLE ARE ELIMINATED
C     TOGETHER. EACH SUPERVARIABLE HAS AS NUMERICAL NAME THAT OF ONE
C     OF ITS VARIABLES (ITS PRINCIPAL VARIABLE).
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. IF
C     SUPERVARIABLE I IS ABSORBED INTO SUPERVARIABLE J THEN IPE(I)=-J.
C     IF SUPERVARIABLE I IS ELIMINATED THEN IPE(I) EITHER POINTS TO THE
C     LIST OF SUPERVARIABLES FOR CREATED ELEMENT I OR IS ZERO IF
C     THE CREATED ELEMENT IS NULL. IF ELEMENT I
C     IS ABSORBED INTO 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.
C     DURING EXECUTION THESE LISTS ARE REVISED AND HOLD
C     LISTS OF ELEMENTS AND SUPERVARIABLES. THE ELEMENTS
C     ALWAYS HEAD THE LISTS. WHEN A SUPERVARIABLE
C     IS ELIMINATED ITS LIST IS REPLACED BY A LIST OF SUPERVARIABLES
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 AND CONTINUES TO HAVE THIS MEANING.
C NV(JS) NEED NOT BE SET. DURING EXECUTION IT IS ZERO IF
C     JS IS NOT A PRINCIPAL VARIABLE AND IF IT IS IT HOLDS
C     THE NUMBER OF VARIABLES IN SUPERVARIABLE JS. FOR ELIMINATED
C     VARIABLES IT IS SET TO THE DEGREE AT THE TIME OF ELIMINATION.
C NXT(JS) NEED NOT BE SET. DURING EXECUTION IT IS THE NEXT
C     SUPERVARIABLE HAVING THE SAME DEGREE AS JS, OR ZERO
C     IF IT IS LAST IN ITS LIST.
C LST(JS) NEED NOT BE SET. DURING EXECUTION IT IS THE
C     LAST SUPERVARIABLE HAVING THE SAME DEGREE AS JS OR
C     -(ITS DEGREE) IF IT IS FIRST IN ITS LIST.
C IPD(ID) NEED NOT BE SET. DURING EXECUTION IT
C     IS THE FIRST SUPERVARIABLE WITH DEGREE ID OR ZERO
C     IF THERE ARE NONE.
C FLAG IS USED AS WORKSPACE FOR ELEMENT AND SUPERVARIABLE FLAGS.
C     WHILE THE CODE IS FINDING THE DEGREE OF SUPERVARIABLE IS
C     FLAG HAS THE FOLLOWING VALUES.
C     A) FOR THE CURRENT PIVOT/NEW ELEMENT ME
C           FLAG(ME)=-1
C     B) FOR VARIABLES JS
C           FLAG(JS)=-1 IF JS IS NOT A PRINCIPAL VARIABLE
C           FLAG(JS)=0 IF JS IS A SUPERVARIABLE IN THE NEW ELEMENT
C           FLAG(JS)=NFLG IF JS IS A SUPERVARIABLE NOT IN THE NEW
C                 ELEMENT THAT HAS BEEN COUNTED IN THE DEGREE
C                 CALCULATION
C           FLAG(JS).GT.NFLG IF JS IS A SUPERVARIABLE NOT IN THE NEW
C                 ELEMENT THAT HAS NOT BEEN COUNTED IN THE DEGREE
C                 CALCULATION
C     C) FOR ELEMENTS IE
C           FLAG(IE)=-1 IF ELEMENT IE HAS BEEN MERGED INTO ANOTHER
C           FLAG(IE)=-NFLG IF ELEMENT IE HAS BEEN USED IN THE DEGREE
C                 CALCULATION FOR IS.
C           FLAG(IE).LT.-NFLG IF ELEMENT IE HAS NOT BEEN USED IN THE
C                 DEGREE CALCULATION FOR IS
C IOVFLO see ICNTL(4) in MA27A/AD.
C NCMPA see INFO(11) in MA27A/AD.
C FRATIO see CNTL(2) in MA27A/AD.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION FRATIO
      INTEGER IWFR,LW,N,IOVFLO,NCMPA
C     ..
C     .. Array Arguments ..
      INTEGER FLAG(N),IPD(N),IPE(N),IW(LW),LST(N),NV(N),NXT(N)
C     ..
C     .. Local Scalars ..
      INTEGER I,ID,IDL,IDN,IE,IP,IS,JP,JP1,JP2,JS,K,K1,K2,KE,KP,KP0,KP1,
     +        KP2,KS,L,LEN,LIMIT,LN,LS,LWFR,MD,ME,ML,MS,NEL,NFLG,NP,
     +        NP0,NS,NVPIV,NVROOT,ROOT
C LIMIT  Limit on number of variables for putting node in root.
C NVROOT Number of variables in the root node
C ROOT   Index of the root node (N+1 if none chosen yet).
C     ..
C     .. External Subroutines ..
      EXTERNAL MA27UD
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,MIN
C     ..
C If a column of the reduced matrix has relative density greater than
C CNTL(2), it is forced into the root. All such columns are taken to
C have sparsity pattern equal to their merged patterns, so the fill
C and operation counts may be overestimated.
C
C IS,JS,KS,LS,MS,NS are used to refer to supervariables.
C IE,JE,KE are used to refer to elements.
C IP,JP,KP,K,NP are used to point to lists of elements
C     or supervariables.
C ID is used for the degree of a supervariable.
C MD is used for the current minimum degree.
C IDN is used for the no. of variables in a newly created element
C NEL is used to hold the no. of variables that have been
C     eliminated.
C ME=MS is the name of the supervariable eliminated and
C     of the element created in the main loop.
C NFLG is used for the current flag value in array FLAG. It starts
C     with the value IOVFLO and is reduced by 1 each time it is used
C     until it has the value 2 when it is reset to the value IOVFLO.
C
C     .. Executable Statements ..
C Initializations
      DO 10 I = 1,N
        IPD(I) = 0
        NV(I) = 1
        FLAG(I) = IOVFLO
   10 CONTINUE
      MD = 1
      NCMPA = 0
      NFLG = IOVFLO
      NEL = 0
      ROOT = N+1