      SUBROUTINE PREP (NEQ, Y, YH, SAVF, EWT, FTEM, IA, JA, 
     1                     WK, IWK, IPPER, F, JAC)
CLLL. OPTIMIZE
      EXTERNAL F,JAC
      INTEGER NEQ, IA, JA, IWK, IPPER
      INTEGER IOWND, IOWNS,
     1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
     2   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
      INTEGER IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
     1   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
     2   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
     3   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
      INTEGER I, IBR, IER, IPIL, IPIU, IPTT1, IPTT2, J, JFOUND, K,
     1   KNEW, KMAX, KMIN, LDIF, LENIGP, LIWK, MAXG, NP1, NZSUT
      DOUBLE PRECISION Y, YH, SAVF, EWT, FTEM, WK 
      DOUBLE PRECISION ROWNS, 
     1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
      DOUBLE PRECISION CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH
      DOUBLE PRECISION DQ, DYJ, ERWT, FAC, YJ
      DIMENSION NEQ(1), Y(1), YH(1), SAVF(1), EWT(1), FTEM(1),
     1   IA(1), JA(1), WK(1), IWK(1)
      COMMON /LS0001/ ROWNS(209),
     2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
     3   IOWND(14), IOWNS(6), 
     4   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
     5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
      COMMON /LSS001/ CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH,
     1   IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
     2   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
     3   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
     4   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
C-----------------------------------------------------------------------
C THIS ROUTINE PERFORMS PREPROCESSING RELATED TO THE SPARSE LINEAR
C SYSTEMS THAT MUST BE SOLVED IF MITER = 1 OR 2.
C THE OPERATIONS THAT ARE PERFORMED HERE ARE..
C  * COMPUTE SPARSENESS STRUCTURE OF JACOBIAN ACCORDING TO MOSS,
C  * COMPUTE GROUPING OF COLUMN INDICES (MITER = 2),
C  * COMPUTE A NEW ORDERING OF ROWS AND COLUMNS OF THE MATRIX,
C  * REORDER JA CORRESPONDING TO THE NEW ORDERING,
C  * PERFORM A SYMBOLIC LU FACTORIZATION OF THE MATRIX, AND 
C  * SET POINTERS FOR SEGMENTS OF THE IWK/WK ARRAY.
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, PREP USES THE
C FOLLOWING FOR COMMUNICATION..
C YH     = THE HISTORY ARRAY.  ONLY THE FIRST COLUMN, CONTAINING THE
C          CURRENT Y VECTOR, IS USED.  USED ONLY IF MOSS .NE. 0.
C SAVF   = A WORK ARRAY OF LENGTH NEQ, USED ONLY IF MOSS .NE. 0.
C EWT    = ARRAY OF LENGTH NEQ CONTAINING (INVERTED) ERROR WEIGHTS.
C          USED ONLY IF MOSS = 2 OR IF ISTATE = MOSS = 1.
C FTEM   = A WORK ARRAY OF LENGTH NEQ, IDENTICAL TO ACOR IN THE DRIVER,
C          USED ONLY IF MOSS = 2.
C WK     = A REAL WORK ARRAY OF LENGTH LENWK, IDENTICAL TO WM IN
C          THE DRIVER.
C IWK    = INTEGER WORK ARRAY, ASSUMED TO OCCUPY THE SAME SPACE AS WK.
C LENWK  = THE LENGTH OF THE WORK ARRAYS WK AND IWK.
C ISTATC = A COPY OF THE DRIVER INPUT ARGUMENT ISTATE (= 1 ON THE
C          FIRST CALL, = 3 ON A CONTINUATION CALL).
C IYS    = FLAG VALUE FROM ODRV OR CDRV.
C IPPER  = OUTPUT ERROR FLAG WITH THE FOLLOWING VALUES AND MEANINGS.. 
C          0  NO ERROR.
C         -1  INSUFFICIENT STORAGE FOR INTERNAL STRUCTURE POINTERS.
C         -2  INSUFFICIENT STORAGE FOR JGROUP.
C         -3  INSUFFICIENT STORAGE FOR ODRV.
C         -4  OTHER ERROR FLAG FROM ODRV (SHOULD NEVER OCCUR).
C         -5  INSUFFICIENT STORAGE FOR CDRV.
C         -6  OTHER ERROR FLAG FROM CDRV.
C-----------------------------------------------------------------------
      IBIAN = LRAT*2
      IPIAN = IBIAN + 1
      NP1 = N + 1
      IPJAN = IPIAN + NP1
      IBJAN = IPJAN - 1
      LIWK = LENWK*LRAT
      IF (IPJAN+N-1 .GT. LIWK) GO TO 210
      IF (MOSS .EQ. 0) GO TO 30
