ma27ad.f

   20 DO 30 I = 1,N
        LSTKI(I) = 1
   30 CONTINUE
C ACCUMULATE NUMBER OF NON-ZEROS WITH INDICES IN RANGE IN NZ1
C     DUPLICATES ON THE DIAGONAL ARE IGNORED BUT NZ1 INCLUDES ANY
C     DIAGONALS NOT PRESENT ON INPUT.
C ACCUMULATE ROW COUNTS IN LSTKI.
      NZ1 = N
      IF (NZ.EQ.0) GO TO 50
      DO 40 I = 1,NZ
        IOLD = IRN(I)
        JOLD = ICN(I)
C JUMP IF INDEX IS OUT OF RANGE.
        IF (IOLD.LT.1 .OR. IOLD.GT.N) GO TO 40
        IF (JOLD.LT.1 .OR. JOLD.GT.N) GO TO 40
        IF (IOLD.EQ.JOLD) GO TO 40
        NZ1 = NZ1 + 1
        IROW = MIN(PERM(IOLD)+0,PERM(JOLD)+0)
        LSTKI(IROW) = LSTKI(IROW) + 1
   40 CONTINUE
   50 NZ2 = NZ1
C ISTKR,ISTKI CURRENT NUMBER OF STACK ENTRIES IN
C     REAL AND INTEGER STORAGE RESPECTIVELY.
C OPS,NRLADU,NIRADU,NIRTOT,NRLTOT,NIRNEC,NRLNEC,NZ2 ARE DEFINED ABOVE.
C NZ2 CURRENT NUMBER OF ORIGINAL MATRIX ENTRIES NOT YET PROCESSED.
C NUMORG CURRENT TOTAL NUMBER OF ROWS ELIMINATED.
C ITOP CURRENT NUMBER OF ELEMENTS ON THE STACK.
   60 ISTKI = 0
      ISTKR = 0
      OPS = 0.0D0
      NRLADU = 0
C     ONE LOCATION IS NEEDED TO RECORD THE NUMBER OF BLOCKS
C     ACTUALLY USED.
      NIRADU = 1
      NIRTOT = NZ1
      NRLTOT = NZ1
      NIRNEC = NZ2
      NRLNEC = NZ2
      NUMORG = 0
      ITOP = 0
C
C EACH PASS THROUGH THIS LOOP PROCESSES A NODE OF THE TREE.
      DO 100 ITREE = 1,NSTEPS
        NELIM = NE(ITREE)
        DELIM = NELIM
        NFR = ND(ITREE)
        NSTK = NA(ITREE)
C ADJUST STORAGE COUNTS ON ASSEMBLY OF CURRENT FRONTAL MATRIX.
        NASSR = NFR* (NFR+1)/2
        IF (NSTK.NE.0) NASSR = NASSR - LSTKR(ITOP) + 1
        NRLTOT = MAX(NRLTOT,NRLADU+NASSR+ISTKR+NZ1)
        NIRTOT = MAX(NIRTOT,NIRADU+NFR+2+ISTKI+NZ1)
        NRLNEC = MAX(NRLNEC,NRLADU+NASSR+ISTKR+NZ2)
        NIRNEC = MAX(NIRNEC,NIRADU+NFR+2+ISTKI+NZ2)
C DECREASE NZ2 BY THE NUMBER OF ENTRIES IN ROWS BEING ELIMINATED AT
C     THIS STAGE.
        DO 70 IORG = 1,NELIM
          JORG = NUMORG + IORG
          NZ2 = NZ2 - LSTKI(JORG)
   70   CONTINUE
        NUMORG = NUMORG + NELIM
C JUMP IF THERE ARE NO STACK ASSEMBLIES AT THIS NODE.
        IF (NSTK.LE.0) GO TO 90
C REMOVE ELEMENTS FROM THE STACK.  THERE ARE ITOP ELEMENTS ON THE
C     STACK WITH THE APPROPRIATE ENTRIES IN LSTKR,LSTKI GIVING
C     THE REAL AND INTEGER STORAGE RESPECTIVELY FOR EACH STACK
C     ELEMENT.
        DO 80 K = 1,NSTK
          LSTK = LSTKR(ITOP)
          ISTKR = ISTKR - LSTK
          LSTK = LSTKI(ITOP)
          ISTKI = ISTKI - LSTK
          ITOP = ITOP - 1
   80   CONTINUE
C ACCUMULATE NON-ZEROS IN FACTORS AND NUMBER OF OPERATIONS.
   90   NRLADU = NRLADU + (NELIM* (2*NFR-NELIM+1))/2
        NIRADU = NIRADU + 2 + NFR
        IF (NELIM.EQ.1) NIRADU = NIRADU - 1
        OPS = OPS + ((NFR*DELIM*(NFR+1)-(2*NFR+1)*DELIM*(DELIM+1)/2+
     +        DELIM* (DELIM+1)* (2*DELIM+1)/6)/2)
        IF (ITREE.EQ.NSTEPS) GO TO 100
C JUMP IF ALL OF FRONTAL MATRIX HAS BEEN ELIMINATED.
        IF (NFR.EQ.NELIM) GO TO 100
C STACK REMAINDER OF ELEMENT.
        ITOP = ITOP + 1
        LSTKR(ITOP) = (NFR-NELIM)* (NFR-NELIM+1)/2
        LSTKI(ITOP) = NFR - NELIM + 1
        ISTKI = ISTKI + LSTKI(ITOP)
        ISTKR = ISTKR + LSTKR(ITOP)
C WE DO NOT NEED TO ADJUST THE COUNTS FOR THE REAL STORAGE BECAUSE
C     THE REMAINDER OF THE FRONTAL MATRIX IS SIMPLY MOVED IN THE
C     STORAGE FROM FACTORS TO STACK AND NO EXTRA STORAGE IS REQUIRED.
        NIRTOT = MAX(NIRTOT,NIRADU+ISTKI+NZ1)
        NIRNEC = MAX(NIRNEC,NIRADU+ISTKI+NZ2)
  100 CONTINUE
C
C ADJUST THE STORAGE COUNTS TO ALLOW FOR THE USE OF THE REAL AND
C     INTEGER STORAGE FOR PURPOSES OTHER THAN PURELY THE
C     FACTORIZATION ITSELF.
C THE SECOND TWO TERMS ARE THE MINUMUM AMOUNT REQUIRED BY MA27N/ND.
      NRLNEC = MAX(NRLNEC,N+MAX(NZ,NZ1))
      NRLTOT = MAX(NRLTOT,N+MAX(NZ,NZ1))
      NRLNEC = MIN(NRLNEC,NRLTOT)
      NIRNEC = MAX(NZ,NIRNEC)
      NIRTOT = MAX(NZ,NIRTOT)
      NIRNEC = MIN(NIRNEC,NIRTOT)

      INFO(3) = NRLTOT
      INFO(4) = NIRTOT
      INFO(5) = NRLNEC
      INFO(6) = NIRNEC
      INFO(7) = NRLADU
      INFO(8) = NIRADU
      RETURN

      END 
      SUBROUTINE MA27ND(N,NZ,NZ1,A,LA,IRN,ICN,IW,LIW,PERM,IW2,ICNTL,
     +                 INFO)
