2
SUBROUTINE DPRJIS (NEQ, Y, YH, NYH, EWT, RTEM, SAVR, S, WK, IWK,
4
EXTERNAL RES, JAC, ADDA
6
DOUBLE PRECISION Y, YH, EWT, RTEM, SAVR, S, WK
7
DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), RTEM(*),
8
1 S(*), SAVR(*), WK(*), IWK(*)
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
20
COMMON /DLS001/ ROWNS(209),
21
1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
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.
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
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).
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-----------------------------------------------------------------------
72
GO TO (100, 200), MITER
74
C If MITER = 1, call RES, then call JAC and ADDA for each column. ------
76
CALL RES (NEQ, TN, Y, S, SAVR, IRES)
78
IF (IRES .GT. 1) GO TO 600
84
CALL JAC (NEQ, TN, Y, S, J, IWK(IPIAN), IWK(IPJAN), RTEM)
86
120 RTEM(I) = RTEM(I)*CON
87
CALL ADDA (NEQ, TN, Y, J, IWK(IPIAN), IWK(IPJAN), RTEM)
96
C If MITER = 2, make NGP + 1 calls to RES to approximate J and P. ------
99
CALL RES (NEQ, TN, Y, S, SAVR, IRES)
101
IF (IRES .GT. 1) GO TO 600
105
JMAX = IWK(IPIGP+NG) - 1
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)
112
IF (IRES .GT. 1) GO TO 600
116
R = MAX(SRUR*ABS(Y(JJ)),0.01D0/EWT(JJ))
119
KMAX = IWK(IBIAN+JJ+1) - 1
122
RTEM(I) = (RTEM(I) - SAVR(I))*FAC
124
CALL ADDA (NEQ, TN, Y, JJ, IWK(IPIAN), IWK(IPJAN), RTEM)
133
CALL RES (NEQ, TN, Y, S, SAVR, IRES)
135
IF (IRES .GT. 1) GO TO 600
137
C Do numerical factorization of P matrix. ------------------------------
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
147
IF (IMUL .EQ. 8) IERPJ = 1
148
IF (IMUL .EQ. 10) IERPJ = -1
150
C Error return for IRES = 2 or IRES = 3 return from RES. ---------------
153
C----------------------- End of Subroutine DPRJIS ----------------------
added
removed
packages/odepack/dprjis.f
