share/odepack/fortran/dpjibt.f

   1 *DECK DPJIBT
   2       SUBROUTINE DPJIBT (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, IER, IIA, IIB, IIC, IPA, IPB, IPC, IRES, J, J1, J2,
  22      1   K, K1, LENP, LBLOX, LPB, LPC, MB, MBSQ, MWID, NB
  23       DOUBLE PRECISION CON, FAC, HL0, R, SRUR
  24 C-----------------------------------------------------------------------
  25 C DPJIBT 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 and r = g(t,y) - A(t,y)*s.  Here J is computed by the user-supplied
  28 C routine JAC if MITER = 1, or by finite differencing if MITER = 2.
  29 C J is stored in WM, rescaled, and ADDA is called to generate P.
  30 C P is then subjected to LU decomposition by DDECBT in preparation
  31 C for later solution of linear systems with P as coefficient matrix.
  32 C
  33 C In addition to variables described previously, communication
  34 C with DPJIBT uses the following:
  35 C Y     = array containing predicted values on entry.
  36 C RTEM  = work array of length N (ACOR in DSTODI).
  37 C SAVR  = array used for output only.  On output it contains the
  38 C         residual evaluated at current values of t and y.
  39 C S     = array containing predicted values of dy/dt (SAVF in DSTODI).
  40 C WM    = real work space for matrices.  On output it contains the
  41 C         LU decomposition of P.
  42 C         Storage of matrix elements starts at WM(3).
  43 C         WM also contains the following matrix-related data:
  44 C         WM(1) = SQRT(UROUND), used in numerical Jacobian increments.
  45 C IWM   = integer work space containing pivot information, starting at
  46 C         IWM(21).  IWM also contains block structure parameters
  47 C         MB = IWM(1) and NB = IWM(2).
  48 C EL0   = EL(1) (input).
  49 C IERPJ = output error flag.
  50 C         = 0 if no trouble occurred,
  51 C         = 1 if the P matrix was found to be unfactorable,
  52 C         = IRES (= 2 or 3) if RES returned IRES = 2 or 3.
  53 C JCUR  = output flag = 1 to indicate that the Jacobian matrix
  54 C         (or approximation) is now current.
  55 C This routine also uses the Common variables EL0, H, TN, UROUND,
  56 C MITER, N, NFE, and NJE.
  57 C-----------------------------------------------------------------------
  58       NJE = NJE + 1
  59       HL0 = H*EL0
  60       IERPJ = 0
  61       JCUR = 1
  62       MB = IWM(1)
  63       NB = IWM(2)
  64       MBSQ = MB*MB
  65       LBLOX = MBSQ*NB
  66       LPB = 3 + LBLOX
  67       LPC = LPB + LBLOX
  68       LENP = 3*LBLOX
  69       GO TO (100, 200), MITER
  70 C If MITER = 1, call RES, then JAC, and multiply by scalar. ------------
  71  100  IRES = 1
  72       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
  73       NFE = NFE + 1
  74       IF (IRES .GT. 1) GO TO 600
  75       DO 110 I = 1,LENP
  76  110    WM(I+2) = 0.0D0
  77       CALL JAC (NEQ, TN, Y, S, MB, NB, WM(3), WM(LPB), WM(LPC))
  78       CON = -HL0
  79       DO 120 I = 1,LENP
  80  120    WM(I+2) = WM(I+2)*CON
  81       GO TO 260
  82 C
  83 C If MITER = 2, make 3*MB + 1 calls to RES to approximate J. -----------
  84  200  CONTINUE
  85       IRES = -1
  86       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
  87       NFE = NFE + 1
  88       IF (IRES .GT. 1) GO TO 600
  89       MWID = 3*MB
  90       SRUR = WM(1)
  91       DO 205 I = 1,LENP
  92  205    WM(2+I) = 0.0D0
  93       DO 250 K = 1,3
  94         DO 240 J = 1,MB
  95 C         Increment Y(I) for group of column indices, and call RES. ----
  96           J1 = J+(K-1)*MB
  97           DO 210 I = J1,N,MWID
  98             R = MAX(SRUR*ABS(Y(I)),0.01D0/EWT(I))
  99             Y(I) = Y(I) + R
 100  210      CONTINUE
 101           CALL RES (NEQ, TN, Y, S, RTEM, IRES)
 102           NFE = NFE + 1
 103           IF (IRES .GT. 1) GO TO 600
 104           DO 215 I = 1,N
 105  215        RTEM(I) = RTEM(I) - SAVR(I)
 106           K1 = K
 107           DO 230 I = J1,N,MWID
 108 C           Get Jacobian elements in column I (block-column K1). -------
 109             Y(I) = YH(I,1)
 110             R = MAX(SRUR*ABS(Y(I)),0.01D0/EWT(I))
 111             FAC = -HL0/R
 112 C           Compute and load elements PA(*,J,K1). ----------------------
 113             IIA = I - J
 114             IPA = 2 + (J-1)*MB + (K1-1)*MBSQ
 115             DO 221 J2 = 1,MB
 116  221          WM(IPA+J2) = RTEM(IIA+J2)*FAC
 117             IF (K1 .LE. 1) GO TO 223
 118 C           Compute and load elements PB(*,J,K1-1). --------------------
 119             IIB = IIA - MB
 120             IPB = IPA + LBLOX - MBSQ
 121             DO 222 J2 = 1,MB
 122  222          WM(IPB+J2) = RTEM(IIB+J2)*FAC
 123  223        CONTINUE
 124             IF (K1 .GE. NB) GO TO 225
 125 C           Compute and load elements PC(*,J,K1+1). --------------------
 126             IIC = IIA + MB
 127             IPC = IPA + 2*LBLOX + MBSQ
 128             DO 224 J2 = 1,MB
 129  224          WM(IPC+J2) = RTEM(IIC+J2)*FAC
 130  225        CONTINUE
 131             IF (K1 .NE. 3) GO TO 227
 132 C           Compute and load elements PC(*,J,1). -----------------------
 133             IPC = IPA - 2*MBSQ + 2*LBLOX
 134             DO 226 J2 = 1,MB
 135  226          WM(IPC+J2) = RTEM(J2)*FAC
 136  227        CONTINUE
 137             IF (K1 .NE. NB-2) GO TO 229
 138 C           Compute and load elements PB(*,J,NB). ----------------------
 139             IIB = N - MB
 140             IPB = IPA + 2*MBSQ + LBLOX
 141             DO 228 J2 = 1,MB
 142  228          WM(IPB+J2) = RTEM(IIB+J2)*FAC
 143  229      K1 = K1 + 3
 144  230      CONTINUE
 145  240    CONTINUE
 146  250  CONTINUE
 147 C RES call for first corrector iteration. ------------------------------
 148       IRES = 1
 149       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
 150       NFE = NFE + 1
 151       IF (IRES .GT. 1) GO TO 600
 152 C Add matrix A. --------------------------------------------------------
 153  260  CONTINUE
 154       CALL ADDA (NEQ, TN, Y, MB, NB, WM(3), WM(LPB), WM(LPC))
 155 C Do LU decomposition on P. --------------------------------------------
 156       CALL DDECBT (MB, NB, WM(3), WM(LPB), WM(LPC), IWM(21), IER)
 157       IF (IER .NE. 0) IERPJ = 1
 158       RETURN
 159 C Error return for IRES = 2 or IRES = 3 return from RES. ---------------
 160  600  IERPJ = IRES
 161       RETURN
 162 C----------------------- End of Subroutine DPJIBT ----------------------
 163       END