C
C SORT PRIOR TO FACTORIZATION USING MA27O/OD.
C
C THIS SUBROUTINE REORDERS THE USER'S INPUT SO THAT THE UPPER TRIANGLE
C     OF THE PERMUTED MATRIX, INCLUDING THE DIAGONAL, IS
C     HELD ORDERED BY ROWS AT THE END OF THE STORAGE FOR A AND IW.
C     IT IGNORES ENTRIES WITH ONE OR BOTH INDICES OUT OF RANGE AND
C     ACCUMULATES DIAGONAL ENTRIES.
C     IT ADDS EXPLICIT ZEROS ON THE DIAGONAL WHERE NECESSARY.
C N      - MUST BE SET TO THE ORDER OF THE MATRIX.
C          IT IS NOT ALTERED BY MA27N/ND.
C NZ     - ON ENTRY NZ MUST BE SET TO THE NUMBER
C          OF NON-ZEROS INPUT.  NOT ALTERED BY MA27N/ND.
C NZ1    - ON EXIT NZ1 WILL BE EQUAL TO THE NUMBER OF ENTRIES IN THE
C          SORTED MATRIX.
C A      - ON ENTRY A(I) MUST
C          HOLD THE VALUE OF THE ORIGINAL MATRIX ELEMENT IN POSITION
C          (IRN(I),ICN(I)),I=1,NZ.  ON EXIT A(LA-NZ1+I),I=1,NZ1 HOLDS
C          THE UPPER TRIANGLE OF THE PERMUTED MATRIX BY ROWS WITH
C          THE DIAGONAL ENTRY FIRST ALTHOUGH THERE IS NO FURTHER
C          ORDERING WITHIN THE ROWS THEMSELVES.
C LA     - LENGTH OF ARRAY A. MUST BE AT LEAST N+MAX(NZ,NZ1).
C          IT IS NOT ALTERED BY MA27N/ND.
C IRN    - IRN(I) MUST BE SET TO
C          HOLD THE ROW INDEX OF ENTRY A(I),I=1,NZ.  IRN WILL BE
C          UNALTERED BY MA27N/ND, UNLESS IT IS EQUIVALENCED WITH IW.
C ICN    - ICN(I) MUST BE SET TO
C          HOLD THE COLUMN INDEX OF ENTRY A(I),I=1,NZ.  ICN WILL BE
C          UNALTERED BY MA27N/ND, UNLESS IT IS EQUIVALENCED WITH IW.
C IW     - USED AS WORKSPACE AND ON
C          EXIT, ENTRIES IW(LIW-NZ1+I),I=1,NZ1 HOLD THE COLUMN INDICES
C          (THE ORIGINAL UNPERMUTED INDICES) OF THE CORRESPONDING ENTRY
C          OF A WITH THE FIRST ENTRY FOR EACH ROW FLAGGED NEGATIVE.
C          IRN(1) MAY BE EQUIVALENCED WITH IW(1) AND ICN(1) MAY BE
C          EQUIVALENCED WITH IW(K) WHERE K.GT.NZ.
C LIW    - LENGTH OF ARRAY IW. MUST BE AT LEAST AS
C          GREAT AS THE MAXIMUM OF NZ AND NZ1.
C          NOT ALTERED BY MA27N/ND.
C PERM   - PERM(I) HOLDS THE
C          POSITION IN THE TENTATIVE PIVOT ORDER OF ROW I IN THE
C          ORIGINAL MATRIX,I=1,N. NOT ALTERED BY MA27N/ND.
C IW2    - USED AS WORKSPACE.
C          SEE COMMENTS IN CODE IMMEDIATELY PRIOR TO
C          EACH USE.
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   INFO(1)  - ON EXIT FROM MA27N/ND, A ZERO VALUE OF
C          INFO(1) INDICATES THAT NO ERROR HAS BEEN DETECTED.
C          POSSIBLE NON-ZERO VALUES ARE ..
C          +1  WARNING.  INDICES OUT OF RANGE.  THESE ARE IGNORED,
C              THEIR NUMBER IS RECORDED IN INFO(2) OF MA27E/ED AND
C              MESSAGES IDENTIFYING THE FIRST TEN ARE OUTPUT ON UNIT
C              ICNTL(2).
C          -3  INTEGER ARRAY IW IS TOO SMALL.
C          -4  DOUBLE PRECISION ARRAY A IS TOO SMALL.
C
C     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
C     ..
C     .. Scalar Arguments ..
      INTEGER LA,LIW,N,NZ,NZ1
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(LA)
      INTEGER ICN(*),IRN(*),IW(LIW),IW2(N),PERM(N),ICNTL(30),INFO(20)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ANEXT,ANOW
      INTEGER I,IA,ICH,II,IIW,INEW,IOLD,IPOS,J1,J2,JJ,JNEW,JOLD,JPOS,K
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC MIN
C     ..
C     .. Executable Statements ..
      INFO(1) = 0
C INITIALIZE WORK ARRAY (IW2) IN PREPARATION FOR
C     COUNTING NUMBERS OF NON-ZEROS IN THE ROWS AND INITIALIZE
C     LAST N ENTRIES IN A WHICH WILL HOLD THE DIAGONAL ENTRIES
      IA = LA
      DO 10 IOLD = 1,N
        IW2(IOLD) = 1
        A(IA) = ZERO
        IA = IA - 1
   10 CONTINUE
C SCAN INPUT COPYING ROW INDICES FROM IRN TO THE FIRST NZ POSITIONS
C     IN IW.  THE NEGATIVE OF THE INDEX IS HELD TO FLAG ENTRIES FOR
C     THE IN-PLACE SORT.  ENTRIES IN IW CORRESPONDING TO DIAGONALS AND
C     ENTRIES WITH OUT-OF-RANGE INDICES ARE SET TO ZERO.
C     FOR DIAGONAL ENTRIES, REALS ARE ACCUMULATED IN THE LAST N
C     LOCATIONS OF A.
C     THE NUMBER OF ENTRIES IN EACH ROW OF THE PERMUTED MATRIX IS
C     ACCUMULATED IN IW2.
C INDICES OUT OF RANGE ARE IGNORED  AFTER BEING COUNTED AND
C     AFTER APPROPRIATE MESSAGES HAVE BEEN PRINTED.
      INFO(2) = 0
C NZ1 IS THE NUMBER OF NON-ZEROS HELD AFTER INDICES OUT OF RANGE HAVE
C     BEEN IGNORED AND DIAGONAL ENTRIES ACCUMULATED.
      NZ1 = N
      IF (NZ.EQ.0) GO TO 80
      DO 70 K = 1,NZ
        IOLD = IRN(K)
        IF (IOLD.GT.N .OR. IOLD.LE.0) GO TO 30
        JOLD = ICN(K)
        IF (JOLD.GT.N .OR. JOLD.LE.0) GO TO 30
        INEW = PERM(IOLD)
        JNEW = PERM(JOLD)
        IF (INEW.NE.JNEW) GO TO 20
        IA = LA - N + IOLD
        A(IA) = A(IA) + A(K)
        GO TO 60

   20   INEW = MIN(INEW,JNEW)
