share/odepack/fortran/dprep.f

   1 *DECK DPREP
   2       SUBROUTINE DPREP (NEQ, Y, YH, SAVF, EWT, FTEM, IA, JA,
   3      1                     WK, IWK, IPPER, F, JAC)
   4       EXTERNAL F,JAC
   5       INTEGER NEQ, IA, JA, IWK, IPPER
   6       DOUBLE PRECISION Y, YH, SAVF, EWT, FTEM, WK
   7       DIMENSION NEQ(*), Y(*), YH(*), SAVF(*), EWT(*), FTEM(*),
   8      1   IA(*), JA(*), WK(*), IWK(*)
   9       INTEGER IOWND, IOWNS,
  10      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  11      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  12      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  13       INTEGER IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
  14      1   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
  15      2   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
  16      3   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
  17       DOUBLE PRECISION ROWNS,
  18      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  19       DOUBLE PRECISION CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH
  20       COMMON /DLS001/ ROWNS(209),
  21      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  22      2   IOWND(6), IOWNS(6),
  23      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  24      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  25      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  26       COMMON /DLSS01/ CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH,
  27      1   IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
  28      2   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
  29      3   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
  30      4   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
  31       INTEGER I, IBR, IER, IPIL, IPIU, IPTT1, IPTT2, J, JFOUND, K,
  32      1   KNEW, KMAX, KMIN, LDIF, LENIGP, LIWK, MAXG, NP1, NZSUT
  33       DOUBLE PRECISION DQ, DYJ, ERWT, FAC, YJ
  34 C-----------------------------------------------------------------------
  35 C This routine performs preprocessing related to the sparse linear
  36 C systems that must be solved if MITER = 1 or 2.
  37 C The operations that are performed here are:
  38 C  * compute sparseness structure of Jacobian according to MOSS,
  39 C  * compute grouping of column indices (MITER = 2),
  40 C  * compute a new ordering of rows and columns of the matrix,
  41 C  * reorder JA corresponding to the new ordering,
  42 C  * perform a symbolic LU factorization of the matrix, and
  43 C  * set pointers for segments of the IWK/WK array.
  44 C In addition to variables described previously, DPREP uses the
  45 C following for communication:
  46 C YH     = the history array.  Only the first column, containing the
  47 C          current Y vector, is used.  Used only if MOSS .ne. 0.
  48 C SAVF   = a work array of length NEQ, used only if MOSS .ne. 0.
  49 C EWT    = array of length NEQ containing (inverted) error weights.
  50 C          Used only if MOSS = 2 or if ISTATE = MOSS = 1.
  51 C FTEM   = a work array of length NEQ, identical to ACOR in the driver,
  52 C          used only if MOSS = 2.
  53 C WK     = a real work array of length LENWK, identical to WM in
  54 C          the driver.
  55 C IWK    = integer work array, assumed to occupy the same space as WK.
  56 C LENWK  = the length of the work arrays WK and IWK.
  57 C ISTATC = a copy of the driver input argument ISTATE (= 1 on the
  58 C          first call, = 3 on a continuation call).
  59 C IYS    = flag value from ODRV or CDRV.
  60 C IPPER  = output error flag with the following values and meanings:
  61 C          0  no error.
  62 C         -1  insufficient storage for internal structure pointers.
  63 C         -2  insufficient storage for JGROUP.
  64 C         -3  insufficient storage for ODRV.
  65 C         -4  other error flag from ODRV (should never occur).
  66 C         -5  insufficient storage for CDRV.
  67 C         -6  other error flag from CDRV.
  68 C-----------------------------------------------------------------------
  69       IBIAN = LRAT*2
  70       IPIAN = IBIAN + 1
  71       NP1 = N + 1
  72       IPJAN = IPIAN + NP1
  73       IBJAN = IPJAN - 1
  74       LIWK = LENWK*LRAT
  75       IF (IPJAN+N-1 .GT. LIWK) GO TO 210
  76       IF (MOSS .EQ. 0) GO TO 30
  77 C
  78       IF (ISTATC .EQ. 3) GO TO 20
  79 C ISTATE = 1 and MOSS .ne. 0.  Perturb Y for structure determination. --
  80       DO 10 I = 1,N
  81         ERWT = 1.0D0/EWT(I)
  82         FAC = 1.0D0 + 1.0D0/(I + 1.0D0)
  83         Y(I) = Y(I) + FAC*SIGN(ERWT,Y(I))
  84  10     CONTINUE
  85       GO TO (70, 100), MOSS
  86 C
  87  20   CONTINUE
  88 C ISTATE = 3 and MOSS .ne. 0.  Load Y from YH(*,1). --------------------
  89       DO 25 I = 1,N
  90  25     Y(I) = YH(I)
  91       GO TO (70, 100), MOSS
  92 C
  93 C MOSS = 0.  Process user's IA,JA.  Add diagonal entries if necessary. -
  94  30   KNEW = IPJAN
  95       KMIN = IA(1)
  96       IWK(IPIAN) = 1
  97       DO 60 J = 1,N
  98         JFOUND = 0
  99         KMAX = IA(J+1) - 1
 100         IF (KMIN .GT. KMAX) GO TO 45
 101         DO 40 K = KMIN,KMAX
 102           I = JA(K)
 103           IF (I .EQ. J) JFOUND = 1
 104           IF (KNEW .GT. LIWK) GO TO 210
 105           IWK(KNEW) = I
 106           KNEW = KNEW + 1
 107  40       CONTINUE
 108         IF (JFOUND .EQ. 1) GO TO 50
 109  45     IF (KNEW .GT. LIWK) GO TO 210
 110         IWK(KNEW) = J
 111         KNEW = KNEW + 1
 112  50     IWK(IPIAN+J) = KNEW + 1 - IPJAN
 113         KMIN = KMAX + 1
 114  60     CONTINUE
 115       GO TO 140
 116 C
 117 C MOSS = 1.  Compute structure from user-supplied Jacobian routine JAC.
 118  70   CONTINUE
 119 C A dummy call to F allows user to create temporaries for use in JAC. --
 120       CALL F (NEQ, TN, Y, SAVF)
 121       K = IPJAN
 122       IWK(IPIAN) = 1
 123       DO 90 J = 1,N
 124         IF (K .GT. LIWK) GO TO 210
 125         IWK(K) = J
 126         K = K + 1
 127         DO 75 I = 1,N
 128  75       SAVF(I) = 0.0D0
 129         CALL JAC (NEQ, TN, Y, J, IWK(IPIAN), IWK(IPJAN), SAVF)
 130         DO 80 I = 1,N
 131           IF (ABS(SAVF(I)) .LE. SETH) GO TO 80
 132           IF (I .EQ. J) GO TO 80
 133           IF (K .GT. LIWK) GO TO 210
 134           IWK(K) = I
 135           K = K + 1
 136  80       CONTINUE
 137         IWK(IPIAN+J) = K + 1 - IPJAN
 138  90     CONTINUE
 139       GO TO 140
 140 C
 141 C MOSS = 2.  Compute structure from results of N + 1 calls to F. -------
 142  100  K = IPJAN
 143       IWK(IPIAN) = 1
 144       CALL F (NEQ, TN, Y, SAVF)
 145       DO 120 J = 1,N
 146         IF (K .GT. LIWK) GO TO 210
 147         IWK(K) = J
 148         K = K + 1
 149         YJ = Y(J)
 150         ERWT = 1.0D0/EWT(J)
 151         DYJ = SIGN(ERWT,YJ)
 152         Y(J) = YJ + DYJ
 153         CALL F (NEQ, TN, Y, FTEM)
 154         Y(J) = YJ
 155         DO 110 I = 1,N
 156           DQ = (FTEM(I) - SAVF(I))/DYJ
 157           IF (ABS(DQ) .LE. SETH) GO TO 110
 158           IF (I .EQ. J) GO TO 110
 159           IF (K .GT. LIWK) GO TO 210
 160           IWK(K) = I
 161           K = K + 1
 162  110      CONTINUE
 163         IWK(IPIAN+J) = K + 1 - IPJAN
 164  120    CONTINUE
 165 C
 166  140  CONTINUE
 167       IF (MOSS .EQ. 0 .OR. ISTATC .NE. 1) GO TO 150
 168 C If ISTATE = 1 and MOSS .ne. 0, restore Y from YH. --------------------
 169       DO 145 I = 1,N
 170  145    Y(I) = YH(I)
 171  150  NNZ = IWK(IPIAN+N) - 1
 172       LENIGP = 0
 173       IPIGP = IPJAN + NNZ
 174       IF (MITER .NE. 2) GO TO 160
 175 C
 176 C Compute grouping of column indices (MITER = 2). ----------------------
 177       MAXG = NP1
 178       IPJGP = IPJAN + NNZ
 179       IBJGP = IPJGP - 1
 180       IPIGP = IPJGP + N
 181       IPTT1 = IPIGP + NP1
 182       IPTT2 = IPTT1 + N
 183       LREQ = IPTT2 + N - 1
 184       IF (LREQ .GT. LIWK) GO TO 220
 185       CALL JGROUP (N, IWK(IPIAN), IWK(IPJAN), MAXG, NGP, IWK(IPIGP),
 186      1   IWK(IPJGP), IWK(IPTT1), IWK(IPTT2), IER)
 187       IF (IER .NE. 0) GO TO 220
 188       LENIGP = NGP + 1
 189 C
 190 C Compute new ordering of rows/columns of Jacobian. --------------------
 191  160  IPR = IPIGP + LENIGP
 192       IPC = IPR
 193       IPIC = IPC + N
 194       IPISP = IPIC + N
 195       IPRSP = (IPISP - 2)/LRAT + 2
 196       IESP = LENWK + 1 - IPRSP
 197       IF (IESP .LT. 0) GO TO 230
 198       IBR = IPR - 1
 199       DO 170 I = 1,N
 200  170    IWK(IBR+I) = I
 201       NSP = LIWK + 1 - IPISP
 202       CALL ODRV (N, IWK(IPIAN), IWK(IPJAN), WK, IWK(IPR), IWK(IPIC),
 203      1   NSP, IWK(IPISP), 1, IYS)
 204       IF (IYS .EQ. 11*N+1) GO TO 240
 205       IF (IYS .NE. 0) GO TO 230
 206 C
 207 C Reorder JAN and do symbolic LU factorization of matrix. --------------
 208       IPA = LENWK + 1 - NNZ
 209       NSP = IPA - IPRSP
 210       LREQ = MAX(12*N/LRAT, 6*N/LRAT+2*N+NNZ) + 3
 211       LREQ = LREQ + IPRSP - 1 + NNZ
 212       IF (LREQ .GT. LENWK) GO TO 250
 213       IBA = IPA - 1
 214       DO 180 I = 1,NNZ
 215  180    WK(IBA+I) = 0.0D0
 216       IPISP = LRAT*(IPRSP - 1) + 1
 217       CALL CDRV (N,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN),
 218      1   WK(IPA),WK(IPA),WK(IPA),NSP,IWK(IPISP),WK(IPRSP),IESP,5,IYS)
 219       LREQ = LENWK - IESP
 220       IF (IYS .EQ. 10*N+1) GO TO 250
 221       IF (IYS .NE. 0) GO TO 260
 222       IPIL = IPISP
 223       IPIU = IPIL + 2*N + 1
 224       NZU = IWK(IPIL+N) - IWK(IPIL)
 225       NZL = IWK(IPIU+N) - IWK(IPIU)
 226       IF (LRAT .GT. 1) GO TO 190
 227       CALL ADJLR (N, IWK(IPISP), LDIF)
 228       LREQ = LREQ + LDIF
 229  190  CONTINUE
 230       IF (LRAT .EQ. 2 .AND. NNZ .EQ. N) LREQ = LREQ + 1
 231       NSP = NSP + LREQ - LENWK
 232       IPA = LREQ + 1 - NNZ
 233       IBA = IPA - 1
 234       IPPER = 0
 235       RETURN
 236 C
 237  210  IPPER = -1
 238       LREQ = 2 + (2*N + 1)/LRAT
 239       LREQ = MAX(LENWK+1,LREQ)
 240       RETURN
 241 C
 242  220  IPPER = -2
 243       LREQ = (LREQ - 1)/LRAT + 1
 244       RETURN
 245 C
 246  230  IPPER = -3
 247       CALL CNTNZU (N, IWK(IPIAN), IWK(IPJAN), NZSUT)
 248       LREQ = LENWK - IESP + (3*N + 4*NZSUT - 1)/LRAT + 1
 249       RETURN
 250 C
 251  240  IPPER = -4
 252       RETURN
 253 C
 254  250  IPPER = -5
 255       RETURN
 256 C
 257  260  IPPER = -6
 258       LREQ = LENWK
 259       RETURN
 260 C----------------------- End of Subroutine DPREP -----------------------
 261       END