C ** Correction report.
C ** Correction 1. 05/09/92: 1 line corrected **
C ** Correction 2. 05/09/92: 2 lines added **
C ** Correction 3. 20/12/93: 8 lines added **
C ** Correction 4. 20/12/93: 7 lines added **
C ** Correction 5. 20/12/93: 14 lines added **
C ** Correction 6. 20/12/93: 4 lines added **
C ** Correction 7. 20/12/93: 2 lines added **
C ** Correction 8. 13/01/94: 1 line added **
C ** Correction 9. 19/08/96: 1 line added, 1 line corrected **
C ** Correction 10. 11/07/97: 1 line corrected **
C ** Correction 11. 05/11/97: 2 lines interchanged **
C ** End of Correction report.
C  THIS VERSION: 05/11/1997 AT 18:00:51 PM.
CS    SUBROUTINE SMCFA(N, NZ, IRN, ICN, A, LA, IW, LIW, IKEEP, NSTEPS,
CD    SUBROUTINE DMCFA(N, NZ, IRN, ICN, A, LA, IW, LIW, IKEEP, NSTEPS,
     * MAXFRT, IW1, IFLAG, DIAG, OMEGA, MLP)
C**
C MLP IS STREAM NUMBER FOR DIAGNOSTIC PRINTING, ZERO SUPPRESSES.
C THIS SUBROUTINE COMPUTES THE FACTORISATION OF THE MATRIX INPUT IN
C     A,IRN,ICN USING INFORMATION (IN IKEEP) FROM MA27A/AD.
C** DIAG HOLDS DIAGONALS OF MATRIX AND FACTORED MATRIX
C** OMEGA HOLDS CHANGES TO DIAGONALS (ABSOLUTE CHANGES)
CS    REAL             A(LA), DIAG(N), OMEGA(N), ADDON
CD    DOUBLE PRECISION A(LA), DIAG(N), OMEGA(N), ADDON
      INTEGER IW1(N)
C
      INTEGER     IRN(*),ICN(*),IW(LIW),IKEEP(N,3)
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 MCFA/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 MCFA/BD.
C IW NEED NOT BE SET ON ENTRY.  USED AS A WORK ARRAY BY MCFA/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 MCFA/BD.
C IKEEP MUST BE UNCHANGED SINCE THE CALL TO MA27A/AD.
C     IT IS NOT ALTERED BY MCFA/BD.
C NSTEPS MUST BE UNCHANGED SINCE THE CALL TO MA27A/AD.
C     IT IS NOT ALTERED BY MCFA/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 MCFA/BD.
C IFLAG NEED NOT BE SET ON ENTRY.  ON EXIT A ZERO VALUE INDICATES
C     NO ERROR HAS BEEN DETECTED.
C SEE BLOCK DATA FOR COMMENTS ON VARIABLES IN COMMON
C
C ** Correction 3. 20/12/93: 8 lines added **
      INTEGER K ,KZ , I , NZ1 , IPHASE , IFRLVL
      INTEGER NTWO , NCMPBI , NEIG , NEMIN , IOVFLO , J2
      INTEGER J1 , JJ , IROWS , LEN , NROWS , IPOS 
      INTEGER KBLK , IAPOS , NCOLS , IBLK , MLP , IFLAG 
      INTEGER LP , LDIAG , MP , MAXFRT , NZ , N 
      INTEGER LA , NSTEPS , LIW , NRLBDU , NIRADU , NIRBDU 
      INTEGER NCMPBR , NCMPA , NRLADU , NRLTOT , IERROR , NIRTOT 
      INTEGER NIRNEC , NRLNEC
C ** Correction 3. 20/12/93: end of correction **
CS    REAL             U, OPS
CD    DOUBLE PRECISION U, OPS
CS    COMMON /MA27D / U, LP, MP, LDIAG
CD    COMMON /MA27DD/ U, LP, MP, LDIAG
CS    COMMON /MA27E / OPS, IERROR, NRLTOT, NIRTOT, NRLNEC, NIRNEC,
CD    COMMON /MA27ED/ OPS, IERROR, NRLTOT, NIRTOT, NRLNEC, NIRNEC,
     * NRLADU, NIRADU, NRLBDU, NIRBDU, NCMPA, NCMPBR, NCMPBI, NTWO,NEIG
CS    COMMON /MA27F / IOVFLO, NEMIN, IFRLVL(20)
CD    COMMON /MA27FD/ IOVFLO, NEMIN, IFRLVL(20)
CS    EXTERNAL SMCFB, SMCFC 
CD    EXTERNAL DMCFB, DMCFC
      IF (LDIAG.LE.0 .OR. MP.LE.0) GO TO 60
C PRINT INPUT PARAMETERS.
      WRITE (MP,10) N, NZ, LA, LIW, NSTEPS, U
   10 FORMAT (//, ' ENTERING MCFAD WITH      N     NZ     LA    LIW ',
     * ' NSTEPS      U', / , 21X, 5I7, F7.2)
      KZ = MIN0(6,NZ)
      IF (LDIAG.GT.1) KZ = NZ
      IF (NZ.GT.0) WRITE (MP,20) (A(K),IRN(K),ICN(K),K=1,KZ)
   20 FORMAT (' MATRIX NON-ZEROS', 2(D16.4, 2I6), /,
     *      (17X, D16.4, 2I6, D16.4, 2I6))
      K = MIN0(9,N)
      IF (LDIAG.GT.1) K = N
      IF (K.GT.0) WRITE (MP,30) (IKEEP(I,1),I=1,K)
   30 FORMAT (' IKEEP(.,1)=', 10I6, /, (12X, 10I6))
      K = MIN0(K,NSTEPS)
      IF (K.GT.0) WRITE (MP,40) (IKEEP(I,2),I=1,K)
   40 FORMAT (' IKEEP(.,2)=', 10I6, /, (12X, 10I6))
      IF (K.GT.0) WRITE (MP,50) (IKEEP(I,3),I=1,K)
   50 FORMAT (' IKEEP(.,3)=', 10I6, /, (12X, 10I6))
   60 IF (N.LT.1 .OR. N.GT.IOVFLO) 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
C** SET PHASE OF CHOLESKY MODIFICATION
      IPHASE = 1
C SORT
CS    CALL SMCFB(N, NZ, NZ1, A, LA, IRN, ICN, IW, LIW, IKEEP, IW1,
CD    CALL DMCFB(N, NZ, NZ1, A, LA, IRN, ICN, IW, LIW, IKEEP, IW1,
     * IFLAG, DIAG, ADDON, MLP)
      IF ( MLP .GT. 0 ) WRITE( MLP, * ) ' ADDON = ', ADDON
      IF (IFLAG.EQ.-3) GO TO 130
      IF (IFLAG.EQ.-4) GO TO 160
C FACTORIZE
CS    CALL SMCFC(N, NZ1, A, LA, IW, LIW, IKEEP, IKEEP(1,3), NSTEPS,
CD    CALL DMCFC(N, NZ1, A, LA, IW, LIW, IKEEP, IKEEP(1,3), NSTEPS,
     * MAXFRT, IKEEP(1,2), IW1, IFLAG, DIAG, OMEGA, ADDON, IPHASE, MLP)
      IF (IFLAG.EQ.-3) GO TO 130
      IF (IFLAG.EQ.-4) GO TO 160
      IF (IFLAG.EQ.-5) GO TO 180
      IF (IFLAG.EQ.-6) GO TO 200
C **** WARNING MESSAGE ****
      IF (IFLAG.EQ.3 .AND. MP.GT.0) WRITE (MP,65) IFLAG, IERROR
   65 FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MCFA *** IFLA',
     * 'G =', I2, /, 5X,'MATRIX IS SINGULAR. RANK=',I5)
      GO TO 220
C **** ERROR RETURNS ****
   70 IFLAG = -1
      IF (LP.GT.0) WRITE (LP,80) IFLAG
   80 FORMAT (' **** ERROR RETURN FROM MCFAD **** IFLAG=', I3)
      IF (LP.GT.0) WRITE (LP,90) N
   90 FORMAT (' VALUE OF N OUT OF RANGE ... =', I10)
      GO TO 220
  100 IFLAG = -2
      IF (LP.GT.0) WRITE (LP,80) IFLAG
      IF (LP.GT.0) WRITE (LP,110) NZ
  110 FORMAT (' VALUE OF NZ OUT OF RANGE .. =', I10)
      GO TO 220
  120 IFLAG = -3
      IERROR = NZ
  130 IF (LP.GT.0) WRITE (LP,80) IFLAG
      IF (LP.GT.0) WRITE (LP,140) LIW, IERROR
  140 FORMAT (' LIW TOO SMALL, MUST BE INCREASED FROM', I10, ' TO',
     * ' AT LEAST', I10)
      GO TO 220
  150 IFLAG = -4
      IERROR = NZ + N
  160 IF (LP.GT.0) WRITE (LP,80) IFLAG
      IF (LP.GT.0) WRITE (LP,170) LA, IERROR
  170 FORMAT (' LA TOO SMALL, MUST BE INCREASED FROM ', I10, ' TO',
     * ' AT LEAST', I10)
      GO TO 220
  180 IF (LP.GT.0) WRITE (LP,80) IFLAG
      IF (LP.GT.0) WRITE (LP,190) IERROR
  190 FORMAT (' ZERO PIVOT AT STAGE', I10, ' WHEN INPUT MATRIX DECLA',
     * 'RED DEFINITE')
      GO TO 220
  200 IF (LP.GT.0) WRITE (LP,80) IFLAG
      IF (LP.GT.0) WRITE (LP,210)
  210 FORMAT (' CHANGE IN SIGN OF PIVOT ENCOUNTERED ',
     *   'WHEN FACTORING ALLEGEDLY DEFINITE MATRIX')
  220 IF (LDIAG.LE.0 .OR. MP.LE.0) GO TO 310
C PRINT OUTPUT PARAMETERS.
      WRITE (MP,230) MAXFRT, IFLAG, NRLBDU, NIRBDU, NCMPBR, NCMPBI,
     * NTWO, IERROR
  230 FORMAT (/' LEAVING MCFAD WITH MAXFRT  IFLAG NRLBDU NIRBDU NCMP',
     * 'BR NCMPBI   NTWO IERROR',/, 20X, 8I7)
      IF (IFLAG.LT.0) GO TO 310
C PRINT OUT MATRIX FACTORS FROM MCFA/BD.
      KBLK = IABS(IW(1)+0)
      IF (KBLK.EQ.0) GO TO 310
      IF (LDIAG.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 (MP,250) IBLK, NROWS, NCOLS
  250   FORMAT (' BLOCK PIVOT =', I8, ' NROWS =', I8, ' NCOLS =', I8)
        J2 = J1 + NCOLS - 1
        IPOS = J2 + 1
        WRITE (MP,260) (IW(JJ),JJ=J1,J2)
  260   FORMAT (' COLUMN INDICES =', 10I6, /, (17X, 10I6))
        WRITE (MP,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 (MP,280) (A(JJ),JJ=J1,J2)
  280     FORMAT (1P,5D16.8)
          LEN = LEN - 1
          IAPOS = J2 + 1
  290   CONTINUE
  300 CONTINUE
  310 RETURN
      END
CS    SUBROUTINE SMCFB(N, NZ, NZ1, A, LA, IRN, ICN, IW, LIW, PERM,
CD    SUBROUTINE DMCFB(N, NZ, NZ1, A, LA, IRN, ICN, IW, LIW, PERM,
     * IW2, IFLAG, DIAG, ADDON, MLP)
C
C SORT PRIOR TO FACTORIZATION USING MCFC.
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.
CS    REAL             A(LA), ZERO, ANEXT, ANOW, DIAG(N)
CD    DOUBLE PRECISION A(LA), ZERO, ANEXT, ANOW, DIAG(N)
      INTEGER IW2(N)
C
      INTEGER     IRN(*),ICN(*),IW(LIW),PERM(N)
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+MAX0(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 IFLAG  - ON EXIT FROM MA27N/ND, A ZERO VALUE OF
C          IFLAG 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 IERROR OF MA27E/ED AND
C              MESSAGES IDENTIFYING THE FIRST TEN ARE OUTPUT ON UNIT MP.
C          -3  INTEGER ARRAY IW IS TOO SMALL.
C          -4  REAL ARRAY A IS TOO SMALL.
C
C ** Correction 4. 20/12/93: 7 lines added **
      INTEGER K , IOLD , INEW , JOLD , IA , NCMPBI 
      INTEGER NCMPBR , NEIG , NTWO , JNEW , J2 , J1 
      INTEGER IIW , JJ , II , ICH , I , IPOS 
      INTEGER JPOS , NCMPA , MLP , IFLAG , MP , LP 
      INTEGER LIW , NZ , N , LA , NZ1 , LDIAG 
      INTEGER NIRADU , NRLADU , NIRBDU , NRLBDU , NIRNEC , NRLTOT
      INTEGER IERROR , NRLNEC , NIRTOT
C ** Correction 4. 20/12/93: end of correction **
CS    REAL             U, OPS, MACHEP, ADDON, MAXDAG, SMACHR
CD    DOUBLE PRECISION U, OPS, MACHEP, ADDON, MAXDAG, DMACHR
C SEE BLOCK DATA FOR COMMENTS ON VARIABLES IN COMMON
CS    COMMON /MA27D / U, LP, MP, LDIAG
CD    COMMON /MA27DD/ U, LP, MP, LDIAG
CS    COMMON /MA27E / OPS, IERROR, NRLTOT, NIRTOT, NRLNEC, NIRNEC,
CD    COMMON /MA27ED/ OPS, IERROR, NRLTOT, NIRTOT, NRLNEC, NIRNEC,
     * NRLADU, NIRADU, NRLBDU, NIRBDU, NCMPA, NCMPBR, NCMPBI, NTWO,NEIG
CS    EXTERNAL SMACHR
CD    EXTERNAL DMACHR
CS    DATA ZERO /0.0E0/
CD    DATA ZERO /0.0D0/
C** OBTAIN MACHEP
CS    MACHEP = SMACHR(1)
CD    MACHEP = DMACHR(1)
      IFLAG = 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.
      IERROR = 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 = MIN0(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   IFLAG = 1
        IERROR = IERROR + 1
        IF (IERROR.LE.1 .AND. MP.GT.0) WRITE (MP,40) IFLAG
   40   FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MCFA *** ',
     *   'IFLAG =', I2)
        IF (IERROR.LE.10 .AND. MP.GT.0) WRITE (MP,50) K, IRN(K), ICN(K)
   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 = MIN0(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 MAXDAG = MACHEP
      DO 190 IOLD=1,N
        INEW = PERM(IOLD)
        JPOS = IW2(INEW) - 1
        IA = LA - N + IOLD
        A(JPOS) = A(IA)
C**  SET DIAG TO VALUE OF DIAGONAL ENTRY (ORIGINAL NUMBERING)
        DIAG(IOLD) = A(IA)
        MAXDAG = MAX(MAXDAG,ABS(A(IA)))
        IW(JPOS) = -IOLD
  190 CONTINUE
C**
C COMPUTE ADDITION TO OFF-DIAGONAL 1-NORM
      ADDON = MAXDAG*MACHEP**0.75
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 IFLAG = -3
      IERROR = NZ1
      GO TO 230
  220 IFLAG = -4
      IERROR = NZ1 + N
C
  230 RETURN
      END
CS    SUBROUTINE SMCFC(N, NZ, A, LA, IW, LIW, PERM, NSTK, NSTEPS,
CD    SUBROUTINE DMCFC(N, NZ, A, LA, IW, LIW, PERM, NSTK, NSTEPS,
     * MAXFRT, NELIM, IW2, IFLAG, DIAG, OMEGA, ADDON, IPHASE, MLP)
C
C FACTORIZATION SUBROUTINE
C
C THIS SUBROUTINE OPERATES ON THE INPUT MATRIX ORDERED BY MFCFB 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 ** Correction 5. 20/12/93: 14 lines added **
      INTEGER IDIAG , IDUMMY , NUMORG , JNEW , JJ , J 
      INTEGER LAELL , LAPOS2 , IFR , IORG , JDUMMY , LTOPST 
      INTEGER LNASS , NUMSTK , JFIRST , NFRONT , JLAST , J1
      INTEGER JNEXT , J2 , IELL , ISWOP , IBEG , IEXCH 
      INTEGER KROW , IPOS , LIELL , KMAX , IOLDPS , IEND 
      INTEGER JCOL , NPIV , KDUMMY , LNPIV , IROW , JJJ 
      INTEGER JAY , KK , IPIV , NPIVP1 , JPIV , NEWEL
      INTEGER NIRTOT , NRLNEC , NRLTOT , LDIAG , IERROR , NRLBDU
      INTEGER NIRBDU , NIRADU , NIRNEC , NRLADU , LIW , MAXFRT
      INTEGER LA , N , NZ , LP , MP , MLP
      INTEGER IFLAG , IPHASE , ISTK , ISTK2 , IWPOS , K
      INTEGER NBLK , IASS , NASS , NUMASS , IINPUT , NTOTPV
      INTEGER NTWO , NEIG , NCMPBI , NCMPA , NCMPBR
      INTEGER I , IX , IY 
C ** Correction 5. 20/12/93: end of correction **
CS    REAL             A(LA), ZERO, AMAX, RMAX, SWOP,
CD    DOUBLE PRECISION A(LA), ZERO, AMAX, RMAX, SWOP,
     * ONE, AMULT, HALF, DIAG(N), OMEGA(N),
     * W1, ONENRM, ADDON
      INTEGER NSTEPS, POSFAC, ASTK, ASTK2, APOS, APOS1, APOS2,
     * AINPUT, PIVSIZ, AZERO, IW2(N)
C
      INTEGER     IW(LIW),PERM(N),NSTK(NSTEPS),NELIM(NSTEPS)
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 MCFA/AD.
C IFLAG  - 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
CS    REAL             U, UU, OPS
CD    DOUBLE PRECISION U, UU, OPS
CS    EXTERNAL SMCFD 
CD    EXTERNAL DMCFD
C SEE BLOCK DATA FOR COMMENTS ON VARIABLES IN COMMON
CS    COMMON /MA27D / U, LP, MP, LDIAG
CD    COMMON /MA27DD/ U, LP, MP, LDIAG
CS    COMMON /MA27E / OPS, IERROR, NRLTOT, NIRTOT, NRLNEC, NIRNEC,
CD    COMMON /MA27ED/ OPS, IERROR, NRLTOT, NIRTOT, NRLNEC, NIRNEC,
     * NRLADU, NIRADU, NRLBDU, NIRBDU, NCMPA, NCMPBR, NCMPBI, NTWO,NEIG
CS    DATA ZERO /0.0E0/, HALF /0.5E0/, ONE /1.0E0/
CD    DATA ZERO /0.0D0/, HALF /0.5D0/, ONE /1.0D0/
C THE FOLLOWING ARITHMETIC FUNCTION GIVES THE DISPLACEMENT FROM
C     THE START OF THE ASSEMBLED MATRIX(OF ORDER IX) OF THE DIAGONAL
C     ENTRY IN ITS ROW IY.
      IDIAG(IX,IY) = ((IY-1)*(2*IX-IY+2))/2
C***
      I = MIN( N, 4 )
      IF (MLP.GT.0) WRITE(MLP,6666) (DIAG(K),K=1,I)
 6666 FORMAT( ' DIAG ', 1P, 2D24.16, /, '      ', 2D24.16 )
      IF (MLP.GT.0) WRITE(MLP,*) 'IPHASE',IPHASE
C INITIALIZATION.
C NBLK IS THE NUMBER OF BLOCK PIVOTS USED.
      NBLK = 0
      NTWO = 0
      NEIG = 0
      NCMPBI = 0
      NCMPBR = 0
      MAXFRT = 0
C A PRIVATE VARIABLE UU IS SET TO U, SO THAT U WILL REMAIN UNALTERED.
      UU = MIN(U,HALF)
      UU = MAX(UU,-HALF)
C** OMEGA IS ARRAY OF ADJUSTMENTS TO DIAGONAL
      DO 10 I=1,N
        IW2(I) = 0
        OMEGA(I) = 0.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.
CS      CALL SMCFD(A, IW, ISTK, ISTK2, IINPUT, 2)
CD      CALL DMCFD(A, IW, ISTK, ISTK2, IINPUT, 2)
        IF (NEWEL+NFRONT.LT.ISTK) GO TO 160
        IERROR = 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 = MAX0(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   CONTINUE
CS      CALL SMCFD(A, IW, ASTK, ASTK2, AINPUT, 1)
CD      CALL DMCFD(A, IW, ASTK, ASTK2, AINPUT, 1)
        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 = MIN0(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)
            APOS = POSFAC + IDIAG(NFRONT,IROW)
            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)
          APOS = POSFAC + IDIAG(NFRONT,IROW)
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
C       WRITE(6,3333) IWPOS, ( IW(IWPOS+J-1),J=1,NFRONT+2 )
C3333   FORMAT( ' IWPOS = ', I6, ' IW(IWPOS),...,IW(IWPOS+NFONT+1) =',
C    *          /, ( 12I6 ) )
        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
            APOS = POSFAC + IDIAG(NFRONT-NPIV,IPIV-NPIV)
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 MP.
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.
C******
C           IF (UU.GT.ZERO) GO TO 320
C           IF (A(APOS).EQ.ZERO) GO TO 790
C JUMP IF THIS IS NOT THE FIRST PIVOT TO BE SELECTED.
C           IF (NTOTPV.GT.0) GO TO 300
C SET ISNPIV.
C           IF (A(APOS).GT.ZERO) ISNPIV = 1
C           IF (A(APOS).LT.ZERO) ISNPIV = -1
C 300       IF (A(APOS).GT.ZERO .AND. ISNPIV.EQ.1) GO TO 560
C           IF (A(APOS).LT.ZERO .AND. ISNPIV.EQ.-1) GO TO 560
C           IF (IFLAG.NE.2) IERROR = 0
C           IERROR = IERROR + 1
C           IFLAG = 2
C           I = NTOTPV + 1
C           IF (MP.GT.0 .AND. IERROR.LE.10) WRITE (MP,310) IFLAG, I
C 310       FORMAT (' *** WARNING MESSAGE FROM SUBROUTINE MCFA *** ',
C    *       'IFLAG =', I2,/,' PIVOT', I6,' HAS DIFFERENT SIGN FROM ',
C    *       'THE PREVIOUS ONE')
C           ISNPIV = -ISNPIV
C           IF (UU.EQ.ZERO) GO TO 560
C           GO TO 800
C** FIRST CHECK IF PIVOT IS POSITIVE
            IF (A(APOS).LE.ZERO) THEN
              IF (MLP.GT.0)
     *        WRITE(MLP,*) ' NEGATIVE PIVOT ENCOUNTERED AT STAGE',NPIV
              IPHASE = 2
            ENDIF
C**
            AMAX = ZERO
            RMAX = AMAX
C** I IS PIVOT ROW
            I = IW(IWPOS+1+IPIV)
            IF (MLP.GT.0) WRITE(MLP,*) 'PIVOT, PIVOT ROW', NPIV,I
            IF (MLP.GT.0) WRITE(MLP,6667) A(APOS)
 6667       FORMAT( ' PIVOT HAS VALUE ', 1P, D24.16 )
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.
C ONENRM IS SET TO ONENORM OF OFF-DIAGONALS IN ROW.
            ONENRM = 0.0
            J1 = APOS + 1
            J2 = APOS + NFRONT - IPIV
            IF (J2.LT.J1) GO TO 360
            DO 330 JJ=J1,J2
C ** Correction 11. 05/11/97: 2 lines interchanged **
              JAY = JJ-J1+1
              IF (IPHASE.EQ.1) THEN
C ** Correction 11. 05/11/97: end of correction **
                J = IW(IWPOS+1+IPIV+JAY)
                W1 = DIAG(J) - A(JJ)*A(JJ)/A(APOS)
C               IF (W1.LE.0.0) THEN
C               IF (W1.LE.1.0D-15) THEN
                IF (W1.LE.ADDON) THEN
                  IF (MLP.GT.0) WRITE(MLP,6668) J, W1
 6668       FORMAT( ' PHASE 2, J, W1 ', I6, 1P, D24.16 )
                  IPHASE = 2
                ENDIF
              ENDIF
              ONENRM = ONENRM + ABS(A(JJ))
              RMAX   = MAX(ABS(A(JJ)),RMAX)
              IF (JAY.LE.NASS-IPIV) THEN
                IF (ABS(A(JJ)).GT.AMAX) THEN
                  AMAX = ABS(A(JJ))
                ENDIF
              ENDIF
  330       CONTINUE
C NOW CALCULATE LARGEST ENTRY IN OTHER PART OF ROW.
  360       APOS1 = APOS
            KK = NFRONT - IPIV
            IF (MLP.GT.0) WRITE(MLP,*) 'NPIV,I,ONENRM', NPIV,I,ONENRM
C**  JUMP IF STILL IN PHASE 1.
            IF (IPHASE.EQ.1) GO TO 450
C**
C CHECK TO SEE IF PIVOT MUST BE INCREASED
            IF (A(APOS).GE.ADDON+ONENRM) GO TO 450
C**
C ADJUST DIAGONAL ENTRY AND RECORD CHANGE IN OMEGA
            OMEGA(I) = ADDON + ONENRM - A(APOS)
            IF (MLP.GT.0) WRITE(MLP,*) 'I,OMEGA(I)', I,OMEGA(I)
            A(APOS)  = ONENRM + ADDON
  450       PIVSIZ = 1
            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 550
              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
  550       CONTINUE
C
C PERFORM THE ELIMINATION USING ENTRY (IPIV,IPIV) AS PIVOT.
C WE STORE U AND DINVERSE.
C 560       CONTINUE
            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)
C**
C UPDATE DIAG ARRAY
                J = IW(IWPOS+2+NPIV+JJ-J1+1)
                DIAG(J) = DIAG(J) + AMULT*A(JJ)
                IF (MLP.GT.0) WRITE(MLP,*) 'J,DIAG(J)', J,DIAG(J)
              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
  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
CS      CALL SMCFD(A, IW, ISTK, ISTK2, IINPUT, 2)
CD      CALL DMCFD(A, IW, ISTK, ISTK2, IINPUT, 2)
  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)
  710   CONTINUE
C WE COPY IN REVERSE DIRECTION TO AVOID OVERWRITE PROBLEMS.
        LAELL = ((LIELL+1)*LIELL)/2
        KK = POSFAC + LAELL
        IF (KK.NE.ASTK) GO TO 720
        ASTK = ASTK - LAELL
        GO TO 740
C THE MOVE AND ZEROING OF ARRAY A IS PERFORMED WITH TWO LOOPS SO
C THAT THE CRAY-1 WILL VECTORIZE THEM SAFELY.
  720   KMAX = KK - 1
C THE FOLLOWING SPECIAL COMMENT FORCES VECTORIZATION ON THE CRAY-1.
CDIR$ IVDEP
        DO 730 K=1,LAELL
          KK = KK - 1
          ASTK = ASTK - 1
          A(ASTK) = A(KK)
  730   CONTINUE
        KMAX = MIN0(KMAX,ASTK-1)
        DO 735 K=KK,KMAX
          A(K) = ZERO
  735   CONTINUE
  740   AZERO = MIN0(AZERO,ASTK-1)
  750   IF (NPIV.EQ.0) IWPOS = IOLDPS
  760 CONTINUE
C
C END OF LOOP ON TREE NODES.
C
      IW(1) = NBLK
      IF (NTWO.GT.0) IW(1) = -NBLK
      NRLBDU = POSFAC - 1
      NIRBDU = IWPOS - 1
      IF (NTOTPV.EQ.N) GO TO 810
      IFLAG = 3
      IERROR = NTOTPV
      GO TO 810
C **** ERROR RETURNS ****
  770 IFLAG = -3
      GO TO 810
  780 IFLAG = -4
      IERROR = LA + MAX0(POSFAC+LNASS,APOS2-LTOPST+2) - ASTK
      GO TO 810
C 790 IFLAG = -5
C     IERROR = NTOTPV + 1
C     GO TO 810
C 800 IFLAG = -6
C     IERROR = NTOTPV + 1
  810 RETURN
      END
CS    SUBROUTINE SMCFD(A, IW, J1, J2, ITOP, IREAL)
CD    SUBROUTINE DMCFD(A, IW, J1, J2, ITOP, IREAL)
C THIS SUBROUTINE PERFORMS A VERY SIMPLE COMPRESS (BLOCK MOVE).
C     ENTRIES J1 TO J2 (INCL.) IN A OR IW AS APPROPRIATE ARE MOVED TO
C     OCCUPY THE POSITIONS IMMEDIATELY PRIOR TO POSITION ITOP.
C A/IW HOLD THE ARRAY BEING COMPRESSED.
C J1/J2 DEFINE THE ENTRIES BEING MOVED.
C ITOP DEFINES THE POSITION IMMEDIATELY AFTER THE POSITIONS TO WHICH
C     J1 TO J2 ARE MOVED.
C IREAL MUST BE SET BY THE USER TO 2 IF THE MOVE IS ON ARRAY IW,
C     ANY OTHER VALUE WILL PERFORM THE MOVE ON A.
C ** Correction 6. 20/12/93: 4 lines added **
      INTEGER NCMPBR , NCMPBI , NCMPA , NRLBDU , NIRBDU , JJJ
      INTEGER JJ , IPOS , NTWO , NEIG , NIRADU , IREAL
      INTEGER IERROR , ITOP , J1 , J2 , NIRNEC , NRLADU
      INTEGER NRLNEC , NRLTOT , NIRTOT
C ** Correction 6. 20/12/93: end of correction **
CS    REAL             A(*)
CD    DOUBLE PRECISION A(*)
C
      INTEGER     IW(*)
C
CS    REAL             OPS
CD    DOUBLE PRECISION OPS
C SEE BLOCK DATA FOR COMMENTS ON VARIABLES IN COMMON
CS    COMMON /MA27E / OPS, IERROR, NRLTOT, NIRTOT, NRLNEC, NIRNEC,
CD    COMMON /MA27ED/ OPS, IERROR, NRLTOT, NIRTOT, NRLNEC, NIRNEC,
     * NRLADU, NIRADU, NRLBDU, NIRBDU, NCMPA, NCMPBR, NCMPBI, NTWO,NEIG
      IPOS = ITOP - 1
      IF (J2.EQ.IPOS) GO TO 50
      IF (IREAL.EQ.2) GO TO 20
      NCMPBR = NCMPBR + 1
      IF (J1.GT.J2) GO TO 40
      DO 10 JJJ=J1,J2
        JJ = J2 - JJJ + J1
        A(IPOS) = A(JJ)
        IPOS = IPOS - 1
   10 CONTINUE
      GO TO 40
   20 NCMPBI = NCMPBI + 1
      IF (J1.GT.J2) GO TO 40
      DO 30 JJJ=J1,J2
        JJ = J2 - JJJ + J1
        IW(IPOS) = IW(JJ)
        IPOS = IPOS - 1
   30 CONTINUE
   40 J2 = ITOP - 1
      J1 = IPOS + 1
   50 RETURN
      END
CS    SUBROUTINE SSYPRC( N, LA, LIW, A, IW, NEG1, NEG2 )
CD    SUBROUTINE DSYPRC( N, LA, LIW, A, IW, NEG1, NEG2 )
      INTEGER          N, LA, LIW, NEG1, NEG2
      INTEGER          IW( LIW )
CS    REAL             A( LA )
CD    DOUBLE PRECISION A( LA )
C
C   THE GILL-MURRAY-PONCELEON-SAUNDERS CODE FOR MODIFYING THE NEGATIVE
C   EIGEN-COMPONENTS OBTAINED WHEN FACTORIZING A SYMMETRIC INDEFINITE 
C   MATRIX USING THE HARWELL CODE MA27. (SEE SOL 90-8, P.19-21)
C
C   NICK GOULD, 20TH JULY 1990.
C   FOR CGT PRODUCTIONS.
C
      INTEGER          ALEN, APOS, IBLK, NBLK, IPOS, NROWS, NCOLS, J, K
CS    REAL             ZERO, ONE, TWO, ALPHA, BETA, GAMMA, TAU, T,
CD    DOUBLE PRECISION ZERO, ONE, TWO, ALPHA, BETA, GAMMA, TAU, T,
     *                 C, S, E1, E2 
      LOGICAL          SINGLE
      INTRINSIC        ABS, SQRT
C
C  COMMON BLOCKS.
C
CS    REAL             EPSMCH, EPSNEG, TINY, BIG
CD    DOUBLE PRECISION EPSMCH, EPSNEG, TINY, BIG
CS    COMMON / SMACHN / EPSMCH, EPSNEG, TINY, BIG
CD    COMMON / DMACHN / EPSMCH, EPSNEG, TINY, BIG
C ** Correction 8. 13/01/94: 2 lines added
CS    SAVE   / SMACHN /
CD    SAVE   / DMACHN /
C ** Correction 8. 13/01/94: end of correction **
C
C  SET DATA.
C
CS    DATA             ZERO, ONE, TWO / 0.0E+0, 1.0E+0, 2.0E+0 /
CD    DATA             ZERO, ONE, TWO / 0.0D+0, 1.0D+0, 2.0D+0 /
C
C  NEG1 AND NEG2 ARE THE NUMBER OF NEGATIVE EIGENVALUES WHICH ARISE
C  FROM NEGATIVE 1X1 AND 2X2 BLOCK PIVOTS.
C
      NEG1 = 0
      NEG2 = 0
      NBLK = ABS( IW( 1 ) )
      IPOS = 2
      APOS = 1
C
C  LOOP OVER ALL THE BLOCK PIVOTS.
C
      DO 100 IBLK = 1, NBLK
         NCOLS    = IW( IPOS )
         IF ( NCOLS .LT. 0 ) THEN
            NROWS = 1
            NCOLS = - NCOLS
         ELSE
            IPOS  = IPOS + 1
            NROWS = IW( IPOS )
         END IF
C
C  PROCESS THE DIAGONALS IN THIS BLOCK.
C
         ALEN   = NCOLS
         SINGLE = .TRUE.
         DO 50 K = IPOS + 1, IPOS + NROWS
            IF ( SINGLE ) THEN
               ALPHA = A( APOS )
               J     = IW( K )
               SINGLE = J .GT. 0
               IF ( SINGLE ) THEN
C
C  NEGATIVE 1X1 BLOCK.
C
                  IF ( ALPHA .LT. ZERO ) THEN
                     NEG1 = NEG1 + 1
                     A( APOS ) = - ALPHA
C     write(6,*) ' 1x1 block modified. eval  = ', alpha
                  ELSE 
                     IF ( ALPHA .GT. ONE / EPSMCH ) THEN
                        NEG1 = NEG1 + 1
                        A( APOS ) = ONE / EPSMCH
C     write(6,*) ' 1x1 block modified. eval  = ', alpha
                     END IF
                  END IF
               ELSE
                  BETA  = A( APOS + 1 )
                  GAMMA = A( APOS + ALEN )
C
C  2X2 BLOCK: ( ALPHA  BETA  ) = ( C  S ) ( E1    ) ( C  S )
C             ( BETA   GAMMA )   ( S -C ) (    E2 ) ( S -C )
C 
                  IF ( ALPHA * GAMMA .LT. BETA ** 2 ) THEN
                     TAU = ( GAMMA - ALPHA ) / ( TWO * BETA )
                     T   = - ONE / ( ABS( TAU ) +
     *                               SQRT( ONE + TAU ** 2 ) )
                     IF ( TAU .LT. ZERO ) T = - T
C ** Correction 10. 11/07/97: 1 line corrected **
                     C = ONE / SQRT( ONE + T ** 2 )
C ** Correction 10. 11/07/97: end of correction **
                     S = T * C
                     E1 = ALPHA + BETA * T
                     E2 = GAMMA - BETA * T
C     write(6,*) ' 2x2 block modified. evals = ', e1, e2
C
C  CHANGE E1 AND E2 TO THEIR ABSOLUTE VALUES AND THEN MULTIPLY THE 
C  THREE 2 * 2 MATRICES TO GET THE MODIFIED ALPHA, BETA AND GAMMA.
C
                     IF ( E1 .LT. ZERO ) THEN
                        NEG2 = NEG2 + 1
                        E1   = - E1
                     END IF
                     IF ( E2 .LT. ZERO ) THEN
                        NEG2 = NEG2 + 1
                        E2   = - E2
                     END IF
                     A( APOS )        = C ** 2 * E1 + S ** 2 * E2
                     A( APOS + 1 )    = C * S * ( E1 - E2 )
                     A( APOS + ALEN ) = S ** 2 * E1 + C ** 2 * E2 
                  END IF
               END IF
            ELSE
               SINGLE = .TRUE.
            END IF
            APOS = APOS + ALEN
            ALEN = ALEN - 1
   50    CONTINUE
         IPOS = IPOS + NCOLS + 1
  100 CONTINUE
      RETURN
      END   
C  THIS VERSION: 13/01/1994 AT 04:39:45 PM.
CS    SUBROUTINE SICCGA( N, NZ, A, INI, INJ, IAI, IAJ, IK, IW, W, C,
CD    SUBROUTINE DICCGA( N, NZ, A, INI, INJ, IAI, IAJ, IK, IW, W, C,
     *                   IFLAG )
C
C  THIS IS A MODIFIED VERSION OF MA31A. IT REQUIRES THE HSL PACKAGE
C  MA31.
C
C  THE PARAMETERS ARE AS FOLLOWS.
C  N        INTEGER  ORDER OF MATRIX. NOT ALTERED BY SUBROUTINE.
C  NZ       INTEGER  NUMBER OF NON-ZEROS IN THE UPPER TRIANGULAR
C      PART OF THE MATRIX INCLUDING THE DIAGONAL.
C  A        DOUBLE PRECISION ARRAY OF LENGTH IAJ. HOLD THE NON-ZEROS
C           OF THE INPUT MATRIX.
C  INI      INTEGER ARRAY OF LENGTH IAI. INI(K) HOLDS THE ROW INDEX
C      OF THE NON-ZERO IN A(K). USED AS WORKSPACE BY ICCGC TO HOLD
C      COLUMN FILE OF THE UPPER TRIANGULAR PART OF THE MATRIX PATTERN.
C  INJ      INTEGER*2 ARRAY OF LENGTH IAJ. INJ(K) HOLDS THE COLUMN
C      INDEX OF A(K).  ON EXIT IT CONTAINS THE COLUMN INDICES OF THE
C      NON-ZEROS OF THE UPPER TRIANGULAR PART OF THE DECOMPOSED MATRIX.
C  IAI      INTEGER  LENGTH OF INI.
C  IAJ      INTEGER  LENGTH OF INJ AND A.
C  IK       INTEGER*2 ARRAY OF LENGTH 4*N.  USED AS WORKSPACE BY ICCGC.
C      IT IS NOT REQUIRED TO BE SET ON ENTRY. ON EXIT IT CONTAINS
C      INFORMATION ABOUT THE DECOMPOSITION. IT SHOULD BE PRESERVED
C      BETWEEN CALL TO ICCGA AND SUBSEQUENT CALLS TO ICCGB.
C      IK(I,1),I=1,2,....,N  HOLDS THE NUMBER OF NON-ZEROS IN THE UPPER
C      TRIANGULAR PART OF ROW I IN THE DECOMPOSED MATRIX.
C      IK(I,2),I=1,2,....,N  HOLDS THE NUMBER OF THE ROW WHICH WAS
C      PIVOT IN THE I'TH PIVOT STEP.
C      IK(I,3),I=1,2,....,N  HOLDS THE THE NUMBER OF NON-ZEROS IN
C      ROW I OF THE LOWER TRIANGULAR PART OF THE INPUTMATRIX.
C  IW       INTEGER ARRAY OF LENGTH 3*N.  IT IS ONLY USED FOR WORKSPACE.
C  W        DOUBLE PRECISION ARRAY OF DIMENSION 3*N USED FOR WORK SPACE
C      DURING THE FACTORIZATION.
C  C        DOUBLE PRECISION VARIABLE. NEW ENTRIES IN ROW I AND COLUMN J
C      ARE IGNORED IN THE FACTORIZATION IF THEIR NUMERICAL VALUE
C      ARE LESS THAN ]C]*DSQRT(D(I)*D(J)).
C  IFLAG    INTEGER.  USED TO EXPRESS POSSIBLE ERRORS WHICH MAY OCCUR
C      DURING THE EXECUTION. A NON-NEGATIVE VALUE ON EXIT INDICATES A
C      SUCCESSFUL RUN. POSSIBLE NEGATIVE VALUES ARE -1 THROUGH -5.
C
      INTEGER           IAI, IAJ, IFLAG, N, NZ
CS    REAL              C
CD    DOUBLE PRECISION  C
      INTEGER           IK( N, 4 ), INI( IAI ), INJ( IAJ ), IW( N, 4 )
CS    REAL              A( IAJ ), W( N, 3 )
CD    DOUBLE PRECISION  A( IAJ ), W( N, 3 )
C
C  LOCAL VARIABLES.
C
      REAL              DD
      INTEGER           IPD, LCOL, LP, LROW, MP, NCP, ND, NUAL, NUCL,
     *                  NURL
      LOGICAL           PHASE2
CS    REAL              ONE, ZERO, ADDON, SMACHR, EPSMCH, DIAMAX, A1, A2
CD    DOUBLE PRECISION  ONE, ZERO, ADDON, DMACHR, EPSMCH, DIAMAX, A1, A2
      INTEGER           I, II, IR, J, K, KI, KJ, KK, KL, KLL, KP, KPP,
     *                  KR, NZ0, NZP1, IC1, IC2, IDUMMY, IR1, IR2
CS    EXTERNAL          SICCGC, SMACHR
CD    EXTERNAL          DICCGC, DMACHR
      INTRINSIC         SQRT, MAX, ABS
C
C  THE COMMON BLOCK VARIABLES ARE AS FOLLOWS.
C  LP,MP    INTEGER  DEFAULT VALUE 6.  UNIT NUMBER FOR ERROR MESSAGES
C      AND DUPLICATE ELEMENT WARNINGS RESPECTIVELY.
C  LROW,LCOL INTEGER   MONITOR THE NUMBER OF ACTIVE LOCATIONS IN ARRAYS
C      INJ AND INI. THE DIFFERENCES NURL-LROW AND NUCL-LCOL ARE EQUAL
C      TO THE NUMBER OF LOCATIONS OCCUPIED BY OLD ROW AND COLUMN COPIES.
C      ON EXIT LROW AND LCOL CONTAINS THE MAXIMUM NUMBER OF ACTIVE
C      ELEMENTS IN THE ROW AND COLUMN FILE RESPECTIVELY AT ANY STAGE
C      OF THE DECOMPOSITION.
C  NCP      INTEGER  MONITOR THE NUMBER OF COMPRESSIONS THAT OCCUR
C      DURING THE DECOMPOSITION.
C  ND       INTEGER HOLD THE NUMBER OF DIAGONAL ENTRIES. IT WILL
C      USUALLY BE EQUAL TO N BUT CAN BE GREATER IN CASE OF DUPLICATE
C      ELEMENTS ON THE DIAGONAL.
C  IPD      INTEGER  HOLD THE PIVOT STEP NUMBER IN WHICH THE ACTIVE
C      PART OF THE MATRIX IS FULL.
C
CS    COMMON            / MA31I  / DD, LP, MP
CD    COMMON            / MA31ID / DD, LP, MP
CS    COMMON            / MA31J  / LROW, LCOL, NCP, ND, IPD
CD    COMMON            / MA31JD / LROW, LCOL, NCP, ND, IPD
CS    COMMON            / MA31K  / NURL, NUCL, NUAL
CD    COMMON            / MA31KD / NURL, NUCL, NUAL
CS    DATA              ZERO, ONE / 0.0E+0, 1.0E+0 /
CD    DATA              ZERO, ONE / 0.0D+0, 1.0D+0 /
C
C CHECK RESTRICTIONS ON INPUT PARAMETERS.
C
      IF ( N .LT. 1 ) THEN
         IF ( LP .GT. 0 ) WRITE ( LP, 2030 )
         IFLAG = - 1
         GO TO 150
      END IF
      IF ( NZ .LT. N ) THEN
         IF ( LP .GT. 0 ) WRITE ( LP, 2040 )
         IFLAG = - 2
C        GO TO 150
      END IF
      IF ( IAI .LT. NZ ) THEN
         IF ( LP .GT. 0 ) WRITE ( LP, 2050 )
         IFLAG = - 3
         GO TO 150
      END IF
      IF ( IAJ .LT. 2 * NZ ) THEN
         IF ( LP .GT. 0 ) WRITE ( LP, 2060 )
         IFLAG = - 4
         GO TO 150
      END IF
C
C INITIALIZE WORK ARRAYS.
C
      DO 20 I = 1, N
         DO 10 J = 1, 3
            W( I, J ) = ZERO
            IK( I, J ) = 0
   10    CONTINUE
   20 CONTINUE
      NUAL = 0
      IFLAG = 0
      ND = 0
      NCP = 0
CS    EPSMCH = SMACHR( 1 )
CD    EPSMCH = DMACHR( 1 )
      DIAMAX = EPSMCH
C
C COUNT NUMBER OF ELEMENTS
C
      DO 30 K = 1, NZ
         I = INI( K )
         J = INJ( K )
         IF ( I .LT. 1 .OR. I .GT. N .OR. J .LT. I .OR. J .GT. N ) THEN
            IF ( LP .GT. 0 ) WRITE ( LP, 2070 ) K, I, J
            IFLAG = - 5
            GO TO 150
         END IF
C
C  CHECK FOR DOUBLE ENTRIES ON THE DIAGONAL AND MOVE DIAGONAL FROM
C  A TO W.
C
         IF ( I .EQ. J ) THEN
            ND = ND + 1
            IF ( W( I, 1 ) .NE. ZERO ) THEN
               IF ( MP .GT. 0 ) WRITE ( MP, 2010 ) I
               IFLAG = 1
            END IF
            W( I, 1 ) = W( I, 1 ) + A( K )
            DIAMAX = MAX( DIAMAX, ABS( W( I, 1 ) ) )
         ELSE
C
C  REMOVE ZEROS.
C
            IF ( A( K ) .EQ. ZERO ) THEN
               ND = ND + 1
            ELSE
               IK( I, 1 ) = IK( I, 1 ) + 1
               IK( J, 2 ) = IK( J, 2 ) + 1
               NUAL = NUAL + 1
               A( NUAL ) = A( K )
               INI( NUAL ) = I
               INJ( NUAL ) = J
            END IF
         END IF
   30 CONTINUE
C
C  NZ0 IS THE NUMBER OF OFF DIAGONAL NON-ZEROS.
C
      NZ0    = NZ - ND
      LCOL   = NZ0
      LROW   = NZ0
      PHASE2 = .FALSE.
      ADDON  = DIAMAX * EPSMCH ** 0.25
C
C TREAT DIAGONAL MATRICES SPECIALLY.
C
      IF ( NZ0 .EQ. 0 ) THEN
         DO 40 I = 1, N
            IK( I, 1 ) = 0
            IK( I, 2 ) = I
            IK( I, 3 ) = 0
C
C  MODIFY THE DIAGONALS IF NECESSARY.
C
            IF ( W( I, 1 ) .LE. ADDON ) THEN
               IFLAG     = 2
               W( I, 1 ) = MAX( ADDON, - W( I, 1 ) )
               IF ( LP .GT. 0 ) WRITE ( MP, 2000 ) I
            END IF
            W( I, 1 ) = ONE / SQRT( W( I, 1 ) )
            W( I, 2 ) = ONE
   40    CONTINUE
         RETURN
      END IF
C
C  NON-DIAGONAL MATRIX. INITIALIZE IW(I,1) AND IW(I,2) TO POINT JUST 
C  BEYOND WHERE THE LAST COMPONENT OF ROW/COLUMN I WILL BE STORED.
C
      KJ = IAI - NZ0 + 1
      KI = 1
      DO 50 I = 1, N
         KI   = KI + IK( I, 1 )
         KJ   = KJ + IK( I, 2 )
         IW( I, 1 ) = KI
         IW( I, 2 ) = KJ
C
C  MODIFY THE DIAGONALS IF NECESSARY.
C
         IF ( W( I, 1 ) .LE. ADDON ) THEN
            IFLAG     = 2
            PHASE2    = .TRUE.
            W( I, 1 ) = MAX( ADDON, - W( I, 1 ) )
            IF ( LP .GT. 0 ) WRITE ( MP, 2000 ) I
         END IF
         W( I, 1 ) = ONE / SQRT( W( I, 1 ) )
   50 CONTINUE
C
C REORDER BY ROWS USING IN-PLACE SORT ALGORITHM.
C
C  IW(I,1) CONTAINS THE NUMBER OF NON-ZEROES IN ROW I OF THE INPUT MATRIX.
C  INITIALIZE IW(I,1) TO POINT JUST BEYOND WHERE THE LAST ELEMENT OF THE
C  ROW WILL BE STORED.
C
C  SAVE CURRENT ENTRY.
C
      DO 53 I = 1, NZ0
         IR1 = INI( I )
C
C  IF IR1.LT.0 THE ELEMENT IS IN PLACE ALREADY.
C
         IF ( IR1 .GE. 0 ) THEN
            IC1 = INJ( I )
            A1  = A( I )
C
C  DETERMINE CORRECT POSITION.
C
            KI = IW( IR1, 1 ) - 1
            DO 51 IDUMMY = 1, NZ0
               IF ( I .EQ. KI ) GO TO 52
C
C  SAVE CONTENTS OF THAT POSITION.
C
               IR2 = INI( KI )
               IC2 = INJ( KI )
C
C  STORE CURRENT ENTRY
C
               INI( KI ) = - IR1
               INJ( KI ) = IC1
               IW( IR1, 1 ) = KI
               IR1 = IR2
               IC1 = IC2
C
C  MAKE CORRESPONDING CHANGES FOR REALS IF REQUIRED.
C
               A2 = A( KI )
               A( KI ) = A1
               A1 = A2
               KI = IW( IR1, 1 ) - 1
   51       CONTINUE
C
C  IF CURRENT ENTRY IS IN PLACE IT IS STORED HERE.
C
   52       CONTINUE
            IF ( IDUMMY .NE. 1 ) THEN
               A( KI ) = A1
               INJ( KI ) = IC1
               INI( KI ) = - IR1
            END IF
            IW( IR1, 1 ) = I
         END IF
   53 CONTINUE
C
C CHECK FOR DOUBLE ENTRIES WHILE USING THE CONSTRUCTED ROW FILE TO SET
C UP THE COLUMN FILE AND COMPRESS THE ROWFILE.
C
      KK = 0
      DO 80 IR = 1, N
         KPP = IW( IR, 1 )
         IW( IR, 1 ) = KK + 1
         KLL = KPP + IK( IR, 1 ) - 1
C
C  LOAD ROW IR INTO W(*,3).
C
         DO 60 K = KPP, KLL
            J = INJ( K )
            IF ( W( J, 3 ) .NE. ZERO ) THEN
               IFLAG = 1
               IF ( MP .GT. 0 ) WRITE ( MP, 2020 ) IR, J
            END IF
            W( J, 3 ) = W( J, 3 ) + A( K )
   60    CONTINUE
C
C  RELOAD ROW IR INTO ARRAYS A AND INJ AND ADJUST INI.
C
         DO 70 K = KPP, KLL
            J = INJ( K )
            IF ( W( J, 3 ) .EQ. ZERO ) THEN
               ND = ND + 1
               LROW = LROW - 1
               LCOL = LCOL - 1
               IK( IR, 1 ) = IK( IR, 1 ) - 1
               IK( J, 2 ) = IK( J, 2 ) - 1
            ELSE
               KK = KK + 1
               A( KK ) = W( J, 3 ) * W( IR, 1 ) * W( J, 1 )
               INJ( KK ) = J
               W( J, 3 ) = ZERO
               KR = IW( J, 2 ) - 1
               IW( J, 2 ) = KR
               INI( KR ) = IR
            END IF
   70    CONTINUE
   80 CONTINUE
      IF ( IFLAG .EQ. 1 ) THEN
C
C  ZERO UNUSED LOCATIONS IN INI.
C
         NZ0 = NZ - ND
         DO 100 I = 1, N - 1
            DO 90 K = IW( I, 2 ) + IK( I, 2 ), IW( I + 1, 2 ) - 1
               INI( K ) = 0
   90       CONTINUE
  100    CONTINUE
      END IF
C
C  STORE INPUT MATRIX.
C
      NUAL = IAJ + 1
      DO 120 II = 1, N
         IK( II, 4 ) = IK( II, 1 )
         I = N - II + 1
         W( I, 2 ) = ONE
         KP = IW( I, 1 )
         KL = KP + IK( I, 1 ) - 1
         DO 110 KK = KP, KL
            K = KP + KL - KK
            NUAL = NUAL - 1
            A( NUAL ) = A( K )
            INJ( NUAL ) = INJ( K )
  110    CONTINUE
         IW( I, 1 ) = NUAL - NZ0
  120 CONTINUE
C
C  SET DIFFERENT PARAMETERS.
C
      NURL = 0
      NZP1 = NZ0 + 1
      NUCL = IW( 1, 2 )
      NUAL = NUAL - NZ0
C
C  ACTIVATE INCOMPLETE FACTORIZATION.
C
CS    CALL SICCGC( N, NZ0, W( 1, 2 ), A( NZP1 ), INI, INJ( NZP1 ), IAI,
CD    CALL DICCGC( N, NZ0, W( 1, 2 ), A( NZP1 ), INI, INJ( NZP1 ), IAI,
     *             IAJ - NZ0, IK, IW, IW( 1, 3 ), W( 1, 3 ), IFLAG, C,
     *             PHASE2 )
C     WRITE( 6, 3030 )
C     WRITE( 6, 3000 ) 1, W( 1, 2 ), IK( 1, 1 ), IK( 1, 2 )
C     WRITE( 6, 3000 ) N, W( N, 2 ), IK( N, 1 ), IK( N, 2 )
C     WRITE( 6, 3040 ) IAJ, NZP1-1
C     WRITE( 6, 3050 ) 1, A( NZP1 - 1 + 1 ), INJ( NZP1 - 1 + 1 )
C
C  THE FACTORIZATION IS TERMINATED.
C
      KP = 1
      DO 140 I = 1, N
         KL = KP + IK( I, 4 ) - 1
         IF ( KP .LE. KL ) THEN
            DO 130 K = KP, KL
               INI( K ) = I
  130       CONTINUE
         END IF
         KP = KL + 1
  140 CONTINUE
      RETURN
C
C  UNSUCCESFUL ENTRY.
C
  150 CONTINUE
      IF ( LP .GT. 0 ) WRITE ( LP, 2080 )
      RETURN
C
C  NON-EXECUTABLE STATEMENTS.
C
 2000 FORMAT ( //, ' Warning from MC31A. Diagonal element', I5, 
     *             ' has been modififed to be positive' )
 2010 FORMAT ( //, '+', 
     *        'Warning: more than one diagonal entry in row',I5 )
 2020 FORMAT ( //, '+', 'Warning: there is more than one entry in row',
     *         I5, ' and column', I5 )
 2030 FORMAT ( //, 34X, ' N  is out of range.' )
 2040 FORMAT ( //, 34X, ' NZ is out of range.' )
 2050 FORMAT ( //, 34X, 'IAI is out of range.' )
 2060 FORMAT ( //, 34X, 'IAJ is out of range.' )
 2070 FORMAT ( //, 34X, 'Element', I7, ' is in row', I5, ' and column',
     *       I5 )
 2080 FORMAT ( '+', 'Error return from ICCGA because' )
C3000 FORMAT( I6, 1P, D12.4, 2I6 )
C3030 FORMAT( ' *** In ICCGA. D and IK ' )
C3040 FORMAT( ' *** In ICCGA IAJ and NZ0 are ', 2I6, ' A and INJ ' )
C3050 FORMAT( I6, 1P, D12.4, I6 )
      END
CS    SUBROUTINE SICCGC( N, NZ, D, A, INI, INJ, IAI, IAJ, IK, IP, IW, W,
CD    SUBROUTINE DICCGC( N, NZ, D, A, INI, INJ, IAI, IAJ, IK, IP, IW, W,
     *                   IFLAG, C, PHASE2 )
      INTEGER           IAI, IAJ, IFLAG, N, NZ
C ** Correction 9. 19/08/96: 1 line added, 1 line corrected **
CS    REAL              C
CD    DOUBLE PRECISION  C
C ** Correction 9. 19/08/96: end of correction **
      LOGICAL           PHASE2
CS    REAL              A( IAJ ), D( N ), W( N )
CD    DOUBLE PRECISION  A( IAJ ), D( N ), W( N )
      INTEGER           IK( N, 3 ), INI( IAI ), INJ( IAJ ), IP( N, 2 ),
     *                  IW( N, 2 )
C
C  THIS IS A MODIFIED VERSION OF MA31C IT REQUIRES THE HSL PACKAGE
C  MA31.
C
C IP(I,1),IP(I,2) POINT TO THE START OF ROW/COLUMN I.
C IK(I,1),IK(I,2) HOLD THE NUMBER OF NONZEROS IN ROW/COLUMN I OF THE
C     LOWER TRIANGULAR PART OF A.
C DURING THE MAIN BODY OF THIS SUBROUTINE THE VECTORS IK(*,3),IW(*,1)
C     IW(*,2) ARE USED TO HOLD DOUBLY LINKED LISTS OF ROWS THAT HAVE
C     NOT BEEN PIVOTAL AND HAVE EQUAL NUMBER OF NONZEROS.
C IK(I,3) HOLD FIRST ROW/COLUMN TO HAVE I NONZEROS OR ZERO IF THERE
C     ARE NONE.
C IW(I,1) HOLD ROW/COLUMN NUMBER OF ROW/COLUMN PRIOR TO ROW I IN ITS
C     LIST OR ZERO IF NONE.
C IW(I,2) HOLD ROW/COLUMN NUMBER OF ROW/COLUMN AFTER ROW I IN ITS LIST
C     OR ZERO IF NONE.
C DURING THE MAIN BODY OF THE SUBROUTINE INI AND INJ KEEP A COLUMN
C     FILE AND A ROW FILE CONTAINING RESPECTIVELY THE ROW NUMBERS OF
C     THE NON-ZEROS OF EACH COLUMN AND THE COLUMN NUMBERS OF THE  NON
C     ZEROS OF EACH ROW. THE IP ARRAYS POINT TO THE START POSITION IN
C     INI AND INJ OF EACH COLUMN AND ROW.
C
      REAL              DD
      INTEGER           IPD, LCOL, LP, LROW, MP, NCP, ND, NUAL, NUCL,
     *                  NURL
      REAL              ALFA, B1, B2, PFILL, PIVT
CS    REAL              SMACHR, AA, AL, CMAX, EPSMCH, ONE, U,
CD    DOUBLE PRECISION  DMACHR, AA, AL, CMAX, EPSMCH, ONE, U,
     *                  ZERO, ADDON
C ** Correction 7. 20/12/93: 2 lines added **
CS    REAL              ONENRM
CD    DOUBLE PRECISION  ONENRM
C ** Correction 7. 20/12/93: end of correction **
      INTEGER           I, II, IIP, IL, IN, IP1, IPDP1, IR, J, J1, JJ,
     *                  JP, K, KC, KK, KL, KL1, KLC, KLJ, KLL, KLR, KP,
     *                  KP2, KPC, KPI, KPJ, KPP, KPR, KR, KRL, KS, L,
     *                  LFULDD, LFULL, MCL, NC, NFILL, NM1, NR, NRJP,
     *                  NZ0, NZC, NZI
      LOGICAL           CHANGE
CS    EXTERNAL          MA31D , SMACHR
CD    EXTERNAL          MA31DD, DMACHR
      INTRINSIC         ABS, MAX, SQRT, FLOAT, IFIX
CS    COMMON            / MA31I  / DD, LP, MP
CD    COMMON            / MA31ID / DD, LP, MP
CS    COMMON            / MA31J  / LROW, LCOL, NCP, ND, IPD
CD    COMMON            / MA31JD / LROW, LCOL, NCP, ND, IPD
CS    COMMON            / MA31K  / NURL, NUCL, NUAL
CD    COMMON            / MA31KD / NURL, NUCL, NUAL
CS    DATA              ZERO, ONE, U, CMAX / 0.0E+0, 1.0E+0,
CS   *                                       1.0E+2, 1.0E+20 /
CD    DATA              ZERO, ONE, U, CMAX / 0.0D+0, 1.0D+0,
CD   *                                       1.0D+2, 1.0D+20 /
C
C INITIALIZE LOCAL VARIABLES.
C
CS    EPSMCH = SMACHR( 1 )
CD    EPSMCH = DMACHR( 1 )
      ADDON  = EPSMCH ** 0.25
      CHANGE = C .GT. ZERO
      NZ0    = NZ
      IPD    = N
      ALFA   = 1.0E+0 / 9.0E-1
      B1     = - 3.0E-2
      B2     = 3.0E-2
      NFILL  = IAJ - NZ0 - N
      MCL    = LCOL
      C      = C ** 2
C
C INITIALIZE IK(*,3).
C
      DO 10 I       = 1, N
         IK( I, 3 ) = 0
   10 CONTINUE
C
C SET UP LINKED LISTS OF ROWS/COLUMNS WITH EQUAL NUMBER OF NON-ZEROS.
C
      DO 20 I = 1, N
         NZI = IK( I, 1 ) + IK( I, 2 ) + 1
         IN = IK( NZI, 3 )
         IK( NZI, 3 ) = I
         IW( I, 2 ) = IN
         IW( I, 1 ) = 0
         IF ( IN .NE. 0 ) IW( IN, 1 ) = I
   20 CONTINUE
C
C START THE ELIMINATION LOOP.
C
      DO 250 IIP = 1, N
C
C SEARCH ROWS WITH NRJP NONZEROS.
C
         DO 30 NRJP = 1, N
            JP = IK( NRJP, 3 )
            IF ( JP .GT. 0 ) GO TO 40
   30    CONTINUE
C
C ROW JP IS USED AS PIVOT.
C
C REMOVE ROWS/COLUMNS INVOLVED IN ELIMINATION FROM ORDERING VECTORS.
C
   40    CONTINUE
         DO 60 L = 1, 2
            KPP = IP( JP, L )
            KLL = IK( JP, L ) + KPP - 1
            DO 50 K = KPP, KLL
               IF ( L .EQ. 1 ) THEN
                  J = INJ( K )
               ELSE
                  J = INI( K )
               END IF
               IL = IW( J, 1 )
               IN = IW( J, 2 )
               IW( J, 2 ) = - 1
               IF ( IN .GE. 0 ) THEN
                  IF ( IL .NE. 0 ) THEN
                     IW( IL, 2 ) = IN
                  ELSE
                     NZ = IK( J, 1 ) + IK( J, 2 ) + 1
                     IK( NZ, 3 ) = IN
                  END IF
                  IF ( IN .GT. 0 ) IW( IN, 1 ) = IL
               END IF
   50       CONTINUE
   60    CONTINUE
C
C REMOVE JP FROM ORDERING VECTORS
C
         IL = IW( JP, 1 )
         IN = IW( JP, 2 )
         IW( JP, 2 ) = - 10
         IF ( IN .GE. 0 ) THEN
            NZ = IK( JP, 1 ) + IK( JP, 2 ) + 1
            IK( NZ, 3 ) = IN
            IF ( IN .GT. 0 ) IW( IN, 1 ) = IL
         END IF
C
C  STORE PIVOT.
C
         IW( JP, 1 ) = - IIP
C
C COMPRESS ROW FILE IF NECESSARY.
C
         IF ( LROW + IK( JP, 1 ) + IK( JP, 2 ) .GT. IAJ - N ) C = CMAX
         IF ( NURL + IK( JP, 1 ) + IK( JP, 2 ) .GE. NUAL ) 
CS   *      CALL MA31D ( A, INJ, IAJ, N, IK, IP, .TRUE. )
CD   *      CALL MA31DD( A, INJ, IAJ, N, IK, IP, .TRUE. )
         KP      = IP( JP, 1 )
         IP( JP, 1 ) = NURL + 1
C
C REMOVE JP FROM COLUMNS CONTAINED IN THE PIVOT ROW.
C
         DO 90 K = KP, IK( JP, 1 ) + KP - 1
            J          = INJ( K )
            KPC        = IP( J, 2 )
            NZ         = IK( J, 2 ) - 1
            IK( J, 2 ) = NZ
            KLC        = KPC + NZ
            IF ( KLC .LE. KPC ) THEN
               INI( KPC ) = 0
            ELSE
               DO 70 KC = KPC, KLC
                  IF ( JP .EQ. INI( KC ) ) GO TO 80
   70          CONTINUE
   80          CONTINUE
               INI( KC )  = INI( KLC )
               INI( KLC ) = 0
            END IF
            LCOL        = LCOL - 1
            NURL        = NURL + 1
            INJ( NURL ) = J
            A( NURL )   = A( K )
            INJ( K )    = 0
   90    CONTINUE
C
C TRANSFORM COLUMN PART OF PIVOT ROW TO THE ROW FILE.
C
         KP2 = IP( JP, 2 )
         DO 120 K = KP2, IK( JP, 2 ) + KP2 - 1
            NURL = NURL + 1
            LCOL = LCOL - 1
            I = INI( K )
            KPR = IP( I, 1 )
            KLR = KPR + IK( I, 1 ) - 1
            DO 100 KR = KPR, KLR
               IF ( JP .EQ. INJ( KR ) ) GO TO 110
  100       CONTINUE
  110       CONTINUE
            INJ( KR ) = INJ( KLR )
            A( NURL ) = A( KR )
            A( KR ) = A( KLR )
            INJ( KLR ) = 0
            IK( I, 1 ) = IK( I, 1 ) - 1
            INJ( NURL ) = I
            INI( K ) = 0
  120    CONTINUE
         NZC = IK( JP, 1 ) + IK( JP, 2 )
         IK( JP, 1 ) = NZC
         IK( JP, 2 ) = 0
C
C  UNPACK PIVOT ROW AND CONTROL DIAGONAL VALUE.
C
         KP = IP( JP, 1 )
         KL = KP + NZC - 1
         ONENRM = ZERO
         DO 130 K = KP, KL
            AA = A( K )
            ONENRM  = ONENRM + ABS( AA )
            J = INJ( K )
            W( J ) = AA
  130    CONTINUE
         IF ( PHASE2 ) THEN
            IF ( D( JP ) .LE. ADDON + ONENRM ) THEN
               IFLAG = 2
               IF ( MP .GT. 0 ) WRITE( MP, 2000 ) JP
               D( JP ) = ADDON + ONENRM
            END IF
         END IF 
         IF ( KP .GT. KL ) GO TO 240
C
C PERFORM ROW OPERATIONS.
C
         DO 220 NC = 1, NZC
            KC = IP( JP, 1 ) + NC - 1
            IR = INJ( KC )
            AL = A( KC ) / D( JP )
C
C COMPRESS ROW FILE IF NECESSARY.
C
            IF ( LROW + IK( IR, 1 ) + IK( JP, 1 ) .GT.
     *          IAJ - N ) C = CMAX
            IF ( NURL + IK( IR, 1 ) + IK( JP, 1 ) .GE. NUAL )
CS   *         CALL MA31D ( A, INJ, IAJ, N, IK, IP, .TRUE. )
CD   *         CALL MA31DD( A, INJ, IAJ, N, IK, IP, .TRUE. )
            KR = IP( IR, 1 )
            KRL = KR + IK( IR, 1 ) - 1
C
C  SCAN THE OTHER ROW AND CHANGE SIGN IN IW FOR EACH COMMON COLUMN NUMBER
C
            DO 140 KS = KR, KRL
               J = INJ( KS )
               IF ( IW( J, 2 ) .EQ. - 1 ) THEN
                  IW( J, 2 ) = 1
                  A( KS ) = A( KS ) - AL * W( J )
               END IF
  140       CONTINUE
C
C  SCAN PIVOT ROW FOR FILLS.
C
            DO 210 KS = KP, KL
               J = INJ( KS )
C
C  ONLY ENTRIES IN THE UPPER TRIANGULAR PART ARE CONSIDERED.
C
               IF ( J .LT. IR ) GO TO 205
               IF ( IW( J, 2 ) .EQ. 1 ) GO TO 205
               AA = - AL * W( J )
               IF ( IR .EQ. J ) THEN
                  D( IR ) = D( IR ) + AA
                  IF ( D( IR ) .LE. ADDON ) PHASE2 = .TRUE.
                  GO TO 205
               END IF
               IF ( AA * AA .LE. C * ABS( D( IR ) * D( J ) ) ) THEN
                  D( J ) = D( J ) + AA
                  D( IR ) = D( IR ) + AA
                  IF ( D( IR ) .LE. ADDON ) PHASE2 = .TRUE.
                  GO TO 205
               END IF
               LROW = LROW + 1
               IK( IR, 1 ) = IK( IR, 1 ) + 1
C
C  IF POSSIBLE PLACE THE NEW ELEMENT NEXT TO THE PRESENT ENTRY.
C
C  SEE IF THERE IS ROOM AT THE END OF THE ENTRY.
C
               IF ( KR .LE. KRL ) THEN
                  IF ( KRL .NE. IAJ ) THEN
                     IF ( INJ( KRL + 1 ) .EQ. 0 ) THEN
                        KRL = KRL + 1
                        INJ( KRL ) = J
                        A( KRL ) = AA
                        GO TO 170
                     END IF
                  END IF
C
C  SEE IF THERE IS ROOM AHEAD OF PRESENT ENTRY.
C
                  IF ( KR .EQ. NUAL ) THEN
                     NUAL = NUAL - 1
                  ELSE
                     IF ( INJ( KR - 1 ) .NE. 0 ) GO TO 150
                  END IF
                  KR = KR - 1
                  IP( IR, 1 ) = KR
                  INJ( KR ) = J
                  A( KR ) = AA
                  GO TO 170
C
C  NEW ENTRY HAS TO BE CREATED.
C
  150             CONTINUE
                  DO 160 KK = KR, KRL
                     NUAL = NUAL - 1
                     INJ( NUAL ) = INJ( KK )
                     A( NUAL ) = A( KK )
                     INJ( KK ) = 0
  160             CONTINUE
               END IF
C
C  ADD THE NEW ELEMENT.
C
               NUAL        = NUAL - 1
               INJ( NUAL ) = J
               A( NUAL )   = AA
               IP( IR, 1 ) = NUAL
               KR          = NUAL
               KRL         = KR + IK( IR, 1 ) - 1
C
C  CREATE FILL IN COLUMN FILE.
C
  170          CONTINUE
               NZ   = IK( J, 2 )
               K    = IP( J, 2 )
               KL1  = K + NZ - 1
               LCOL = LCOL + 1
C
C  IF POSSIBLE PLACE NEW ELEMENT AT THE END OF PRESENT ENTRY.
C
               IF ( NZ .NE. 0 ) THEN
                  IF ( KL1 .NE. IAI ) THEN
                     IF ( INI( KL1 + 1 ) .EQ. 0 ) THEN
                        INI( KL1 + 1 ) = IR
                        GO TO 200
                     END IF
                  END IF
C
C  IF POSSIBLE PLACE ELEMENT AHEAD OF PRESENT ENTRY.
C
                  IF ( K .EQ. NUCL ) THEN
                     IF ( NUCL .EQ. 1 ) GO TO 180
                     NUCL = NUCL - 1
                  ELSE
                     IF ( INI( K - 1 ) .NE. 0 ) GO TO 180
                  END IF
                  K = K - 1
                  INI( K ) = IR
                  IP( J, 2 ) = K
                  GO TO 200
               END IF
C
C  NEW ENTRY HAS TO BE CREATED.
C
  180          CONTINUE
               IF ( NZ + 1 .GE. NUCL ) THEN
C
C  COMPRESS COLUMN FILE IF THERE IS NOT ROOM FOR NEW ENTRY.
C
                  IF ( LCOL + NZ + 2 .GE. IAI ) C = CMAX
CS                CALL MA31D ( A, INI, IAI, N, IK( 1, 2 ), IP( 1, 2 ),
CD                CALL MA31DD( A, INI, IAI, N, IK( 1, 2 ), IP( 1, 2 ),
     *                         .FALSE. )
                  K = IP( J, 2 )
                  KL1 = K + NZ - 1
               END IF
C
C  TRANSFER OLD ENTRY INTO NEW.
C
               DO 190 KK = K, KL1
                  NUCL = NUCL - 1
                  INI( NUCL ) = INI( KK )
                  INI( KK ) = 0
  190          CONTINUE
C
C  ADD THE NEW ELEMENT.
C
               NUCL = NUCL - 1
               INI( NUCL ) = IR
               IP( J, 2 ) = NUCL
  200          CONTINUE
               IK( J, 2 ) = NZ + 1
  205          CONTINUE
               IW( J, 2 ) = - 1
  210       CONTINUE
  220    CONTINUE
C
C  UPDATE ORDERING ARRAYS.
C
         DO 230 K       = KP, KL
            J           = INJ( K )
            W( J )      = ZERO
            A( K )      = A( K ) / D( JP )
            NZ          = IK( J, 1 ) + IK( J, 2 ) + 1
            IN          = IK( NZ, 3 )
            IW( J, 2 )  = IN
            IW( J, 1 )  = 0
            IK( NZ, 3 ) = J
            IF ( IN .NE. 0 ) IW( IN, 1 ) = J
  230    CONTINUE
         MCL = MAX( MCL, LCOL )
         PIVT = FLOAT( IIP ) / FLOAT( N )
C
C  GIVE WARNING IF AVAILABLE SPACE IS USED TOO EARLY.
C
         IF ( C .EQ. CMAX ) THEN
            IF ( IPD .LT. IIP ) GO TO 240
            IPD = IIP
            IF ( PIVT .GT. .9 ) GO TO 240
            IFLAG = 4
            IF ( MP .GT. 0 ) WRITE ( MP, 2010 ) IIP
            GO TO 240
         ELSE
C
C  CHANGE C IF NECESSARY.
C
            IF ( .NOT. CHANGE ) GO TO 240
            PFILL = FLOAT( LROW - NZ0 ) / FLOAT( NFILL )
            IF ( PIVT .GE. 9.0E-1 ) GO TO 240
            IF ( PFILL .GT. ALFA * PIVT + B1 ) THEN
               IF ( PFILL .LT. ALFA * PIVT + B2 ) GO TO 240
               C = 2.25D+0 * C
            END IF
            ALFA = ( 1.0E+0 - PFILL ) / ( 9.0E-1 - PIVT )
            B1   = PFILL - PIVT * ALFA - 3.0E-2
            B2   = B1 + 6.0E-2
         END IF
C
C  IF THE MATRIX IS FULL, STOP THE SPARSE ANALYZE.
C
  240    CONTINUE
         NR     = N - IIP
         LFULL  = NR * ( NR - 1 ) / 2
         LFULDD = IFIX( DD * FLOAT( LFULL ) )
         IF ( LCOL .GE. LFULDD .AND. NURL + LFULL .LT. IAJ ) GO TO 260
  250 CONTINUE
C
C  SPARSE ELIMINATION LOOP TERMINATES. 
C  FACTORIZE THE REMAINING FULL MATRIX.
C
  260 CONTINUE
      IPD  = IIP
      C    = SQRT( C )
      LCOL = MCL
      IF ( .NOT. CHANGE ) C = - C
C
C  THE ORDER OF THE FULL MATRIX IS NR.
C  LOOP THROUGH ROWS IN THE ACTIVE MATRIX AND STORE ROW NUMBERS IN INI.
C
      KK       = 0
      DO 280 I = 1, NR
         JP    = IK( I, 3 )
  270    CONTINUE
         IF ( JP .GT. 0 ) THEN
            KK        = KK + 1
            INI( KK ) = JP
            JP        = IW( JP, 2 )
            GO TO 270
         END IF
         IF ( KK .EQ. NR ) GO TO 290
  280 CONTINUE
C
C  MAKE A SORT OF THE ROW NUMBERS IN INI.
C
  290 CONTINUE
      DO 310 I = 1, NR - 1
         J1 = I + 1
         DO 300 J = J1, NR
            IF ( INI( J ) .LE. INI( I ) ) THEN
               JJ       = INI( I )
               INI( I ) = INI( J )
               INI( J ) = JJ
            END IF
  300    CONTINUE
  310 CONTINUE
      DO 320 I       = 1, NR
         II          = INI( I )
         IW( II, 1 ) = - IPD - I
  320 CONTINUE
C
C  MAKE AN ORDERED LIST OF THE PIVOTS.
C
      DO 330 I       = 1, N
         IR          = - IW( I, 1 )
         IK( IR, 2 ) = I
  330 CONTINUE
C
C  MOVE FULL MATRIX TO THE FRONT AND ORDER.
C
      IPDP1 = IPD + 1
      NM1   = N - 1
C ** Correction 1. 05/09/92: 1 line corrected **
      IF ( IPDP1 .GT. NM1 ) GO TO 400
C ** Correction 1. 05/09/92: end of correction **
      DO 360 IIP = IPDP1, NM1
         JP      = IK( IIP, 2 )
         KP      = IP( JP, 1 )
         KL      = KP + IK( JP, 1 ) - 1
C
C   MOVE ROW JP TO W.
C
         DO 340 K    = KP, KL
            J        = INJ( K )
            INJ( K ) = 0
            W( J )   = A( K )
  340    CONTINUE
C
C  COMPRESS FILE IF NECESSARY.
C
CS       IF ( NURL + N - IIP .GE. NUAL ) CALL MA31D ( A, INJ, IAJ, N,
CD       IF ( NURL + N - IIP .GE. NUAL ) CALL MA31DD( A, INJ, IAJ, N,
     *       IK, IP, .TRUE. )
         IP( JP, 1 ) = NURL + 1
         IK( JP, 1 ) = N - IIP
C
C  MOVE ROWS AND COLUMN INDICES INTO PIVOTAL ORDER.
C
         DO 350 I = IIP + 1, N
            J = IK( I, 2 )
            NURL = NURL + 1
            A( NURL ) = W( J )
            INJ( NURL ) = J
            W( J ) = ZERO
  350    CONTINUE
  360 CONTINUE
      LROW = NURL
C
C  FACTORIZE THE FULL MATRIX.
C
      DO 390 IIP   = IPDP1, NM1
         JP        = IK( IIP, 2 )
         KPI       = IP( JP, 1 )
C
C  IF NECESSARY, MODIFY THE DIAGONAL TO ENSURE THAT IT IS POSITIVE.
C
         ONENRM = ZERO
         DO 365 K  = KPI, KPI + IK( JP, 1 ) - 1
            ONENRM  = ONENRM + ABS( A( K ) )
  365    CONTINUE
         IF ( PHASE2 ) THEN
            IF ( D( JP ) .LE. ADDON + ONENRM ) THEN
               D( JP ) = ADDON + ONENRM
               IFLAG   = 2
               IF ( MP .GT. 0 ) WRITE( MP, 2000 ) JP
            END IF
         END IF
C
C  LOOP THROUGH THE OTHER ROW
C
         IP1 = IIP + 1
         IF ( IP1 .NE. N ) THEN
            DO 380 J     = IP1, NM1
               JJ        = IK( J, 2 )
               KPJ       = IP( JJ, 1 )
               KLJ       = KPJ + IK( JJ, 1 ) - 1
               AL        = A( KPI ) / D( JP )
               D( JJ )   = D( JJ ) - AL * A( KPI )
               IF ( D( JJ ) .LE. ADDON ) PHASE2 = .TRUE.
               KK        = KPI + 1
               DO 370 K  = KPJ, KLJ
                  A( K ) = A( K ) - AL * A( KK )
                  KK     = KK + 1
  370          CONTINUE
C
C  STORE FACTOR AND PROCEED TO NEXT ROW.
C
               A( KPI ) = AL
               KPI      = KPI + 1
  380       CONTINUE
         END IF
C
C  MODIFY LAST DIAGONAL ENTRY
C
         JJ       = IK( N, 2 )
         AL       = A( KPI ) / D( JP )
         D( JJ )  = D( JJ ) - AL * A( KPI )
         IF ( D( JJ ) .LE. ADDON ) PHASE2 = .TRUE.
         A( KPI ) = AL
  390 CONTINUE
C
C  IF NECESSARY, MODIFY THE DIAGONAL TO ENSURE THAT IT IS POSITIVE.
C
C ** Correction 2. 05/09/92: 2 lines added **
  400 CONTINUE
      JJ = IK( N, 2 )
C ** Correction 2. 05/09/92: end of correction **
      IF ( D( JJ ) .LE. ADDON ) THEN
         D( JJ ) = ADDON
         IFLAG   = 2
         PHASE2 = .TRUE.
         IF ( MP .GT. 0 ) WRITE ( MP, 2000 ) JJ
      END IF
      RETURN
C
C  NON-EXECUTABLE STATEMENTS.
C
 2000 FORMAT ( //, ' Warning modification of zero or negative',
     *        ' diagonal entry has been performed in location', I7 )
 2010 FORMAT ( //, ' Warning available space used at pivot step', I7 )
      END
CS    SUBROUTINE SICCGB( N, A, INJ, IAJ, D, SCALE, IK, B )
CD    SUBROUTINE DICCGB( N, A, INJ, IAJ, D, SCALE, IK, B )
C
C  THIS SUBROUTINE COMPUTES THE PRODUCT VECTOR OF THE INVERSE OF THE
C  PRECONDITIONING MATRIX AND THE VECTOR B. THE RESULT IS PLACED IN B.
C
      INTEGER           IAJ, N
CS    REAL              A( IAJ ), B( N ), D( N ), SCALE( N )
CD    DOUBLE PRECISION  A( IAJ ), B( N ), D( N ), SCALE( N )
      INTEGER           IK( N, 2 ), INJ( IAJ )
C
C  LOCAL VARIABLES.
C
      INTEGER           I, IPD, LCOL, LROW, NCP, ND
CS    REAL              BIC, BIR
CD    DOUBLE PRECISION  BIC, BIR
      INTEGER           IC, IIP, IPI, IR, K, KL, KP
CS    COMMON            / MA31J  / LROW, LCOL, NCP, ND, IPD
CD    COMMON            / MA31JD / LROW, LCOL, NCP, ND, IPD
C
C  SCALE THE RHS.
C
      DO 10 I = 1, N
         B( I ) = B( I ) * SCALE( I )
   10 CONTINUE
C
C  PRECONDITION THE RHS.
C
      KP = 1
C
C  PERFORM THE FORWARD SUBSTITUTION.
C
      DO 30 IIP = 1, N
         IC = IK( IIP, 2 )
         KL = KP + IK( IC, 1 ) - 1
         BIC = B( IC )
         DO 20 K = KP, KL
            IR = INJ( K )
            B( IR ) = B( IR ) - A( K ) * BIC
   20    CONTINUE
         KP = KL + 1
   30 CONTINUE
      KL = LROW
C
C  PERFORM THE BACK SUBSTITUTION.
C
      DO 50 IPI = 1, N
         IIP = N + 1 - IPI
         IR = IK( IIP, 2 )
         BIR = 0.0
         KP = KL - IK( IR, 1 ) + 1
         DO 40 K = KP, KL
            IC = INJ( K )
            BIR = BIR - A( K ) * B( IC )
   40    CONTINUE
         B( IR ) = B( IR ) / D( IR ) + BIR
         KL = KP - 1
   50 CONTINUE
C
C  RESCALE THE RHS.
C
      DO 60 I = 1, N
         B( I ) = B( I ) * SCALE( I )
   60 CONTINUE
      RETURN
      END