C INCREMENT NUMBER OF ENTRIES IN ROW INEW.
        IW2(INEW) = IW2(INEW) + 1
        IW(K) = -IOLD
        NZ1 = NZ1 + 1
        GO TO 70
C ENTRY OUT OF RANGE.  IT WILL BE IGNORED AND A FLAG SET.
   30   INFO(1) = 1
        INFO(2) = INFO(2) + 1
        IF (INFO(2).LE.1 .AND. ICNTL(2).GT.0) THEN
          WRITE (ICNTL(2),FMT=40) INFO(1)
        ENDIF

   40   FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MA27BD',
     +          '  *** INFO(1) =',I2)

        IF (INFO(2).LE.10 .AND. ICNTL(2).GT.0) THEN
          WRITE (ICNTL(2),FMT=50) K,IRN(K),ICN(K)
        END IF

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

   60   IW(K) = 0
   70 CONTINUE
C CALCULATE POINTERS (IN IW2) TO THE POSITION IMMEDIATELY AFTER THE END
C     OF EACH ROW.
   80 IF (NZ.LT.NZ1 .AND. NZ1.NE.N) GO TO 100
C ROOM IS INCLUDED FOR THE DIAGONALS.
      K = 1
      DO 90 I = 1,N
        K = K + IW2(I)
        IW2(I) = K
   90 CONTINUE
      GO TO 120
C ROOM IS NOT INCLUDED FOR THE DIAGONALS.
  100 K = 1
      DO 110 I = 1,N
        K = K + IW2(I) - 1
        IW2(I) = K
  110 CONTINUE
C FAIL IF INSUFFICIENT SPACE IN ARRAYS A OR IW.
  120 IF (NZ1.GT.LIW) GO TO 210
      IF (NZ1+N.GT.LA) GO TO 220
C NOW RUN THROUGH NON-ZEROS IN ORDER PLACING THEM IN THEIR NEW
C POSITION AND DECREMENTING APPROPRIATE IW2 ENTRY.  IF WE ARE
C ABOUT TO OVERWRITE AN ENTRY NOT YET MOVED, WE MUST DEAL WITH
C THIS AT THIS TIME.
      IF (NZ1.EQ.N) GO TO 180
      DO 140 K = 1,NZ
        IOLD = -IW(K)
        IF (IOLD.LE.0) GO TO 140
        JOLD = ICN(K)
        ANOW = A(K)
        IW(K) = 0
        DO 130 ICH = 1,NZ
          INEW = PERM(IOLD)
          JNEW = PERM(JOLD)
          INEW = MIN(INEW,JNEW)
          IF (INEW.EQ.PERM(JOLD)) JOLD = IOLD
          JPOS = IW2(INEW) - 1
          IOLD = -IW(JPOS)
          ANEXT = A(JPOS)
          A(JPOS) = ANOW
          IW(JPOS) = JOLD
          IW2(INEW) = JPOS
          IF (IOLD.EQ.0) GO TO 140
          ANOW = ANEXT
          JOLD = ICN(JPOS)
  130   CONTINUE
  140 CONTINUE
      IF (NZ.GE.NZ1) GO TO 180
C MOVE UP ENTRIES TO ALLOW FOR DIAGONALS.
      IPOS = NZ1
      JPOS = NZ1 - N
      DO 170 II = 1,N
        I = N - II + 1
        J1 = IW2(I)
        J2 = JPOS
        IF (J1.GT.JPOS) GO TO 160
        DO 150 JJ = J1,J2
          IW(IPOS) = IW(JPOS)
          A(IPOS) = A(JPOS)
          IPOS = IPOS - 1
          JPOS = JPOS - 1
  150   CONTINUE
  160   IW2(I) = IPOS + 1
        IPOS = IPOS - 1
  170 CONTINUE
C RUN THROUGH ROWS INSERTING DIAGONAL ENTRIES AND FLAGGING BEGINNING
C     OF EACH ROW BY NEGATING FIRST COLUMN INDEX.
  180 DO 190 IOLD = 1,N
        INEW = PERM(IOLD)
        JPOS = IW2(INEW) - 1
        IA = LA - N + IOLD
        A(JPOS) = A(IA)
        IW(JPOS) = -IOLD
  190 CONTINUE
C MOVE SORTED MATRIX TO THE END OF THE ARRAYS.
      IPOS = NZ1
      IA = LA
      IIW = LIW
      DO 200 I = 1,NZ1
        A(IA) = A(IPOS)
        IW(IIW) = IW(IPOS)
        IPOS = IPOS - 1
        IA = IA - 1
        IIW = IIW - 1
  200 CONTINUE
      GO TO 230
C **** ERROR RETURN ****
  210 INFO(1) = -3
      INFO(2) = NZ1
      GO TO 230

  220 INFO(1) = -4
      INFO(2) = NZ1 + N
C
  230 RETURN

      END
      SUBROUTINE MA27OD(N,NZ,A,LA,IW,LIW,PERM,NSTK,NSTEPS,MAXFRT,NELIM,
     +                 IW2,ICNTL,CNTL,INFO)
