share/odepack/fortran/dprepj.f

   1 *DECK DPREPJ
   2       SUBROUTINE DPREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
   3      1   F, JAC)
   4 C***BEGIN PROLOGUE  DPREPJ
   5 C***SUBSIDIARY
   6 C***PURPOSE  Compute and process Newton iteration matrix.
   7 C***TYPE      DOUBLE PRECISION (SPREPJ-S, DPREPJ-D)
   8 C***AUTHOR  Hindmarsh, Alan C., (LLNL)
   9 C***DESCRIPTION
  10 C
  11 C  DPREPJ is called by DSTODE to compute and process the matrix
  12 C  P = I - h*el(1)*J , where J is an approximation to the Jacobian.
  13 C  Here J is computed by the user-supplied routine JAC if
  14 C  MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5.
  15 C  If MITER = 3, a diagonal approximation to J is used.
  16 C  J is stored in WM and replaced by P.  If MITER .ne. 3, P is then
  17 C  subjected to LU decomposition in preparation for later solution
  18 C  of linear systems with P as coefficient matrix.  This is done
  19 C  by DGEFA if MITER = 1 or 2, and by DGBFA if MITER = 4 or 5.
  20 C
  21 C  In addition to variables described in DSTODE and DLSODE prologues,
  22 C  communication with DPREPJ uses the following:
  23 C  Y     = array containing predicted values on entry.
  24 C  FTEM  = work array of length N (ACOR in DSTODE).
  25 C  SAVF  = array containing f evaluated at predicted y.
  26 C  WM    = real work space for matrices.  On output it contains the
  27 C          inverse diagonal matrix if MITER = 3 and the LU decomposition
  28 C          of P if MITER is 1, 2 , 4, or 5.
  29 C          Storage of matrix elements starts at WM(3).
  30 C          WM also contains the following matrix-related data:
  31 C          WM(1) = SQRT(UROUND), used in numerical Jacobian increments.
  32 C          WM(2) = H*EL0, saved for later use if MITER = 3.
  33 C  IWM   = integer work space containing pivot information, starting at
  34 C          IWM(21), if MITER is 1, 2, 4, or 5.  IWM also contains band
  35 C          parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
  36 C  EL0   = EL(1) (input).
  37 C  IERPJ = output error flag,  = 0 if no trouble, .gt. 0 if
  38 C          P matrix found to be singular.
  39 C  JCUR  = output flag = 1 to indicate that the Jacobian matrix
  40 C          (or approximation) is now current.
  41 C  This routine also uses the COMMON variables EL0, H, TN, UROUND,
  42 C  MITER, N, NFE, and NJE.
  43 C
  44 C***SEE ALSO  DLSODE
  45 C***ROUTINES CALLED  DGBFA, DGEFA, DVNORM
  46 C***COMMON BLOCKS    DLS001
  47 C***REVISION HISTORY  (YYMMDD)
  48 C   791129  DATE WRITTEN
  49 C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
  50 C   890504  Minor cosmetic changes.  (FNF)
  51 C   930809  Renamed to allow single/double precision versions. (ACH)
  52 C   010418  Reduced size of Common block /DLS001/. (ACH)
  53 C   031105  Restored 'own' variables to Common block /DLS001/, to
  54 C           enable interrupt/restart feature. (ACH)
  55 C***END PROLOGUE  DPREPJ
  56 C**End
  57       EXTERNAL F, JAC
  58       INTEGER NEQ, NYH, IWM
  59       DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WM
  60       DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*),
  61      1   WM(*), IWM(*)
  62       INTEGER IOWND, IOWNS,
  63      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  64      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  65      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  66       DOUBLE PRECISION ROWNS,
  67      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  68       COMMON /DLS001/ ROWNS(209),
  69      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  70      2   IOWND(6), IOWNS(6),
  71      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  72      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  73      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  74       INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP,
  75      1   MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1
  76       DOUBLE PRECISION CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ,
  77      1   DVNORM
  78 C
  79 C***FIRST EXECUTABLE STATEMENT  DPREPJ
  80       NJE = NJE + 1
  81       IERPJ = 0
  82       JCUR = 1
  83       HL0 = H*EL0
  84       GO TO (100, 200, 300, 400, 500), MITER
  85 C If MITER = 1, call JAC and multiply by scalar. -----------------------
  86  100  LENP = N*N
  87       DO 110 I = 1,LENP
  88  110    WM(I+2) = 0.0D0
  89       CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N)
  90       CON = -HL0
  91       DO 120 I = 1,LENP
  92  120    WM(I+2) = WM(I+2)*CON
  93       GO TO 240
  94 C If MITER = 2, make N calls to F to approximate J. --------------------
  95  200  FAC = DVNORM (N, SAVF, EWT)
  96       R0 = 1000.0D0*ABS(H)*UROUND*N*FAC
  97       IF (R0 .EQ. 0.0D0) R0 = 1.0D0
  98       SRUR = WM(1)
  99       J1 = 2
 100       DO 230 J = 1,N
 101         YJ = Y(J)
 102         R = MAX(SRUR*ABS(YJ),R0/EWT(J))
 103         Y(J) = Y(J) + R
 104         FAC = -HL0/R
 105         CALL F (NEQ, TN, Y, FTEM)
 106         DO 220 I = 1,N
 107  220      WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
 108         Y(J) = YJ
 109         J1 = J1 + N
 110  230    CONTINUE
 111       NFE = NFE + N
 112 C Add identity matrix. -------------------------------------------------
 113  240  J = 3
 114       NP1 = N + 1
 115       DO 250 I = 1,N
 116         WM(J) = WM(J) + 1.0D0
 117  250    J = J + NP1
 118 C Do LU decomposition on P. --------------------------------------------
 119       CALL DGEFA (WM(3), N, N, IWM(21), IER)
 120       IF (IER .NE. 0) IERPJ = 1
 121       RETURN
 122 C If MITER = 3, construct a diagonal approximation to J and P. ---------
 123  300  WM(2) = HL0
 124       R = EL0*0.1D0
 125       DO 310 I = 1,N
 126  310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
 127       CALL F (NEQ, TN, Y, WM(3))
 128       NFE = NFE + 1
 129       DO 320 I = 1,N
 130         R0 = H*SAVF(I) - YH(I,2)
 131         DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I))
 132         WM(I+2) = 1.0D0
 133         IF (ABS(R0) .LT. UROUND/EWT(I)) GO TO 320
 134         IF (ABS(DI) .EQ. 0.0D0) GO TO 330
 135         WM(I+2) = 0.1D0*R0/DI
 136  320    CONTINUE
 137       RETURN
 138  330  IERPJ = 1
 139       RETURN
 140 C If MITER = 4, call JAC and multiply by scalar. -----------------------
 141  400  ML = IWM(1)
 142       MU = IWM(2)
 143       ML3 = ML + 3
 144       MBAND = ML + MU + 1
 145       MEBAND = MBAND + ML
 146       LENP = MEBAND*N
 147       DO 410 I = 1,LENP
 148  410    WM(I+2) = 0.0D0
 149       CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND)
 150       CON = -HL0
 151       DO 420 I = 1,LENP
 152  420    WM(I+2) = WM(I+2)*CON
 153       GO TO 570
 154 C If MITER = 5, make MBAND calls to F to approximate J. ----------------
 155  500  ML = IWM(1)
 156       MU = IWM(2)
 157       MBAND = ML + MU + 1
 158       MBA = MIN(MBAND,N)
 159       MEBAND = MBAND + ML
 160       MEB1 = MEBAND - 1
 161       SRUR = WM(1)
 162       FAC = DVNORM (N, SAVF, EWT)
 163       R0 = 1000.0D0*ABS(H)*UROUND*N*FAC
 164       IF (R0 .EQ. 0.0D0) R0 = 1.0D0
 165       DO 560 J = 1,MBA
 166         DO 530 I = J,N,MBAND
 167           YI = Y(I)
 168           R = MAX(SRUR*ABS(YI),R0/EWT(I))
 169  530      Y(I) = Y(I) + R
 170         CALL F (NEQ, TN, Y, FTEM)
 171         DO 550 JJ = J,N,MBAND
 172           Y(JJ) = YH(JJ,1)
 173           YJJ = Y(JJ)
 174           R = MAX(SRUR*ABS(YJJ),R0/EWT(JJ))
 175           FAC = -HL0/R
 176           I1 = MAX(JJ-MU,1)
 177           I2 = MIN(JJ+ML,N)
 178           II = JJ*MEB1 - ML + 2
 179           DO 540 I = I1,I2
 180  540        WM(II+I) = (FTEM(I) - SAVF(I))*FAC
 181  550      CONTINUE
 182  560    CONTINUE
 183       NFE = NFE + MBA
 184 C Add identity matrix. -------------------------------------------------
 185  570  II = MBAND + 2
 186       DO 580 I = 1,N
 187         WM(II) = WM(II) + 1.0D0
 188  580    II = II + MEBAND
 189 C Do LU decomposition of P. --------------------------------------------
 190       CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
 191       IF (IER .NE. 0) IERPJ = 1
 192       RETURN
 193 C----------------------- END OF SUBROUTINE DPREPJ ----------------------
 194       END