1C/MEMBR ADD NAME=PREPJI,SSI=0
2      subroutine prepji (neq, y, yh, nyh, ewt, rtem, savr, s, wm, iwm,
3     1   res, jac, adda)
4clll. optimize
5
6      external res, jac, adda
7      integer neq, nyh, iwm
8      integer iownd, iowns,
9     1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
10     2   maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu
11      integer i, i1, i2, ier, ii, ires, j, j1, jj, lenp,
12     1   mba, mband, meb1, meband, ml, ml3, mu
13      double precision y, yh, ewt, rtem, savr, s, wm
14      double precision rownd, rowns,
15     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
16      double precision con, fac, hl0, r, srur, yi, yj, yjj
17      dimension neq(*), y(*), yh(nyh,*), ewt(*), rtem(*),
18     1   s(*), savr(*), wm(*), iwm(*)
19cDEC$ ATTRIBUTES DLLIMPORT:: /ls0001/
20      common /ls0001/ rownd, rowns(209),
21     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
22     3   iownd(14), iowns(6),
23     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
24     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu
25c-----------------------------------------------------------------------
26c%purpose
27c prepji is called by stodi to compute and process the matrix
28c p = a - h*el(1)*j , where j is an approximation to the jacobian dr/dy,
29c where r = g(t,y) - a(t,y)*s. here j is computed by the user-supplied
30c routine jac if miter = 1 or 4, or by finite differencing if miter =
31c 2 or 5. j is stored in wm, rescaled, and adda is called to generate
32c p. p is then subjected to lu decomposition in preparation
33c for later solution of linear systems with p as coefficient
34c matrix.  this is done by dgefa if miter = 1 or 2, and by
35c dgbfa if miter = 4 or 5.
36c
37c%additional parameters
38c in addition to variables described previously, communication
39c with prepji uses the following..
40c y     = array containing predicted values on entry.
41c rtem  = work array of length n (acor in stodi).
42c savr  = array used for output only.  on output it contains the
43c         residual evaluated at current values of t and y.
44c s     = array containing predicted values of dy/dt (savf in stodi).
45c wm    = real work space for matrices.  on output it contains the
46c         lu decomposition of p.
47c         storage of matrix elements starts at wm(3).
48c         wm also contains the following matrix-related data..
49c         wm(1) = sqrt(uround), used in numerical jacobian increments.
50c iwm   = integer work space containing pivot information, starting at
51c         iwm(21).  iwm also contains the band parameters
52c         ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
53c el0   = el(1) (input).
54c ierpj = output error flag.
55c         = 0 if no trouble occurred,
56c         = 1 if the p matrix was found to be singular,
57c         = ires (= 2 or 3) if res returned ires = 2 or 3.
58c jcur  = output flag = 1 to indicate that the jacobian matrix
59c         (or approximation) is now current.
60c this routine also uses the common variables el0, h, tn, uround,
61c miter, n, nre, and nje.
62c!
63c-----------------------------------------------------------------------
64      nje = nje + 1
65      hl0 = h*el0
66      ierpj = 0
67      jcur = 1
68      go to (100, 200, 300, 400, 500), miter
69c 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(ierror.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(ierror.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
84c 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(ierror.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(ierror.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(ierror.gt.0) return
110      nre = nre + 1
111      if (ires .gt. 1) go to 600
112c add matrix a. --------------------------------------------------------
113 240  continue
114      call adda(neq, tn, y, 0, 0, wm(3), n)
115      if(ierror.gt.0) return
116c do lu decomposition on p. --------------------------------------------
117      call dgefa (wm(3), n, n, iwm(21), ier)
118      if (ier .ne. 0) ierpj = 1
119      return
120c dummy section for miter = 3
121 300  return
122c 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(ierror.gt.0) return
126      nre = nre + 1
127      if (ires .gt. 1) go to 600
128      ml = iwm(1)
129      mu = iwm(2)
130cc mod 06-01-89
131cc      ml3 = ml + 3
132      ml3 = 3
133cc 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(ierror.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
145c 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(ierror.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(ierror.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(ierror.gt.0) return
184      nre = nre + 1
185      if (ires .gt. 1) go to 600
186c add matrix a. --------------------------------------------------------
187 570  continue
188      call adda(neq, tn, y, ml, mu, wm(ml3), meband)
189      if(ierror.gt.0) return
190c 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
194c error return for ires = 2 or ires = 3 return from res. ---------------
195 600  ierpj = ires
196      return
197c----------------------- end of subroutine prepji ----------------------
198      end
199