C
      IF (ISTATC .EQ. 3) GO TO 20
C ISTATE = 1 AND MOSS .NE. 0.  PERTURB Y FOR STRUCTURE DETERMINATION. --
      DO 10 I = 1,N 
        ERWT = 1.0D0/EWT(I)
        FAC = 1.0D0 + 1.0D0/(DBLE(I)+1.0D0)     
        Y(I) = Y(I) + FAC*DSIGN(ERWT,Y(I))
 10     CONTINUE
      GO TO (70, 100), MOSS
C
 20   CONTINUE
C ISTATE = 3 AND MOSS .NE. 0.  LOAD Y FROM YH(*,1). --------------------
      DO 25 I = 1,N 
 25     Y(I) = YH(I)
      GO TO (70, 100), MOSS
C
C MOSS = 0.  PROCESS USER-S IA,JA.  ADD DIAGONAL ENTRIES IF NECESSARY. -
 30   KNEW = IPJAN
      KMIN = IA(1)
      IWK(IPIAN) = 1
      DO 60 J = 1,N 
        JFOUND = 0
        KMAX = IA(J+1) - 1
        IF (KMIN .GT. KMAX) GO TO 45
        DO 40 K = KMIN,KMAX
          I = JA(K) 
          IF (I .EQ. J) JFOUND = 1
          IF (KNEW .GT. LIWK) GO TO 210 
          IWK(KNEW) = I
          KNEW = KNEW + 1
 40       CONTINUE
        IF (JFOUND .EQ. 1) GO TO 50
 45     IF (KNEW .GT. LIWK) GO TO 210
        IWK(KNEW) = J
        KNEW = KNEW + 1
 50     IWK(IPIAN+J) = KNEW + 1 - IPJAN 
        KMIN = KMAX + 1
 60     CONTINUE
      GO TO 140
C
C MOSS = 1.  COMPUTE STRUCTURE FROM USER-SUPPLIED JACOBIAN ROUTINE JAC.
 70   CONTINUE
C A DUMMY CALL TO F ALLOWS USER TO CREATE TEMPORARIES FOR USE IN JAC. --
      CALL F (NEQ, TN, Y, SAVF)
      K = IPJAN
      IWK(IPIAN) = 1
      DO 90 J = 1,N 
        IF (K .GT. LIWK) GO TO 210
        IWK(K) = J
        K = K + 1
        DO 75 I = 1,N
 75       SAVF(I) = 0.0D0
        CALL JAC (NEQ, TN, Y, J, IWK(IPIAN), IWK(IPJAN), SAVF)
        DO 80 I = 1,N
          IF (DABS(SAVF(I)) .LE. SETH) GO TO 80
          IF (I .EQ. J) GO TO 80
          IF (K .GT. LIWK) GO TO 210
          IWK(K) = I
          K = K + 1 
 80       CONTINUE
        IWK(IPIAN+J) = K + 1 - IPJAN
 90     CONTINUE
      GO TO 140
C
C MOSS = 2.  COMPUTE STRUCTURE FROM RESULTS OF N + 1 CALLS TO F. -------
 100  K = IPJAN
      IWK(IPIAN) = 1
      CALL F (NEQ, TN, Y, SAVF)
      DO 120 J = 1,N
        IF (K .GT. LIWK) GO TO 210
        IWK(K) = J
        K = K + 1
        YJ = Y(J)
        ERWT = 1.0D0/EWT(J)
        DYJ = DSIGN(ERWT,YJ)
        Y(J) = YJ + DYJ
        CALL F (NEQ, TN, Y, FTEM)
        Y(J) = YJ
        DO 110 I = 1,N
          DQ = (FTEM(I) - SAVF(I))/DYJ
          IF (DABS(DQ) .LE. SETH) GO TO 110
          IF (I .EQ. J) GO TO 110
          IF (K .GT. LIWK) GO TO 210
          IWK(K) = I
          K = K + 1 
 110      CONTINUE
        IWK(IPIAN+J) = K + 1 - IPJAN
 120    CONTINUE
C
 140  CONTINUE
      IF (MOSS .EQ. 0 .OR. ISTATC .NE. 1) GO TO 150
