acfa0dac12e07ea3d1383fd37ef1cdea98d47f34
[scilab.git] / scilab / modules / differential_equations / src / fortran / prepji.f
1 C/MEMBR ADD NAME=PREPJI,SSI=0
2       subroutine prepji (neq, y, yh, nyh, ewt, rtem, savr, s, wm, iwm,
3      1   res, jac, adda)
4 clll. optimize
5       external res, jac, adda
6       integer neq, nyh, iwm
7       integer iownd, iowns,
8      1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
9      2   maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu
10       integer i, i1, i2, ier, ii, ires, j, j1, jj, lenp,
11      1   mba, mband, meb1, meband, ml, ml3, mu
12       double precision y, yh, ewt, rtem, savr, s, wm
13       double precision rownd, rowns,
14      1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
15       double precision con, fac, hl0, r, srur, yi, yj, yjj
16       dimension neq(*), y(*), yh(nyh,*), ewt(*), rtem(*),
17      1   s(*), savr(*), wm(*), iwm(*)
18       common /ls0001/ rownd, rowns(209),
19      2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
20      3   iownd(14), iowns(6),
21      4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
22      5   maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu
23       integer         iero
24       common /ierode/ iero
25 c-----------------------------------------------------------------------
26 c%purpose
27 c prepji is called by stodi to compute and process the matrix
28 c p = a - h*el(1)*j , where j is an approximation to the jacobian dr/dy,
29 c where r = g(t,y) - a(t,y)*s. here j is computed by the user-supplied
30 c routine jac if miter = 1 or 4, or by finite differencing if miter =
31 c 2 or 5. j is stored in wm, rescaled, and adda is called to generate
32 c p. p is then subjected to lu decomposition in preparation
33 c for later solution of linear systems with p as coefficient
34 c matrix.  this is done by dgefa if miter = 1 or 2, and by
35 c dgbfa if miter = 4 or 5.
36 c
37 c%additional parameters
38 c in addition to variables described previously, communication
39 c with prepji uses the following..
40 c y     = array containing predicted values on entry.
41 c rtem  = work array of length n (acor in stodi).
42 c savr  = array used for output only.  on output it contains the
43 c         residual evaluated at current values of t and y.
44 c s     = array containing predicted values of dy/dt (savf in stodi).
45 c wm    = real work space for matrices.  on output it contains the
46 c         lu decomposition of p.
47 c         storage of matrix elements starts at wm(3).
48 c         wm also contains the following matrix-related data..
49 c         wm(1) = sqrt(uround), used in numerical jacobian increments.
50 c iwm   = integer work space containing pivot information, starting at
51 c         iwm(21).  iwm also contains the band parameters
52 c         ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
53 c el0   = el(1) (input).
54 c ierpj = output error flag.
55 c         = 0 if no trouble occurred,
56 c         = 1 if the p matrix was found to be singular,
57 c         = ires (= 2 or 3) if res returned ires = 2 or 3.
58 c jcur  = output flag = 1 to indicate that the jacobian matrix
59 c         (or approximation) is now current.
60 c this routine also uses the common variables el0, h, tn, uround,
61 c miter, n, nre, and nje.
62 c!
63 c-----------------------------------------------------------------------
64       nje = nje + 1
65       hl0 = h*el0
66       ierpj = 0
67       jcur = 1
68       go to (100, 200, 300, 400, 500), miter
69 c if miter = 1, call res, then jac, and multiply by scalar. ------------
70  100  ires = 1
71       call res (neq, tn, y, s, savr, ires)
72       if(iero.gt.0) return
73       nre = nre + 1
74       if (ires .gt. 1) go to 600
75       lenp = n*n
76       do 110 i = 1,lenp
77  110    wm(i+2) = 0.0d+0
78       call jac ( neq, tn, y, s, 0, 0, wm(3), n )
79       if(iero.gt.0) return
80       con = -hl0
81       do 120 i = 1,lenp
82  120    wm(i+2) = wm(i+2)*con
83       go to 240
84 c if miter = 2, make n + 1 calls to res to approximate j. --------------
85  200  continue
86       ires = -1
87       call res (neq, tn, y, s, savr, ires)
88       if(iero.gt.0) return
89       nre = nre + 1
90       if (ires .gt. 1) go to 600
91       srur = wm(1)
92       j1 = 2
93       do 230 j = 1,n
94         yj = y(j)
95         r = max(srur*abs(yj),0.010d+0/ewt(j))
96         y(j) = y(j) + r
97         fac = -hl0/r
98         call res ( neq, tn, y, s, rtem, ires )
99         if(iero.gt.0) return
100         nre = nre + 1
101         if (ires .gt. 1) go to 600
102         do 220 i = 1,n
103  220      wm(i+j1) = (rtem(i) - savr(i))*fac
104         y(j) = yj
105         j1 = j1 + n
106  230    continue
107       ires = 1
108       call res (neq, tn, y, s, savr, ires)
109       if(iero.gt.0) return
110       nre = nre + 1
111       if (ires .gt. 1) go to 600
112 c add matrix a. --------------------------------------------------------
113  240  continue
114       call adda(neq, tn, y, 0, 0, wm(3), n)
115       if(iero.gt.0) return
116 c do lu decomposition on p. --------------------------------------------
117       call dgefa (wm(3), n, n, iwm(21), ier)
118       if (ier .ne. 0) ierpj = 1
119       return
120 c dummy section for miter = 3
121  300  return
122 c if miter = 4, call res, then jac, and multiply by scalar. ------------
123  400  ires = 1
124       call res (neq, tn, y, s, savr, ires)
125       if(iero.gt.0) return
126       nre = nre + 1
127       if (ires .gt. 1) go to 600
128       ml = iwm(1)
129       mu = iwm(2)
130 cc mod 06-01-89
131 cc      ml3 = ml + 3
132       ml3 = 3
133 cc fin
134       mband = ml + mu + 1
135       meband = mband + ml
136       lenp = meband*n
137       do 410 i = 1,lenp
138  410    wm(i+2) = 0.0d+0
139       call jac ( neq, tn, y, s, ml, mu, wm(ml3), meband)
140       if(iero.gt.0) return
141       con = -hl0
142       do 420 i = 1,lenp
143  420    wm(i+2) = wm(i+2)*con
144       go to 570
145 c if miter = 5, make ml + mu + 2 calls to res to approximate j. --------
146  500  continue
147       ires = -1
148       call res (neq, tn, y, s, savr, ires)
149       if(iero.gt.0) return
150       nre = nre + 1
151       if (ires .gt. 1) go to 600
152       ml = iwm(1)
153       mu = iwm(2)
154       ml3 = ml + 3
155       mband = ml + mu + 1
156       mba = min(mband,n)
157       meband = mband + ml
158       meb1 = meband - 1
159       srur = wm(1)
160       do 560 j = 1,mba
161         do 530 i = j,n,mband
162           yi = y(i)
163           r = max(srur*abs(yi),0.010d+0/ewt(i))
164  530      y(i) = y(i) + r
165         call res ( neq, tn, y, s, rtem, ires)
166         if(iero.gt.0) return
167         nre = nre + 1
168         if (ires .gt. 1) go to 600
169         do 550 jj = j,n,mband
170           y(jj) = yh(jj,1)
171           yjj = y(jj)
172           r = max(srur*abs(yjj),0.010d+0/ewt(jj))
173           fac = -hl0/r
174           i1 = max(jj-mu,1)
175           i2 = min(jj+ml,n)
176           ii = jj*meb1 - ml + 2
177           do 540 i = i1,i2
178  540        wm(ii+i) = (rtem(i) - savr(i))*fac
179  550      continue
180  560    continue
181       ires = 1
182       call res (neq, tn, y, s, savr, ires)
183       if(iero.gt.0) return
184       nre = nre + 1
185       if (ires .gt. 1) go to 600
186 c add matrix a. --------------------------------------------------------
187  570  continue
188       call adda(neq, tn, y, ml, mu, wm(ml3), meband)
189       if(iero.gt.0) return
190 c do lu decomposition of p. --------------------------------------------
191       call dgbfa (wm(3), meband, n, ml, mu, iwm(21), ier)
192       if (ier .ne. 0) ierpj = 1
193       return
194 c error return for ires = 2 or ires = 3 return from res. ---------------
195  600  ierpj = ires
196       return
197 c----------------------- end of subroutine prepji ----------------------
198       end