share/odepack/fortran/dprjis.f

   1 *DECK DPRJIS
   2       SUBROUTINE DPRJIS (NEQ, Y, YH, NYH, EWT, RTEM, SAVR, S, WK, IWK,
   3      1   RES, JAC, ADDA)
   4       EXTERNAL RES, JAC, ADDA
   5       INTEGER NEQ, NYH, IWK
   6       DOUBLE PRECISION Y, YH, EWT, RTEM, SAVR, S, WK
   7       DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), RTEM(*),
   8      1   S(*), SAVR(*), 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 RLSS
  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/ RLSS(6),
  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, IMUL, IRES, J, JJ, JMAX, JMIN, K, KMAX, KMIN, NG
  32       DOUBLE PRECISION CON, FAC, HL0, R, SRUR
  33 C-----------------------------------------------------------------------
  34 C DPRJIS is called to compute and process the matrix
  35 C P = A - H*EL(1)*J, where J is an approximation to the Jacobian dr/dy,
  36 C where r = g(t,y) - A(t,y)*s.  J is computed by columns, either by
  37 C the user-supplied routine JAC if MITER = 1, or by finite differencing
  38 C if MITER = 2.  J is stored in WK, rescaled, and ADDA is called to
  39 C generate P.  The matrix P is subjected to LU decomposition in CDRV.
  40 C P and its LU decomposition are stored separately in WK.
  41 C
  42 C In addition to variables described previously, communication
  43 C with DPRJIS uses the following:
  44 C Y     = array containing predicted values on entry.
  45 C RTEM  = work array of length N (ACOR in DSTODI).
  46 C SAVR  = array containing r evaluated at predicted y. On output it
  47 C         contains the residual evaluated at current values of t and y.
  48 C S     = array containing predicted values of dy/dt (SAVF in DSTODI).
  49 C WK    = real work space for matrices.  On output it contains P and
  50 C         its sparse LU decomposition.  Storage of matrix elements
  51 C         starts at WK(3).
  52 C         WK also contains the following matrix-related data.
  53 C         WK(1) = SQRT(UROUND), used in numerical Jacobian increments.
  54 C IWK   = integer work space for matrix-related data, assumed to be
  55 C         equivalenced to WK.  In addition,  WK(IPRSP) and IWK(IPISP)
  56 C         are assumed to have identical locations.
  57 C EL0   = EL(1) (input).
  58 C IERPJ = output error flag (in COMMON).
  59 C         =  0 if no error.
  60 C         =  1 if zero pivot found in CDRV.
  61 C         = IRES (= 2 or 3) if RES returned IRES = 2 or 3.
  62 C         = -1 if insufficient storage for CDRV (should not occur).
  63 C         = -2 if other error found in CDRV (should not occur here).
  64 C JCUR  = output flag = 1 to indicate that the Jacobian matrix
  65 C         (or approximation) is now current.
  66 C This routine also uses other variables in Common.
  67 C-----------------------------------------------------------------------
  68       HL0 = H*EL0
  69       CON = -HL0
  70       JCUR = 1
  71       NJE = NJE + 1
  72       GO TO (100, 200), MITER
  73 C
  74 C If MITER = 1, call RES, then call JAC and ADDA for each column. ------
  75  100  IRES = 1
  76       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
  77       NFE = NFE + 1
  78       IF (IRES .GT. 1) GO TO 600
  79       KMIN = IWK(IPIAN)
  80       DO 130 J = 1,N
  81         KMAX = IWK(IPIAN+J)-1
  82         DO 110 I = 1,N
  83  110      RTEM(I) = 0.0D0
  84         CALL JAC (NEQ, TN, Y, S, J, IWK(IPIAN), IWK(IPJAN), RTEM)
  85         DO 120 I = 1,N
  86  120      RTEM(I) = RTEM(I)*CON
  87         CALL ADDA (NEQ, TN, Y, J, IWK(IPIAN), IWK(IPJAN), RTEM)
  88         DO 125 K = KMIN,KMAX
  89           I = IWK(IBJAN+K)
  90           WK(IBA+K) = RTEM(I)
  91  125      CONTINUE
  92         KMIN = KMAX + 1
  93  130    CONTINUE
  94       GO TO 290
  95 C
  96 C If MITER = 2, make NGP + 1 calls to RES to approximate J and P. ------
  97  200  CONTINUE
  98       IRES = -1
  99       CALL RES (NEQ, TN, Y, S, SAVR, IRES)
 100       NFE = NFE + 1
 101       IF (IRES .GT. 1) GO TO 600
 102       SRUR = WK(1)
 103       JMIN = IWK(IPIGP)
 104       DO 240 NG = 1,NGP
 105         JMAX = IWK(IPIGP+NG) - 1
 106         DO 210 J = JMIN,JMAX
 107           JJ = IWK(IBJGP+J)
 108           R = MAX(SRUR*ABS(Y(JJ)),0.01D0/EWT(JJ))
 109  210      Y(JJ) = Y(JJ) + R
 110         CALL RES (NEQ,TN,Y,S,RTEM,IRES)
 111         NFE = NFE + 1
 112         IF (IRES .GT. 1) GO TO 600
 113         DO 230 J = JMIN,JMAX
 114           JJ = IWK(IBJGP+J)
 115           Y(JJ) = YH(JJ,1)
 116           R = MAX(SRUR*ABS(Y(JJ)),0.01D0/EWT(JJ))
 117           FAC = -HL0/R
 118           KMIN = IWK(IBIAN+JJ)
 119           KMAX = IWK(IBIAN+JJ+1) - 1
 120           DO 220 K = KMIN,KMAX
 121             I = IWK(IBJAN+K)
 122             RTEM(I) = (RTEM(I) - SAVR(I))*FAC
 123  220        CONTINUE
 124         CALL ADDA (NEQ, TN, Y, JJ, IWK(IPIAN), IWK(IPJAN), RTEM)
 125         DO 225 K = KMIN,KMAX
 126           I = IWK(IBJAN+K)
 127           WK(IBA+K) = RTEM(I)
 128  225      CONTINUE
 129  230      CONTINUE
 130         JMIN = JMAX + 1
 131  240    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 C
 137 C Do numerical factorization of P matrix. ------------------------------
 138  290  NLU = NLU + 1
 139       IERPJ = 0
 140       DO 295 I = 1,N
 141  295    RTEM(I) = 0.0D0
 142       CALL CDRV (N,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN),
 143      1  WK(IPA),RTEM,RTEM,NSP,IWK(IPISP),WK(IPRSP),IESP,2,IYS)
 144       IF (IYS .EQ. 0) RETURN
 145       IMUL = (IYS - 1)/N
 146       IERPJ = -2
 147       IF (IMUL .EQ. 8) IERPJ = 1
 148       IF (IMUL .EQ. 10) IERPJ = -1
 149       RETURN
 150 C Error return for IRES = 2 or IRES = 3 return from RES. ---------------
 151  600  IERPJ = IRES
 152       RETURN
 153 C----------------------- End of Subroutine DPRJIS ----------------------
 154       END