C
C FACTORIZATION SUBROUTINE
C
C THIS SUBROUTINE OPERATES ON THE INPUT MATRIX ORDERED BY MA27N/ND AND
C     PRODUCES THE FACTORS OF U AND D ('A'=UTRANSPOSE*D*U) FOR USE IN
C     SUBSEQUENT SOLUTIONS.  GAUSSIAN ELIMINATION IS USED WITH PIVOTS
C     CHOSEN FROM THE DIAGONAL.  TO ENSURE STABILITY, BLOCK PIVOTS OF
C     ORDER 2 WILL BE USED IF THE DIAGONAL ENTRY IS NOT LARGE ENOUGH.
C
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 IN UPPER TRIANGLE OF
C          PERMUTED MATRIX.  NOT ALTERED BY MA27O/OD.
C A      - MUST BE SET ON INPUT TO MATRIX HELD BY ROWS REORDERED BY
C          PERMUTATION FROM MA27A/AD IN A(LA-NZ+I),I=1,NZ.   ON
C          EXIT FROM MA27O/OD, THE FACTORS OF U AND D ARE HELD IN
C          POSITIONS 1 TO POSFAC-1.
C LA     - LENGTH OF ARRAY A.  A VALUE FOR LA
C          SUFFICIENT FOR DEFINITE SYSTEMS
C          WILL HAVE BEEN PROVIDED BY MA27A/AD. NOT ALTERED BY MA27O/OD.
C IW     - MUST BE SET SO THAT,ON INPUT, IW(LIW-NZ+I),I=1,NZ
C          HOLDS THE COLUMN INDEX OF THE ENTRY IN A(LA-NZ+I).  ON EXIT,
C          IW HOLDS INTEGER INDEXING INFORMATION ON THE FACTORS.
C          THE ABSOLUTE VALUE OF THE FIRST ENTRY IN IW WILL BE SET TO
C          THE NUMBER OF BLOCK PIVOTS ACTUALLY USED.  THIS MAY BE
C          DIFFERENT FROM NSTEPS SINCE NUMERICAL CONSIDERATIONS
C          MAY PREVENT US CHOOSING A PIVOT AT EACH STAGE.  IF THIS ENTRY
C          IN IW IS NEGATIVE, THEN AT LEAST ONE TWO BY TWO
C          PIVOT HAS BEEN USED DURING THE DECOMPOSITION.
C          INTEGER INFORMATION ON EACH BLOCK PIVOT ROW FOLLOWS.  FOR
C          EACH BLOCK PIVOT ROW THE COLUMN INDICES ARE PRECEDED BY A
C          COUNT OF THE NUMBER OF ROWS AND COLUMNS IN THE BLOCK PIVOT
C          WHERE, IF ONLY ONE ROW IS PRESENT, ONLY THE NUMBER OF
C          COLUMNS TOGETHER WITH A NEGATIVE FLAG IS HELD.  THE FIRST
C          COLUMN INDEX FOR A TWO BY TWO PIVOT IS FLAGGED NEGATIVE.
C LIW    - LENGTH OF ARRAY IW.  A VALUE FOR LIW SUFFICIENT FOR
C          DEFINITE SYSTEMS
C          WILL HAVE BEEN PROVIDED BY MA27A/AD.  NOT ALTERED BY MA27O/OD
C PERM   - PERM(I) MUST BE SET TO POSITION OF ROW/COLUMN I IN THE
C          TENTATIVE PIVOT ORDER GENERATED BY MA27A/AD.
C          IT IS NOT ALTERED BY MA27O/OD.
C NSTK   - MUST BE LEFT UNCHANGED SINCE OUTPUT FROM MA27A/AD. NSTK(I)
C          GIVES THE NUMBER OF GENERATED STACK ELEMENTS ASSEMBLED AT
C          STAGE I.  IT IS NOT ALTERED BY MA27O/OD.
C NSTEPS - LENGTH OF ARRAYS NSTK AND NELIM. VALUE IS GIVEN ON OUTPUT
C          FROM MA27A/AD (WILL NEVER EXCEED N). IT IS NOT ALTERED BY
C          MA27O/OD.
C MAXFRT - NEED NOT BE SET ON INPUT.  ON OUTPUT
C          MAXFRT WILL BE SET TO THE MAXIMUM FRONT SIZE ENCOUNTERED
C          DURING THE DECOMPOSITION.
C NELIM  - MUST BE UNCHANGED SINCE OUTPUT FROM MA27A/AD. NELIM(I)
C          GIVES THE NUMBER OF ORIGINAL ROWS ASSEMBLED AT STAGE I.
C          IT IS NOT ALTERED BY MA27O/OD.
C IW2    - INTEGER ARRAY OF LENGTH N. USED AS WORKSPACE BY MA27O/OD.
C          ALTHOUGH WE COULD HAVE USED A SHORT WORD INTEGER IN THE IBM
C          VERSION, WE HAVE NOT DONE SO BECAUSE THERE IS A SPARE
C          FULL INTEGER ARRAY (USED IN THE SORT BEFORE MA27O/OD)
C          AVAILABLE WHEN MA27O/OD IS CALLED FROM 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   INFO(1)  - INTEGER VARIABLE.  DIAGNOSTIC FLAG.  A ZERO VALUE ON EXIT
C          INDICATES SUCCESS.  POSSIBLE NEGATIVE VALUES ARE ...
C          -3  INSUFFICIENT STORAGE FOR IW.
C          -4  INSUFFICIENT STORAGE FOR A.
C          -5  ZERO PIVOT FOUND IN FACTORIZATION OF DEFINITE MATRIX.
C
C     .. Parameters ..
      DOUBLE PRECISION ZERO,HALF,ONE
      PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0)
C     ..
C     .. Scalar Arguments ..
      INTEGER LA,LIW,MAXFRT,N,NSTEPS,NZ
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(LA),CNTL(5)
      INTEGER IW(LIW),IW2(N),NELIM(NSTEPS),NSTK(NSTEPS),PERM(N)
      INTEGER ICNTL(30),INFO(20)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION AMAX,AMULT,AMULT1,AMULT2,DETPIV,RMAX,SWOP,
     +        THRESH,TMAX,UU
      INTEGER AINPUT,APOS,APOS1,APOS2,APOS3,ASTK,ASTK2,AZERO,I,IASS,
     +        IBEG,IDUMMY,IELL,IEND,IEXCH,IFR,IINPUT,IOLDPS,IORG,IPIV,
     +        IPMNP,IPOS,IROW,ISNPIV,ISTK,ISTK2,ISWOP,IWPOS,IX,IY,J,J1,
     +        J2,JCOL,JDUMMY,JFIRST,JJ,JJ1,JJJ,JLAST,JMAX,JMXMIP,JNEW,
     +        JNEXT,JPIV,JPOS,K,K1,K2,KDUMMY,KK,KMAX,KROW,LAELL,LAPOS2,
     +        LIELL,LNASS,LNPIV,LT,LTOPST,NASS,NBLK,NEWEL,NFRONT,NPIV,
     +        NPIVP1,NTOTPV,NUMASS,NUMORG,NUMSTK,PIVSIZ,POSFAC,POSPV1,
     +        POSPV2
      INTEGER IDIAG
C IDIAG IS A TEMPORARY FOR THE DISPLACEMENT FROM THE START OF THE 
C     ASSEMBLED MATRIX (OF ORDER IX) OF THE DIAGONAL ENTRY IN ITS ROW IY.
      INTEGER NTWO,NEIG,NCMPBI,NCMPBR,NRLBDU,NIRBDU
C     ..
C     .. External Subroutines ..
      EXTERNAL MA27PD
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,MAX,MIN
C     ..
C     .. Executable Statements ..
C INITIALIZATION.
C NBLK IS THE NUMBER OF BLOCK PIVOTS USED.
      NBLK = 0
      NTWO = 0
      NEIG = 0
      NCMPBI = 0
      NCMPBR = 0
      MAXFRT = 0
      NRLBDU = 0
      NIRBDU = 0
C A PRIVATE VARIABLE UU IS SET TO CNTL(1), SO THAT CNTL(1) WILL REMAIN
C UNALTERED.
      UU = MIN(CNTL(1),HALF)
      UU = MAX(UU,-HALF)
      DO 10 I = 1,N
        IW2(I) = 0
   10 CONTINUE
