share/odepack/fortran/dprepji.f

   1 *DECK DPREPJI
   2       SUBROUTINE DPREPJI (NEQ, Y, YH, NYH, EWT, RTEM, SAVR, S, WM, IWM,
   3      1   RES, JAC, ADDA)
   4       EXTERNAL RES, JAC, ADDA
   5       INTEGER NEQ, NYH, IWM
   6       DOUBLE PRECISION Y, YH, EWT, RTEM, SAVR, S, WM
   7       DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), RTEM(*),
   8      1   S(*), SAVR(*), WM(*), IWM(*)
   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       DOUBLE PRECISION ROWNS,
  14      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  15       COMMON /DLS001/ ROWNS(209),
  16      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  17      2   IOWND(6), IOWNS(6),
  18      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  19      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  20      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  21       INTEGER I, I1, I2, IER, II, IRES, J, J1, JJ, LENP,
  22      1   MBA, MBAND, MEB1, MEBAND, ML, ML3, MU
  23       DOUBLE PRECISION CON, FAC, HL0, R, SRUR, YI, YJ, YJJ
  24 C-----------------------------------------------------------------------
  25 C DPREPJI is called by DSTODI to compute and process the matrix
  26 C P = A - H*EL(1)*J , where J is an approximation to the Jacobian dr/dy,
  27 C where r = g(t,y) - A(t,y)*s.  Here J is computed by the user-supplied
  28 C routine JAC if MITER = 1 or 4, or by finite differencing if MITER =
  29 C 2 or 5.  J is stored in WM, rescaled, and ADDA is called to generate
  30 C P. P is then subjected to LU decomposition in preparation
  31 C for later solution of linear systems with P as coefficient
  32 C matrix.  This is done by DGEFA if MITER = 1 or 2, and by
  33 C DGBFA if MITER = 4 or 5.
  34 C
  35 C In addition to variables described previously, communication
  36 C with DPREPJI uses the following:
  37 C Y     = array containing predicted values on entry.
  38 C RTEM  = work array of length N (ACOR in DSTODI).
  39 C SAVR  = array used for output only.  On output it contains the
  40 C         residual evaluated at current values of t and y.
  41 C S     = array containing predicted values of dy/dt (SAVF in DSTODI).
  42 C WM    = real work space for matrices.  On output it contains the
  43 C         LU decomposition of P.
  44 C         Storage of matrix elements starts at WM(3).
  45 C         WM also contains the following matrix-related data:
  46 C         WM(1) = SQRT(UROUND), used in numerical Jacobian increments.
  47 C IWM   = integer work space containing pivot information, starting at
  48 C         IWM(21).  IWM also contains the band parameters
  49 C         ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
  50 C EL0   = el(1) (input).
  51 C IERPJ = output error flag.
  52 C         = 0 if no trouble occurred,
  53 C         = 1 if the P matrix was found to be singular,
  54 C         = IRES (= 2 or 3) if RES returned IRES = 2 or 3.
  55 C JCUR  = output flag = 1 to indicate that the Jacobian matrix
  56 C         (or approximation) is now current.
  57 C This routine also uses the Common variables EL0, H, TN, UROUND,
  58 C MITER, N, NFE, and NJE.
  59 C-----------------------------------------------------------------------
  60       NJE = NJE + 1
  61       HL0 = H*EL0
  62       IERPJ = 0
  63       JCUR = 1
  64       GO TO (100, 200, 300, 400, 500), MITER
  65 C If MITER = 1, call RES, then JAC, and multiply by scalar. ------------
  66  100  IRES = 1
  67       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
  68       NFE = NFE + 1
  69       IF (IRES .GT. 1) GO TO 600
  70       LENP = N*N
  71       DO 110 I = 1,LENP
  72  110    WM(I+2) = 0.0D0
  73       CALL JAC ( NEQ, TN, Y, S, 0, 0, WM(3), N )
  74       CON = -HL0
  75       DO 120 I = 1,LENP
  76  120    WM(I+2) = WM(I+2)*CON
  77       GO TO 240
  78 C If MITER = 2, make N + 1 calls to RES to approximate J. --------------
  79  200  CONTINUE
  80       IRES = -1
  81       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
  82       NFE = NFE + 1
  83       IF (IRES .GT. 1) GO TO 600
  84       SRUR = WM(1)
  85       J1 = 2
  86       DO 230 J = 1,N
  87         YJ = Y(J)
  88         R = MAX(SRUR*ABS(YJ),0.01D0/EWT(J))
  89         Y(J) = Y(J) + R
  90         FAC = -HL0/R
  91         CALL RES ( NEQ, TN, Y, S, RTEM, IRES )
  92         NFE = NFE + 1
  93         IF (IRES .GT. 1) GO TO 600
  94         DO 220 I = 1,N
  95  220      WM(I+J1) = (RTEM(I) - SAVR(I))*FAC
  96         Y(J) = YJ
  97         J1 = J1 + N
  98  230    CONTINUE
  99       IRES = 1
 100       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
 101       NFE = NFE + 1
 102       IF (IRES .GT. 1) GO TO 600
 103 C Add matrix A. --------------------------------------------------------
 104  240  CONTINUE
 105       CALL ADDA(NEQ, TN, Y, 0, 0, WM(3), N)
 106 C Do LU decomposition on P. --------------------------------------------
 107       CALL DGEFA (WM(3), N, N, IWM(21), IER)
 108       IF (IER .NE. 0) IERPJ = 1
 109       RETURN
 110 C Dummy section for MITER = 3
 111  300  RETURN
 112 C If MITER = 4, call RES, then JAC, and multiply by scalar. ------------
 113  400  IRES = 1
 114       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
 115       NFE = NFE + 1
 116       IF (IRES .GT. 1) GO TO 600
 117       ML = IWM(1)
 118       MU = IWM(2)
 119       ML3 = ML + 3
 120       MBAND = ML + MU + 1
 121       MEBAND = MBAND + ML
 122       LENP = MEBAND*N
 123       DO 410 I = 1,LENP
 124  410    WM(I+2) = 0.0D0
 125       CALL JAC ( NEQ, TN, Y, S, ML, MU, WM(ML3), MEBAND)
 126       CON = -HL0
 127       DO 420 I = 1,LENP
 128  420    WM(I+2) = WM(I+2)*CON
 129       GO TO 570
 130 C If MITER = 5, make ML + MU + 2 calls to RES to approximate J. --------
 131  500  CONTINUE
 132       IRES = -1
 133       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
 134       NFE = NFE + 1
 135       IF (IRES .GT. 1) GO TO 600
 136       ML = IWM(1)
 137       MU = IWM(2)
 138       ML3 = ML + 3
 139       MBAND = ML + MU + 1
 140       MBA = MIN(MBAND,N)
 141       MEBAND = MBAND + ML
 142       MEB1 = MEBAND - 1
 143       SRUR = WM(1)
 144       DO 560 J = 1,MBA
 145         DO 530 I = J,N,MBAND
 146           YI = Y(I)
 147           R = MAX(SRUR*ABS(YI),0.01D0/EWT(I))
 148  530      Y(I) = Y(I) + R
 149         CALL RES ( NEQ, TN, Y, S, RTEM, IRES)
 150         NFE = NFE + 1
 151         IF (IRES .GT. 1) GO TO 600
 152         DO 550 JJ = J,N,MBAND
 153           Y(JJ) = YH(JJ,1)
 154           YJJ = Y(JJ)
 155           R = MAX(SRUR*ABS(YJJ),0.01D0/EWT(JJ))
 156           FAC = -HL0/R
 157           I1 = MAX(JJ-MU,1)
 158           I2 = MIN(JJ+ML,N)
 159           II = JJ*MEB1 - ML + 2
 160           DO 540 I = I1,I2
 161  540        WM(II+I) = (RTEM(I) - SAVR(I))*FAC
 162  550      CONTINUE
 163  560    CONTINUE
 164       IRES = 1
 165       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
 166       NFE = NFE + 1
 167       IF (IRES .GT. 1) GO TO 600
 168 C Add matrix A. --------------------------------------------------------
 169   570 CONTINUE
 170       CALL ADDA(NEQ, TN, Y, ML, MU, WM(ML3), MEBAND)
 171 C Do LU decomposition of P. --------------------------------------------
 172       CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
 173       IF (IER .NE. 0) IERPJ = 1
 174       RETURN
 175 C Error return for IRES = 2 or IRES = 3 return from RES. ---------------
 176  600  IERPJ = IRES
 177       RETURN
 178 C----------------------- End of Subroutine DPREPJI ---------------------
 179       END