C IF ISTATE = 1 AND MOSS .NE. 0, RESTORE Y FROM YH. --------------------
      DO 145 I = 1,N
 145    Y(I) = YH(I)
 150  NNZ = IWK(IPIAN+N) - 1
      LENIGP = 0
      IPIGP = IPJAN + NNZ
      IF (MITER .NE. 2) GO TO 160
C
C COMPUTE GROUPING OF COLUMN INDICES (MITER = 2). ----------------------
      MAXG = NP1
      IPJGP = IPJAN + NNZ
      IBJGP = IPJGP - 1
      IPIGP = IPJGP + N
      IPTT1 = IPIGP + NP1
      IPTT2 = IPTT1 + N
      LREQ = IPTT2 + N - 1
      IF (LREQ .GT. LIWK) GO TO 220
      CALL JGROUP (N, IWK(IPIAN), IWK(IPJAN), MAXG, NGP, IWK(IPIGP),
     1   IWK(IPJGP), IWK(IPTT1), IWK(IPTT2), IER) 
      IF (IER .NE. 0) GO TO 220
      LENIGP = NGP + 1
C
C COMPUTE NEW ORDERING OF ROWS/COLUMNS OF JACOBIAN. --------------------
 160  IPR = IPIGP + LENIGP
      IPC = IPR
      IPIC = IPC + N
      IPISP = IPIC + N
      IPRSP = (IPISP - 2)/LRAT + 2
      IESP = LENWK + 1 - IPRSP
      IF (IESP .LT. 0) GO TO 230
      IBR = IPR - 1 
      DO 170 I = 1,N
 170    IWK(IBR+I) = I
      NSP = LIWK + 1 - IPISP
      CALL ODRV (N, IWK(IPIAN), IWK(IPJAN), WK, IWK(IPR), IWK(IPIC),
     1   NSP, IWK(IPISP), 1, IYS)
      IF (IYS .EQ. 11*N+1) GO TO 240
      IF (IYS .NE. 0) GO TO 230
C
C REORDER JAN AND DO SYMBOLIC LU FACTORIZATION OF MATRIX. --------------
      IPA = LENWK + 1 - NNZ
      NSP = IPA - IPRSP
      LREQ = MAX0(12*N/LRAT, 6*N/LRAT+2*N+NNZ) + 3
      LREQ = LREQ + IPRSP - 1 + NNZ
      IF (LREQ .GT. LENWK) GO TO 250
      IBA = IPA - 1 
      DO 180 I = 1,NNZ
 180    WK(IBA+I) = 0.0D0
      IPISP = LRAT*(IPRSP - 1) + 1
      CALL CDRV (N,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN), 
     1   WK(IPA),WK(IPA),WK(IPA),NSP,IWK(IPISP),WK(IPRSP),IESP,5,IYS) 
      LREQ = LENWK - IESP
      IF (IYS .EQ. 10*N+1) GO TO 250
      IF (IYS .NE. 0) GO TO 260
      IPIL = IPISP
      IPIU = IPIL + 2*N + 1
      NZU = IWK(IPIL+N) - IWK(IPIL)
      NZL = IWK(IPIU+N) - IWK(IPIU)
      IF (LRAT .GT. 1) GO TO 190
      CALL ADJLR (N, IWK(IPISP), LDIF)
      LREQ = LREQ + LDIF
 190  CONTINUE
      IF (LRAT .EQ. 2 .AND. NNZ .EQ. N) LREQ = LREQ + 1
      NSP = NSP + LREQ - LENWK
      IPA = LREQ + 1 - NNZ
      IBA = IPA - 1 
      IPPER = 0
      RETURN
C
 210  IPPER = -1
      LREQ = 2 + (2*N + 1)/LRAT
      LREQ = MAX0(LENWK+1,LREQ)
      RETURN
C
 220  IPPER = -2
      LREQ = (LREQ - 1)/LRAT + 1
      RETURN
C
 230  IPPER = -3
      CALL CNTNZU (N, IWK(IPIAN), IWK(IPJAN), NZSUT)
      LREQ = LENWK - IESP + (3*N + 4*NZSUT - 1)/LRAT + 1
      RETURN
C
 240  IPPER = -4
      RETURN
C
 250  IPPER = -5
      RETURN
C
 260  IPPER = -6
      LREQ = LENWK
      RETURN
C----------------------- END OF SUBROUTINE PREP ------------------------
      END 