C IWPOS IS POINTER TO FIRST FREE POSITION FOR FACTORS IN IW.
C POSFAC IS POINTER FOR FACTORS IN A. AT EACH PASS THROUGH THE
C     MAJOR LOOP POSFAC INITIALLY POINTS TO THE FIRST FREE LOCATION
C     IN A AND THEN IS SET TO THE POSITION OF THE CURRENT PIVOT IN A.
C ISTK IS POINTER TO TOP OF STACK IN IW.
C ISTK2 IS POINTER TO BOTTOM OF STACK IN IW (NEEDED BY COMPRESS).
C ASTK IS POINTER TO TOP OF STACK IN A.
C ASTK2 IS POINTER TO BOTTOM OF STACK IN A (NEEDED BY COMPRESS).
C IINPUT IS POINTER TO CURRENT POSITION IN ORIGINAL ROWS IN IW.
C AINPUT IS POINTER TO CURRENT POSITION IN ORIGINAL ROWS IN A.
C AZERO IS POINTER TO LAST POSITION ZEROED IN A.
C NTOTPV IS THE TOTAL NUMBER OF PIVOTS SELECTED. THIS IS USED
C     TO DETERMINE WHETHER THE MATRIX IS SINGULAR.
      IWPOS = 2
      POSFAC = 1
      ISTK = LIW - NZ + 1
      ISTK2 = ISTK - 1
      ASTK = LA - NZ + 1
      ASTK2 = ASTK - 1
      IINPUT = ISTK
      AINPUT = ASTK
      AZERO = 0
      NTOTPV = 0
C NUMASS IS THE ACCUMULATED NUMBER OF ROWS ASSEMBLED SO FAR.
      NUMASS = 0
C
C EACH PASS THROUGH THIS MAIN LOOP PERFORMS ALL THE OPERATIONS
C     ASSOCIATED WITH ONE SET OF ASSEMBLY/ELIMINATIONS.
      DO 760 IASS = 1,NSTEPS
C NASS WILL BE SET TO THE NUMBER OF FULLY ASSEMBLED VARIABLES IN
C     CURRENT NEWLY CREATED ELEMENT.
        NASS = NELIM(IASS)
C NEWEL IS A POINTER INTO IW TO CONTROL OUTPUT OF INTEGER INFORMATION
C     FOR NEWLY CREATED ELEMENT.
        NEWEL = IWPOS + 1
C SYMBOLICALLY ASSEMBLE INCOMING ROWS AND GENERATED STACK ELEMENTS
C     ORDERING THE RESULTANT ELEMENT ACCORDING TO PERMUTATION PERM.  WE
C     ASSEMBLE THE STACK ELEMENTS FIRST BECAUSE THESE WILL ALREADY BE
C     ORDERED.
C SET HEADER POINTER FOR MERGE OF INDEX LISTS.
        JFIRST = N + 1
C INITIALIZE NUMBER OF VARIABLES IN CURRENT FRONT.
        NFRONT = 0
        NUMSTK = NSTK(IASS)
        LTOPST = 1
        LNASS = 0
C JUMP IF NO STACK ELEMENTS ARE BEING ASSEMBLED AT THIS STAGE.
        IF (NUMSTK.EQ.0) GO TO 80
        J2 = ISTK - 1
        LNASS = NASS
        LTOPST = ((IW(ISTK)+1)*IW(ISTK))/2
        DO 70 IELL = 1,NUMSTK
C ASSEMBLE ELEMENT IELL PLACING
C     THE INDICES INTO A LINKED LIST IN IW2 ORDERED
C     ACCORDING TO PERM.
          JNEXT = JFIRST
          JLAST = N + 1
          J1 = J2 + 2
          J2 = J1 - 1 + IW(J1-1)
C RUN THROUGH INDEX LIST OF STACK ELEMENT IELL.
          DO 60 JJ = J1,J2
            J = IW(JJ)
C JUMP IF ALREADY ASSEMBLED
            IF (IW2(J).GT.0) GO TO 60
            JNEW = PERM(J)
C IF VARIABLE WAS PREVIOUSLY FULLY SUMMED BUT WAS NOT PIVOTED ON
C     EARLIER BECAUSE OF NUMERICAL TEST, INCREMENT NUMBER OF FULLY
C     SUMMED ROWS/COLUMNS IN FRONT.
            IF (JNEW.LE.NUMASS) NASS = NASS + 1
C FIND POSITION IN LINKED LIST FOR NEW VARIABLE.  NOTE THAT WE START
C     FROM WHERE WE LEFT OFF AFTER ASSEMBLY OF PREVIOUS VARIABLE.
            DO 20 IDUMMY = 1,N
              IF (JNEXT.EQ.N+1) GO TO 30
              IF (PERM(JNEXT).GT.JNEW) GO TO 30
              JLAST = JNEXT
              JNEXT = IW2(JLAST)
   20       CONTINUE
   30       IF (JLAST.NE.N+1) GO TO 40
            JFIRST = J
            GO TO 50

   40       IW2(JLAST) = J
   50       IW2(J) = JNEXT
            JLAST = J
C INCREMENT NUMBER OF VARIABLES IN THE FRONT.
            NFRONT = NFRONT + 1
   60     CONTINUE
   70   CONTINUE
        LNASS = NASS - LNASS
C NOW INCORPORATE ORIGINAL ROWS.  NOTE THAT THE COLUMNS IN THESE ROWS
C     NEED NOT BE IN ORDER. WE ALSO PERFORM
C     A SWOP SO THAT THE DIAGONAL ENTRY IS THE FIRST IN ITS
C     ROW.  THIS ALLOWS US TO AVOID STORING THE INVERSE OF ARRAY PERM.
   80   NUMORG = NELIM(IASS)
        J1 = IINPUT
        DO 150 IORG = 1,NUMORG
          J = -IW(J1)
          DO 140 IDUMMY = 1,LIW
            JNEW = PERM(J)
C JUMP IF VARIABLE ALREADY INCLUDED.
            IF (IW2(J).GT.0) GO TO 130
C HERE WE MUST ALWAYS START OUR SEARCH AT THE BEGINNING.
            JLAST = N + 1
            JNEXT = JFIRST
            DO 90 JDUMMY = 1,N
              IF (JNEXT.EQ.N+1) GO TO 100
              IF (PERM(JNEXT).GT.JNEW) GO TO 100
              JLAST = JNEXT
              JNEXT = IW2(JLAST)
   90       CONTINUE
  100       IF (JLAST.NE.N+1) GO TO 110
            JFIRST = J
            GO TO 120

  110       IW2(JLAST) = J
  120       IW2(J) = JNEXT
C INCREMENT NUMBER OF VARIABLES IN FRONT.
            NFRONT = NFRONT + 1
  130       J1 = J1 + 1
            IF (J1.GT.LIW) GO TO 150
            J = IW(J1)
            IF (J.LT.0) GO TO 150
  140     CONTINUE
  150   CONTINUE
C NOW RUN THROUGH LINKED LIST IW2 PUTTING INDICES OF VARIABLES IN NEW
C     ELEMENT INTO IW AND SETTING IW2 ENTRY TO POINT TO THE RELATIVE
C     POSITION OF THE VARIABLE IN THE NEW ELEMENT.
        IF (NEWEL+NFRONT.LT.ISTK) GO TO 160
