-
Notifications
You must be signed in to change notification settings - Fork 0
/
UnitLegendMin.pas
125 lines (125 loc) · 4.7 KB
/
UnitLegendMin.pas
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
meband, ml, ml3, mu, np1
double precision y, yh, ewt, ftem, savf, wm
double precision rowns,
1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
double precision con, di, fac, hl0, r, r0, srur, yi, yj, yjj,
1 vnorm
dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
1 wm(*), iwm(*)
common /ls0001/ rowns(209),
2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
3 iownd(14), iowns(6),
4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
c-----------------------------------------------------------------------
c prepj is called by stode to compute and process the matrix
c p = i - h*el(1)*j , where j is an approximation to the jacobian.
c here j is computed by the user-supplied routine jac if
c miter = 1 or 4, or by finite differencing if miter = 2, 3, or 5.
c if miter = 3, a diagonal approximation to j is used.
c j is stored in wm and replaced by p. if miter .ne. 3, p is then
c subjected to lu decomposition in preparation for later solution
c of linear systems with p as coefficient matrix. this is done
c by dgefa if miter = 1 or 2, and by dgbfa if miter = 4 or 5.
c
c in addition to variables described previously, communication
c with prepj uses the following..
c y = array containing predicted values on entry.
c ftem = work array of length n (acor in stode).
c savf = array containing f evaluated at predicted y.
c wm = real work space for matrices. on output it contains the
c inverse diagonal matrix if miter = 3 and the lu decomposition
c of p if miter is 1, 2 , 4, or 5.
c storage of matrix elements starts at wm(3).
c wm also contains the following matrix-related data..
c wm(1) = sqrt(uround), used in numerical jacobian increments.
c wm(2) = h*el0, saved for later use if miter = 3.
c iwm = integer work space containing pivot information, starting at
c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band
c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
c el0 = el(1) (input).
c ierpj = output error flag, = 0 if no trouble, .gt. 0 if
c p matrix found to be singular.
c jcur = output flag = 1 to indicate that the jacobian matrix
c (or approximation) is now current.
c this routine also uses the common variables el0, h, tn, uround,
c miter, n, nfe, and nje.
c-----------------------------------------------------------------------
nje = nje + 1
ierpj = 0
jcur = 1
hl0 = h*el0
go to (100, 200, 300, 400, 500), miter
c if miter = 1, call jac and multiply by scalar. -----------------------
100 lenp = n*n
do 110 i = 1,lenp
110 wm(i+2) = 0.0d0
call jac (neq, tn, y, 0, 0, wm(3), n)
con = -hl0
do 120 i = 1,lenp
120 wm(i+2) = wm(i+2)*con
go to 240
c if miter = 2, make n calls to f to approximate j. --------------------
200 fac = vnorm (n, savf, ewt)
r0 = 1000.0d0*dabs(h)*uround*dfloat(n)*fac
if (r0 .eq. 0.0d0) r0 = 1.0d0
srur = wm(1)
j1 = 2
do 230 j = 1,n
yj = y(j)
r = dmax1(srur*dabs(yj),r0/ewt(j))
y(j) = y(j) + r
fac = -hl0/r
call f (neq, tn, y, ftem)
do 220 i = 1,n
220 wm(i+j1) = (ftem(i) - savf(i))*fac
y(j) = yj
j1 = j1 + n
230 continue
nfe = nfe + n
c add identity matrix. -------------------------------------------------
240 j = 3
np1 = n + 1
do 250 i = 1,n
wm(j) = wm(j) + 1.0d0
250 j = j + np1
c do lu decomposition on p. --------------------------------------------
call dgefa (wm(3), n, n, iwm(21), ier)
if (ier .ne. 0) ierpj = 1
return
c if miter = 3, construct a diagonal approximation to j and p. ---------
300 wm(2) = hl0
r = el0*0.1d0
do 310 i = 1,n
310 y(i) = y(i) + r*(h*savf(i) - yh(i,2))
call f (neq, tn, y, wm(3))
nfe = nfe + 1
do 320 i = 1,n
r0 = h*savf(i) - yh(i,2)
di = 0.1d0*r0 - h*(wm(i+2) - savf(i))
wm(i+2) = 1.0d0
if (dabs(r0) .lt. uround/ewt(i)) go to 320
if (dabs(di) .eq. 0.0d0) go to 330
wm(i+2) = 0.1d0*r0/di
320 continue
return
330 ierpj = 1
return
c if miter = 4, call jac and multiply by scalar. -----------------------
400 ml = iwm(1)
mu = iwm(2)
ml3 = ml + 3
mband = ml + mu + 1
meband = mband + ml
lenp = meband*n
do 410 i = 1,lenp
410 wm(i+2) = 0.0d0
call jac (neq, tn, y, ml, mu, wm(ml3), meband)
con = -hl0
do 420 i = 1,lenp
420 wm(i+2) = wm(i+2)*con
go to 570
c if miter = 5, make mband calls to f to approximate j. ----------------
500 ml = iwm(1)
mu = iwm(2)
mband