share/odepack/fortran/dprjs.f

   1 *DECK DPRJS
   2       SUBROUTINE DPRJS (NEQ,Y,YH,NYH,EWT,FTEM,SAVF,WK,IWK,F,JAC)
   3       EXTERNAL F,JAC
   4       INTEGER NEQ, NYH, IWK
   5       DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WK
   6       DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*),
   7      1   WK(*), IWK(*)
   8       INTEGER IOWND, IOWNS,
   9      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  10      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  11      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  12       INTEGER IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
  13      1   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
  14      2   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
  15      3   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
  16       DOUBLE PRECISION ROWNS,
  17      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  18       DOUBLE PRECISION CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH
  19       COMMON /DLS001/ ROWNS(209),
  20      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  21      2   IOWND(6), IOWNS(6),
  22      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  23      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  24      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  25       COMMON /DLSS01/ CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH,
  26      1   IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
  27      2   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
  28      3   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
  29      4   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
  30       INTEGER I, IMUL, J, JJ, JOK, JMAX, JMIN, K, KMAX, KMIN, NG
  31       DOUBLE PRECISION CON, DI, FAC, HL0, PIJ, R, R0, RCON, RCONT,
  32      1   SRUR, DVNORM
  33 C-----------------------------------------------------------------------
  34 C DPRJS is called to compute and process the matrix
  35 C P = I - H*EL(1)*J , where J is an approximation to the Jacobian.
  36 C J is computed by columns, either by the user-supplied routine JAC
  37 C if MITER = 1, or by finite differencing if MITER = 2.
  38 C if MITER = 3, a diagonal approximation to J is used.
  39 C if MITER = 1 or 2, and if the existing value of the Jacobian
  40 C (as contained in P) is considered acceptable, then a new value of
  41 C P is reconstructed from the old value.  In any case, when MITER
  42 C is 1 or 2, the P matrix is subjected to LU decomposition in CDRV.
  43 C P and its LU decomposition are stored (separately) in WK.
  44 C
  45 C In addition to variables described previously, communication
  46 C with DPRJS uses the following:
  47 C Y     = array containing predicted values on entry.
  48 C FTEM  = work array of length N (ACOR in DSTODE).
  49 C SAVF  = array containing f evaluated at predicted y.
  50 C WK    = real work space for matrices.  On output it contains the
  51 C         inverse diagonal matrix if MITER = 3, and P and its sparse
  52 C         LU decomposition if MITER is 1 or 2.
  53 C         Storage of matrix elements starts at WK(3).
  54 C         WK also contains the following matrix-related data:
  55 C         WK(1) = SQRT(UROUND), used in numerical Jacobian increments.
  56 C         WK(2) = H*EL0, saved for later use if MITER = 3.
  57 C IWK   = integer work space for matrix-related data, assumed to
  58 C         be equivalenced to WK.  In addition, WK(IPRSP) and IWK(IPISP)
  59 C         are assumed to have identical locations.
  60 C EL0   = EL(1) (input).
  61 C IERPJ = output error flag (in Common).
  62 C       = 0 if no error.
  63 C       = 1  if zero pivot found in CDRV.
  64 C       = 2  if a singular matrix arose with MITER = 3.
  65 C       = -1 if insufficient storage for CDRV (should not occur here).
  66 C       = -2 if other error found in CDRV (should not occur here).
  67 C JCUR  = output flag showing status of (approximate) Jacobian matrix:
  68 C          = 1 to indicate that the Jacobian is now current, or
  69 C          = 0 to indicate that a saved value was used.
  70 C This routine also uses other variables in Common.
  71 C-----------------------------------------------------------------------
  72       HL0 = H*EL0
  73       CON = -HL0
  74       IF (MITER .EQ. 3) GO TO 300
  75 C See whether J should be reevaluated (JOK = 0) or not (JOK = 1). ------
  76       JOK = 1
  77       IF (NST .EQ. 0 .OR. NST .GE. NSLJ+MSBJ) JOK = 0
  78       IF (ICF .EQ. 1 .AND. ABS(RC - 1.0D0) .LT. CCMXJ) JOK = 0
  79       IF (ICF .EQ. 2) JOK = 0
  80       IF (JOK .EQ. 1) GO TO 250
  81 C
  82 C MITER = 1 or 2, and the Jacobian is to be reevaluated. ---------------
  83  20   JCUR = 1
  84       NJE = NJE + 1
  85       NSLJ = NST
  86       IPLOST = 0
  87       CONMIN = ABS(CON)
  88       GO TO (100, 200), MITER
  89 C
  90 C If MITER = 1, call JAC, multiply by scalar, and add identity. --------
  91  100  CONTINUE
  92       KMIN = IWK(IPIAN)
  93       DO 130 J = 1, N
  94         KMAX = IWK(IPIAN+J) - 1
  95         DO 110 I = 1,N
  96  110      FTEM(I) = 0.0D0
  97         CALL JAC (NEQ, TN, Y, J, IWK(IPIAN), IWK(IPJAN), FTEM)
  98         DO 120 K = KMIN, KMAX
  99           I = IWK(IBJAN+K)
 100           WK(IBA+K) = FTEM(I)*CON
 101           IF (I .EQ. J) WK(IBA+K) = WK(IBA+K) + 1.0D0
 102  120      CONTINUE
 103         KMIN = KMAX + 1
 104  130    CONTINUE
 105       GO TO 290
 106 C
 107 C If MITER = 2, make NGP calls to F to approximate J and P. ------------
 108  200  CONTINUE
 109       FAC = DVNORM(N, SAVF, EWT)
 110       R0 = 1000.0D0 * ABS(H) * UROUND * N * FAC
 111       IF (R0 .EQ. 0.0D0) R0 = 1.0D0
 112       SRUR = WK(1)
 113       JMIN = IWK(IPIGP)
 114       DO 240 NG = 1,NGP
 115         JMAX = IWK(IPIGP+NG) - 1
 116         DO 210 J = JMIN,JMAX
 117           JJ = IWK(IBJGP+J)
 118           R = MAX(SRUR*ABS(Y(JJ)),R0/EWT(JJ))
 119  210      Y(JJ) = Y(JJ) + R
 120         CALL F (NEQ, TN, Y, FTEM)
 121         DO 230 J = JMIN,JMAX
 122           JJ = IWK(IBJGP+J)
 123           Y(JJ) = YH(JJ,1)
 124           R = MAX(SRUR*ABS(Y(JJ)),R0/EWT(JJ))
 125           FAC = -HL0/R
 126           KMIN =IWK(IBIAN+JJ)
 127           KMAX =IWK(IBIAN+JJ+1) - 1
 128           DO 220 K = KMIN,KMAX
 129             I = IWK(IBJAN+K)
 130             WK(IBA+K) = (FTEM(I) - SAVF(I))*FAC
 131             IF (I .EQ. JJ) WK(IBA+K) = WK(IBA+K) + 1.0D0
 132  220        CONTINUE
 133  230      CONTINUE
 134         JMIN = JMAX + 1
 135  240    CONTINUE
 136       NFE = NFE + NGP
 137       GO TO 290
 138 C
 139 C If JOK = 1, reconstruct new P from old P. ----------------------------
 140  250  JCUR = 0
 141       RCON = CON/CON0
 142       RCONT = ABS(CON)/CONMIN
 143       IF (RCONT .GT. RBIG .AND. IPLOST .EQ. 1) GO TO 20
 144       KMIN = IWK(IPIAN)
 145       DO 275 J = 1,N
 146         KMAX = IWK(IPIAN+J) - 1
 147         DO 270 K = KMIN,KMAX
 148           I = IWK(IBJAN+K)
 149           PIJ = WK(IBA+K)
 150           IF (I .NE. J) GO TO 260
 151           PIJ = PIJ - 1.0D0
 152           IF (ABS(PIJ) .GE. PSMALL) GO TO 260
 153             IPLOST = 1
 154             CONMIN = MIN(ABS(CON0),CONMIN)
 155  260      PIJ = PIJ*RCON
 156           IF (I .EQ. J) PIJ = PIJ + 1.0D0
 157           WK(IBA+K) = PIJ
 158  270      CONTINUE
 159         KMIN = KMAX + 1
 160  275    CONTINUE
 161 C
 162 C Do numerical factorization of P matrix. ------------------------------
 163  290  NLU = NLU + 1
 164       CON0 = CON
 165       IERPJ = 0
 166       DO 295 I = 1,N
 167  295    FTEM(I) = 0.0D0
 168       CALL CDRV (N,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN),
 169      1   WK(IPA),FTEM,FTEM,NSP,IWK(IPISP),WK(IPRSP),IESP,2,IYS)
 170       IF (IYS .EQ. 0) RETURN
 171       IMUL = (IYS - 1)/N
 172       IERPJ = -2
 173       IF (IMUL .EQ. 8) IERPJ = 1
 174       IF (IMUL .EQ. 10) IERPJ = -1
 175       RETURN
 176 C
 177 C If MITER = 3, construct a diagonal approximation to J and P. ---------
 178  300  CONTINUE
 179       JCUR = 1
 180       NJE = NJE + 1
 181       WK(2) = HL0
 182       IERPJ = 0
 183       R = EL0*0.1D0
 184       DO 310 I = 1,N
 185  310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
 186       CALL F (NEQ, TN, Y, WK(3))
 187       NFE = NFE + 1
 188       DO 320 I = 1,N
 189         R0 = H*SAVF(I) - YH(I,2)
 190         DI = 0.1D0*R0 - H*(WK(I+2) - SAVF(I))
 191         WK(I+2) = 1.0D0
 192         IF (ABS(R0) .LT. UROUND/EWT(I)) GO TO 320
 193         IF (ABS(DI) .EQ. 0.0D0) GO TO 330
 194         WK(I+2) = 0.1D0*R0/DI
 195  320    CONTINUE
 196       RETURN
 197  330  IERPJ = 2
 198       RETURN
 199 C----------------------- End of Subroutine DPRJS -----------------------
 200       END