C COMPRESS IW.
        CALL MA27PD(A,IW,ISTK,ISTK2,IINPUT,2,NCMPBR,NCMPBI)
        IF (NEWEL+NFRONT.LT.ISTK) GO TO 160
        INFO(2) = LIW + 1 + NEWEL + NFRONT - ISTK
        GO TO 770

  160   J = JFIRST
        DO 170 IFR = 1,NFRONT
          NEWEL = NEWEL + 1
          IW(NEWEL) = J
          JNEXT = IW2(J)
          IW2(J) = NEWEL - (IWPOS+1)
          J = JNEXT
  170   CONTINUE
C
C ASSEMBLE REALS INTO FRONTAL MATRIX.
        MAXFRT = MAX(MAXFRT,NFRONT)
        IW(IWPOS) = NFRONT
C FIRST ZERO OUT FRONTAL MATRIX AS APPROPRIATE FIRST CHECKING TO SEE
C     IF THERE IS SUFFICIENT SPACE.
        LAELL = ((NFRONT+1)*NFRONT)/2
        APOS2 = POSFAC + LAELL - 1
        IF (NUMSTK.NE.0) LNASS = LNASS* (2*NFRONT-LNASS+1)/2
        IF (POSFAC+LNASS-1.GE.ASTK) GO TO 180
        IF (APOS2.LT.ASTK+LTOPST-1) GO TO 190
C COMPRESS A.
  180   CALL MA27PD(A,IW,ASTK,ASTK2,AINPUT,1,NCMPBR,NCMPBI)
        IF (POSFAC+LNASS-1.GE.ASTK) GO TO 780
        IF (APOS2.GE.ASTK+LTOPST-1) GO TO 780
  190   IF (APOS2.LE.AZERO) GO TO 220
        APOS = AZERO + 1
        LAPOS2 = MIN(APOS2,ASTK-1)
        IF (LAPOS2.LT.APOS) GO TO 210
        DO 200 K = APOS,LAPOS2
          A(K) = ZERO
  200   CONTINUE
  210   AZERO = APOS2
C JUMP IF THERE ARE NO STACK ELEMENTS TO ASSEMBLE.
  220   IF (NUMSTK.EQ.0) GO TO 260
C PLACE REALS CORRESPONDING TO STACK ELEMENTS IN CORRECT POSITIONS IN A.
        DO 250 IELL = 1,NUMSTK
          J1 = ISTK + 1
          J2 = ISTK + IW(ISTK)
          DO 240 JJ = J1,J2
            IROW = IW(JJ)
            IROW = IW2(IROW)
            IX = NFRONT
            IY = IROW
            IDIAG = ((IY-1)* (2*IX-IY+2))/2
            APOS = POSFAC + IDIAG
            DO 230 JJJ = JJ,J2
              J = IW(JJJ)
              APOS2 = APOS + IW2(J) - IROW
              A(APOS2) = A(APOS2) + A(ASTK)
              A(ASTK) = ZERO
              ASTK = ASTK + 1
  230       CONTINUE
  240     CONTINUE
C INCREMENT STACK POINTER.
          ISTK = J2 + 1
  250   CONTINUE
C INCORPORATE REALS FROM ORIGINAL ROWS.
  260   DO 280 IORG = 1,NUMORG
          J = -IW(IINPUT)
C WE CAN DO THIS BECAUSE THE DIAGONAL IS NOW THE FIRST ENTRY.
          IROW = IW2(J)
          IX = NFRONT
          IY = IROW
          IDIAG = ((IY-1)* (2*IX-IY+2))/2
          APOS = POSFAC + IDIAG
C THE FOLLOWING LOOP GOES FROM 1 TO NZ BECAUSE THERE MAY BE DUPLICATES.
          DO 270 IDUMMY = 1,NZ
            APOS2 = APOS + IW2(J) - IROW
            A(APOS2) = A(APOS2) + A(AINPUT)
            AINPUT = AINPUT + 1
            IINPUT = IINPUT + 1
            IF (IINPUT.GT.LIW) GO TO 280
            J = IW(IINPUT)
            IF (J.LT.0) GO TO 280
  270     CONTINUE
  280   CONTINUE
C RESET IW2 AND NUMASS.
        NUMASS = NUMASS + NUMORG
        J1 = IWPOS + 2
        J2 = IWPOS + NFRONT + 1
        DO 290 K = J1,J2
          J = IW(K)
          IW2(J) = 0
  290   CONTINUE
C PERFORM PIVOTING ON ASSEMBLED ELEMENT.
C NPIV IS THE NUMBER OF PIVOTS SO FAR SELECTED.
C LNPIV IS THE NUMBER OF PIVOTS SELECTED AFTER THE LAST PASS THROUGH
C     THE THE FOLLOWING LOOP.
        LNPIV = -1
        NPIV = 0
        DO 650 KDUMMY = 1,NASS
          IF (NPIV.EQ.NASS) GO TO 660
          IF (NPIV.EQ.LNPIV) GO TO 660
          LNPIV = NPIV
          NPIVP1 = NPIV + 1
C JPIV IS USED AS A FLAG TO INDICATE WHEN 2 BY 2 PIVOTING HAS OCCURRED
C     SO THAT IPIV IS INCREMENTED CORRECTLY.
          JPIV = 1
C NASS IS MAXIMUM POSSIBLE NUMBER OF PIVOTS.
C WE EITHER TAKE THE DIAGONAL ENTRY OR THE 2 BY 2 PIVOT WITH THE
C     LARGEST OFF-DIAGONAL AT EACH STAGE.
C EACH PASS THROUGH THIS LOOP TRIES TO CHOOSE ONE PIVOT.
          DO 640 IPIV = NPIVP1,NASS
            JPIV = JPIV - 1
C JUMP IF WE HAVE JUST PROCESSED A 2 BY 2 PIVOT.
            IF (JPIV.EQ.1) GO TO 640
            IX = NFRONT-NPIV
            IY = IPIV-NPIV
            IDIAG = ((IY-1)* (2*IX-IY+2))/2
            APOS = POSFAC + IDIAG
C IF THE USER HAS INDICATED THAT THE MATRIX IS DEFINITE, WE
C     DO NOT NEED TO TEST FOR STABILITY BUT WE DO CHECK TO SEE IF THE
C     PIVOT IS NON-ZERO OR HAS CHANGED SIGN.
C     IF IT IS ZERO, WE EXIT WITH AN ERROR. IF IT HAS CHANGED SIGN
C     AND U WAS SET NEGATIVE, THEN WE AGAIN EXIT IMMEDIATELY.  IF THE
C     PIVOT CHANGES SIGN AND U WAS ZERO, WE CONTINUE WITH THE
C     FACTORIZATION BUT PRINT A WARNING MESSAGE ON UNIT ICNTL(2).
C ISNPIV HOLDS A FLAG FOR THE SIGN OF THE PIVOTS TO DATE SO THAT
C     A SIGN CHANGE WHEN DECOMPOSING AN ALLEGEDLY DEFINITE MATRIX CAN
C     BE DETECTED.
            IF (UU.GT.ZERO) GO TO 320
            IF (ABS(A(APOS)).LE.CNTL(3)) GO TO 790
