1
subroutine prepji (neq, y, yh, nyh, ewt, rtem, savr, s, wm, iwm,
4
external res, jac, adda
7
1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
8
2 maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu
9
integer i, i1, i2, ier, ii, ires, j, j1, jj, lenp,
10
1 mba, mband, meb1, meband, ml, ml3, mu
11
double precision y, yh, ewt, rtem, savr, s, wm
12
double precision rowns,
13
1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
14
double precision con, fac, hl0, r, srur, yi, yj, yjj
15
dimension neq(1), y(1), yh(nyh,*), ewt(1), rtem(1),
16
1 s(1), savr(1), wm(*), iwm(*)
17
common /ls0001/ rowns(209),
18
2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
19
3 iownd(14), iowns(6),
20
4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
21
5 maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu
22
c-----------------------------------------------------------------------
23
c prepji is called by stodi to compute and process the matrix
24
c p = a - h*el(1)*j , where j is an approximation to the jacobian dr/dy,
25
c where r = g(t,y) - a(t,y)*s. here j is computed by the user-supplied
26
c routine jac if miter = 1 or 4, or by finite differencing if miter =
27
c 2 or 5. j is stored in wm, rescaled, and adda is called to generate
28
c p. p is then subjected to lu decomposition in preparation
29
c for later solution of linear systems with p as coefficient
30
c matrix. this is done by dgefa if miter = 1 or 2, and by
31
c dgbfa if miter = 4 or 5.
33
c in addition to variables described previously, communication
34
c with prepji uses the following..
35
c y = array containing predicted values on entry.
36
c rtem = work array of length n (acor in stodi).
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 stodi).
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 the band parameters
47
c ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
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 singular,
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, nre, and nje.
57
c-----------------------------------------------------------------------
62
go to (100, 200, 300, 400, 500), miter
63
c if miter = 1, call res, then jac, and multiply by scalar. ------------
65
call res (neq, tn, y, s, savr, ires)
67
if (ires .gt. 1) go to 600
71
call jac ( neq, tn, y, s, 0, 0, wm(3), n )
74
120 wm(i+2) = wm(i+2)*con
76
c if miter = 2, make n + 1 calls to res to approximate j. --------------
79
call res (neq, tn, y, s, savr, ires)
81
if (ires .gt. 1) go to 600
86
r = dmax1(srur*dabs(yj),0.01d0/ewt(j))
89
call res ( neq, tn, y, s, rtem, ires )
91
if (ires .gt. 1) go to 600
93
220 wm(i+j1) = (rtem(i) - savr(i))*fac
98
call res (neq, tn, y, s, savr, ires)
100
if (ires .gt. 1) go to 600
101
c add matrix a. --------------------------------------------------------
103
call adda(neq, tn, y, 0, 0, wm(3), n)
104
c do lu decomposition on p. --------------------------------------------
105
call dgefa (wm(3), n, n, iwm(21), ier)
106
if (ier .ne. 0) ierpj = 1
108
c dummy section for miter = 3
110
c if miter = 4, call res, then jac, and multiply by scalar. ------------
112
call res (neq, tn, y, s, savr, ires)
114
if (ires .gt. 1) go to 600
123
call jac ( neq, tn, y, s, ml, mu, wm(ml3), meband)
126
420 wm(i+2) = wm(i+2)*con
128
c if miter = 5, make ml + mu + 2 calls to res to approximate j. --------
131
call res (neq, tn, y, s, savr, ires)
133
if (ires .gt. 1) go to 600
145
r = dmax1(srur*dabs(yi),0.01d0/ewt(i))
147
call res ( neq, tn, y, s, rtem, ires)
149
if (ires .gt. 1) go to 600
150
do 550 jj = j,n,mband
153
r = dmax1(srur*dabs(yjj),0.01d0/ewt(jj))
157
ii = jj*meb1 - ml + 2
159
540 wm(ii+i) = (rtem(i) - savr(i))*fac
163
call res (neq, tn, y, s, savr, ires)
165
if (ires .gt. 1) go to 600
166
c add matrix a. --------------------------------------------------------
168
call adda(neq, tn, y, ml, mu, wm(ml3), meband)
169
c do lu decomposition of p. --------------------------------------------
170
call dgbfa (wm(3), meband, n, ml, mu, iwm(21), ier)
171
if (ier .ne. 0) ierpj = 1
173
c error return for ires = 2 or ires = 3 return from res. ---------------
176
c----------------------- end of subroutine prepji ----------------------
added
removed
Lib/integrate/odepack/prepji.f