C JUMP IF THIS IS NOT THE FIRST PIVOT TO BE SELECTED.
            IF (NTOTPV.GT.0) GO TO 300
C SET ISNPIV.
            IF (A(APOS).GT.ZERO) ISNPIV = 1
            IF (A(APOS).LT.ZERO) ISNPIV = -1
  300       IF (A(APOS).GT.ZERO .AND. ISNPIV.EQ.1) GO TO 560
            IF (A(APOS).LT.ZERO .AND. ISNPIV.EQ.-1) GO TO 560
            IF (INFO(1).NE.2) INFO(2) = 0
            INFO(2) = INFO(2) + 1
            INFO(1) = 2
            I = NTOTPV + 1
            IF (ICNTL(2).GT.0 .AND. INFO(2).LE.10) THEN
              WRITE (ICNTL(2),FMT=310) INFO(1),I
            END IF

  310       FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MA27BD',
     +              '  *** INFO(1) =',I2,/,' PIVOT',I6,
     +             ' HAS DIFFERENT SIGN FROM THE PREVIOUS ONE')

            ISNPIV = -ISNPIV
            IF (UU.EQ.ZERO) GO TO 560
            GO TO 800

  320       AMAX = ZERO
            TMAX = AMAX
C FIND LARGEST ENTRY TO RIGHT OF DIAGONAL IN ROW OF PROSPECTIVE PIVOT
C     IN THE FULLY-SUMMED PART.  ALSO RECORD COLUMN OF THIS LARGEST
C     ENTRY.
            J1 = APOS + 1
            J2 = APOS + NASS - IPIV
            IF (J2.LT.J1) GO TO 340
            DO 330 JJ = J1,J2
              IF (ABS(A(JJ)).LE.AMAX) GO TO 330
              JMAX = IPIV + JJ - J1 + 1
              AMAX = ABS(A(JJ))
  330       CONTINUE
C DO SAME AS ABOVE FOR NON-FULLY-SUMMED PART ONLY HERE WE DO NOT NEED
C     TO RECORD COLUMN SO LOOP IS SIMPLER.
  340       J1 = J2 + 1
            J2 = APOS + NFRONT - IPIV
            IF (J2.LT.J1) GO TO 360
            DO 350 JJ = J1,J2
              TMAX = MAX(ABS(A(JJ)),TMAX)
  350       CONTINUE
C NOW CALCULATE LARGEST ENTRY IN OTHER PART OF ROW.
  360       RMAX = MAX(TMAX,AMAX)
            APOS1 = APOS
            KK = NFRONT - IPIV
            LT = IPIV - (NPIV+1)
            IF (LT.EQ.0) GO TO 380
            DO 370 K = 1,LT
              KK = KK + 1
              APOS1 = APOS1 - KK
              RMAX = MAX(RMAX,ABS(A(APOS1)))
  370       CONTINUE
C JUMP IF STABILITY TEST SATISFIED.
  380       IF (ABS(A(APOS)).GT.MAX(CNTL(3),UU*RMAX)) GO TO 450
C CHECK BLOCK PIVOT OF ORDER 2 FOR STABILITY.
            IF (ABS(AMAX).LE.CNTL(3)) GO TO 640
            IX = NFRONT-NPIV
            IY = JMAX-NPIV
            IDIAG = ((IY-1)* (2*IX-IY+2))/2
            APOS2 = POSFAC + IDIAG
            DETPIV = A(APOS)*A(APOS2) - AMAX*AMAX
            THRESH = ABS(DETPIV)
C SET THRESH TO U TIMES THE RECIPROCAL OF THE MAX-NORM OF THE INVERSE
C     OF THE PROSPECTIVE BLOCK.
            THRESH = THRESH/ (UU*MAX(ABS(A(APOS))+AMAX,
     +               ABS(A(APOS2))+AMAX))
C CHECK 2 BY 2 PIVOT FOR STABILITY.
C FIRST CHECK AGAINST ROW IPIV.
            IF (THRESH.LE.RMAX) GO TO 640
C FIND LARGEST ENTRY IN ROW JMAX.
C FIND MAXIMUM TO THE RIGHT OF THE DIAGONAL.
            RMAX = ZERO
            J1 = APOS2 + 1
            J2 = APOS2 + NFRONT - JMAX
            IF (J2.LT.J1) GO TO 400
            DO 390 JJ = J1,J2
              RMAX = MAX(RMAX,ABS(A(JJ)))
  390       CONTINUE
C NOW CHECK TO THE LEFT OF THE DIAGONAL.
C WE USE TWO LOOPS TO AVOID TESTING FOR ROW IPIV INSIDE THE LOOP.
  400       KK = NFRONT - JMAX + 1
            APOS3 = APOS2
            JMXMIP = JMAX - IPIV - 1
            IF (JMXMIP.EQ.0) GO TO 420
            DO 410 K = 1,JMXMIP
              APOS2 = APOS2 - KK
              KK = KK + 1
              RMAX = MAX(RMAX,ABS(A(APOS2)))
  410       CONTINUE
  420       IPMNP = IPIV - NPIV - 1
            IF (IPMNP.EQ.0) GO TO 440
            APOS2 = APOS2 - KK
            KK = KK + 1
            DO 430 K = 1,IPMNP
              APOS2 = APOS2 - KK
              KK = KK + 1
              RMAX = MAX(RMAX,ABS(A(APOS2)))
  430       CONTINUE
  440       IF (THRESH.LE.RMAX) GO TO 640
            PIVSIZ = 2
            GO TO 460

  450       PIVSIZ = 1
  460       IROW = IPIV - NPIV
C
C PIVOT HAS BEEN CHOSEN.  IF BLOCK PIVOT OF ORDER 2, PIVSIZ IS EQUAL TO
C     TWO OTHERWISE PIVSIZ EQUALS ONE..
C THE FOLLOWING LOOP MOVES THE PIVOT BLOCK TO THE TOP LEFT HAND CORNER
C     OF THE FRONTAL MATRIX.
            DO 550 KROW = 1,PIVSIZ
C WE JUMP IF SWOP IS NOT NECESSARY.
              IF (IROW.EQ.1) GO TO 530
              J1 = POSFAC + IROW
              J2 = POSFAC + NFRONT - (NPIV+1)
              IF (J2.LT.J1) GO TO 480
              APOS2 = APOS + 1
C SWOP PORTION OF ROWS WHOSE COLUMN INDICES ARE GREATER THAN LATER ROW.
              DO 470 JJ = J1,J2
                SWOP = A(APOS2)
                A(APOS2) = A(JJ)
                A(JJ) = SWOP
                APOS2 = APOS2 + 1
  470         CONTINUE
  480         J1 = POSFAC + 1
              J2 = POSFAC + IROW - 2
              APOS2 = APOS
              KK = NFRONT - (IROW+NPIV)
              IF (J2.LT.J1) GO TO 500
C SWOP PORTION OF ROWS/COLUMNS WHOSE INDICES LIE BETWEEN THE TWO ROWS.
              DO 490 JJJ = J1,J2
                JJ = J2 - JJJ + J1
                KK = KK + 1
                APOS2 = APOS2 - KK
                SWOP = A(APOS2)
                A(APOS2) = A(JJ)
                A(JJ) = SWOP
  490         CONTINUE
  500         IF (NPIV.EQ.0) GO TO 520
              APOS1 = POSFAC
              KK = KK + 1
              APOS2 = APOS2 - KK
C SWOP PORTION OF COLUMNS WHOSE INDICES ARE LESS THAN EARLIER ROW.
              DO 510 JJ = 1,NPIV
                KK = KK + 1
                APOS1 = APOS1 - KK
                APOS2 = APOS2 - KK
                SWOP = A(APOS2)
                A(APOS2) = A(APOS1)
                A(APOS1) = SWOP
  510         CONTINUE
C SWOP DIAGONALS AND INTEGER INDEXING INFORMATION
  520         SWOP = A(APOS)
              A(APOS) = A(POSFAC)
              A(POSFAC) = SWOP
              IPOS = IWPOS + NPIV + 2
              IEXCH = IWPOS + IROW + NPIV + 1
              ISWOP = IW(IPOS)
              IW(IPOS) = IW(IEXCH)
              IW(IEXCH) = ISWOP
  530         IF (PIVSIZ.EQ.1) GO TO 550
C SET VARIABLES FOR THE SWOP OF SECOND ROW OF BLOCK PIVOT.
              IF (KROW.EQ.2) GO TO 540
              IROW = JMAX - (NPIV+1)
              JPOS = POSFAC
              POSFAC = POSFAC + NFRONT - NPIV
              NPIV = NPIV + 1
              APOS = APOS3
              GO TO 550
C RESET VARIABLES PREVIOUSLY SET FOR SECOND PASS.
  540         NPIV = NPIV - 1
              POSFAC = JPOS
  550       CONTINUE
C
            IF (PIVSIZ.EQ.2) GO TO 600
C PERFORM THE ELIMINATION USING ENTRY (IPIV,IPIV) AS PIVOT.
C WE STORE U AND DINVERSE.
  560       A(POSFAC) = ONE/A(POSFAC)
            IF (A(POSFAC).LT.ZERO) NEIG = NEIG + 1
            J1 = POSFAC + 1
            J2 = POSFAC + NFRONT - (NPIV+1)
            IF (J2.LT.J1) GO TO 590
            IBEG = J2 + 1
            DO 580 JJ = J1,J2
              AMULT = -A(JJ)*A(POSFAC)
              IEND = IBEG + NFRONT - (NPIV+JJ-J1+2)
C THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1.
CDIR$ IVDEP
              DO 570 IROW = IBEG,IEND
                JCOL = JJ + IROW - IBEG
                A(IROW) = A(IROW) + AMULT*A(JCOL)
  570         CONTINUE
              IBEG = IEND + 1
              A(JJ) = AMULT
  580       CONTINUE
  590       NPIV = NPIV + 1
            NTOTPV = NTOTPV + 1
            JPIV = 1
            POSFAC = POSFAC + NFRONT - NPIV + 1
            GO TO 640
C PERFORM ELIMINATION USING BLOCK PIVOT OF ORDER TWO.
C REPLACE BLOCK PIVOT BY ITS INVERSE.
C SET FLAG TO INDICATE USE OF 2 BY 2 PIVOT IN IW.
  600       IPOS = IWPOS + NPIV + 2
            NTWO = NTWO + 1
            IW(IPOS) = -IW(IPOS)
            POSPV1 = POSFAC
            POSPV2 = POSFAC + NFRONT - NPIV
            SWOP = A(POSPV2)
            IF (DETPIV.LT.ZERO) NEIG = NEIG + 1
            IF (DETPIV.GT.ZERO .AND. SWOP.LT.ZERO) NEIG = NEIG + 2
            A(POSPV2) = A(POSPV1)/DETPIV
            A(POSPV1) = SWOP/DETPIV
            A(POSPV1+1) = -A(POSPV1+1)/DETPIV
            J1 = POSPV1 + 2
            J2 = POSPV1 + NFRONT - (NPIV+1)
            IF (J2.LT.J1) GO TO 630
            JJ1 = POSPV2
            IBEG = POSPV2 + NFRONT - (NPIV+1)
            DO 620 JJ = J1,J2
              JJ1 = JJ1 + 1
              AMULT1 = - (A(POSPV1)*A(JJ)+A(POSPV1+1)*A(JJ1))
              AMULT2 = - (A(POSPV1+1)*A(JJ)+A(POSPV2)*A(JJ1))
              IEND = IBEG + NFRONT - (NPIV+JJ-J1+3)
C THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1.
CDIR$ IVDEP
              DO 610 IROW = IBEG,IEND
                K1 = JJ + IROW - IBEG
                K2 = JJ1 + IROW - IBEG
                A(IROW) = A(IROW) + AMULT1*A(K1) + AMULT2*A(K2)
  610         CONTINUE
              IBEG = IEND + 1
              A(JJ) = AMULT1
              A(JJ1) = AMULT2
  620       CONTINUE
  630       NPIV = NPIV + 2
            NTOTPV = NTOTPV + 2
            JPIV = 2
            POSFAC = POSPV2 + NFRONT - NPIV + 1
  640     CONTINUE
  650   CONTINUE
C END OF MAIN ELIMINATION LOOP.
C
  660   IF (NPIV.NE.0) NBLK = NBLK + 1
        IOLDPS = IWPOS
        IWPOS = IWPOS + NFRONT + 2
        IF (NPIV.EQ.0) GO TO 690
        IF (NPIV.GT.1) GO TO 680
        IW(IOLDPS) = -IW(IOLDPS)
        DO 670 K = 1,NFRONT
          J1 = IOLDPS + K
          IW(J1) = IW(J1+1)
  670   CONTINUE
        IWPOS = IWPOS - 1
        GO TO 690

  680   IW(IOLDPS+1) = NPIV
C COPY REMAINDER OF ELEMENT TO TOP OF STACK
  690   LIELL = NFRONT - NPIV
        IF (LIELL.EQ.0 .OR. IASS.EQ.NSTEPS) GO TO 750
        IF (IWPOS+LIELL.LT.ISTK) GO TO 700
        CALL MA27PD(A,IW,ISTK,ISTK2,IINPUT,2,NCMPBR,NCMPBI)
  700   ISTK = ISTK - LIELL - 1
        IW(ISTK) = LIELL
        J1 = ISTK
        KK = IWPOS - LIELL - 1
C THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1.
CDIR$ IVDEP
        DO 710 K = 1,LIELL
          J1 = J1 + 1
          KK = KK + 1
          IW(J1) = IW(KK)