module_cbmz_lsodes_solver.F
References to this file elsewhere.
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! CBMZ module: see module_cbmz.F for information and terms of use
8 !**********************************************************************************
9 module module_cbmz_lsodes_solver
10
11 contains
12
13
14 !ZZ
15 !
16 ! Obtained Oct 11, 1994 from ODEPACK in NETLIB by RDS
17 subroutine lsodes_solver ( &
18 f, neq, y, t, tout, itol, rtol, atol, itask, &
19 istate, iopt, rwork, lrw, iwork, liw, jac, mf, &
20 ruserpar, nruserpar, iuserpar, niuserpar )
21 external f, jac
22 integer neq, itol, itask, istate, iopt, lrw, iwork, liw, mf
23 integer nruserpar, iuserpar, niuserpar
24 real y, t, tout, rtol, atol, rwork
25 real ruserpar
26 !jdf dimension neq(1), y(1), rtol(1), atol(1), rwork(lrw), iwork(liw)
27 dimension neq(*), y(*), rtol(*), atol(*), rwork(lrw), iwork(liw)
28 dimension ruserpar(nruserpar), iuserpar(niuserpar)
29 !-----------------------------------------------------------------------
30 ! this is the march 30, 1987 version of
31 ! lsodes.. livermore solver for ordinary differential equations
32 ! with general sparse jacobian matrices.
33 ! this version is in single precision.
34 !
35 ! lsodes solves the initial value problem for stiff or nonstiff
36 ! systems of first order ode-s,
37 ! dy/dt = f(t,y) , or, in component form,
38 ! dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(neq)) (i = 1,...,neq).
39 ! lsodes is a variant of the lsode package, and is intended for
40 ! problems in which the jacobian matrix df/dy has an arbitrary
41 ! sparse structure (when the problem is stiff).
42 !
43 ! authors.. alan c. hindmarsh,
44 ! computing and mathematics research division, l-316
45 ! lawrence livermore national laboratory
46 ! livermore, ca 94550.
47 !
48 ! and andrew h. sherman
49 ! j. s. nolen and associates
50 ! houston, tx 77084
51 !-----------------------------------------------------------------------
52 ! references..
53 ! 1. alan c. hindmarsh, odepack, a systematized collection of ode
54 ! solvers, in scientific computing, r. s. stepleman et al. (eds.),
55 ! north-holland, amsterdam, 1983, pp. 55-64.
56 !
57 ! 2. s. c. eisenstat, m. c. gursky, m. h. schultz, and a. h. sherman,
58 ! yale sparse matrix package.. i. the symmetric codes,
59 ! int. j. num. meth. eng., 18 (1982), pp. 1145-1151.
60 !
61 ! 3. s. c. eisenstat, m. c. gursky, m. h. schultz, and a. h. sherman,
62 ! yale sparse matrix package.. ii. the nonsymmetric codes,
63 ! research report no. 114, dept. of computer sciences, yale
64 ! university, 1977.
65 !-----------------------------------------------------------------------
66 ! summary of usage.
67 !
68 ! communication between the user and the lsodes package, for normal
69 ! situations, is summarized here. this summary describes only a subset
70 ! of the full set of options available. see the full description for
71 ! details, including optional communication, nonstandard options,
72 ! and instructions for special situations. see also the example
73 ! problem (with program and output) following this summary.
74 !
75 ! a. first provide a subroutine of the form..
76 ! subroutine f (neq, t, y, ydot)
77 ! dimension y(neq), ydot(neq)
78 ! which supplies the vector function f by loading ydot(i) with f(i).
79 !
80 ! b. next determine (or guess) whether or not the problem is stiff.
81 ! stiffness occurs when the jacobian matrix df/dy has an eigenvalue
82 ! whose real part is negative and large in magnitude, compared to the
83 ! reciprocal of the t span of interest. if the problem is nonstiff,
84 ! use a method flag mf = 10. if it is stiff, there are two standard
85 ! for the method flag, mf = 121 and mf = 222. in both cases, lsodes
86 ! requires the jacobian matrix in some form, and it treats this matrix
87 ! in general sparse form, with sparsity structure determined internally.
88 ! (for options where the user supplies the sparsity structure, see
89 ! the full description of mf below.)
90 !
91 ! c. if the problem is stiff, you are encouraged to supply the jacobian
92 ! directly (mf = 121), but if this is not feasible, lsodes will
93 ! compute it internally by difference quotients (mf = 222).
94 ! if you are supplying the jacobian, provide a subroutine of the form..
95 ! subroutine jac (neq, t, y, j, ian, jan, pdj)
96 ! dimension y(1), ian(1), jan(1), pdj(1)
97 ! here neq, t, y, and j are input arguments, and the jac routine is to
98 ! load the array pdj (of length neq) with the j-th column of df/dy.
99 ! i.e., load pdj(i) with df(i)/dy(j) for all relevant values of i.
100 ! the arguments ian and jan should be ignored for normal situations.
101 ! lsodes will call the jac routine with j = 1,2,...,neq.
102 ! only nonzero elements need be loaded. usually, a crude approximation
103 ! to df/dy, possibly with fewer nonzero elements, will suffice.
104 !
105 ! d. write a main program which calls subroutine lsodes once for
106 ! each point at which answers are desired. this should also provide
107 ! for possible use of logical unit 6 for output of error messages
108 ! by lsodes. on the first call to lsodes, supply arguments as follows..
109 ! f = name of subroutine for right-hand side vector f.
110 ! this name must be declared external in calling program.
111 ! neq = number of first order ode-s.
112 ! y = array of initial values, of length neq.
113 ! t = the initial value of the independent variable.
114 ! tout = first point where output is desired (.ne. t).
115 ! itol = 1 or 2 according as atol (below) is a scalar or array.
116 ! rtol = relative tolerance parameter (scalar).
117 ! atol = absolute tolerance parameter (scalar or array).
118 ! the estimated local error in y(i) will be controlled so as
119 ! to be roughly less (in magnitude) than
120 ! ewt(i) = rtol*abs(y(i)) + atol if itol = 1, or
121 ! ewt(i) = rtol*abs(y(i)) + atol(i) if itol = 2.
122 ! thus the local error test passes if, in each component,
123 ! either the absolute error is less than atol (or atol(i)),
124 ! or the relative error is less than rtol.
125 ! use rtol = 0.0 for pure absolute error control, and
126 ! use atol = 0.0 (or atol(i) = 0.0) for pure relative error
127 ! control. caution.. actual (global) errors may exceed these
128 ! local tolerances, so choose them conservatively.
129 ! itask = 1 for normal computation of output values of y at t = tout.
130 ! istate = integer flag (input and output). set istate = 1.
131 ! iopt = 0 to indicate no optional inputs used.
132 ! rwork = real work array of length at least..
133 ! 20 + 16*neq for mf = 10,
134 ! 20 + (2 + 1./lenrat)*nnz + (11 + 9./lenrat)*neq
135 ! for mf = 121 or 222,
136 ! where..
137 ! nnz = the number of nonzero elements in the sparse
138 ! jacobian (if this is unknown, use an estimate), and
139 ! lenrat = the real to integer wordlength ratio (usually 1 in
140 ! single precision and 2 in double precision).
141 ! in any case, the required size of rwork cannot generally
142 ! be predicted in advance if mf = 121 or 222, and the value
143 ! above is a rough estimate of a crude lower bound. some
144 ! experimentation with this size may be necessary.
145 ! (when known, the correct required length is an optional
146 ! output, available in iwork(17).)
147 ! lrw = declared length of rwork (in user-s dimension).
148 ! iwork = integer work array of length at least 30.
149 ! liw = declared length of iwork (in user-s dimension).
150 ! jac = name of subroutine for jacobian matrix (mf = 121).
151 ! if used, this name must be declared external in calling
152 ! program. if not used, pass a dummy name.
153 ! mf = method flag. standard values are..
154 ! 10 for nonstiff (adams) method, no jacobian used.
155 ! 121 for stiff (bdf) method, user-supplied sparse jacobian.
156 ! 222 for stiff method, internally generated sparse jacobian.
157 ! note that the main program must declare arrays y, rwork, iwork,
158 ! and possibly atol.
159 !
160 ! e. the output from the first call (or any call) is..
161 ! y = array of computed values of y(t) vector.
162 ! t = corresponding value of independent variable (normally tout).
163 ! istate = 2 if lsodes was successful, negative otherwise.
164 ! -1 means excess work done on this call (perhaps wrong mf).
165 ! -2 means excess accuracy requested (tolerances too small).
166 ! -3 means illegal input detected (see printed message).
167 ! -4 means repeated error test failures (check all inputs).
168 ! -5 means repeated convergence failures (perhaps bad jacobian
169 ! supplied or wrong choice of mf or tolerances).
170 ! -6 means error weight became zero during problem. (solution
171 ! component i vanished, and atol or atol(i) = 0.)
172 ! -7 means a fatal error return flag came from the sparse
173 ! solver cdrv by way of prjs or slss. should never happen.
174 ! a return with istate = -1, -4, or -5 may result from using
175 ! an inappropriate sparsity structure, one that is quite
176 ! different from the initial structure. consider calling
177 ! lsodes again with istate = 3 to force the structure to be
178 ! reevaluated. see the full description of istate below.
179 !
180 ! f. to continue the integration after a successful return, simply
181 ! reset tout and call lsodes again. no other parameters need be reset.
182 !
183 !-----------------------------------------------------------------------
184 ! example problem.
185 !
186 ! the following is a simple example problem, with the coding
187 ! needed for its solution by lsodes. the problem is from chemical
188 ! kinetics, and consists of the following 12 rate equations..
189 ! dy1/dt = -rk1*y1
190 ! dy2/dt = rk1*y1 + rk11*rk14*y4 + rk19*rk14*y5
191 ! - rk3*y2*y3 - rk15*y2*y12 - rk2*y2
192 ! dy3/dt = rk2*y2 - rk5*y3 - rk3*y2*y3 - rk7*y10*y3
193 ! + rk11*rk14*y4 + rk12*rk14*y6
194 ! dy4/dt = rk3*y2*y3 - rk11*rk14*y4 - rk4*y4
195 ! dy5/dt = rk15*y2*y12 - rk19*rk14*y5 - rk16*y5
196 ! dy6/dt = rk7*y10*y3 - rk12*rk14*y6 - rk8*y6
197 ! dy7/dt = rk17*y10*y12 - rk20*rk14*y7 - rk18*y7
198 ! dy8/dt = rk9*y10 - rk13*rk14*y8 - rk10*y8
199 ! dy9/dt = rk4*y4 + rk16*y5 + rk8*y6 + rk18*y7
200 ! dy10/dt = rk5*y3 + rk12*rk14*y6 + rk20*rk14*y7
201 ! + rk13*rk14*y8 - rk7*y10*y3 - rk17*y10*y12
202 ! - rk6*y10 - rk9*y10
203 ! dy11/dt = rk10*y8
204 ! dy12/dt = rk6*y10 + rk19*rk14*y5 + rk20*rk14*y7
205 ! - rk15*y2*y12 - rk17*y10*y12
206 !
207 ! with rk1 = rk5 = 0.1, rk4 = rk8 = rk16 = rk18 = 2.5,
208 ! rk10 = 5.0, rk2 = rk6 = 10.0, rk14 = 30.0,
209 ! rk3 = rk7 = rk9 = rk11 = rk12 = rk13 = rk19 = rk20 = 50.0,
210 ! rk15 = rk17 = 100.0.
211 !
212 ! the t interval is from 0 to 1000, and the initial conditions
213 ! are y1 = 1, y2 = y3 = ... = y12 = 0. the problem is stiff.
214 !
215 ! the following coding solves this problem with lsodes, using mf = 121
216 ! and printing results at t = .1, 1., 10., 100., 1000. it uses
217 ! itol = 1 and mixed relative/absolute tolerance controls.
218 ! during the run and at the end, statistical quantities of interest
219 ! are printed (see optional outputs in the full description below).
220 !
221 ! external fex, jex
222 ! dimension y(12), rwork(500), iwork(30)
223 ! data lrw/500/, liw/30/
224 ! neq = 12
225 ! do 10 i = 1,neq
226 ! 10 y(i) = 0.0e0
227 ! y(1) = 1.0e0
228 ! t = 0.0e0
229 ! tout = 0.1e0
230 ! itol = 1
231 ! rtol = 1.0e-4
232 ! atol = 1.0e-6
233 ! itask = 1
234 ! istate = 1
235 ! iopt = 0
236 ! mf = 121
237 ! do 40 iout = 1,5
238 ! call lsodes (fex, neq, y, t, tout, itol, rtol, atol,
239 ! 1 itask, istate, iopt, rwork, lrw, iwork, liw, jex, mf)
240 ! write(6,30)t,iwork(11),rwork(11),(y(i),i=1,neq)
241 ! 30 format(//7h at t =,e11.3,4x,
242 ! 1 12h no. steps =,i5,4x,12h last step =,e11.3/
243 ! 2 13h y array = ,4e14.5/13x,4e14.5/13x,4e14.5)
244 ! if (istate .lt. 0) go to 80
245 ! tout = tout*10.0e0
246 ! 40 continue
247 ! lenrw = iwork(17)
248 ! leniw = iwork(18)
249 ! nst = iwork(11)
250 ! nfe = iwork(12)
251 ! nje = iwork(13)
252 ! nlu = iwork(21)
253 ! nnz = iwork(19)
254 ! nnzlu = iwork(25) + iwork(26) + neq
255 ! write (6,70) lenrw,leniw,nst,nfe,nje,nlu,nnz,nnzlu
256 ! 70 format(//22h required rwork size =,i4,15h iwork size =,i4/
257 ! 1 12h no. steps =,i4,12h no. f-s =,i4,12h no. j-s =,i4,
258 ! 2 13h no. lu-s =,i4/23h no. of nonzeros in j =,i5,
259 ! 3 26h no. of nonzeros in lu =,i5)
260 ! stop
261 ! 80 write(6,90)istate
262 ! 90 format(///22h error halt.. istate =,i3)
263 ! stop
264 ! end
265 !
266 ! subroutine fex (neq, t, y, ydot)
267 ! real t, y, ydot
268 ! real rk1, rk2, rk3, rk4, rk5, rk6, rk7, rk8, rk9,
269 ! 1 rk10, rk11, rk12, rk13, rk14, rk15, rk16, rk17
270 ! dimension y(12), ydot(12)
271 ! data rk1/0.1e0/, rk2/10.0e0/, rk3/50.0e0/, rk4/2.5e0/, rk5/0.1e0/,
272 ! 1 rk6/10.0e0/, rk7/50.0e0/, rk8/2.5e0/, rk9/50.0e0/, rk10/5.0e0/,
273 ! 2 rk11/50.0e0/, rk12/50.0e0/, rk13/50.0e0/, rk14/30.0e0/,
274 ! 3 rk15/100.0e0/, rk16/2.5e0/, rk17/100.0e0/, rk18/2.5e0/,
275 ! 4 rk19/50.0e0/, rk20/50.0e0/
276 ! ydot(1) = -rk1*y(1)
277 ! ydot(2) = rk1*y(1) + rk11*rk14*y(4) + rk19*rk14*y(5)
278 ! 1 - rk3*y(2)*y(3) - rk15*y(2)*y(12) - rk2*y(2)
279 ! ydot(3) = rk2*y(2) - rk5*y(3) - rk3*y(2)*y(3) - rk7*y(10)*y(3)
280 ! 1 + rk11*rk14*y(4) + rk12*rk14*y(6)
281 ! ydot(4) = rk3*y(2)*y(3) - rk11*rk14*y(4) - rk4*y(4)
282 ! ydot(5) = rk15*y(2)*y(12) - rk19*rk14*y(5) - rk16*y(5)
283 ! ydot(6) = rk7*y(10)*y(3) - rk12*rk14*y(6) - rk8*y(6)
284 ! ydot(7) = rk17*y(10)*y(12) - rk20*rk14*y(7) - rk18*y(7)
285 ! ydot(8) = rk9*y(10) - rk13*rk14*y(8) - rk10*y(8)
286 ! ydot(9) = rk4*y(4) + rk16*y(5) + rk8*y(6) + rk18*y(7)
287 ! ydot(10) = rk5*y(3) + rk12*rk14*y(6) + rk20*rk14*y(7)
288 ! 1 + rk13*rk14*y(8) - rk7*y(10)*y(3) - rk17*y(10)*y(12)
289 ! 2 - rk6*y(10) - rk9*y(10)
290 ! ydot(11) = rk10*y(8)
291 ! ydot(12) = rk6*y(10) + rk19*rk14*y(5) + rk20*rk14*y(7)
292 ! 1 - rk15*y(2)*y(12) - rk17*y(10)*y(12)
293 ! return
294 ! end
295 !
296 ! subroutine jex (neq, t, y, j, ia, ja, pdj)
297 ! real t, y, pdj
298 ! real rk1, rk2, rk3, rk4, rk5, rk6, rk7, rk8, rk9,
299 ! 1 rk10, rk11, rk12, rk13, rk14, rk15, rk16, rk17
300 ! dimension y(1), ia(1), ja(1), pdj(1)
301 ! data rk1/0.1e0/, rk2/10.0e0/, rk3/50.0e0/, rk4/2.5e0/, rk5/0.1e0/,
302 ! 1 rk6/10.0e0/, rk7/50.0e0/, rk8/2.5e0/, rk9/50.0e0/, rk10/5.0e0/,
303 ! 2 rk11/50.0e0/, rk12/50.0e0/, rk13/50.0e0/, rk14/30.0e0/,
304 ! 3 rk15/100.0e0/, rk16/2.5e0/, rk17/100.0e0/, rk18/2.5e0/,
305 ! 4 rk19/50.0e0/, rk20/50.0e0/
306 ! go to (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), j
307 ! 1 pdj(1) = -rk1
308 ! pdj(2) = rk1
309 ! return
310 ! 2 pdj(2) = -rk3*y(3) - rk15*y(12) - rk2
311 ! pdj(3) = rk2 - rk3*y(3)
312 ! pdj(4) = rk3*y(3)
313 ! pdj(5) = rk15*y(12)
314 ! pdj(12) = -rk15*y(12)
315 ! return
316 ! 3 pdj(2) = -rk3*y(2)
317 ! pdj(3) = -rk5 - rk3*y(2) - rk7*y(10)
318 ! pdj(4) = rk3*y(2)
319 ! pdj(6) = rk7*y(10)
320 ! pdj(10) = rk5 - rk7*y(10)
321 ! return
322 ! 4 pdj(2) = rk11*rk14
323 ! pdj(3) = rk11*rk14
324 ! pdj(4) = -rk11*rk14 - rk4
325 ! pdj(9) = rk4
326 ! return
327 ! 5 pdj(2) = rk19*rk14
328 ! pdj(5) = -rk19*rk14 - rk16
329 ! pdj(9) = rk16
330 ! pdj(12) = rk19*rk14
331 ! return
332 ! 6 pdj(3) = rk12*rk14
333 ! pdj(6) = -rk12*rk14 - rk8
334 ! pdj(9) = rk8
335 ! pdj(10) = rk12*rk14
336 ! return
337 ! 7 pdj(7) = -rk20*rk14 - rk18
338 ! pdj(9) = rk18
339 ! pdj(10) = rk20*rk14
340 ! pdj(12) = rk20*rk14
341 ! return
342 ! 8 pdj(8) = -rk13*rk14 - rk10
343 ! pdj(10) = rk13*rk14
344 ! pdj(11) = rk10
345 ! 9 return
346 ! 10 pdj(3) = -rk7*y(3)
347 ! pdj(6) = rk7*y(3)
348 ! pdj(7) = rk17*y(12)
349 ! pdj(8) = rk9
350 ! pdj(10) = -rk7*y(3) - rk17*y(12) - rk6 - rk9
351 ! pdj(12) = rk6 - rk17*y(12)
352 ! 11 return
353 ! 12 pdj(2) = -rk15*y(2)
354 ! pdj(5) = rk15*y(2)
355 ! pdj(7) = rk17*y(10)
356 ! pdj(10) = -rk17*y(10)
357 ! pdj(12) = -rk15*y(2) - rk17*y(10)
358 ! return
359 ! end
360 !
361 ! the output of this program (on a cray-1 in single precision)
362 ! is as follows..
363 !
364 !
365 ! at t = 1.000e-01 no. steps = 12 last step = 1.515e-02
366 ! y array = 9.90050e-01 6.28228e-03 3.65313e-03 7.51934e-07
367 ! 1.12167e-09 1.18458e-09 1.77291e-12 3.26476e-07
368 ! 5.46720e-08 9.99500e-06 4.48483e-08 2.76398e-06
369 !
370 !
371 ! at t = 1.000e+00 no. steps = 33 last step = 7.880e-02
372 ! y array = 9.04837e-01 9.13105e-03 8.20622e-02 2.49177e-05
373 ! 1.85055e-06 1.96797e-06 1.46157e-07 2.39557e-05
374 ! 3.26306e-05 7.21621e-04 5.06433e-05 3.05010e-03
375 !
376 !
377 ! at t = 1.000e+01 no. steps = 48 last step = 1.239e+00
378 ! y array = 3.67876e-01 3.68958e-03 3.65133e-01 4.48325e-05
379 ! 6.10798e-05 4.33148e-05 5.90211e-05 1.18449e-04
380 ! 3.15235e-03 3.56531e-03 4.15520e-03 2.48741e-01
381 !
382 !
383 ! at t = 1.000e+02 no. steps = 91 last step = 3.764e+00
384 ! y array = 4.44981e-05 4.42666e-07 4.47273e-04 -3.53257e-11
385 ! 2.81577e-08 -9.67741e-11 2.77615e-07 1.45322e-07
386 ! 1.56230e-02 4.37394e-06 1.60104e-02 9.52246e-01
387 !
388 !
389 ! at t = 1.000e+03 no. steps = 111 last step = 4.156e+02
390 ! y array = -2.65492e-13 2.60539e-14 -8.59563e-12 6.29355e-14
391 ! -1.78066e-13 5.71471e-13 -1.47561e-12 4.58078e-15
392 ! 1.56314e-02 1.37878e-13 1.60184e-02 9.52719e-01
393 !
394 !
395 ! required rwork size = 442 iwork size = 30
396 ! no. steps = 111 no. f-s = 142 no. j-s = 2 no. lu-s = 20
397 ! no. of nonzeros in j = 44 no. of nonzeros in lu = 50
398 !-----------------------------------------------------------------------
399 ! full description of user interface to lsodes.
400 !
401 ! the user interface to lsodes consists of the following parts.
402 !
403 ! i. the call sequence to subroutine lsodes, which is a driver
404 ! routine for the solver. this includes descriptions of both
405 ! the call sequence arguments and of user-supplied routines.
406 ! following these descriptions is a description of
407 ! optional inputs available through the call sequence, and then
408 ! a description of optional outputs (in the work arrays).
409 !
410 ! ii. descriptions of other routines in the lsodes package that may be
411 ! (optionally) called by the user. these provide the ability to
412 ! alter error message handling, save and restore the internal
413 ! common, and obtain specified derivatives of the solution y(t).
414 !
415 ! iii. descriptions of common blocks to be declared in overlay
416 ! or similar environments, or to be saved when doing an interrupt
417 ! of the problem and continued solution later.
418 !
419 ! iv. description of two routines in the lsodes package, either of
420 ! which the user may replace with his own version, if desired.
421 ! these relate to the measurement of errors.
422 !
423 !-----------------------------------------------------------------------
424 ! part i. call sequence.
425 !
426 ! the call sequence parameters used for input only are
427 ! f, neq, tout, itol, rtol, atol, itask, iopt, lrw, liw, jac, mf,
428 ! and those used for both input and output are
429 ! y, t, istate.
430 ! the work arrays rwork and iwork are also used for conditional and
431 ! optional inputs and optional outputs. (the term output here refers
432 ! to the return from subroutine lsodes to the user-s calling program.)
433 !
434 ! the legality of input parameters will be thoroughly checked on the
435 ! initial call for the problem, but not checked thereafter unless a
436 ! change in input parameters is flagged by istate = 3 on input.
437 !
438 ! the descriptions of the call arguments are as follows.
439 !
440 ! f = the name of the user-supplied subroutine defining the
441 ! ode system. the system must be put in the first-order
442 ! form dy/dt = f(t,y), where f is a vector-valued function
443 ! of the scalar t and the vector y. subroutine f is to
444 ! compute the function f. it is to have the form
445 ! subroutine f (neq, t, y, ydot)
446 ! dimension y(1), ydot(1)
447 ! where neq, t, and y are input, and the array ydot = f(t,y)
448 ! is output. y and ydot are arrays of length neq.
449 ! (in the dimension statement above, 1 is a dummy
450 ! dimension.. it can be replaced by any value.)
451 ! subroutine f should not alter y(1),...,y(neq).
452 ! f must be declared external in the calling program.
453 !
454 ! subroutine f may access user-defined quantities in
455 ! neq(2),... and/or in y(neq(1)+1),... if neq is an array
456 ! (dimensioned in f) and/or y has length exceeding neq(1).
457 ! see the descriptions of neq and y below.
458 !
459 ! if quantities computed in the f routine are needed
460 ! externally to lsodes, an extra call to f should be made
461 ! for this purpose, for consistent and accurate results.
462 ! if only the derivative dy/dt is needed, use intdy instead.
463 !
464 ! neq = the size of the ode system (number of first order
465 ! ordinary differential equations). used only for input.
466 ! neq may be decreased, but not increased, during the problem.
467 ! if neq is decreased (with istate = 3 on input), the
468 ! remaining components of y should be left undisturbed, if
469 ! these are to be accessed in f and/or jac.
470 !
471 ! normally, neq is a scalar, and it is generally referred to
472 ! as a scalar in this user interface description. however,
473 ! neq may be an array, with neq(1) set to the system size.
474 ! (the lsodes package accesses only neq(1).) in either case,
475 ! this parameter is passed as the neq argument in all calls
476 ! to f and jac. hence, if it is an array, locations
477 ! neq(2),... may be used to store other integer data and pass
478 ! it to f and/or jac. subroutines f and/or jac must include
479 ! neq in a dimension statement in that case.
480 !
481 ! y = a real array for the vector of dependent variables, of
482 ! length neq or more. used for both input and output on the
483 ! first call (istate = 1), and only for output on other calls.
484 ! on the first call, y must contain the vector of initial
485 ! values. on output, y contains the computed solution vector,
486 ! evaluated at t. if desired, the y array may be used
487 ! for other purposes between calls to the solver.
488 !
489 ! this array is passed as the y argument in all calls to
490 ! f and jac. hence its length may exceed neq, and locations
491 ! y(neq+1),... may be used to store other real data and
492 ! pass it to f and/or jac. (the lsodes package accesses only
493 ! y(1),...,y(neq).)
494 !
495 ! t = the independent variable. on input, t is used only on the
496 ! first call, as the initial point of the integration.
497 ! on output, after each call, t is the value at which a
498 ! computed solution y is evaluated (usually the same as tout).
499 ! on an error return, t is the farthest point reached.
500 !
501 ! tout = the next value of t at which a computed solution is desired.
502 ! used only for input.
503 !
504 ! when starting the problem (istate = 1), tout may be equal
505 ! to t for one call, then should .ne. t for the next call.
506 ! for the initial t, an input value of tout .ne. t is used
507 ! in order to determine the direction of the integration
508 ! (i.e. the algebraic sign of the step sizes) and the rough
509 ! scale of the problem. integration in either direction
510 ! (forward or backward in t) is permitted.
511 !
512 ! if itask = 2 or 5 (one-step modes), tout is ignored after
513 ! the first call (i.e. the first call with tout .ne. t).
514 ! otherwise, tout is required on every call.
515 !
516 ! if itask = 1, 3, or 4, the values of tout need not be
517 ! monotone, but a value of tout which backs up is limited
518 ! to the current internal t interval, whose endpoints are
519 ! tcur - hu and tcur (see optional outputs, below, for
520 ! tcur and hu).
521 !
522 ! itol = an indicator for the type of error control. see
523 ! description below under atol. used only for input.
524 !
525 ! rtol = a relative error tolerance parameter, either a scalar or
526 ! an array of length neq. see description below under atol.
527 ! input only.
528 !
529 ! atol = an absolute error tolerance parameter, either a scalar or
530 ! an array of length neq. input only.
531 !
532 ! the input parameters itol, rtol, and atol determine
533 ! the error control performed by the solver. the solver will
534 ! control the vector e = (e(i)) of estimated local errors
535 ! in y, according to an inequality of the form
536 ! rms-norm of ( e(i)/ewt(i) ) .le. 1,
537 ! where ewt(i) = rtol(i)*abs(y(i)) + atol(i),
538 ! and the rms-norm (root-mean-square norm) here is
539 ! rms-norm(v) = sqrt(sum v(i)**2 / neq). here ewt = (ewt(i))
540 ! is a vector of weights which must always be positive, and
541 ! the values of rtol and atol should all be non-negative.
542 ! the following table gives the types (scalar/array) of
543 ! rtol and atol, and the corresponding form of ewt(i).
544 !
545 ! itol rtol atol ewt(i)
546 ! 1 scalar scalar rtol*abs(y(i)) + atol
547 ! 2 scalar array rtol*abs(y(i)) + atol(i)
548 ! 3 array scalar rtol(i)*abs(y(i)) + atol
549 ! 4 array array rtol(i)*abs(y(i)) + atol(i)
550 !
551 ! when either of these parameters is a scalar, it need not
552 ! be dimensioned in the user-s calling program.
553 !
554 ! if none of the above choices (with itol, rtol, and atol
555 ! fixed throughout the problem) is suitable, more general
556 ! error controls can be obtained by substituting
557 ! user-supplied routines for the setting of ewt and/or for
558 ! the norm calculation. see part iv below.
559 !
560 ! if global errors are to be estimated by making a repeated
561 ! run on the same problem with smaller tolerances, then all
562 ! components of rtol and atol (i.e. of ewt) should be scaled
563 ! down uniformly.
564 !
565 ! itask = an index specifying the task to be performed.
566 ! input only. itask has the following values and meanings.
567 ! 1 means normal computation of output values of y(t) at
568 ! t = tout (by overshooting and interpolating).
569 ! 2 means take one step only and return.
570 ! 3 means stop at the first internal mesh point at or
571 ! beyond t = tout and return.
572 ! 4 means normal computation of output values of y(t) at
573 ! t = tout but without overshooting t = tcrit.
574 ! tcrit must be input as rwork(1). tcrit may be equal to
575 ! or beyond tout, but not behind it in the direction of
576 ! integration. this option is useful if the problem
577 ! has a singularity at or beyond t = tcrit.
578 ! 5 means take one step, without passing tcrit, and return.
579 ! tcrit must be input as rwork(1).
580 !
581 ! note.. if itask = 4 or 5 and the solver reaches tcrit
582 ! (within roundoff), it will return t = tcrit (exactly) to
583 ! indicate this (unless itask = 4 and tout comes before tcrit,
584 ! in which case answers at t = tout are returned first).
585 !
586 ! istate = an index used for input and output to specify the
587 ! the state of the calculation.
588 !
589 ! on input, the values of istate are as follows.
590 ! 1 means this is the first call for the problem
591 ! (initializations will be done). see note below.
592 ! 2 means this is not the first call, and the calculation
593 ! is to continue normally, with no change in any input
594 ! parameters except possibly tout and itask.
595 ! (if itol, rtol, and/or atol are changed between calls
596 ! with istate = 2, the new values will be used but not
597 ! tested for legality.)
598 ! 3 means this is not the first call, and the
599 ! calculation is to continue normally, but with
600 ! a change in input parameters other than
601 ! tout and itask. changes are allowed in
602 ! neq, itol, rtol, atol, iopt, lrw, liw, mf,
603 ! the conditional inputs ia and ja,
604 ! and any of the optional inputs except h0.
605 ! in particular, if miter = 1 or 2, a call with istate = 3
606 ! will cause the sparsity structure of the problem to be
607 ! recomputed (or reread from ia and ja if moss = 0).
608 ! note.. a preliminary call with tout = t is not counted
609 ! as a first call here, as no initialization or checking of
610 ! input is done. (such a call is sometimes useful for the
611 ! purpose of outputting the initial conditions.)
612 ! thus the first call for which tout .ne. t requires
613 ! istate = 1 on input.
614 !
615 ! on output, istate has the following values and meanings.
616 ! 1 means nothing was done, as tout was equal to t with
617 ! istate = 1 on input. (however, an internal counter was
618 ! set to detect and prevent repeated calls of this type.)
619 ! 2 means the integration was performed successfully.
620 ! -1 means an excessive amount of work (more than mxstep
621 ! steps) was done on this call, before completing the
622 ! requested task, but the integration was otherwise
623 ! successful as far as t. (mxstep is an optional input
624 ! and is normally 500.) to continue, the user may
625 ! simply reset istate to a value .gt. 1 and call again
626 ! (the excess work step counter will be reset to 0).
627 ! in addition, the user may increase mxstep to avoid
628 ! this error return (see below on optional inputs).
629 ! -2 means too much accuracy was requested for the precision
630 ! of the machine being used. this was detected before
631 ! completing the requested task, but the integration
632 ! was successful as far as t. to continue, the tolerance
633 ! parameters must be reset, and istate must be set
634 ! to 3. the optional output tolsf may be used for this
635 ! purpose. (note.. if this condition is detected before
636 ! taking any steps, then an illegal input return
637 ! (istate = -3) occurs instead.)
638 ! -3 means illegal input was detected, before taking any
639 ! integration steps. see written message for details.
640 ! note.. if the solver detects an infinite loop of calls
641 ! to the solver with illegal input, it will cause
642 ! the run to stop.
643 ! -4 means there were repeated error test failures on
644 ! one attempted step, before completing the requested
645 ! task, but the integration was successful as far as t.
646 ! the problem may have a singularity, or the input
647 ! may be inappropriate.
648 ! -5 means there were repeated convergence test failures on
649 ! one attempted step, before completing the requested
650 ! task, but the integration was successful as far as t.
651 ! this may be caused by an inaccurate jacobian matrix,
652 ! if one is being used.
653 ! -6 means ewt(i) became zero for some i during the
654 ! integration. pure relative error control (atol(i)=0.0)
655 ! was requested on a variable which has now vanished.
656 ! the integration was successful as far as t.
657 ! -7 means a fatal error return flag came from the sparse
658 ! solver cdrv by way of prjs or slss (numerical
659 ! factorization or backsolve). this should never happen.
660 ! the integration was successful as far as t.
661 !
662 ! note.. an error return with istate = -1, -4, or -5 and with
663 ! miter = 1 or 2 may mean that the sparsity structure of the
664 ! problem has changed significantly since it was last
665 ! determined (or input). in that case, one can attempt to
666 ! complete the integration by setting istate = 3 on the next
667 ! call, so that a new structure determination is done.
668 !
669 ! note.. since the normal output value of istate is 2,
670 ! it does not need to be reset for normal continuation.
671 ! also, since a negative input value of istate will be
672 ! regarded as illegal, a negative output value requires the
673 ! user to change it, and possibly other inputs, before
674 ! calling the solver again.
675 !
676 ! iopt = an integer flag to specify whether or not any optional
677 ! inputs are being used on this call. input only.
678 ! the optional inputs are listed separately below.
679 ! iopt = 0 means no optional inputs are being used.
680 ! default values will be used in all cases.
681 ! iopt = 1 means one or more optional inputs are being used.
682 !
683 ! rwork = a work array used for a mixture of real (single precision)
684 ! and integer work space.
685 ! the length of rwork (in real words) must be at least
686 ! 20 + nyh*(maxord + 1) + 3*neq + lwm where
687 ! nyh = the initial value of neq,
688 ! maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
689 ! smaller value is given as an optional input),
690 ! lwm = 0 if miter = 0,
691 ! lwm = 2*nnz + 2*neq + (nnz+9*neq)/lenrat if miter = 1,
692 ! lwm = 2*nnz + 2*neq + (nnz+10*neq)/lenrat if miter = 2,
693 ! lwm = neq + 2 if miter = 3.
694 ! in the above formulas,
695 ! nnz = number of nonzero elements in the jacobian matrix.
696 ! lenrat = the real to integer wordlength ratio (usually 1 in
697 ! single precision and 2 in double precision).
698 ! (see the mf description for meth and miter.)
699 ! thus if maxord has its default value and neq is constant,
700 ! the minimum length of rwork is..
701 ! 20 + 16*neq for mf = 10,
702 ! 20 + 16*neq + lwm for mf = 11, 111, 211, 12, 112, 212,
703 ! 22 + 17*neq for mf = 13,
704 ! 20 + 9*neq for mf = 20,
705 ! 20 + 9*neq + lwm for mf = 21, 121, 221, 22, 122, 222,
706 ! 22 + 10*neq for mf = 23.
707 ! if miter = 1 or 2, the above formula for lwm is only a
708 ! crude lower bound. the required length of rwork cannot
709 ! be readily predicted in general, as it depends on the
710 ! sparsity structure of the problem. some experimentation
711 ! may be necessary.
712 !
713 ! the first 20 words of rwork are reserved for conditional
714 ! and optional inputs and optional outputs.
715 !
716 ! the following word in rwork is a conditional input..
717 ! rwork(1) = tcrit = critical value of t which the solver
718 ! is not to overshoot. required if itask is
719 ! 4 or 5, and ignored otherwise. (see itask.)
720 !
721 ! lrw = the length of the array rwork, as declared by the user.
722 ! (this will be checked by the solver.)
723 !
724 ! iwork = an integer work array. the length of iwork must be at least
725 ! 31 + neq + nnz if moss = 0 and miter = 1 or 2, or
726 ! 30 otherwise.
727 ! (nnz is the number of nonzero elements in df/dy.)
728 !
729 ! in lsodes, iwork is used only for conditional and
730 ! optional inputs and optional outputs.
731 !
732 ! the following two blocks of words in iwork are conditional
733 ! inputs, required if moss = 0 and miter = 1 or 2, but not
734 ! otherwise (see the description of mf for moss).
735 ! iwork(30+j) = ia(j) (j=1,...,neq+1)
736 ! iwork(31+neq+k) = ja(k) (k=1,...,nnz)
737 ! the two arrays ia and ja describe the sparsity structure
738 ! to be assumed for the jacobian matrix. ja contains the row
739 ! indices where nonzero elements occur, reading in columnwise
740 ! order, and ia contains the starting locations in ja of the
741 ! descriptions of columns 1,...,neq, in that order, with
742 ! ia(1) = 1. thus, for each column index j = 1,...,neq, the
743 ! values of the row index i in column j where a nonzero
744 ! element may occur are given by
745 ! i = ja(k), where ia(j) .le. k .lt. ia(j+1).
746 ! if nnz is the total number of nonzero locations assumed,
747 ! then the length of the ja array is nnz, and ia(neq+1) must
748 ! be nnz + 1. duplicate entries are not allowed.
749 !
750 ! liw = the length of the array iwork, as declared by the user.
751 ! (this will be checked by the solver.)
752 !
753 ! note.. the work arrays must not be altered between calls to lsodes
754 ! for the same problem, except possibly for the conditional and
755 ! optional inputs, and except for the last 3*neq words of rwork.
756 ! the latter space is used for internal scratch space, and so is
757 ! available for use by the user outside lsodes between calls, if
758 ! desired (but not for use by f or jac).
759 !
760 ! jac = name of user-supplied routine (miter = 1 or moss = 1) to
761 ! compute the jacobian matrix, df/dy, as a function of
762 ! the scalar t and the vector y. it is to have the form
763 ! subroutine jac (neq, t, y, j, ian, jan, pdj)
764 ! dimension y(1), ian(1), jan(1), pdj(1)
765 ! where neq, t, y, j, ian, and jan are input, and the array
766 ! pdj, of length neq, is to be loaded with column j
767 ! of the jacobian on output. thus df(i)/dy(j) is to be
768 ! loaded into pdj(i) for all relevant values of i.
769 ! here t and y have the same meaning as in subroutine f,
770 ! and j is a column index (1 to neq). ian and jan are
771 ! undefined in calls to jac for structure determination
772 ! (moss = 1). otherwise, ian and jan are structure
773 ! descriptors, as defined under optional outputs below, and
774 ! so can be used to determine the relevant row indices i, if
775 ! desired. (in the dimension statement above, 1 is a
776 ! dummy dimension.. it can be replaced by any value.)
777 ! jac need not provide df/dy exactly. a crude
778 ! approximation (possibly with greater sparsity) will do.
779 ! in any case, pdj is preset to zero by the solver,
780 ! so that only the nonzero elements need be loaded by jac.
781 ! calls to jac are made with j = 1,...,neq, in that order, and
782 ! each such set of calls is preceded by a call to f with the
783 ! same arguments neq, t, and y. thus to gain some efficiency,
784 ! intermediate quantities shared by both calculations may be
785 ! saved in a user common block by f and not recomputed by jac,
786 ! if desired. jac must not alter its input arguments.
787 ! jac must be declared external in the calling program.
788 ! subroutine jac may access user-defined quantities in
789 ! neq(2),... and/or in y(neq(1)+1),... if neq is an array
790 ! (dimensioned in jac) and/or y has length exceeding neq(1).
791 ! see the descriptions of neq and y above.
792 !
793 ! mf = the method flag. used only for input.
794 ! mf has three decimal digits-- moss, meth, miter--
795 ! mf = 100*moss + 10*meth + miter.
796 ! moss indicates the method to be used to obtain the sparsity
797 ! structure of the jacobian matrix if miter = 1 or 2..
798 ! moss = 0 means the user has supplied ia and ja
799 ! (see descriptions under iwork above).
800 ! moss = 1 means the user has supplied jac (see below)
801 ! and the structure will be obtained from neq
802 ! initial calls to jac.
803 ! moss = 2 means the structure will be obtained from neq+1
804 ! initial calls to f.
805 ! meth indicates the basic linear multistep method..
806 ! meth = 1 means the implicit adams method.
807 ! meth = 2 means the method based on backward
808 ! differentiation formulas (bdf-s).
809 ! miter indicates the corrector iteration method..
810 ! miter = 0 means functional iteration (no jacobian matrix
811 ! is involved).
812 ! miter = 1 means chord iteration with a user-supplied
813 ! sparse jacobian, given by subroutine jac.
814 ! miter = 2 means chord iteration with an internally
815 ! generated (difference quotient) sparse jacobian
816 ! (using ngp extra calls to f per df/dy value,
817 ! where ngp is an optional output described below.)
818 ! miter = 3 means chord iteration with an internally
819 ! generated diagonal jacobian approximation.
820 ! (using 1 extra call to f per df/dy evaluation).
821 ! if miter = 1 or moss = 1, the user must supply a subroutine
822 ! jac (the name is arbitrary) as described above under jac.
823 ! otherwise, a dummy argument can be used.
824 !
825 ! the standard choices for mf are..
826 ! mf = 10 for a nonstiff problem,
827 ! mf = 21 or 22 for a stiff problem with ia/ja supplied
828 ! (21 if jac is supplied, 22 if not),
829 ! mf = 121 for a stiff problem with jac supplied,
830 ! but not ia/ja,
831 ! mf = 222 for a stiff problem with neither ia/ja nor
832 ! jac supplied.
833 ! the sparseness structure can be changed during the
834 ! problem by making a call to lsodes with istate = 3.
835 !-----------------------------------------------------------------------
836 ! optional inputs.
837 !
838 ! the following is a list of the optional inputs provided for in the
839 ! call sequence. (see also part ii.) for each such input variable,
840 ! this table lists its name as used in this documentation, its
841 ! location in the call sequence, its meaning, and the default value.
842 ! the use of any of these inputs requires iopt = 1, and in that
843 ! case all of these inputs are examined. a value of zero for any
844 ! of these optional inputs will cause the default value to be used.
845 ! thus to use a subset of the optional inputs, simply preload
846 ! locations 5 to 10 in rwork and iwork to 0.0 and 0 respectively, and
847 ! then set those of interest to nonzero values.
848 !
849 ! name location meaning and default value
850 !
851 ! h0 rwork(5) the step size to be attempted on the first step.
852 ! the default value is determined by the solver.
853 !
854 ! hmax rwork(6) the maximum absolute step size allowed.
855 ! the default value is infinite.
856 !
857 ! hmin rwork(7) the minimum absolute step size allowed.
858 ! the default value is 0. (this lower bound is not
859 ! enforced on the final step before reaching tcrit
860 ! when itask = 4 or 5.)
861 !
862 ! seth rwork(8) the element threshhold for sparsity determination
863 ! when moss = 1 or 2. if the absolute value of
864 ! an estimated jacobian element is .le. seth, it
865 ! will be assumed to be absent in the structure.
866 ! the default value of seth is 0.
867 !
868 ! maxord iwork(5) the maximum order to be allowed. the default
869 ! value is 12 if meth = 1, and 5 if meth = 2.
870 ! if maxord exceeds the default value, it will
871 ! be reduced to the default value.
872 ! if maxord is changed during the problem, it may
873 ! cause the current order to be reduced.
874 !
875 ! mxstep iwork(6) maximum number of (internally defined) steps
876 ! allowed during one call to the solver.
877 ! the default value is 500.
878 !
879 ! mxhnil iwork(7) maximum number of messages printed (per problem)
880 ! warning that t + h = t on a step (h = step size).
881 ! this must be positive to result in a non-default
882 ! value. the default value is 10.
883 !-----------------------------------------------------------------------
884 ! optional outputs.
885 !
886 ! as optional additional output from lsodes, the variables listed
887 ! below are quantities related to the performance of lsodes
888 ! which are available to the user. these are communicated by way of
889 ! the work arrays, but also have internal mnemonic names as shown.
890 ! except where stated otherwise, all of these outputs are defined
891 ! on any successful return from lsodes, and on any return with
892 ! istate = -1, -2, -4, -5, or -6. on an illegal input return
893 ! (istate = -3), they will be unchanged from their existing values
894 ! (if any), except possibly for tolsf, lenrw, and leniw.
895 ! on any error return, outputs relevant to the error will be defined,
896 ! as noted below.
897 !
898 ! name location meaning
899 !
900 ! hu rwork(11) the step size in t last used (successfully).
901 !
902 ! hcur rwork(12) the step size to be attempted on the next step.
903 !
904 ! tcur rwork(13) the current value of the independent variable
905 ! which the solver has actually reached, i.e. the
906 ! current internal mesh point in t. on output, tcur
907 ! will always be at least as far as the argument
908 ! t, but may be farther (if interpolation was done).
909 !
910 ! tolsf rwork(14) a tolerance scale factor, greater than 1.0,
911 ! computed when a request for too much accuracy was
912 ! detected (istate = -3 if detected at the start of
913 ! the problem, istate = -2 otherwise). if itol is
914 ! left unaltered but rtol and atol are uniformly
915 ! scaled up by a factor of tolsf for the next call,
916 ! then the solver is deemed likely to succeed.
917 ! (the user may also ignore tolsf and alter the
918 ! tolerance parameters in any other way appropriate.)
919 !
920 ! nst iwork(11) the number of steps taken for the problem so far.
921 !
922 ! nfe iwork(12) the number of f evaluations for the problem so far,
923 ! excluding those for structure determination
924 ! (moss = 2).
925 !
926 ! nje iwork(13) the number of jacobian evaluations for the problem
927 ! so far, excluding those for structure determination
928 ! (moss = 1).
929 !
930 ! nqu iwork(14) the method order last used (successfully).
931 !
932 ! nqcur iwork(15) the order to be attempted on the next step.
933 !
934 ! imxer iwork(16) the index of the component of largest magnitude in
935 ! the weighted local error vector ( e(i)/ewt(i) ),
936 ! on an error return with istate = -4 or -5.
937 !
938 ! lenrw iwork(17) the length of rwork actually required.
939 ! this is defined on normal returns and on an illegal
940 ! input return for insufficient storage.
941 !
942 ! leniw iwork(18) the length of iwork actually required.
943 ! this is defined on normal returns and on an illegal
944 ! input return for insufficient storage.
945 !
946 ! nnz iwork(19) the number of nonzero elements in the jacobian
947 ! matrix, including the diagonal (miter = 1 or 2).
948 ! (this may differ from that given by ia(neq+1)-1
949 ! if moss = 0, because of added diagonal entries.)
950 !
951 ! ngp iwork(20) the number of groups of column indices, used in
952 ! difference quotient jacobian aproximations if
953 ! miter = 2. this is also the number of extra f
954 ! evaluations needed for each jacobian evaluation.
955 !
956 ! nlu iwork(21) the number of sparse lu decompositions for the
957 ! problem so far.
958 !
959 ! lyh iwork(22) the base address in rwork of the history array yh,
960 ! described below in this list.
961 !
962 ! ipian iwork(23) the base address of the structure descriptor array
963 ! ian, described below in this list.
964 !
965 ! ipjan iwork(24) the base address of the structure descriptor array
966 ! jan, described below in this list.
967 !
968 ! nzl iwork(25) the number of nonzero elements in the strict lower
969 ! triangle of the lu factorization used in the chord
970 ! iteration (miter = 1 or 2).
971 !
972 ! nzu iwork(26) the number of nonzero elements in the strict upper
973 ! triangle of the lu factorization used in the chord
974 ! iteration (miter = 1 or 2).
975 ! the total number of nonzeros in the factorization
976 ! is therefore nzl + nzu + neq.
977 !
978 ! the following four arrays are segments of the rwork array which
979 ! may also be of interest to the user as optional outputs.
980 ! for each array, the table below gives its internal name,
981 ! its base address, and its description.
982 ! for yh and acor, the base addresses are in rwork (a real array).
983 ! the integer arrays ian and jan are to be obtained by declaring an
984 ! integer array iwk and identifying iwk(1) with rwork(21), using either
985 ! an equivalence statement or a subroutine call. then the base
986 ! addresses ipian (of ian) and ipjan (of jan) in iwk are to be obtained
987 ! as optional outputs iwork(23) and iwork(24), respectively.
988 ! thus ian(1) is iwk(ipian), etc.
989 !
990 ! name base address description
991 !
992 ! ian ipian (in iwk) structure descriptor array of size neq + 1.
993 ! jan ipjan (in iwk) structure descriptor array of size nnz.
994 ! (see above) ian and jan together describe the sparsity
995 ! structure of the jacobian matrix, as used by
996 ! lsodes when miter = 1 or 2.
997 ! jan contains the row indices of the nonzero
998 ! locations, reading in columnwise order, and
999 ! ian contains the starting locations in jan of
1000 ! the descriptions of columns 1,...,neq, in
1001 ! that order, with ian(1) = 1. thus for each
1002 ! j = 1,...,neq, the row indices i of the
1003 ! nonzero locations in column j are
1004 ! i = jan(k), ian(j) .le. k .lt. ian(j+1).
1005 ! note that ian(neq+1) = nnz + 1.
1006 ! (if moss = 0, ian/jan may differ from the
1007 ! input ia/ja because of a different ordering
1008 ! in each column, and added diagonal entries.)
1009 !
1010 ! yh lyh the nordsieck history array, of size nyh by
1011 ! (optional (nqcur + 1), where nyh is the initial value
1012 ! output) of neq. for j = 0,1,...,nqcur, column j+1
1013 ! of yh contains hcur**j/factorial(j) times
1014 ! the j-th derivative of the interpolating
1015 ! polynomial currently representing the solution,
1016 ! evaluated at t = tcur. the base address lyh
1017 ! is another optional output, listed above.
1018 !
1019 ! acor lenrw-neq+1 array of size neq used for the accumulated
1020 ! corrections on each step, scaled on output
1021 ! to represent the estimated local error in y
1022 ! on the last step. this is the vector e in
1023 ! the description of the error control. it is
1024 ! defined only on a successful return from
1025 ! lsodes.
1026 !
1027 !-----------------------------------------------------------------------
1028 ! part ii. other routines callable.
1029 !
1030 ! the following are optional calls which the user may make to
1031 ! gain additional capabilities in conjunction with lsodes.
1032 ! (the routines xsetun and xsetf are designed to conform to the
1033 ! slatec error handling package.)
1034 !
1035 ! form of call function
1036 ! call xsetun(lun) set the logical unit number, lun, for
1037 ! output of messages from lsodes, if
1038 ! the default is not desired.
1039 ! the default value of lun is 6.
1040 !
1041 ! call xsetf(mflag) set a flag to control the printing of
1042 ! messages by lsodes.
1043 ! mflag = 0 means do not print. (danger..
1044 ! this risks losing valuable information.)
1045 ! mflag = 1 means print (the default).
1046 !
1047 ! either of the above calls may be made at
1048 ! any time and will take effect immediately.
1049 !
1050 ! call srcms(rsav,isav,job) saves and restores the contents of
1051 ! the internal common blocks used by
1052 ! lsodes (see part iii below).
1053 ! rsav must be a real array of length 224
1054 ! or more, and isav must be an integer
1055 ! array of length 75 or more.
1056 ! job=1 means save common into rsav/isav.
1057 ! job=2 means restore common from rsav/isav.
1058 ! srcms is useful if one is
1059 ! interrupting a run and restarting
1060 ! later, or alternating between two or
1061 ! more problems solved with lsodes.
1062 !
1063 ! call intdy(,,,,,) provide derivatives of y, of various
1064 ! (see below) orders, at a specified point t, if
1065 ! desired. it may be called only after
1066 ! a successful return from lsodes.
1067 !
1068 ! the detailed instructions for using intdy are as follows.
1069 ! the form of the call is..
1070 !
1071 ! lyh = iwork(22)
1072 ! call intdy (t, k, rwork(lyh), nyh, dky, iflag)
1073 !
1074 ! the input parameters are..
1075 !
1076 ! t = value of independent variable where answers are desired
1077 ! (normally the same as the t last returned by lsodes).
1078 ! for valid results, t must lie between tcur - hu and tcur.
1079 ! (see optional outputs for tcur and hu.)
1080 ! k = integer order of the derivative desired. k must satisfy
1081 ! 0 .le. k .le. nqcur, where nqcur is the current order
1082 ! (see optional outputs). the capability corresponding
1083 ! to k = 0, i.e. computing y(t), is already provided
1084 ! by lsodes directly. since nqcur .ge. 1, the first
1085 ! derivative dy/dt is always available with intdy.
1086 ! lyh = the base address of the history array yh, obtained
1087 ! as an optional output as shown above.
1088 ! nyh = column length of yh, equal to the initial value of neq.
1089 !
1090 ! the output parameters are..
1091 !
1092 ! dky = a real array of length neq containing the computed value
1093 ! of the k-th derivative of y(t).
1094 ! iflag = integer flag, returned as 0 if k and t were legal,
1095 ! -1 if k was illegal, and -2 if t was illegal.
1096 ! on an error return, a message is also written.
1097 !-----------------------------------------------------------------------
1098 ! part iii. common blocks.
1099 !
1100 ! if lsodes is to be used in an overlay situation, the user
1101 ! must declare, in the primary overlay, the variables in..
1102 ! (1) the call sequence to lsodes,
1103 ! (2) the three internal common blocks
1104 ! /ls0001/ of length 257 (218 single precision words
1105 ! followed by 39 integer words),
1106 ! /lss001/ of length 40 ( 6 single precision words
1107 ! followed by 34 integer words),
1108 ! /eh0001/ of length 2 (integer words).
1109 !
1110 ! if lsodes is used on a system in which the contents of internal
1111 ! common blocks are not preserved between calls, the user should
1112 ! declare the above three common blocks in his main program to insure
1113 ! that their contents are preserved.
1114 !
1115 ! if the solution of a given problem by lsodes is to be interrupted
1116 ! and then later continued, such as when restarting an interrupted run
1117 ! or alternating between two or more problems, the user should save,
1118 ! following the return from the last lsodes call prior to the
1119 ! interruption, the contents of the call sequence variables and the
1120 ! internal common blocks, and later restore these values before the
1121 ! next lsodes call for that problem. to save and restore the common
1122 ! blocks, use subroutine srcms (see part ii above).
1123 !
1124 ! note.. in this version of lsodes, there are two data statements,
1125 ! in subroutines lsodes and xerrwv, which load variables into these
1126 ! labeled common blocks. on some systems, it may be necessary to
1127 ! move these to a separate block data subprogram.
1128 !
1129 !-----------------------------------------------------------------------
1130 ! part iv. optionally replaceable solver routines.
1131 !
1132 ! below are descriptions of two routines in the lsodes package which
1133 ! relate to the measurement of errors. either routine can be
1134 ! replaced by a user-supplied version, if desired. however, since such
1135 ! a replacement may have a major impact on performance, it should be
1136 ! done only when absolutely necessary, and only with great caution.
1137 ! (note.. the means by which the package version of a routine is
1138 ! superseded by the user-s version may be system-dependent.)
1139 !
1140 ! (a) ewset.
1141 ! the following subroutine is called just before each internal
1142 ! integration step, and sets the array of error weights, ewt, as
1143 ! described under itol/rtol/atol above..
1144 ! subroutine ewset (neq, itol, rtol, atol, ycur, ewt)
1145 ! where neq, itol, rtol, and atol are as in the lsodes call sequence,
1146 ! ycur contains the current dependent variable vector, and
1147 ! ewt is the array of weights set by ewset.
1148 !
1149 ! if the user supplies this subroutine, it must return in ewt(i)
1150 ! (i = 1,...,neq) a positive quantity suitable for comparing errors
1151 ! in y(i) to. the ewt array returned by ewset is passed to the
1152 ! vnorm routine (see below), and also used by lsodes in the computation
1153 ! of the optional output imxer, the diagonal jacobian approximation,
1154 ! and the increments for difference quotient jacobians.
1155 !
1156 ! in the user-supplied version of ewset, it may be desirable to use
1157 ! the current values of derivatives of y. derivatives up to order nq
1158 ! are available from the history array yh, described above under
1159 ! optional outputs. in ewset, yh is identical to the ycur array,
1160 ! extended to nq + 1 columns with a column length of nyh and scale
1161 ! factors of h**j/factorial(j). on the first call for the problem,
1162 ! given by nst = 0, nq is 1 and h is temporarily set to 1.0.
1163 ! the quantities nq, nyh, h, and nst can be obtained by including
1164 ! in ewset the statements..
1165 ! common /ls0001/ rls(218),ils(39)
1166 ! nq = ils(35)
1167 ! nyh = ils(14)
1168 ! nst = ils(36)
1169 ! h = rls(212)
1170 ! thus, for example, the current value of dy/dt can be obtained as
1171 ! ycur(nyh+i)/h (i=1,...,neq) (and the division by h is
1172 ! unnecessary when nst = 0).
1173 !
1174 ! (b) vnorm.
1175 ! the following is a real function routine which computes the weighted
1176 ! root-mean-square norm of a vector v..
1177 ! d = vnorm (n, v, w)
1178 ! where..
1179 ! n = the length of the vector,
1180 ! v = real array of length n containing the vector,
1181 ! w = real array of length n containing weights,
1182 ! d = sqrt( (1/n) * sum(v(i)*w(i))**2 ).
1183 ! vnorm is called with n = neq and with w(i) = 1.0/ewt(i), where
1184 ! ewt is as set by subroutine ewset.
1185 !
1186 ! if the user supplies this function, it should return a non-negative
1187 ! value of vnorm suitable for use in the error control in lsodes.
1188 ! none of the arguments should be altered by vnorm.
1189 ! for example, a user-supplied vnorm routine might..
1190 ! -substitute a max-norm of (v(i)*w(i)) for the rms-norm, or
1191 ! -ignore some components of v in the norm, with the effect of
1192 ! suppressing the error control on those components of y.
1193 !-----------------------------------------------------------------------
1194 !-----------------------------------------------------------------------
1195 ! other routines in the lsodes package.
1196 !
1197 ! in addition to subroutine lsodes, the lsodes package includes the
1198 ! following subroutines and function routines..
1199 ! iprep acts as an iterface between lsodes and prep, and also does
1200 ! adjusting of work space pointers and work arrays.
1201 ! prep is called by iprep to compute sparsity and do sparse matrix
1202 ! preprocessing if miter = 1 or 2.
1203 ! jgroup is called by prep to compute groups of jacobian column
1204 ! indices for use when miter = 2.
1205 ! adjlr adjusts the length of required sparse matrix work space.
1206 ! it is called by prep.
1207 ! cntnzu is called by prep and counts the nonzero elements in the
1208 ! strict upper triangle of j + j-transpose, where j = df/dy.
1209 ! intdy computes an interpolated value of the y vector at t = tout.
1210 ! stode is the core integrator, which does one step of the
1211 ! integration and the associated error control.
1212 ! cfode sets all method coefficients and test constants.
1213 ! prjs computes and preprocesses the jacobian matrix j = df/dy
1214 ! and the newton iteration matrix p = i - h*l0*j.
1215 ! slss manages solution of linear system in chord iteration.
1216 ! ewset sets the error weight vector ewt before each step.
1217 ! vnorm computes the weighted r.m.s. norm of a vector.
1218 ! srcms is a user-callable routine to save and restore
1219 ! the contents of the internal common blocks.
1220 ! odrv constructs a reordering of the rows and columns of
1221 ! a matrix by the minimum degree algorithm. odrv is a
1222 ! driver routine which calls subroutines md, mdi, mdm,
1223 ! mdp, mdu, and sro. see ref. 2 for details. (the odrv
1224 ! module has been modified since ref. 2, however.)
1225 ! cdrv performs reordering, symbolic factorization, numerical
1226 ! factorization, or linear system solution operations,
1227 ! depending on a path argument ipath. cdrv is a
1228 ! driver routine which calls subroutines nroc, nsfc,
1229 ! nnfc, nnsc, and nntc. see ref. 3 for details.
1230 ! lsodes uses cdrv to solve linear systems in which the
1231 ! coefficient matrix is p = i - con*j, where i is the
1232 ! identity, con is a scalar, and j is an approximation to
1233 ! the jacobian df/dy. because cdrv deals with rowwise
1234 ! sparsity descriptions, cdrv works with p-transpose, not p.
1235 ! r1mach computes the unit roundoff in a machine-independent manner.
1236 ! xerrwv, xsetun, and xsetf handle the printing of all error
1237 ! messages and warnings. xerrwv is machine-dependent.
1238 ! note.. vnorm and r1mach are function routines.
1239 ! all the others are subroutines.
1240 !
1241 ! the intrinsic and external routines used by lsodes are..
1242 ! abs, amax1, amin1, float, max0, min0, mod, sign, sqrt, and write.
1243 !
1244 !-----------------------------------------------------------------------
1245 ! the following card is for optimized compilation on lll compilers.
1246 !lll. optimize
1247 !-----------------------------------------------------------------------
1248 !rce external prjs, slss
1249 integer illin, init, lyh, lewt, lacor, lsavf, lwm, liwm, &
1250 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns
1251 integer icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
1252 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1253 integer iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp, &
1254 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa, &
1255 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj, &
1256 nslj, ngp, nlu, nnz, nsp, nzl, nzu
1257 integer i, i1, i2, iflag, imax, imul, imxer, ipflag, ipgo, irem, &
1258 j, kgo, lenrat, lenyht, leniw, lenrw, lf0, lia, lja, &
1259 lrtem, lwtem, lyhd, lyhn, mf1, mord, mxhnl0, mxstp0, ncolm
1260 real rowns, &
1261 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1262 real con0, conmin, ccmxj, psmall, rbig, seth
1263 !rce real atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli, &
1264 !rce tcrit, tdist, tnext, tol, tolsf, tp, size, sum, w0, &
1265 !rce r1mach, vnorm
1266 real atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli, &
1267 tcrit, tdist, tnext, tol, tolsf, tp, size, sum, w0
1268 dimension mord(2)
1269 logical ihit
1270 !-----------------------------------------------------------------------
1271 ! the following two internal common blocks contain
1272 ! (a) variables which are local to any subroutine but whose values must
1273 ! be preserved between calls to the routine (own variables), and
1274 ! (b) variables which are communicated between subroutines.
1275 ! the structure of each block is as follows.. all real variables are
1276 ! listed first, followed by all integers. within each type, the
1277 ! variables are grouped with those local to subroutine lsodes first,
1278 ! then those local to subroutine stode or subroutine prjs
1279 ! (no other routines have own variables), and finally those used
1280 ! for communication. the block ls0001 is declared in subroutines
1281 ! lsodes, iprep, prep, intdy, stode, prjs, and slss. the block lss001
1282 ! is declared in subroutines lsodes, iprep, prep, prjs, and slss.
1283 ! groups of variables are replaced by dummy arrays in the common
1284 ! declarations in routines where those variables are not used.
1285 !-----------------------------------------------------------------------
1286 common /ls0001/ rowns(209), &
1287 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, &
1288 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm, &
1289 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6), &
1290 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
1291 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1292 !
1293 common /lss001/ con0, conmin, ccmxj, psmall, rbig, seth, &
1294 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp, &
1295 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa, &
1296 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj, &
1297 nslj, ngp, nlu, nnz, nsp, nzl, nzu
1298 !
1299 data mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
1300 !raz data illin/0/, ntrep/0/
1301 !-----------------------------------------------------------------------
1302 ! in the data statement below, set lenrat equal to the ratio of
1303 ! the wordlength for a real number to that for an integer. usually,
1304 ! lenrat = 1 for single precision and 2 for double precision. if the
1305 ! true ratio is not an integer, use the next smaller integer (.ge. 1).
1306 !-----------------------------------------------------------------------
1307 data lenrat/1/
1308 !-----------------------------------------------------------------------
1309 ! block a.
1310 ! this code block is executed on every call.
1311 ! it tests istate and itask for legality and branches appropriately.
1312 ! if istate .gt. 1 but the flag init shows that initialization has
1313 ! not yet been done, an error return occurs.
1314 ! if istate = 1 and tout = t, jump to block g and return immediately.
1315 !-----------------------------------------------------------------------
1316 if (istate .lt. 1 .or. istate .gt. 3) go to 601
1317 if (itask .lt. 1 .or. itask .gt. 5) go to 602
1318 if (istate .eq. 1) go to 10
1319 if (init .eq. 0) go to 603
1320 if (istate .eq. 2) go to 200
1321 go to 20
1322 10 init = 0
1323 if (tout .eq. t) go to 430
1324 20 ntrep = 0
1325 !-----------------------------------------------------------------------
1326 ! block b.
1327 ! the next code block is executed for the initial call (istate = 1),
1328 ! or for a continuation call with parameter changes (istate = 3).
1329 ! it contains checking of all inputs and various initializations.
1330 ! if istate = 1, the final setting of work space pointers, the matrix
1331 ! preprocessing, and other initializations are done in block c.
1332 !
1333 ! first check legality of the non-optional inputs neq, itol, iopt,
1334 ! mf, ml, and mu.
1335 !-----------------------------------------------------------------------
1336 if (neq(1) .le. 0) go to 604
1337 if (istate .eq. 1) go to 25
1338 if (neq(1) .gt. n) go to 605
1339 25 n = neq(1)
1340 if (itol .lt. 1 .or. itol .gt. 4) go to 606
1341 if (iopt .lt. 0 .or. iopt .gt. 1) go to 607
1342 moss = mf/100
1343 mf1 = mf - 100*moss
1344 meth = mf1/10
1345 miter = mf1 - 10*meth
1346 if (moss .lt. 0 .or. moss .gt. 2) go to 608
1347 if (meth .lt. 1 .or. meth .gt. 2) go to 608
1348 if (miter .lt. 0 .or. miter .gt. 3) go to 608
1349 if (miter .eq. 0 .or. miter .eq. 3) moss = 0
1350 ! next process and check the optional inputs. --------------------------
1351 if (iopt .eq. 1) go to 40
1352 maxord = mord(meth)
1353 mxstep = mxstp0
1354 mxhnil = mxhnl0
1355 if (istate .eq. 1) h0 = 0.0e0
1356 hmxi = 0.0e0
1357 hmin = 0.0e0
1358 seth = 0.0e0
1359 go to 60
1360 40 maxord = iwork(5)
1361 if (maxord .lt. 0) go to 611
1362 if (maxord .eq. 0) maxord = 100
1363 maxord = min0(maxord,mord(meth))
1364 mxstep = iwork(6)
1365 if (mxstep .lt. 0) go to 612
1366 if (mxstep .eq. 0) mxstep = mxstp0
1367 mxhnil = iwork(7)
1368 if (mxhnil .lt. 0) go to 613
1369 if (mxhnil .eq. 0) mxhnil = mxhnl0
1370 if (istate .ne. 1) go to 50
1371 h0 = rwork(5)
1372 if ((tout - t)*h0 .lt. 0.0e0) go to 614
1373 50 hmax = rwork(6)
1374 if (hmax .lt. 0.0e0) go to 615
1375 hmxi = 0.0e0
1376 if (hmax .gt. 0.0e0) hmxi = 1.0e0/hmax
1377 hmin = rwork(7)
1378 if (hmin .lt. 0.0e0) go to 616
1379 seth = rwork(8)
1380 if (seth .lt. 0.0e0) go to 609
1381 ! check rtol and atol for legality. ------------------------------------
1382 60 rtoli = rtol(1)
1383 atoli = atol(1)
1384 do 65 i = 1,n
1385 if (itol .ge. 3) rtoli = rtol(i)
1386 if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
1387 if (rtoli .lt. 0.0e0) go to 619
1388 if (atoli .lt. 0.0e0) go to 620
1389 65 continue
1390 !-----------------------------------------------------------------------
1391 ! compute required work array lengths, as far as possible, and test
1392 ! these against lrw and liw. then set tentative pointers for work
1393 ! arrays. pointers to rwork/iwork segments are named by prefixing l to
1394 ! the name of the segment. e.g., the segment yh starts at rwork(lyh).
1395 ! segments of rwork (in order) are denoted wm, yh, savf, ewt, acor.
1396 ! if miter = 1 or 2, the required length of the matrix work space wm
1397 ! is not yet known, and so a crude minimum value is used for the
1398 ! initial tests of lrw and liw, and yh is temporarily stored as far
1399 ! to the right in rwork as possible, to leave the maximum amount
1400 ! of space for wm for matrix preprocessing. thus if miter = 1 or 2
1401 ! and moss .ne. 2, some of the segments of rwork are temporarily
1402 ! omitted, as they are not needed in the preprocessing. these
1403 ! omitted segments are.. acor if istate = 1, ewt and acor if istate = 3
1404 ! and moss = 1, and savf, ewt, and acor if istate = 3 and moss = 0.
1405 !-----------------------------------------------------------------------
1406 lrat = lenrat
1407 if (istate .eq. 1) nyh = n
1408 lwmin = 0
1409 if (miter .eq. 1) lwmin = 4*n + 10*n/lrat
1410 if (miter .eq. 2) lwmin = 4*n + 11*n/lrat
1411 if (miter .eq. 3) lwmin = n + 2
1412 lenyh = (maxord+1)*nyh
1413 lrest = lenyh + 3*n
1414 lenrw = 20 + lwmin + lrest
1415 iwork(17) = lenrw
1416 leniw = 30
1417 if (moss .eq. 0 .and. miter .ne. 0 .and. miter .ne. 3) &
1418 leniw = leniw + n + 1
1419 iwork(18) = leniw
1420 if (lenrw .gt. lrw) go to 617
1421 if (leniw .gt. liw) go to 618
1422 lia = 31
1423 if (moss .eq. 0 .and. miter .ne. 0 .and. miter .ne. 3) &
1424 leniw = leniw + iwork(lia+n) - 1
1425 iwork(18) = leniw
1426 if (leniw .gt. liw) go to 618
1427 lja = lia + n + 1
1428 lia = min0(lia,liw)
1429 lja = min0(lja,liw)
1430 lwm = 21
1431 if (istate .eq. 1) nq = 1
1432 ncolm = min0(nq+1,maxord+2)
1433 lenyhm = ncolm*nyh
1434 lenyht = lenyh
1435 if (miter .eq. 1 .or. miter .eq. 2) lenyht = lenyhm
1436 imul = 2
1437 if (istate .eq. 3) imul = moss
1438 if (moss .eq. 2) imul = 3
1439 lrtem = lenyht + imul*n
1440 lwtem = lwmin
1441 if (miter .eq. 1 .or. miter .eq. 2) lwtem = lrw - 20 - lrtem
1442 lenwk = lwtem
1443 lyhn = lwm + lwtem
1444 lsavf = lyhn + lenyht
1445 lewt = lsavf + n
1446 lacor = lewt + n
1447 istatc = istate
1448 if (istate .eq. 1) go to 100
1449 !-----------------------------------------------------------------------
1450 ! istate = 3. move yh to its new location.
1451 ! note that only the part of yh needed for the next step, namely
1452 ! min(nq+1,maxord+2) columns, is actually moved.
1453 ! a temporary error weight array ewt is loaded if moss = 2.
1454 ! sparse matrix processing is done in iprep/prep if miter = 1 or 2.
1455 ! if maxord was reduced below nq, then the pointers are finally set
1456 ! so that savf is identical to yh(*,maxord+2).
1457 !-----------------------------------------------------------------------
1458 lyhd = lyh - lyhn
1459 imax = lyhn - 1 + lenyhm
1460 ! move yh. branch for move right, no move, or move left. --------------
1461 if (lyhd) 70,80,74
1462 70 do 72 i = lyhn,imax
1463 j = imax + lyhn - i
1464 72 rwork(j) = rwork(j+lyhd)
1465 go to 80
1466 74 do 76 i = lyhn,imax
1467 76 rwork(i) = rwork(i+lyhd)
1468 80 lyh = lyhn
1469 iwork(22) = lyh
1470 if (miter .eq. 0 .or. miter .eq. 3) go to 92
1471 if (moss .ne. 2) go to 85
1472 ! temporarily load ewt if miter = 1 or 2 and moss = 2. -----------------
1473 call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1474 do 82 i = 1,n
1475 if (rwork(i+lewt-1) .le. 0.0e0) go to 621
1476 82 rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
1477 85 continue
1478 ! iprep and prep do sparse matrix preprocessing if miter = 1 or 2. -----
1479 lsavf = min0(lsavf,lrw)
1480 lewt = min0(lewt,lrw)
1481 lacor = min0(lacor,lrw)
1482 call iprep (neq, y, rwork, iwork(lia), iwork(lja), ipflag, f, jac, &
1483 ruserpar, nruserpar, iuserpar, niuserpar)
1484 lenrw = lwm - 1 + lenwk + lrest
1485 iwork(17) = lenrw
1486 if (ipflag .ne. -1) iwork(23) = ipian
1487 if (ipflag .ne. -1) iwork(24) = ipjan
1488 ipgo = -ipflag + 1
1489 go to (90, 628, 629, 630, 631, 632, 633), ipgo
1490 90 iwork(22) = lyh
1491 if (lenrw .gt. lrw) go to 617
1492 ! set flag to signal parameter changes to stode. -----------------------
1493 92 jstart = -1
1494 if (n .eq. nyh) go to 200
1495 ! neq was reduced. zero part of yh to avoid undefined references. -----
1496 i1 = lyh + l*nyh
1497 i2 = lyh + (maxord + 1)*nyh - 1
1498 if (i1 .gt. i2) go to 200
1499 do 95 i = i1,i2
1500 95 rwork(i) = 0.0e0
1501 go to 200
1502 !-----------------------------------------------------------------------
1503 ! block c.
1504 ! the next block is for the initial call only (istate = 1).
1505 ! it contains all remaining initializations, the initial call to f,
1506 ! the sparse matrix preprocessing (miter = 1 or 2), and the
1507 ! calculation of the initial step size.
1508 ! the error weights in ewt are inverted after being loaded.
1509 !-----------------------------------------------------------------------
1510 100 continue
1511 lyh = lyhn
1512 iwork(22) = lyh
1513 tn = t
1514 nst = 0
1515 h = 1.0e0
1516 nnz = 0
1517 ngp = 0
1518 nzl = 0
1519 nzu = 0
1520 ! load the initial value vector in yh. ---------------------------------
1521 do 105 i = 1,n
1522 105 rwork(i+lyh-1) = y(i)
1523 ! initial call to f. (lf0 points to yh(*,2).) -------------------------
1524 lf0 = lyh + nyh
1525 call f (neq, t, y, rwork(lf0), &
1526 ruserpar, nruserpar, iuserpar, niuserpar)
1527 nfe = 1
1528 ! load and invert the ewt array. (h is temporarily set to 1.0.) -------
1529 call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1530 do 110 i = 1,n
1531 if (rwork(i+lewt-1) .le. 0.0e0) go to 621
1532 110 rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
1533 if (miter .eq. 0 .or. miter .eq. 3) go to 120
1534 ! iprep and prep do sparse matrix preprocessing if miter = 1 or 2. -----
1535 lacor = min0(lacor,lrw)
1536 call iprep (neq, y, rwork, iwork(lia), iwork(lja), ipflag, f, jac, &
1537 ruserpar, nruserpar, iuserpar, niuserpar)
1538 lenrw = lwm - 1 + lenwk + lrest
1539 iwork(17) = lenrw
1540 if (ipflag .ne. -1) iwork(23) = ipian
1541 if (ipflag .ne. -1) iwork(24) = ipjan
1542 ipgo = -ipflag + 1
1543 go to (115, 628, 629, 630, 631, 632, 633), ipgo
1544 115 iwork(22) = lyh
1545 if (lenrw .gt. lrw) go to 617
1546 ! check tcrit for legality (itask = 4 or 5). ---------------------------
1547 120 continue
1548 if (itask .ne. 4 .and. itask .ne. 5) go to 125
1549 tcrit = rwork(1)
1550 if ((tcrit - tout)*(tout - t) .lt. 0.0e0) go to 625
1551 if (h0 .ne. 0.0e0 .and. (t + h0 - tcrit)*h0 .gt. 0.0e0) &
1552 h0 = tcrit - t
1553 ! initialize all remaining parameters. ---------------------------------
1554 125 uround = r1mach(4)
1555 jstart = 0
1556 if (miter .ne. 0) rwork(lwm) = sqrt(uround)
1557 msbj = 50
1558 nslj = 0
1559 ccmxj = 0.2e0
1560 psmall = 1000.0e0*uround
1561 rbig = 0.01e0/psmall
1562 nhnil = 0
1563 nje = 0
1564 nlu = 0
1565 nslast = 0
1566 hu = 0.0e0
1567 nqu = 0
1568 ccmax = 0.3e0
1569 maxcor = 3
1570 msbp = 20
1571 mxncf = 10
1572 !-----------------------------------------------------------------------
1573 ! the coding below computes the step size, h0, to be attempted on the
1574 ! first step, unless the user has supplied a value for this.
1575 ! first check that tout - t differs significantly from zero.
1576 ! a scalar tolerance quantity tol is computed, as max(rtol(i))
1577 ! if this is positive, or max(atol(i)/abs(y(i))) otherwise, adjusted
1578 ! so as to be between 100*uround and 1.0e-3.
1579 ! then the computed value h0 is given by..
1580 ! neq
1581 ! h0**2 = tol / ( w0**-2 + (1/neq) * sum ( f(i)/ywt(i) )**2 )
1582 ! 1
1583 ! where w0 = max ( abs(t), abs(tout) ),
1584 ! f(i) = i-th component of initial value of f,
1585 ! ywt(i) = ewt(i)/tol (a weight for y(i)).
1586 ! the sign of h0 is inferred from the initial values of tout and t.
1587 !-----------------------------------------------------------------------
1588 lf0 = lyh + nyh
1589 if (h0 .ne. 0.0e0) go to 180
1590 tdist = abs(tout - t)
1591 w0 = amax1(abs(t),abs(tout))
1592 if (tdist .lt. 2.0e0*uround*w0) go to 622
1593 tol = rtol(1)
1594 if (itol .le. 2) go to 140
1595 do 130 i = 1,n
1596 130 tol = amax1(tol,rtol(i))
1597 140 if (tol .gt. 0.0e0) go to 160
1598 atoli = atol(1)
1599 do 150 i = 1,n
1600 if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
1601 ayi = abs(y(i))
1602 if (ayi .ne. 0.0e0) tol = amax1(tol,atoli/ayi)
1603 150 continue
1604 160 tol = amax1(tol,100.0e0*uround)
1605 tol = amin1(tol,0.001e0)
1606 sum = vnorm (n, rwork(lf0), rwork(lewt))
1607 sum = 1.0e0/(tol*w0*w0) + tol*sum**2
1608 h0 = 1.0e0/sqrt(sum)
1609 h0 = amin1(h0,tdist)
1610 h0 = sign(h0,tout-t)
1611 ! adjust h0 if necessary to meet hmax bound. ---------------------------
1612 180 rh = abs(h0)*hmxi
1613 if (rh .gt. 1.0e0) h0 = h0/rh
1614 ! load h with h0 and scale yh(*,2) by h0. ------------------------------
1615 h = h0
1616 do 190 i = 1,n
1617 190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
1618 go to 270
1619 !-----------------------------------------------------------------------
1620 ! block d.
1621 ! the next code block is for continuation calls only (istate = 2 or 3)
1622 ! and is to check stop conditions before taking a step.
1623 !-----------------------------------------------------------------------
1624 200 nslast = nst
1625 go to (210, 250, 220, 230, 240), itask
1626 210 if ((tn - tout)*h .lt. 0.0e0) go to 250
1627 call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1628 if (iflag .ne. 0) go to 627
1629 t = tout
1630 go to 420
1631 220 tp = tn - hu*(1.0e0 + 100.0e0*uround)
1632 if ((tp - tout)*h .gt. 0.0e0) go to 623
1633 if ((tn - tout)*h .lt. 0.0e0) go to 250
1634 go to 400
1635 230 tcrit = rwork(1)
1636 if ((tn - tcrit)*h .gt. 0.0e0) go to 624
1637 if ((tcrit - tout)*h .lt. 0.0e0) go to 625
1638 if ((tn - tout)*h .lt. 0.0e0) go to 245
1639 call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1640 if (iflag .ne. 0) go to 627
1641 t = tout
1642 go to 420
1643 240 tcrit = rwork(1)
1644 if ((tn - tcrit)*h .gt. 0.0e0) go to 624
1645 245 hmx = abs(tn) + abs(h)
1646 ihit = abs(tn - tcrit) .le. 100.0e0*uround*hmx
1647 if (ihit) go to 400
1648 tnext = tn + h*(1.0e0 + 4.0e0*uround)
1649 if ((tnext - tcrit)*h .le. 0.0e0) go to 250
1650 h = (tcrit - tn)*(1.0e0 - 4.0e0*uround)
1651 if (istate .eq. 2) jstart = -2
1652 !-----------------------------------------------------------------------
1653 ! block e.
1654 ! the next block is normally executed for all calls and contains
1655 ! the call to the one-step core integrator stode.
1656 !
1657 ! this is a looping point for the integration steps.
1658 !
1659 ! first check for too many steps being taken, update ewt (if not at
1660 ! start of problem), check for too much accuracy being requested, and
1661 ! check for h below the roundoff level in t.
1662 !-----------------------------------------------------------------------
1663 250 continue
1664 if ((nst-nslast) .ge. mxstep) go to 500
1665 call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1666 do 260 i = 1,n
1667 if (rwork(i+lewt-1) .le. 0.0e0) go to 510
1668 260 rwork(i+lewt-1) = 1.0e0/rwork(i+lewt-1)
1669 270 tolsf = uround*vnorm (n, rwork(lyh), rwork(lewt))
1670 if (tolsf .le. 1.0e0) go to 280
1671 ! diagnostic dump
1672 tolsf = tolsf*2.0e0
1673 if (nst .eq. 0) go to 626
1674 go to 520
1675 280 if ((tn + h) .ne. tn) go to 290
1676 nhnil = nhnil + 1
1677 if (nhnil .gt. mxhnil) go to 290
1678 call xerrwv('lsodes-- warning..internal t (=r1) and h (=r2) are', &
1679 50, 101, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1680 call xerrwv( &
1681 ' such that in the machine, t + h = t on the next step ', &
1682 60, 101, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1683 call xerrwv(' (h = step size). solver will continue anyway', &
1684 50, 101, 0, 0, 0, 0, 2, tn, h)
1685 if (nhnil .lt. mxhnil) go to 290
1686 call xerrwv('lsodes-- above warning has been issued i1 times. ', &
1687 50, 102, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1688 call xerrwv(' it will not be issued again for this problem', &
1689 50, 102, 0, 1, mxhnil, 0, 0, 0.0e0, 0.0e0)
1690 290 continue
1691 !-----------------------------------------------------------------------
1692 ! call stode(neq,y,yh,nyh,yh,ewt,savf,acor,wm,wm,f,jac,prjs,slss)
1693 !-----------------------------------------------------------------------
1694 call stode_lsodes (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt), &
1695 rwork(lsavf), rwork(lacor), rwork(lwm), rwork(lwm), &
1696 f, jac, prjs, slss, &
1697 ruserpar, nruserpar, iuserpar, niuserpar )
1698 kgo = 1 - kflag
1699 go to (300, 530, 540, 550), kgo
1700 !-----------------------------------------------------------------------
1701 ! block f.
1702 ! the following block handles the case of a successful return from the
1703 ! core integrator (kflag = 0). test for stop conditions.
1704 !-----------------------------------------------------------------------
1705 300 init = 1
1706 go to (310, 400, 330, 340, 350), itask
1707 ! itask = 1. if tout has been reached, interpolate. -------------------
1708 310 if ((tn - tout)*h .lt. 0.0e0) go to 250
1709 call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1710 t = tout
1711 go to 420
1712 ! itask = 3. jump to exit if tout was reached. ------------------------
1713 330 if ((tn - tout)*h .ge. 0.0e0) go to 400
1714 go to 250
1715 ! itask = 4. see if tout or tcrit was reached. adjust h if necessary.
1716 340 if ((tn - tout)*h .lt. 0.0e0) go to 345
1717 call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1718 t = tout
1719 go to 420
1720 345 hmx = abs(tn) + abs(h)
1721 ihit = abs(tn - tcrit) .le. 100.0e0*uround*hmx
1722 if (ihit) go to 400
1723 tnext = tn + h*(1.0e0 + 4.0e0*uround)
1724 if ((tnext - tcrit)*h .le. 0.0e0) go to 250
1725 h = (tcrit - tn)*(1.0e0 - 4.0e0*uround)
1726 jstart = -2
1727 go to 250
1728 ! itask = 5. see if tcrit was reached and jump to exit. ---------------
1729 350 hmx = abs(tn) + abs(h)
1730 ihit = abs(tn - tcrit) .le. 100.0e0*uround*hmx
1731 !-----------------------------------------------------------------------
1732 ! block g.
1733 ! the following block handles all successful returns from lsodes.
1734 ! if itask .ne. 1, y is loaded from yh and t is set accordingly.
1735 ! istate is set to 2, the illegal input counter is zeroed, and the
1736 ! optional outputs are loaded into the work arrays before returning.
1737 ! if istate = 1 and tout = t, there is a return with no action taken,
1738 ! except that if this has happened repeatedly, the run is terminated.
1739 !-----------------------------------------------------------------------
1740 400 do 410 i = 1,n
1741 410 y(i) = rwork(i+lyh-1)
1742 t = tn
1743 if (itask .ne. 4 .and. itask .ne. 5) go to 420
1744 if (ihit) t = tcrit
1745 420 istate = 2
1746 illin = 0
1747 rwork(11) = hu
1748 rwork(12) = h
1749 rwork(13) = tn
1750 iwork(11) = nst
1751 iwork(12) = nfe
1752 iwork(13) = nje
1753 iwork(14) = nqu
1754 iwork(15) = nq
1755 iwork(19) = nnz
1756 iwork(20) = ngp
1757 iwork(21) = nlu
1758 iwork(25) = nzl
1759 iwork(26) = nzu
1760 return
1761 !
1762 430 ntrep = ntrep + 1
1763 if (ntrep .lt. 5) return
1764 call xerrwv( &
1765 'lsodes-- repeated calls with istate = 1 and tout = t (=r1) ', &
1766 60, 301, 0, 0, 0, 0, 1, t, 0.0e0)
1767 go to 800
1768 !-----------------------------------------------------------------------
1769 ! block h.
1770 ! the following block handles all unsuccessful returns other than
1771 ! those for illegal input. first the error message routine is called.
1772 ! if there was an error test or convergence test failure, imxer is set.
1773 ! then y is loaded from yh, t is set to tn, and the illegal input
1774 ! counter illin is set to 0. the optional outputs are loaded into
1775 ! the work arrays before returning.
1776 !-----------------------------------------------------------------------
1777 ! the maximum number of steps was taken before reaching tout. ----------
1778 500 call xerrwv('lsodes-- at current t (=r1), mxstep (=i1) steps ', &
1779 50, 201, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1780 call xerrwv(' taken on this call before reaching tout ', &
1781 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0e0)
1782 istate = -1
1783 go to 580
1784 ! ewt(i) .le. 0.0 for some i (not at start of problem). ----------------
1785 510 ewti = rwork(lewt+i-1)
1786 call xerrwv('lsodes-- at t (=r1), ewt(i1) has become r2 .le. 0.', &
1787 50, 202, 0, 1, i, 0, 2, tn, ewti)
1788 istate = -6
1789 go to 580
1790 ! too much accuracy requested for machine precision. -------------------
1791 520 call xerrwv('lsodes-- at t (=r1), too much accuracy requested ', &
1792 50, 203, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1793 call xerrwv(' for precision of machine.. see tolsf (=r2) ', &
1794 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
1795 rwork(14) = tolsf
1796 istate = -2
1797 go to 580
1798 ! kflag = -1. error test failed repeatedly or with abs(h) = hmin. -----
1799 530 call xerrwv('lsodes-- at t(=r1) and step size h(=r2), the error', &
1800 50, 204, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1801 call xerrwv(' test failed repeatedly or with abs(h) = hmin', &
1802 50, 204, 0, 0, 0, 0, 2, tn, h)
1803 istate = -4
1804 go to 560
1805 ! kflag = -2. convergence failed repeatedly or with abs(h) = hmin. ----
1806 540 call xerrwv('lsodes-- at t (=r1) and step size h (=r2), the ', &
1807 50, 205, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1808 call xerrwv(' corrector convergence failed repeatedly ', &
1809 50, 205, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1810 call xerrwv(' or with abs(h) = hmin ', &
1811 30, 205, 0, 0, 0, 0, 2, tn, h)
1812 istate = -5
1813 go to 560
1814 ! kflag = -3. fatal error flag returned by prjs or slss (cdrv). -------
1815 550 call xerrwv('lsodes-- at t (=r1) and step size h (=r2), a fatal', &
1816 50, 207, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1817 call xerrwv(' error flag was returned by cdrv (by way of ', &
1818 50, 207, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1819 call xerrwv(' subroutine prjs or slss)', &
1820 30, 207, 0, 0, 0, 0, 2, tn, h)
1821 istate = -7
1822 go to 580
1823 ! compute imxer if relevant. -------------------------------------------
1824 560 big = 0.0e0
1825 imxer = 1
1826 do 570 i = 1,n
1827 size = abs(rwork(i+lacor-1)*rwork(i+lewt-1))
1828 if (big .ge. size) go to 570
1829 big = size
1830 imxer = i
1831 570 continue
1832 iwork(16) = imxer
1833 ! set y vector, t, illin, and optional outputs. ------------------------
1834 580 do 590 i = 1,n
1835 590 y(i) = rwork(i+lyh-1)
1836 t = tn
1837 illin = 0
1838 rwork(11) = hu
1839 rwork(12) = h
1840 rwork(13) = tn
1841 iwork(11) = nst
1842 iwork(12) = nfe
1843 iwork(13) = nje
1844 iwork(14) = nqu
1845 iwork(15) = nq
1846 iwork(19) = nnz
1847 iwork(20) = ngp
1848 iwork(21) = nlu
1849 iwork(25) = nzl
1850 iwork(26) = nzu
1851 return
1852 !-----------------------------------------------------------------------
1853 ! block i.
1854 ! the following block handles all error returns due to illegal input
1855 ! (istate = -3), as detected before calling the core integrator.
1856 ! first the error message routine is called. then if there have been
1857 ! 5 consecutive such returns just before this call to the solver,
1858 ! the run is halted.
1859 !-----------------------------------------------------------------------
1860 601 call xerrwv('lsodes-- istate (=i1) illegal ', &
1861 30, 1, 0, 1, istate, 0, 0, 0.0e0, 0.0e0)
1862 go to 700
1863 602 call xerrwv('lsodes-- itask (=i1) illegal ', &
1864 30, 2, 0, 1, itask, 0, 0, 0.0e0, 0.0e0)
1865 go to 700
1866 603 call xerrwv('lsodes-- istate .gt. 1 but lsodes not initialized ', &
1867 50, 3, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1868 go to 700
1869 604 call xerrwv('lsodes-- neq (=i1) .lt. 1 ', &
1870 30, 4, 0, 1, neq(1), 0, 0, 0.0e0, 0.0e0)
1871 go to 700
1872 605 call xerrwv('lsodes-- istate = 3 and neq increased (i1 to i2) ', &
1873 50, 5, 0, 2, n, neq(1), 0, 0.0e0, 0.0e0)
1874 go to 700
1875 606 call xerrwv('lsodes-- itol (=i1) illegal ', &
1876 30, 6, 0, 1, itol, 0, 0, 0.0e0, 0.0e0)
1877 go to 700
1878 607 call xerrwv('lsodes-- iopt (=i1) illegal ', &
1879 30, 7, 0, 1, iopt, 0, 0, 0.0e0, 0.0e0)
1880 go to 700
1881 608 call xerrwv('lsodes-- mf (=i1) illegal ', &
1882 30, 8, 0, 1, mf, 0, 0, 0.0e0, 0.0e0)
1883 go to 700
1884 609 call xerrwv('lsodes-- seth (=r1) .lt. 0.0 ', &
1885 30, 9, 0, 0, 0, 0, 1, seth, 0.0e0)
1886 go to 700
1887 611 call xerrwv('lsodes-- maxord (=i1) .lt. 0 ', &
1888 30, 11, 0, 1, maxord, 0, 0, 0.0e0, 0.0e0)
1889 go to 700
1890 612 call xerrwv('lsodes-- mxstep (=i1) .lt. 0 ', &
1891 30, 12, 0, 1, mxstep, 0, 0, 0.0e0, 0.0e0)
1892 go to 700
1893 613 call xerrwv('lsodes-- mxhnil (=i1) .lt. 0 ', &
1894 30, 13, 0, 1, mxhnil, 0, 0, 0.0e0, 0.0e0)
1895 go to 700
1896 614 call xerrwv('lsodes-- tout (=r1) behind t (=r2) ', &
1897 40, 14, 0, 0, 0, 0, 2, tout, t)
1898 call xerrwv(' integration direction is given by h0 (=r1) ', &
1899 50, 14, 0, 0, 0, 0, 1, h0, 0.0e0)
1900 go to 700
1901 615 call xerrwv('lsodes-- hmax (=r1) .lt. 0.0 ', &
1902 30, 15, 0, 0, 0, 0, 1, hmax, 0.0e0)
1903 go to 700
1904 616 call xerrwv('lsodes-- hmin (=r1) .lt. 0.0 ', &
1905 30, 16, 0, 0, 0, 0, 1, hmin, 0.0e0)
1906 go to 700
1907 617 call xerrwv('lsodes-- rwork length is insufficient to proceed. ', &
1908 50, 17, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1909 call xerrwv( &
1910 ' length needed is .ge. lenrw (=i1), exceeds lrw (=i2)', &
1911 60, 17, 0, 2, lenrw, lrw, 0, 0.0e0, 0.0e0)
1912 go to 700
1913 618 call xerrwv('lsodes-- iwork length is insufficient to proceed. ', &
1914 50, 18, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1915 call xerrwv( &
1916 ' length needed is .ge. leniw (=i1), exceeds liw (=i2)', &
1917 60, 18, 0, 2, leniw, liw, 0, 0.0e0, 0.0e0)
1918 go to 700
1919 619 call xerrwv('lsodes-- rtol(i1) is r1 .lt. 0.0 ', &
1920 40, 19, 0, 1, i, 0, 1, rtoli, 0.0e0)
1921 go to 700
1922 620 call xerrwv('lsodes-- atol(i1) is r1 .lt. 0.0 ', &
1923 40, 20, 0, 1, i, 0, 1, atoli, 0.0e0)
1924 go to 700
1925 621 ewti = rwork(lewt+i-1)
1926 call xerrwv('lsodes-- ewt(i1) is r1 .le. 0.0 ', &
1927 40, 21, 0, 1, i, 0, 1, ewti, 0.0e0)
1928 go to 700
1929 622 call xerrwv( &
1930 'lsodes-- tout (=r1) too close to t(=r2) to start integration', &
1931 60, 22, 0, 0, 0, 0, 2, tout, t)
1932 go to 700
1933 623 call xerrwv( &
1934 'lsodes-- itask = i1 and tout (=r1) behind tcur - hu (= r2) ', &
1935 60, 23, 0, 1, itask, 0, 2, tout, tp)
1936 go to 700
1937 624 call xerrwv( &
1938 'lsodes-- itask = 4 or 5 and tcrit (=r1) behind tcur (=r2) ', &
1939 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
1940 go to 700
1941 625 call xerrwv( &
1942 'lsodes-- itask = 4 or 5 and tcrit (=r1) behind tout (=r2) ', &
1943 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
1944 go to 700
1945 626 call xerrwv('lsodes-- at start of problem, too much accuracy ', &
1946 50, 26, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1947 call xerrwv( &
1948 ' requested for precision of machine.. see tolsf (=r1) ', &
1949 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0e0)
1950 rwork(14) = tolsf
1951 go to 700
1952 627 call xerrwv('lsodes-- trouble from intdy. itask = i1, tout = r1', &
1953 50, 27, 0, 1, itask, 0, 1, tout, 0.0e0)
1954 go to 700
1955 628 call xerrwv( &
1956 'lsodes-- rwork length insufficient (for subroutine prep). ', &
1957 60, 28, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1958 call xerrwv( &
1959 ' length needed is .ge. lenrw (=i1), exceeds lrw (=i2)', &
1960 60, 28, 0, 2, lenrw, lrw, 0, 0.0e0, 0.0e0)
1961 go to 700
1962 629 call xerrwv( &
1963 'lsodes-- rwork length insufficient (for subroutine jgroup). ', &
1964 60, 29, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1965 call xerrwv( &
1966 ' length needed is .ge. lenrw (=i1), exceeds lrw (=i2)', &
1967 60, 29, 0, 2, lenrw, lrw, 0, 0.0e0, 0.0e0)
1968 go to 700
1969 630 call xerrwv( &
1970 'lsodes-- rwork length insufficient (for subroutine odrv). ', &
1971 60, 30, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1972 call xerrwv( &
1973 ' length needed is .ge. lenrw (=i1), exceeds lrw (=i2)', &
1974 60, 30, 0, 2, lenrw, lrw, 0, 0.0e0, 0.0e0)
1975 go to 700
1976 631 call xerrwv( &
1977 'lsodes-- error from odrv in yale sparse matrix package ', &
1978 60, 31, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1979 imul = (iys - 1)/n
1980 irem = iys - imul*n
1981 call xerrwv( &
1982 ' at t (=r1), odrv returned error flag = i1*neq + i2. ', &
1983 60, 31, 0, 2, imul, irem, 1, tn, 0.0e0)
1984 go to 700
1985 632 call xerrwv( &
1986 'lsodes-- rwork length insufficient (for subroutine cdrv). ', &
1987 60, 32, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1988 call xerrwv( &
1989 ' length needed is .ge. lenrw (=i1), exceeds lrw (=i2)', &
1990 60, 32, 0, 2, lenrw, lrw, 0, 0.0e0, 0.0e0)
1991 go to 700
1992 633 call xerrwv( &
1993 'lsodes-- error from cdrv in yale sparse matrix package ', &
1994 60, 33, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
1995 imul = (iys - 1)/n
1996 irem = iys - imul*n
1997 call xerrwv( &
1998 ' at t (=r1), cdrv returned error flag = i1*neq + i2. ', &
1999 60, 33, 0, 2, imul, irem, 1, tn, 0.0e0)
2000 if (imul .eq. 2) call xerrwv( &
2001 ' duplicate entry in sparsity structure descriptors ', &
2002 60, 33, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
2003 if (imul .eq. 3 .or. imul .eq. 6) call xerrwv( &
2004 ' insufficient storage for nsfc (called by cdrv) ', &
2005 60, 33, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
2006 !
2007 700 if (illin .eq. 5) go to 710
2008 illin = illin + 1
2009 istate = -3
2010 return
2011 710 call xerrwv('lsodes-- repeated occurrences of illegal input ', &
2012 50, 302, 0, 0, 0, 0, 0, 0.0e0, 0.0e0)
2013 !
2014 800 call xerrwv('lsodes-- run aborted.. apparent infinite loop ', &
2015 50, 303, 2, 0, 0, 0, 0, 0.0e0, 0.0e0)
2016 return
2017 !----------------------- end of subroutine lsodes ----------------------
2018 end subroutine lsodes_solver
2019 subroutine adjlr (n, isp, ldif)
2020 integer n, isp, ldif
2021 !jdf dimension isp(1)
2022 dimension isp(*)
2023 !-----------------------------------------------------------------------
2024 ! this routine computes an adjustment, ldif, to the required
2025 ! integer storage space in iwk (sparse matrix work space).
2026 ! it is called only if the word length ratio is lrat = 1.
2027 ! this is to account for the possibility that the symbolic lu phase
2028 ! may require more storage than the numerical lu and solution phases.
2029 !-----------------------------------------------------------------------
2030 integer ip, jlmax, jumax, lnfc, lsfc, nzlu
2031 !
2032 ip = 2*n + 1
2033 ! get jlmax = ijl(n) and jumax = iju(n) (sizes of jl and ju). ----------
2034 jlmax = isp(ip)
2035 jumax = isp(ip+ip)
2036 ! nzlu = (size of l) + (size of u) = (il(n+1)-il(1)) + (iu(n+1)-iu(1)).
2037 nzlu = isp(n+1) - isp(1) + isp(ip+n+1) - isp(ip+1)
2038 lsfc = 12*n + 3 + 2*max0(jlmax,jumax)
2039 lnfc = 9*n + 2 + jlmax + jumax + nzlu
2040 ldif = max0(0, lsfc - lnfc)
2041 return
2042 !----------------------- end of subroutine adjlr -----------------------
2043 end subroutine adjlr
2044 subroutine cdrv &
2045 (n, r,c,ic, ia,ja,a, b, z, nsp,isp,rsp,esp, path, flag)
2046 !lll. optimize
2047 !*** subroutine cdrv
2048 !*** driver for subroutines for solving sparse nonsymmetric systems of
2049 ! linear equations (compressed pointer storage)
2050 !
2051 !
2052 ! parameters
2053 ! class abbreviations are--
2054 ! n - integer variable
2055 ! f - real variable
2056 ! v - supplies a value to the driver
2057 ! r - returns a result from the driver
2058 ! i - used internally by the driver
2059 ! a - array
2060 !
2061 ! class - parameter
2062 ! ------+----------
2063 ! -
2064 ! the nonzero entries of the coefficient matrix m are stored
2065 ! row-by-row in the array a. to identify the individual nonzero
2066 ! entries in each row, we need to know in which column each entry
2067 ! lies. the column indices which correspond to the nonzero entries
2068 ! of m are stored in the array ja. i.e., if a(k) = m(i,j), then
2069 ! ja(k) = j. in addition, we need to know where each row starts and
2070 ! how long it is. the index positions in ja and a where the rows of
2071 ! m begin are stored in the array ia. i.e., if m(i,j) is the first
2072 ! nonzero entry (stored) in the i-th row and a(k) = m(i,j), then
2073 ! ia(i) = k. moreover, the index in ja and a of the first location
2074 ! following the last element in the last row is stored in ia(n+1).
2075 ! thus, the number of entries in the i-th row is given by
2076 ! ia(i+1) - ia(i), the nonzero entries of the i-th row are stored
2077 ! consecutively in
2078 ! a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1),
2079 ! and the corresponding column indices are stored consecutively in
2080 ! ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1).
2081 ! for example, the 5 by 5 matrix
2082 ! ( 1. 0. 2. 0. 0.)
2083 ! ( 0. 3. 0. 0. 0.)
2084 ! m = ( 0. 4. 5. 6. 0.)
2085 ! ( 0. 0. 0. 7. 0.)
2086 ! ( 0. 0. 0. 8. 9.)
2087 ! would be stored as
2088 ! - 1 2 3 4 5 6 7 8 9
2089 ! ---+--------------------------
2090 ! ia - 1 3 4 7 8 10
2091 ! ja - 1 3 2 2 3 4 4 4 5
2092 ! a - 1. 2. 3. 4. 5. 6. 7. 8. 9. .
2093 !
2094 ! nv - n - number of variables/equations.
2095 ! fva - a - nonzero entries of the coefficient matrix m, stored
2096 ! - by rows.
2097 ! - size = number of nonzero entries in m.
2098 ! nva - ia - pointers to delimit the rows in a.
2099 ! - size = n+1.
2100 ! nva - ja - column numbers corresponding to the elements of a.
2101 ! - size = size of a.
2102 ! fva - b - right-hand side b. b and z can the same array.
2103 ! - size = n.
2104 ! fra - z - solution x. b and z can be the same array.
2105 ! - size = n.
2106 !
2107 ! the rows and columns of the original matrix m can be
2108 ! reordered (e.g., to reduce fillin or ensure numerical stability)
2109 ! before calling the driver. if no reordering is done, then set
2110 ! r(i) = c(i) = ic(i) = i for i=1,...,n. the solution z is returned
2111 ! in the original order.
2112 ! if the columns have been reordered (i.e., c(i).ne.i for some
2113 ! i), then the driver will call a subroutine (nroc) which rearranges
2114 ! each row of ja and a, leaving the rows in the original order, but
2115 ! placing the elements of each row in increasing order with respect
2116 ! to the new ordering. if path.ne.1, then nroc is assumed to have
2117 ! been called already.
2118 !
2119 ! nva - r - ordering of the rows of m.
2120 ! - size = n.
2121 ! nva - c - ordering of the columns of m.
2122 ! - size = n.
2123 ! nva - ic - inverse of the ordering of the columns of m. i.e.,
2124 ! - ic(c(i)) = i for i=1,...,n.
2125 ! - size = n.
2126 !
2127 ! the solution of the system of linear equations is divided into
2128 ! three stages --
2129 ! nsfc -- the matrix m is processed symbolically to determine where
2130 ! fillin will occur during the numeric factorization.
2131 ! nnfc -- the matrix m is factored numerically into the product ldu
2132 ! of a unit lower triangular matrix l, a diagonal matrix
2133 ! d, and a unit upper triangular matrix u, and the system
2134 ! mx = b is solved.
2135 ! nnsc -- the linear system mx = b is solved using the ldu
2136 ! or factorization from nnfc.
2137 ! nntc -- the transposed linear system mt x = b is solved using
2138 ! the ldu factorization from nnf.
2139 ! for several systems whose coefficient matrices have the same
2140 ! nonzero structure, nsfc need be done only once (for the first
2141 ! system). then nnfc is done once for each additional system. for
2142 ! several systems with the same coefficient matrix, nsfc and nnfc
2143 ! need be done only once (for the first system). then nnsc or nntc
2144 ! is done once for each additional right-hand side.
2145 !
2146 ! nv - path - path specification. values and their meanings are --
2147 ! - 1 perform nroc, nsfc, and nnfc.
2148 ! - 2 perform nnfc only (nsfc is assumed to have been
2149 ! - done in a manner compatible with the storage
2150 ! - allocation used in the driver).
2151 ! - 3 perform nnsc only (nsfc and nnfc are assumed to
2152 ! - have been done in a manner compatible with the
2153 ! - storage allocation used in the driver).
2154 ! - 4 perform nntc only (nsfc and nnfc are assumed to
2155 ! - have been done in a manner compatible with the
2156 ! - storage allocation used in the driver).
2157 ! - 5 perform nroc and nsfc.
2158 !
2159 ! various errors are detected by the driver and the individual
2160 ! subroutines.
2161 !
2162 ! nr - flag - error flag. values and their meanings are --
2163 ! - 0 no errors detected
2164 ! - n+k null row in a -- row = k
2165 ! - 2n+k duplicate entry in a -- row = k
2166 ! - 3n+k insufficient storage in nsfc -- row = k
2167 ! - 4n+1 insufficient storage in nnfc
2168 ! - 5n+k null pivot -- row = k
2169 ! - 6n+k insufficient storage in nsfc -- row = k
2170 ! - 7n+1 insufficient storage in nnfc
2171 ! - 8n+k zero pivot -- row = k
2172 ! - 10n+1 insufficient storage in cdrv
2173 ! - 11n+1 illegal path specification
2174 !
2175 ! working storage is needed for the factored form of the matrix
2176 ! m plus various temporary vectors. the arrays isp and rsp should be
2177 ! equivalenced. integer storage is allocated from the beginning of
2178 ! isp and real storage from the end of rsp.
2179 !
2180 ! nv - nsp - declared dimension of rsp. nsp generally must
2181 ! - be larger than 8n+2 + 2k (where k = (number of
2182 ! - nonzero entries in m)).
2183 ! nvira - isp - integer working storage divided up into various arrays
2184 ! - needed by the subroutines. isp and rsp should be
2185 ! - equivalenced.
2186 ! - size = lratio*nsp.
2187 ! fvira - rsp - real working storage divided up into various arrays
2188 ! - needed by the subroutines. isp and rsp should be
2189 ! - equivalenced.
2190 ! - size = nsp.
2191 ! nr - esp - if sufficient storage was available to perform the
2192 ! - symbolic factorization (nsfc), then esp is set to
2193 ! - the amount of excess storage provided (negative if
2194 ! - insufficient storage was available to perform the
2195 ! - numeric factorization (nnfc)).
2196 !
2197 !
2198 ! conversion to double precision
2199 !
2200 ! to convert these routines for double precision arrays..
2201 ! (1) use the double precision declarations in place of the real
2202 ! declarations in each subprogram, as given in comment cards.
2203 ! (2) change the data-loaded value of the integer lratio
2204 ! in subroutine cdrv, as indicated below.
2205 ! (3) change e0 to d0 in the constants in statement number 10
2206 ! in subroutine nnfc and the line following that.
2207 !
2208 !jdf integer r(1), c(1), ic(1), ia(1), ja(1), isp(1), esp, path,
2209 !jdf * flag, d, u, q, row, tmp, ar, umax
2210 !jdf real a(1), b(1), z(1), rsp(1)
2211 integer r(*), c(*), ic(*), ia(*), ja(*), isp(*), esp, path, &
2212 flag, d, u, q, row, tmp, ar, umax
2213 real a(*), b(*), z(*), rsp(*)
2214 ! double precision a(1), b(1), z(1), rsp(1)
2215 !
2216 ! set lratio equal to the ratio between the length of floating point
2217 ! and integer array data. e. g., lratio = 1 for (real, integer),
2218 ! lratio = 2 for (double precision, integer)
2219 !
2220 data lratio/1/
2221 !
2222 if (path.lt.1 .or. 5.lt.path) go to 111
2223 !******initialize and divide up temporary storage *******************
2224 il = 1
2225 ijl = il + (n+1)
2226 iu = ijl + n
2227 iju = iu + (n+1)
2228 irl = iju + n
2229 jrl = irl + n
2230 jl = jrl + n
2231 !
2232 ! ****** reorder a if necessary, call nsfc if flag is set ***********
2233 if ((path-1) * (path-5) .ne. 0) go to 5
2234 max = (lratio*nsp + 1 - jl) - (n+1) - 5*n
2235 jlmax = max/2
2236 q = jl + jlmax
2237 ira = q + (n+1)
2238 jra = ira + n
2239 irac = jra + n
2240 iru = irac + n
2241 jru = iru + n
2242 jutmp = jru + n
2243 jumax = lratio*nsp + 1 - jutmp
2244 esp = max/lratio
2245 if (jlmax.le.0 .or. jumax.le.0) go to 110
2246 !
2247 do 1 i=1,n
2248 if (c(i).ne.i) go to 2
2249 1 continue
2250 go to 3
2251 2 ar = nsp + 1 - n
2252 call nroc &
2253 (n, ic, ia,ja,a, isp(il), rsp(ar), isp(iu), flag)
2254 if (flag.ne.0) go to 100
2255 !
2256 3 call nsfc &
2257 (n, r, ic, ia,ja, &
2258 jlmax, isp(il), isp(jl), isp(ijl), &
2259 jumax, isp(iu), isp(jutmp), isp(iju), &
2260 isp(q), isp(ira), isp(jra), isp(irac), &
2261 isp(irl), isp(jrl), isp(iru), isp(jru), flag)
2262 if(flag .ne. 0) go to 100
2263 ! ****** move ju next to jl *****************************************
2264 jlmax = isp(ijl+n-1)
2265 ju = jl + jlmax
2266 jumax = isp(iju+n-1)
2267 if (jumax.le.0) go to 5
2268 do 4 j=1,jumax
2269 4 isp(ju+j-1) = isp(jutmp+j-1)
2270 !
2271 ! ****** call remaining subroutines *********************************
2272 5 jlmax = isp(ijl+n-1)
2273 ju = jl + jlmax
2274 jumax = isp(iju+n-1)
2275 l = (ju + jumax - 2 + lratio) / lratio + 1
2276 lmax = isp(il+n) - 1
2277 d = l + lmax
2278 u = d + n
2279 row = nsp + 1 - n
2280 tmp = row - n
2281 umax = tmp - u
2282 esp = umax - (isp(iu+n) - 1)
2283 !
2284 if ((path-1) * (path-2) .ne. 0) go to 6
2285 if (umax.lt.0) go to 110
2286 call nnfc &
2287 (n, r, c, ic, ia, ja, a, z, b, &
2288 lmax, isp(il), isp(jl), isp(ijl), rsp(l), rsp(d), &
2289 umax, isp(iu), isp(ju), isp(iju), rsp(u), &
2290 rsp(row), rsp(tmp), isp(irl), isp(jrl), flag)
2291 if(flag .ne. 0) go to 100
2292 !
2293 6 if ((path-3) .ne. 0) go to 7
2294 call nnsc &
2295 (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l), &
2296 rsp(d), isp(iu), isp(ju), isp(iju), rsp(u), &
2297 z, b, rsp(tmp))
2298 !
2299 7 if ((path-4) .ne. 0) go to 8
2300 call nntc &
2301 (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l), &
2302 rsp(d), isp(iu), isp(ju), isp(iju), rsp(u), &
2303 z, b, rsp(tmp))
2304 8 return
2305 !
2306 ! ** error.. error detected in nroc, nsfc, nnfc, or nnsc
2307 100 return
2308 ! ** error.. insufficient storage
2309 110 flag = 10*n + 1
2310 return
2311 ! ** error.. illegal path specification
2312 111 flag = 11*n + 1
2313 return
2314 end subroutine cdrv
2315 subroutine cfode (meth, elco, tesco)
2316 !lll. optimize
2317 integer meth
2318 integer i, ib, nq, nqm1, nqp1
2319 real elco, tesco
2320 real agamq, fnq, fnqm1, pc, pint, ragq, &
2321 rqfac, rq1fac, tsign, xpin
2322 dimension elco(13,12), tesco(3,12)
2323 !-----------------------------------------------------------------------
2324 ! cfode is called by the integrator routine to set coefficients
2325 ! needed there. the coefficients for the current method, as
2326 ! given by the value of meth, are set for all orders and saved.
2327 ! the maximum order assumed here is 12 if meth = 1 and 5 if meth = 2.
2328 ! (a smaller value of the maximum order is also allowed.)
2329 ! cfode is called once at the beginning of the problem,
2330 ! and is not called again unless and until meth is changed.
2331 !
2332 ! the elco array contains the basic method coefficients.
2333 ! the coefficients el(i), 1 .le. i .le. nq+1, for the method of
2334 ! order nq are stored in elco(i,nq). they are given by a genetrating
2335 ! polynomial, i.e.,
2336 ! l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq.
2337 ! for the implicit adams methods, l(x) is given by
2338 ! dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0.
2339 ! for the bdf methods, l(x) is given by
2340 ! l(x) = (x+1)*(x+2)* ... *(x+nq)/k,
2341 ! where k = factorial(nq)*(1 + 1/2 + ... + 1/nq).
2342 !
2343 ! the tesco array contains test constants used for the
2344 ! local error test and the selection of step size and/or order.
2345 ! at order nq, tesco(k,nq) is used for the selection of step
2346 ! size at order nq - 1 if k = 1, at order nq if k = 2, and at order
2347 ! nq + 1 if k = 3.
2348 !-----------------------------------------------------------------------
2349 dimension pc(12)
2350 !
2351 go to (100, 200), meth
2352 !
2353 100 elco(1,1) = 1.0e0
2354 elco(2,1) = 1.0e0
2355 tesco(1,1) = 0.0e0
2356 tesco(2,1) = 2.0e0
2357 tesco(1,2) = 1.0e0
2358 tesco(3,12) = 0.0e0
2359 pc(1) = 1.0e0
2360 rqfac = 1.0e0
2361 do 140 nq = 2,12
2362 !-----------------------------------------------------------------------
2363 ! the pc array will contain the coefficients of the polynomial
2364 ! p(x) = (x+1)*(x+2)*...*(x+nq-1).
2365 ! initially, p(x) = 1.
2366 !-----------------------------------------------------------------------
2367 rq1fac = rqfac
2368 rqfac = rqfac/float(nq)
2369 nqm1 = nq - 1
2370 fnqm1 = float(nqm1)
2371 nqp1 = nq + 1
2372 ! form coefficients of p(x)*(x+nq-1). ----------------------------------
2373 pc(nq) = 0.0e0
2374 do 110 ib = 1,nqm1
2375 i = nqp1 - ib
2376 110 pc(i) = pc(i-1) + fnqm1*pc(i)
2377 pc(1) = fnqm1*pc(1)
2378 ! compute integral, -1 to 0, of p(x) and x*p(x). -----------------------
2379 pint = pc(1)
2380 xpin = pc(1)/2.0e0
2381 tsign = 1.0e0
2382 do 120 i = 2,nq
2383 tsign = -tsign
2384 pint = pint + tsign*pc(i)/float(i)
2385 120 xpin = xpin + tsign*pc(i)/float(i+1)
2386 ! store coefficients in elco and tesco. --------------------------------
2387 elco(1,nq) = pint*rq1fac
2388 elco(2,nq) = 1.0e0
2389 do 130 i = 2,nq
2390 130 elco(i+1,nq) = rq1fac*pc(i)/float(i)
2391 agamq = rqfac*xpin
2392 ragq = 1.0e0/agamq
2393 tesco(2,nq) = ragq
2394 if (nq .lt. 12) tesco(1,nqp1) = ragq*rqfac/float(nqp1)
2395 tesco(3,nqm1) = ragq
2396 140 continue
2397 return
2398 !
2399 200 pc(1) = 1.0e0
2400 rq1fac = 1.0e0
2401 do 230 nq = 1,5
2402 !-----------------------------------------------------------------------
2403 ! the pc array will contain the coefficients of the polynomial
2404 ! p(x) = (x+1)*(x+2)*...*(x+nq).
2405 ! initially, p(x) = 1.
2406 !-----------------------------------------------------------------------
2407 fnq = float(nq)
2408 nqp1 = nq + 1
2409 ! form coefficients of p(x)*(x+nq). ------------------------------------
2410 pc(nqp1) = 0.0e0
2411 do 210 ib = 1,nq
2412 i = nq + 2 - ib
2413 210 pc(i) = pc(i-1) + fnq*pc(i)
2414 pc(1) = fnq*pc(1)
2415 ! store coefficients in elco and tesco. --------------------------------
2416 do 220 i = 1,nqp1
2417 220 elco(i,nq) = pc(i)/pc(2)
2418 elco(2,nq) = 1.0e0
2419 tesco(1,nq) = rq1fac
2420 tesco(2,nq) = float(nqp1)/elco(1,nq)
2421 tesco(3,nq) = float(nq+2)/elco(1,nq)
2422 rq1fac = rq1fac/fnq
2423 230 continue
2424 return
2425 !----------------------- end of subroutine cfode -----------------------
2426 end subroutine cfode
2427 subroutine cntnzu (n, ia, ja, nzsut)
2428 integer n, ia, ja, nzsut
2429 !jdf dimension ia(1), ja(1)
2430 dimension ia(*), ja(*)
2431 !-----------------------------------------------------------------------
2432 ! this routine counts the number of nonzero elements in the strict
2433 ! upper triangle of the matrix m + m(transpose), where the sparsity
2434 ! structure of m is given by pointer arrays ia and ja.
2435 ! this is needed to compute the storage requirements for the
2436 ! sparse matrix reordering operation in odrv.
2437 !-----------------------------------------------------------------------
2438 integer ii, jj, j, jmin, jmax, k, kmin, kmax, num
2439 !
2440 num = 0
2441 do 50 ii = 1,n
2442 jmin = ia(ii)
2443 jmax = ia(ii+1) - 1
2444 if (jmin .gt. jmax) go to 50
2445 do 40 j = jmin,jmax
2446 if (ja(j) - ii) 10, 40, 30
2447 10 jj =ja(j)
2448 kmin = ia(jj)
2449 kmax = ia(jj+1) - 1
2450 if (kmin .gt. kmax) go to 30
2451 do 20 k = kmin,kmax
2452 if (ja(k) .eq. ii) go to 40
2453 20 continue
2454 30 num = num + 1
2455 40 continue
2456 50 continue
2457 nzsut = num
2458 return
2459 !----------------------- end of subroutine cntnzu ----------------------
2460 end subroutine cntnzu
2461 subroutine ewset (n, itol, rtol, atol, ycur, ewt)
2462 !lll. optimize
2463 !-----------------------------------------------------------------------
2464 ! this subroutine sets the error weight vector ewt according to
2465 ! ewt(i) = rtol(i)*abs(ycur(i)) + atol(i), i = 1,...,n,
2466 ! with the subscript on rtol and/or atol possibly replaced by 1 above,
2467 ! depending on the value of itol.
2468 !-----------------------------------------------------------------------
2469 integer n, itol
2470 integer i
2471 real rtol, atol, ycur, ewt
2472 !jdf dimension rtol(1), atol(1), ycur(n), ewt(n)
2473 dimension rtol(*), atol(*), ycur(n), ewt(n)
2474 !
2475 go to (10, 20, 30, 40), itol
2476 10 continue
2477 do 15 i = 1,n
2478 15 ewt(i) = rtol(1)*abs(ycur(i)) + atol(1)
2479 return
2480 20 continue
2481 do 25 i = 1,n
2482 25 ewt(i) = rtol(1)*abs(ycur(i)) + atol(i)
2483 return
2484 30 continue
2485 do 35 i = 1,n
2486 35 ewt(i) = rtol(i)*abs(ycur(i)) + atol(1)
2487 return
2488 40 continue
2489 do 45 i = 1,n
2490 45 ewt(i) = rtol(i)*abs(ycur(i)) + atol(i)
2491 return
2492 !----------------------- end of subroutine ewset -----------------------
2493 end subroutine ewset
2494 subroutine intdy (t, k, yh, nyh, dky, iflag)
2495 !lll. optimize
2496 integer k, nyh, iflag
2497 integer iownd, iowns, &
2498 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
2499 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
2500 integer i, ic, j, jb, jb2, jj, jj1, jp1
2501 real t, yh, dky
2502 real rowns, &
2503 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
2504 real c, r, s, tp
2505 !jdf dimension yh(nyh,1), dky(1)
2506 dimension yh(nyh,*), dky(*)
2507 common /ls0001/ rowns(209), &
2508 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, &
2509 iownd(14), iowns(6), &
2510 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
2511 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
2512 !-----------------------------------------------------------------------
2513 ! intdy computes interpolated values of the k-th derivative of the
2514 ! dependent variable vector y, and stores it in dky. this routine
2515 ! is called within the package with k = 0 and t = tout, but may
2516 ! also be called by the user for any k up to the current order.
2517 ! (see detailed instructions in the usage documentation.)
2518 !-----------------------------------------------------------------------
2519 ! the computed values in dky are gotten by interpolation using the
2520 ! nordsieck history array yh. this array corresponds uniquely to a
2521 ! vector-valued polynomial of degree nqcur or less, and dky is set
2522 ! to the k-th derivative of this polynomial at t.
2523 ! the formula for dky is..
2524 ! q
2525 ! dky(i) = sum c(j,k) * (t - tn)**(j-k) * h**(-j) * yh(i,j+1)
2526 ! j=k
2527 ! where c(j,k) = j*(j-1)*...*(j-k+1), q = nqcur, tn = tcur, h = hcur.
2528 ! the quantities nq = nqcur, l = nq+1, n = neq, tn, and h are
2529 ! communicated by common. the above sum is done in reverse order.
2530 ! iflag is returned negative if either k or t is out of bounds.
2531 !-----------------------------------------------------------------------
2532 iflag = 0
2533 if (k .lt. 0 .or. k .gt. nq) go to 80
2534 tp = tn - hu - 100.0e0*uround*(tn + hu)
2535 if ((t-tp)*(t-tn) .gt. 0.0e0) go to 90
2536 !
2537 s = (t - tn)/h
2538 ic = 1
2539 if (k .eq. 0) go to 15
2540 jj1 = l - k
2541 do 10 jj = jj1,nq
2542 10 ic = ic*jj
2543 15 c = float(ic)
2544 do 20 i = 1,n
2545 20 dky(i) = c*yh(i,l)
2546 if (k .eq. nq) go to 55
2547 jb2 = nq - k
2548 do 50 jb = 1,jb2
2549 j = nq - jb
2550 jp1 = j + 1
2551 ic = 1
2552 if (k .eq. 0) go to 35
2553 jj1 = jp1 - k
2554 do 30 jj = jj1,j
2555 30 ic = ic*jj
2556 35 c = float(ic)
2557 do 40 i = 1,n
2558 40 dky(i) = c*yh(i,jp1) + s*dky(i)
2559 50 continue
2560 if (k .eq. 0) return
2561 55 r = h**(-k)
2562 do 60 i = 1,n
2563 60 dky(i) = r*dky(i)
2564 return
2565 !
2566 80 call xerrwv('intdy-- k (=i1) illegal ', &
2567 30, 51, 0, 1, k, 0, 0, 0.0e0, 0.0e0)
2568 iflag = -1
2569 return
2570 90 call xerrwv('intdy-- t (=r1) illegal ', &
2571 30, 52, 0, 0, 0, 0, 1, t, 0.0e0)
2572 call xerrwv( &
2573 ' t not in interval tcur - hu (= r1) to tcur (=r2) ', &
2574 60, 52, 0, 0, 0, 0, 2, tp, tn)
2575 iflag = -2
2576 return
2577 !----------------------- end of subroutine intdy -----------------------
2578 end subroutine intdy
2579 subroutine iprep (neq, y, rwork, ia, ja, ipflag, f, jac, &
2580 ruserpar, nruserpar, iuserpar, niuserpar )
2581 !lll. optimize
2582 external f, jac
2583 integer neq, ia, ja, ipflag
2584 integer illin, init, lyh, lewt, lacor, lsavf, lwm, liwm, &
2585 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns
2586 integer icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
2587 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
2588 integer iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp, &
2589 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa, &
2590 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj, &
2591 nslj, ngp, nlu, nnz, nsp, nzl, nzu
2592 integer i, imax, lewtn, lyhd, lyhn
2593 integer nruserpar, iuserpar, niuserpar
2594 real y, rwork
2595 real rowns, &
2596 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
2597 real rlss
2598 real ruserpar
2599 !jdf dimension neq(1), y(1), rwork(1), ia(1), ja(1)
2600 dimension neq(*), y(*), rwork(*), ia(*), ja(*)
2601 dimension ruserpar(nruserpar), iuserpar(niuserpar)
2602 common /ls0001/ rowns(209), &
2603 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, &
2604 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm, &
2605 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6), &
2606 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
2607 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
2608 common /lss001/ rlss(6), &
2609 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp, &
2610 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa, &
2611 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj, &
2612 nslj, ngp, nlu, nnz, nsp, nzl, nzu
2613 !-----------------------------------------------------------------------
2614 ! this routine serves as an interface between the driver and
2615 ! subroutine prep. it is called only if miter is 1 or 2.
2616 ! tasks performed here are..
2617 ! * call prep,
2618 ! * reset the required wm segment length lenwk,
2619 ! * move yh back to its final location (following wm in rwork),
2620 ! * reset pointers for yh, savf, ewt, and acor, and
2621 ! * move ewt to its new position if istate = 1.
2622 ! ipflag is an output error indication flag. ipflag = 0 if there was
2623 ! no trouble, and ipflag is the value of the prep error flag ipper
2624 ! if there was trouble in subroutine prep.
2625 !-----------------------------------------------------------------------
2626 ipflag = 0
2627 ! call prep to do matrix preprocessing operations. ---------------------
2628 call prep_lsodes (neq, y, rwork(lyh), rwork(lsavf), rwork(lewt), &
2629 rwork(lacor), ia, ja, rwork(lwm), rwork(lwm), ipflag, f, jac, &
2630 ruserpar, nruserpar, iuserpar, niuserpar )
2631 lenwk = max0(lreq,lwmin)
2632 if (ipflag .lt. 0) return
2633 ! if prep was successful, move yh to end of required space for wm. -----
2634 lyhn = lwm + lenwk
2635 if (lyhn .gt. lyh) return
2636 lyhd = lyh - lyhn
2637 if (lyhd .eq. 0) go to 20
2638 imax = lyhn - 1 + lenyhm
2639 do 10 i = lyhn,imax
2640 10 rwork(i) = rwork(i+lyhd)
2641 lyh = lyhn
2642 ! reset pointers for savf, ewt, and acor. ------------------------------
2643 20 lsavf = lyh + lenyh
2644 lewtn = lsavf + n
2645 lacor = lewtn + n
2646 if (istatc .eq. 3) go to 40
2647 ! if istate = 1, move ewt (left) to its new position. ------------------
2648 if (lewtn .gt. lewt) return
2649 do 30 i = 1,n
2650 30 rwork(i+lewtn-1) = rwork(i+lewt-1)
2651 40 lewt = lewtn
2652 return
2653 !----------------------- end of subroutine iprep -----------------------
2654 end subroutine iprep
2655 subroutine jgroup (n,ia,ja,maxg,ngrp,igp,jgp,incl,jdone,ier)
2656 !lll. optimize
2657 integer n, ia, ja, maxg, ngrp, igp, jgp, incl, jdone, ier
2658 !jdf dimension ia(1), ja(1), igp(1), jgp(n), incl(n), jdone(n)
2659 dimension ia(*), ja(*), igp(*), jgp(n), incl(n), jdone(n)
2660 !-----------------------------------------------------------------------
2661 ! this subroutine constructs groupings of the column indices of
2662 ! the jacobian matrix, used in the numerical evaluation of the
2663 ! jacobian by finite differences.
2664 !
2665 ! input..
2666 ! n = the order of the matrix.
2667 ! ia,ja = sparse structure descriptors of the matrix by rows.
2668 ! maxg = length of available storate in the igp array.
2669 !
2670 ! output..
2671 ! ngrp = number of groups.
2672 ! jgp = array of length n containing the column indices by groups.
2673 ! igp = pointer array of length ngrp + 1 to the locations in jgp
2674 ! of the beginning of each group.
2675 ! ier = error indicator. ier = 0 if no error occurred, or 1 if
2676 ! maxg was insufficient.
2677 !
2678 ! incl and jdone are working arrays of length n.
2679 !-----------------------------------------------------------------------
2680 integer i, j, k, kmin, kmax, ncol, ng
2681 !
2682 ier = 0
2683 do 10 j = 1,n
2684 10 jdone(j) = 0
2685 ncol = 1
2686 do 60 ng = 1,maxg
2687 igp(ng) = ncol
2688 do 20 i = 1,n
2689 20 incl(i) = 0
2690 do 50 j = 1,n
2691 ! reject column j if it is already in a group.--------------------------
2692 if (jdone(j) .eq. 1) go to 50
2693 kmin = ia(j)
2694 kmax = ia(j+1) - 1
2695 do 30 k = kmin,kmax
2696 ! reject column j if it overlaps any column already in this group.------
2697 i = ja(k)
2698 if (incl(i) .eq. 1) go to 50
2699 30 continue
2700 ! accept column j into group ng.----------------------------------------
2701 jgp(ncol) = j
2702 ncol = ncol + 1
2703 jdone(j) = 1
2704 do 40 k = kmin,kmax
2705 i = ja(k)
2706 40 incl(i) = 1
2707 50 continue
2708 ! stop if this group is empty (grouping is complete).-------------------
2709 if (ncol .eq. igp(ng)) go to 70
2710 60 continue
2711 ! error return if not all columns were chosen (maxg too small).---------
2712 if (ncol .le. n) go to 80
2713 ng = maxg
2714 70 ngrp = ng - 1
2715 return
2716 80 ier = 1
2717 return
2718 !----------------------- end of subroutine jgroup ----------------------
2719 end subroutine jgroup
2720 subroutine md &
2721 (n, ia,ja, max, v,l, head,last,next, mark, flag)
2722 !lll. optimize
2723 !***********************************************************************
2724 ! md -- minimum degree algorithm (based on element model)
2725 !***********************************************************************
2726 !
2727 ! description
2728 !
2729 ! md finds a minimum degree ordering of the rows and columns of a
2730 ! general sparse matrix m stored in (ia,ja,a) format.
2731 ! when the structure of m is nonsymmetric, the ordering is that
2732 ! obtained for the symmetric matrix m + m-transpose.
2733 !
2734 !
2735 ! additional parameters
2736 !
2737 ! max - declared dimension of the one-dimensional arrays v and l.
2738 ! max must be at least n+2k, where k is the number of
2739 ! nonzeroes in the strict upper triangle of m + m-transpose
2740 !
2741 ! v - integer one-dimensional work array. dimension = max
2742 !
2743 ! l - integer one-dimensional work array. dimension = max
2744 !
2745 ! head - integer one-dimensional work array. dimension = n
2746 !
2747 ! last - integer one-dimensional array used to return the permutation
2748 ! of the rows and columns of m corresponding to the minimum
2749 ! degree ordering. dimension = n
2750 !
2751 ! next - integer one-dimensional array used to return the inverse of
2752 ! the permutation returned in last. dimension = n
2753 !
2754 ! mark - integer one-dimensional work array (may be the same as v).
2755 ! dimension = n
2756 !
2757 ! flag - integer error flag. values and their meanings are -
2758 ! 0 no errors detected
2759 ! 9n+k insufficient storage in md
2760 !
2761 !
2762 ! definitions of internal parameters
2763 !
2764 ! ---------+---------------------------------------------------------
2765 ! v(s) - value field of list entry
2766 ! ---------+---------------------------------------------------------
2767 ! l(s) - link field of list entry (0 =) end of list)
2768 ! ---------+---------------------------------------------------------
2769 ! l(vi) - pointer to element list of uneliminated vertex vi
2770 ! ---------+---------------------------------------------------------
2771 ! l(ej) - pointer to boundary list of active element ej
2772 ! ---------+---------------------------------------------------------
2773 ! head(d) - vj =) vj head of d-list d
2774 ! - 0 =) no vertex in d-list d
2775 !
2776 !
2777 ! - vi uneliminated vertex
2778 ! - vi in ek - vi not in ek
2779 ! ---------+-----------------------------+---------------------------
2780 ! next(vi) - undefined but nonnegative - vj =) vj next in d-list
2781 ! - - 0 =) vi tail of d-list
2782 ! ---------+-----------------------------+---------------------------
2783 ! last(vi) - (not set until mdp) - -d =) vi head of d-list d
2784 ! --vk =) compute degree - vj =) vj last in d-list
2785 ! - ej =) vi prototype of ej - 0 =) vi not in any d-list
2786 ! - 0 =) do not compute degree -
2787 ! ---------+-----------------------------+---------------------------
2788 ! mark(vi) - mark(vk) - nonneg. tag .lt. mark(vk)
2789 !
2790 !
2791 ! - vi eliminated vertex
2792 ! - ei active element - otherwise
2793 ! ---------+-----------------------------+---------------------------
2794 ! next(vi) - -j =) vi was j-th vertex - -j =) vi was j-th vertex
2795 ! - to be eliminated - to be eliminated
2796 ! ---------+-----------------------------+---------------------------
2797 ! last(vi) - m =) size of ei = m - undefined
2798 ! ---------+-----------------------------+---------------------------
2799 ! mark(vi) - -m =) overlap count of ei - undefined
2800 ! - with ek = m -
2801 ! - otherwise nonnegative tag -
2802 ! - .lt. mark(vk) -
2803 !
2804 !-----------------------------------------------------------------------
2805 !
2806 !jdf integer ia(1), ja(1), v(1), l(1), head(1), last(1), next(1),
2807 !jdf * mark(1), flag, tag, dmin, vk,ek, tail
2808 integer ia(*), ja(*), v(*), l(*), head(*), last(*), next(*), &
2809 mark(*), flag, tag, dmin, vk,ek, tail
2810 equivalence (vk,ek)
2811 !
2812 !----initialization
2813 tag = 0
2814 call mdi &
2815 (n, ia,ja, max,v,l, head,last,next, mark,tag, flag)
2816 if (flag.ne.0) return
2817 !
2818 k = 0
2819 dmin = 1
2820 !
2821 !----while k .lt. n do
2822 1 if (k.ge.n) go to 4
2823 !
2824 !------search for vertex of minimum degree
2825 2 if (head(dmin).gt.0) go to 3
2826 dmin = dmin + 1
2827 go to 2
2828 !
2829 !------remove vertex vk of minimum degree from degree list
2830 3 vk = head(dmin)
2831 head(dmin) = next(vk)
2832 if (head(dmin).gt.0) last(head(dmin)) = -dmin
2833 !
2834 !------number vertex vk, adjust tag, and tag vk
2835 k = k+1
2836 next(vk) = -k
2837 last(ek) = dmin - 1
2838 tag = tag + last(ek)
2839 mark(vk) = tag
2840 !
2841 !------form element ek from uneliminated neighbors of vk
2842 call mdm &
2843 (vk,tail, v,l, last,next, mark)
2844 !
2845 !------purge inactive elements and do mass elimination
2846 call mdp &
2847 (k,ek,tail, v,l, head,last,next, mark)
2848 !
2849 !------update degrees of uneliminated vertices in ek
2850 call mdu &
2851 (ek,dmin, v,l, head,last,next, mark)
2852 !
2853 go to 1
2854 !
2855 !----generate inverse permutation from permutation
2856 4 do 5 k=1,n
2857 next(k) = -next(k)
2858 5 last(next(k)) = k
2859 !
2860 return
2861 end subroutine md
2862 subroutine mdi &
2863 (n, ia,ja, max,v,l, head,last,next, mark,tag, flag)
2864 !lll. optimize
2865 !***********************************************************************
2866 ! mdi -- initialization
2867 !***********************************************************************
2868 !jdf integer ia(1), ja(1), v(1), l(1), head(1), last(1), next(1),
2869 !jdf * mark(1), tag, flag, sfs, vi,dvi, vj
2870 integer ia(*), ja(*), v(*), l(*), head(*), last(*), next(*), &
2871 mark(*), tag, flag, sfs, vi,dvi, vj
2872 !
2873 !----initialize degrees, element lists, and degree lists
2874 do 1 vi=1,n
2875 mark(vi) = 1
2876 l(vi) = 0
2877 1 head(vi) = 0
2878 sfs = n+1
2879 !
2880 !----create nonzero structure
2881 !----for each nonzero entry a(vi,vj)
2882 do 6 vi=1,n
2883 jmin = ia(vi)
2884 jmax = ia(vi+1) - 1
2885 if (jmin.gt.jmax) go to 6
2886 do 5 j=jmin,jmax
2887 vj = ja(j)
2888 if (vj-vi) 2, 5, 4
2889 !
2890 !------if a(vi,vj) is in strict lower triangle
2891 !------check for previous occurrence of a(vj,vi)
2892 2 lvk = vi
2893 kmax = mark(vi) - 1
2894 if (kmax .eq. 0) go to 4
2895 do 3 k=1,kmax
2896 lvk = l(lvk)
2897 if (v(lvk).eq.vj) go to 5
2898 3 continue
2899 !----for unentered entries a(vi,vj)
2900 4 if (sfs.ge.max) go to 101
2901 !
2902 !------enter vj in element list for vi
2903 mark(vi) = mark(vi) + 1
2904 v(sfs) = vj
2905 l(sfs) = l(vi)
2906 l(vi) = sfs
2907 sfs = sfs+1
2908 !
2909 !------enter vi in element list for vj
2910 mark(vj) = mark(vj) + 1
2911 v(sfs) = vi
2912 l(sfs) = l(vj)
2913 l(vj) = sfs
2914 sfs = sfs+1
2915 5 continue
2916 6 continue
2917 !
2918 !----create degree lists and initialize mark vector
2919 do 7 vi=1,n
2920 dvi = mark(vi)
2921 next(vi) = head(dvi)
2922 head(dvi) = vi
2923 last(vi) = -dvi
2924 nextvi = next(vi)
2925 if (nextvi.gt.0) last(nextvi) = vi
2926 7 mark(vi) = tag
2927 !
2928 return
2929 !
2930 ! ** error- insufficient storage
2931 101 flag = 9*n + vi
2932 return
2933 end subroutine mdi
2934 subroutine mdm &
2935 (vk,tail, v,l, last,next, mark)
2936 !lll. optimize
2937 !***********************************************************************
2938 ! mdm -- form element from uneliminated neighbors of vk
2939 !***********************************************************************
2940 !jdf integer vk, tail, v(1), l(1), last(1), next(1), mark(1),
2941 !jdf * tag, s,ls,vs,es, b,lb,vb, blp,blpmax
2942 integer vk, tail, v(*), l(*), last(*), next(*), mark(*), &
2943 tag, s,ls,vs,es, b,lb,vb, blp,blpmax
2944 equivalence (vs, es)
2945 !
2946 !----initialize tag and list of uneliminated neighbors
2947 tag = mark(vk)
2948 tail = vk
2949 !
2950 !----for each vertex/element vs/es in element list of vk
2951 ls = l(vk)
2952 1 s = ls
2953 if (s.eq.0) go to 5
2954 ls = l(s)
2955 vs = v(s)
2956 if (next(vs).lt.0) go to 2
2957 !
2958 !------if vs is uneliminated vertex, then tag and append to list of
2959 !------uneliminated neighbors
2960 mark(vs) = tag
2961 l(tail) = s
2962 tail = s
2963 go to 4
2964 !
2965 !------if es is active element, then ...
2966 !--------for each vertex vb in boundary list of element es
2967 2 lb = l(es)
2968 blpmax = last(es)
2969 do 3 blp=1,blpmax
2970 b = lb
2971 lb = l(b)
2972 vb = v(b)
2973 !
2974 !----------if vb is untagged vertex, then tag and append to list of
2975 !----------uneliminated neighbors
2976 if (mark(vb).ge.tag) go to 3
2977 mark(vb) = tag
2978 l(tail) = b
2979 tail = b
2980 3 continue
2981 !
2982 !--------mark es inactive
2983 mark(es) = tag
2984 !
2985 4 go to 1
2986 !
2987 !----terminate list of uneliminated neighbors
2988 5 l(tail) = 0
2989 !
2990 return
2991 end subroutine mdm
2992 subroutine mdp &
2993 (k,ek,tail, v,l, head,last,next, mark)
2994 !lll. optimize
2995 !***********************************************************************
2996 ! mdp -- purge inactive elements and do mass elimination
2997 !***********************************************************************
2998 !jdf integer ek, tail, v(1), l(1), head(1), last(1), next(1),
2999 !jdf * mark(1), tag, free, li,vi,lvi,evi, s,ls,es, ilp,ilpmax
3000 integer ek, tail, v(*), l(*), head(*), last(*), next(*), &
3001 mark(*), tag, free, li,vi,lvi,evi, s,ls,es, ilp,ilpmax
3002 !
3003 !----initialize tag
3004 tag = mark(ek)
3005 !
3006 !----for each vertex vi in ek
3007 li = ek
3008 ilpmax = last(ek)
3009 if (ilpmax.le.0) go to 12
3010 do 11 ilp=1,ilpmax
3011 i = li
3012 li = l(i)
3013 vi = v(li)
3014 !
3015 !------remove vi from degree list
3016 if (last(vi).eq.0) go to 3
3017 if (last(vi).gt.0) go to 1
3018 head(-last(vi)) = next(vi)
3019 go to 2
3020 1 next(last(vi)) = next(vi)
3021 2 if (next(vi).gt.0) last(next(vi)) = last(vi)
3022 !
3023 !------remove inactive items from element list of vi
3024 3 ls = vi
3025 4 s = ls
3026 ls = l(s)
3027 if (ls.eq.0) go to 6
3028 es = v(ls)
3029 if (mark(es).lt.tag) go to 5
3030 free = ls
3031 l(s) = l(ls)
3032 ls = s
3033 5 go to 4
3034 !
3035 !------if vi is interior vertex, then remove from list and eliminate
3036 6 lvi = l(vi)
3037 if (lvi.ne.0) go to 7
3038 l(i) = l(li)
3039 li = i
3040 !
3041 k = k+1
3042 next(vi) = -k
3043 last(ek) = last(ek) - 1
3044 go to 11
3045 !
3046 !------else ...
3047 !--------classify vertex vi
3048 7 if (l(lvi).ne.0) go to 9
3049 evi = v(lvi)
3050 if (next(evi).ge.0) go to 9
3051 if (mark(evi).lt.0) go to 8
3052 !
3053 !----------if vi is prototype vertex, then mark as such, initialize
3054 !----------overlap count for corresponding element, and move vi to end
3055 !----------of boundary list
3056 last(vi) = evi
3057 mark(evi) = -1
3058 l(tail) = li
3059 tail = li
3060 l(i) = l(li)
3061 li = i
3062 go to 10
3063 !
3064 !----------else if vi is duplicate vertex, then mark as such and adjust
3065 !----------overlap count for corresponding element
3066 8 last(vi) = 0
3067 mark(evi) = mark(evi) - 1
3068 go to 10
3069 !
3070 !----------else mark vi to compute degree
3071 9 last(vi) = -ek
3072 !
3073 !--------insert ek in element list of vi
3074 10 v(free) = ek
3075 l(free) = l(vi)
3076 l(vi) = free
3077 11 continue
3078 !
3079 !----terminate boundary list
3080 12 l(tail) = 0
3081 !
3082 return
3083 end subroutine mdp
3084 subroutine mdu &
3085 (ek,dmin, v,l, head,last,next, mark)
3086 !lll. optimize
3087 !***********************************************************************
3088 ! mdu -- update degrees of uneliminated vertices in ek
3089 !***********************************************************************
3090 !jdf integer ek, dmin, v(1), l(1), head(1), last(1), next(1),
3091 !jdf * mark(1), tag, vi,evi,dvi, s,vs,es, b,vb, ilp,ilpmax,
3092 !jdf * blp,blpmax
3093 integer ek, dmin, v(*), l(*), head(*), last(*), next(*), &
3094 mark(*), tag, vi,evi,dvi, s,vs,es, b,vb, ilp,ilpmax, &
3095 blp,blpmax
3096 equivalence (vs, es)
3097 !
3098 !----initialize tag
3099 tag = mark(ek) - last(ek)
3100 !
3101 !----for each vertex vi in ek
3102 i = ek
3103 ilpmax = last(ek)
3104 if (ilpmax.le.0) go to 11
3105 do 10 ilp=1,ilpmax
3106 i = l(i)
3107 vi = v(i)
3108 if (last(vi)) 1, 10, 8
3109 !
3110 !------if vi neither prototype nor duplicate vertex, then merge elements
3111 !------to compute degree
3112 1 tag = tag + 1
3113 dvi = last(ek)
3114 !
3115 !--------for each vertex/element vs/es in element list of vi
3116 s = l(vi)
3117 2 s = l(s)
3118 if (s.eq.0) go to 9
3119 vs = v(s)
3120 if (next(vs).lt.0) go to 3
3121 !
3122 !----------if vs is uneliminated vertex, then tag and adjust degree
3123 mark(vs) = tag
3124 dvi = dvi + 1
3125 go to 5
3126 !
3127 !----------if es is active element, then expand
3128 !------------check for outmatched vertex
3129 3 if (mark(es).lt.0) go to 6
3130 !
3131 !------------for each vertex vb in es
3132 b = es
3133 blpmax = last(es)
3134 do 4 blp=1,blpmax
3135 b = l(b)
3136 vb = v(b)
3137 !
3138 !--------------if vb is untagged, then tag and adjust degree
3139 if (mark(vb).ge.tag) go to 4
3140 mark(vb) = tag
3141 dvi = dvi + 1
3142 4 continue
3143 !
3144 5 go to 2
3145 !
3146 !------else if vi is outmatched vertex, then adjust overlaps but do not
3147 !------compute degree
3148 6 last(vi) = 0
3149 mark(es) = mark(es) - 1
3150 7 s = l(s)
3151 if (s.eq.0) go to 10
3152 es = v(s)
3153 if (mark(es).lt.0) mark(es) = mark(es) - 1
3154 go to 7
3155 !
3156 !------else if vi is prototype vertex, then calculate degree by
3157 !------inclusion/exclusion and reset overlap count
3158 8 evi = last(vi)
3159 dvi = last(ek) + last(evi) + mark(evi)
3160 mark(evi) = 0
3161 !
3162 !------insert vi in appropriate degree list
3163 9 next(vi) = head(dvi)
3164 head(dvi) = vi
3165 last(vi) = -dvi
3166 if (next(vi).gt.0) last(next(vi)) = vi
3167 if (dvi.lt.dmin) dmin = dvi
3168 !
3169 10 continue
3170 !
3171 11 return
3172 end subroutine mdu
3173 subroutine nnfc &
3174 (n, r,c,ic, ia,ja,a, z, b, &
3175 lmax,il,jl,ijl,l, d, umax,iu,ju,iju,u, &
3176 row, tmp, irl,jrl, flag)
3177 !lll. optimize
3178 !*** subroutine nnfc
3179 !*** numerical ldu-factorization of sparse nonsymmetric matrix and
3180 ! solution of system of linear equations (compressed pointer
3181 ! storage)
3182 !
3183 !
3184 ! input variables.. n, r, c, ic, ia, ja, a, b,
3185 ! il, jl, ijl, lmax, iu, ju, iju, umax
3186 ! output variables.. z, l, d, u, flag
3187 !
3188 ! parameters used internally..
3189 ! nia - irl, - vectors used to find the rows of l. at the kth step
3190 ! nia - jrl of the factorization, jrl(k) points to the head
3191 ! - of a linked list in jrl of column indices j
3192 ! - such j .lt. k and l(k,j) is nonzero. zero
3193 ! - indicates the end of the list. irl(j) (j.lt.k)
3194 ! - points to the smallest i such that i .ge. k and
3195 ! - l(i,j) is nonzero.
3196 ! - size of each = n.
3197 ! fia - row - holds intermediate values in calculation of u and l.
3198 ! - size = n.
3199 ! fia - tmp - holds new right-hand side b* for solution of the
3200 ! - equation ux = b*.
3201 ! - size = n.
3202 !
3203 ! internal variables..
3204 ! jmin, jmax - indices of the first and last positions in a row to
3205 ! be examined.
3206 ! sum - used in calculating tmp.
3207 !
3208 integer rk,umax
3209 !jdf integer r(1), c(1), ic(1), ia(1), ja(1), il(1), jl(1), ijl(1)
3210 !jdf integer iu(1), ju(1), iju(1), irl(1), jrl(1), flag
3211 !jdf real a(1), l(1), d(1), u(1), z(1), b(1), row(1)
3212 !jdf real tmp(1), lki, sum, dk
3213 integer r(*), c(*), ic(*), ia(*), ja(*), il(*), jl(*), ijl(*)
3214 integer iu(*), ju(*), iju(*), irl(*), jrl(*), flag
3215 real a(*), l(*), d(*), u(*), z(*), b(*), row(*)
3216 real tmp(*), lki, sum, dk
3217 ! double precision a(1), l(1), d(1), u(1), z(1), b(1), row(1)
3218 ! double precision tmp(1), lki, sum, dk
3219 !
3220 ! ****** initialize pointers and test storage ***********************
3221 if(il(n+1)-1 .gt. lmax) go to 104
3222 if(iu(n+1)-1 .gt. umax) go to 107
3223 do 1 k=1,n
3224 irl(k) = il(k)
3225 jrl(k) = 0
3226 1 continue
3227 !
3228 ! ****** for each row ***********************************************
3229 do 19 k=1,n
3230 ! ****** reverse jrl and zero row where kth row of l will fill in ***
3231 row(k) = 0
3232 i1 = 0
3233 if (jrl(k) .eq. 0) go to 3
3234 i = jrl(k)
3235 2 i2 = jrl(i)
3236 jrl(i) = i1
3237 i1 = i
3238 row(i) = 0
3239 i = i2
3240 if (i .ne. 0) go to 2
3241 ! ****** set row to zero where u will fill in ***********************
3242 3 jmin = iju(k)
3243 jmax = jmin + iu(k+1) - iu(k) - 1
3244 if (jmin .gt. jmax) go to 5
3245 do 4 j=jmin,jmax
3246 4 row(ju(j)) = 0
3247 ! ****** place kth row of a in row **********************************
3248 5 rk = r(k)
3249 jmin = ia(rk)
3250 jmax = ia(rk+1) - 1
3251 do 6 j=jmin,jmax
3252 row(ic(ja(j))) = a(j)
3253 6 continue
3254 ! ****** initialize sum, and link through jrl ***********************
3255 sum = b(rk)
3256 i = i1
3257 if (i .eq. 0) go to 10
3258 ! ****** assign the kth row of l and adjust row, sum ****************
3259 7 lki = -row(i)
3260 ! ****** if l is not required, then comment out the following line **
3261 l(irl(i)) = -lki
3262 sum = sum + lki * tmp(i)
3263 jmin = iu(i)
3264 jmax = iu(i+1) - 1
3265 if (jmin .gt. jmax) go to 9
3266 mu = iju(i) - jmin
3267 do 8 j=jmin,jmax
3268 8 row(ju(mu+j)) = row(ju(mu+j)) + lki * u(j)
3269 9 i = jrl(i)
3270 if (i .ne. 0) go to 7
3271 !
3272 ! ****** assign kth row of u and diagonal d, set tmp(k) *************
3273 10 if (row(k) .eq. 0.0e0) go to 108
3274 dk = 1.0e0 / row(k)
3275 d(k) = dk
3276 tmp(k) = sum * dk
3277 if (k .eq. n) go to 19
3278 jmin = iu(k)
3279 jmax = iu(k+1) - 1
3280 if (jmin .gt. jmax) go to 12
3281 mu = iju(k) - jmin
3282 do 11 j=jmin,jmax
3283 11 u(j) = row(ju(mu+j)) * dk
3284 12 continue
3285 !
3286 ! ****** update irl and jrl, keeping jrl in decreasing order ********
3287 i = i1
3288 if (i .eq. 0) go to 18
3289 14 irl(i) = irl(i) + 1
3290 i1 = jrl(i)
3291 if (irl(i) .ge. il(i+1)) go to 17
3292 ijlb = irl(i) - il(i) + ijl(i)
3293 j = jl(ijlb)
3294 15 if (i .gt. jrl(j)) go to 16
3295 j = jrl(j)
3296 go to 15
3297 16 jrl(i) = jrl(j)
3298 jrl(j) = i
3299 17 i = i1
3300 if (i .ne. 0) go to 14
3301 18 if (irl(k) .ge. il(k+1)) go to 19
3302 j = jl(ijl(k))
3303 jrl(k) = jrl(j)
3304 jrl(j) = k
3305 19 continue
3306 !
3307 ! ****** solve ux = tmp by back substitution **********************
3308 k = n
3309 do 22 i=1,n
3310 sum = tmp(k)
3311 jmin = iu(k)
3312 jmax = iu(k+1) - 1
3313 if (jmin .gt. jmax) go to 21
3314 mu = iju(k) - jmin
3315 do 20 j=jmin,jmax
3316 20 sum = sum - u(j) * tmp(ju(mu+j))
3317 21 tmp(k) = sum
3318 z(c(k)) = sum
3319 22 k = k-1
3320 flag = 0
3321 return
3322 !
3323 ! ** error.. insufficient storage for l
3324 104 flag = 4*n + 1
3325 return
3326 ! ** error.. insufficient storage for u
3327 107 flag = 7*n + 1
3328 return
3329 ! ** error.. zero pivot
3330 108 flag = 8*n + k
3331 return
3332 end subroutine nnfc
3333 subroutine nnsc &
3334 (n, r, c, il, jl, ijl, l, d, iu, ju, iju, u, z, b, tmp)
3335 !lll. optimize
3336 !*** subroutine nnsc
3337 !*** numerical solution of sparse nonsymmetric system of linear
3338 ! equations given ldu-factorization (compressed pointer storage)
3339 !
3340 !
3341 ! input variables.. n, r, c, il, jl, ijl, l, d, iu, ju, iju, u, b
3342 ! output variables.. z
3343 !
3344 ! parameters used internally..
3345 ! fia - tmp - temporary vector which gets result of solving ly = b.
3346 ! - size = n.
3347 !
3348 ! internal variables..
3349 ! jmin, jmax - indices of the first and last positions in a row of
3350 ! u or l to be used.
3351 !
3352 !jdf integer r(1), c(1), il(1), jl(1), ijl(1), iu(1), ju(1), iju(1)
3353 !jdf real l(1), d(1), u(1), b(1), z(1), tmp(1), tmpk, sum
3354 integer r(*), c(*), il(*), jl(*), ijl(*), iu(*), ju(*), iju(*)
3355 real l(*), d(*), u(*), b(*), z(*), tmp(*), tmpk, sum
3356 ! double precision l(1), d(1), u(1), b(1), z(1), tmp(1), tmpk,sum
3357 !
3358 ! ****** set tmp to reordered b *************************************
3359 do 1 k=1,n
3360 1 tmp(k) = b(r(k))
3361 ! ****** solve ly = b by forward substitution *********************
3362 do 3 k=1,n
3363 jmin = il(k)
3364 jmax = il(k+1) - 1
3365 tmpk = -d(k) * tmp(k)
3366 tmp(k) = -tmpk
3367 if (jmin .gt. jmax) go to 3
3368 ml = ijl(k) - jmin
3369 do 2 j=jmin,jmax
3370 2 tmp(jl(ml+j)) = tmp(jl(ml+j)) + tmpk * l(j)
3371 3 continue
3372 ! ****** solve ux = y by back substitution ************************
3373 k = n
3374 do 6 i=1,n
3375 sum = -tmp(k)
3376 jmin = iu(k)
3377 jmax = iu(k+1) - 1
3378 if (jmin .gt. jmax) go to 5
3379 mu = iju(k) - jmin
3380 do 4 j=jmin,jmax
3381 4 sum = sum + u(j) * tmp(ju(mu+j))
3382 5 tmp(k) = -sum
3383 z(c(k)) = -sum
3384 k = k - 1
3385 6 continue
3386 return
3387 end subroutine nnsc
3388 subroutine nntc &
3389 (n, r, c, il, jl, ijl, l, d, iu, ju, iju, u, z, b, tmp)
3390 !lll. optimize
3391 !*** subroutine nntc
3392 !*** numeric solution of the transpose of a sparse nonsymmetric system
3393 ! of linear equations given lu-factorization (compressed pointer
3394 ! storage)
3395 !
3396 !
3397 ! input variables.. n, r, c, il, jl, ijl, l, d, iu, ju, iju, u, b
3398 ! output variables.. z
3399 !
3400 ! parameters used internally..
3401 ! fia - tmp - temporary vector which gets result of solving ut y = b
3402 ! - size = n.
3403 !
3404 ! internal variables..
3405 ! jmin, jmax - indices of the first and last positions in a row of
3406 ! u or l to be used.
3407 !
3408 !jdf integer r(1), c(1), il(1), jl(1), ijl(1), iu(1), ju(1), iju(1)
3409 !jdf real l(1), d(1), u(1), b(1), z(1), tmp(1), tmpk,sum
3410 integer r(*), c(*), il(*), jl(*), ijl(*), iu(*), ju(*), iju(*)
3411 real l(*), d(*), u(*), b(*), z(*), tmp(*), tmpk,sum
3412 ! double precision l(1), d(1), u(1), b(1), z(1), tmp(1), tmpk,sum
3413 !
3414 ! ****** set tmp to reordered b *************************************
3415 do 1 k=1,n
3416 1 tmp(k) = b(c(k))
3417 ! ****** solve ut y = b by forward substitution *******************
3418 do 3 k=1,n
3419 jmin = iu(k)
3420 jmax = iu(k+1) - 1
3421 tmpk = -tmp(k)
3422 if (jmin .gt. jmax) go to 3
3423 mu = iju(k) - jmin
3424 do 2 j=jmin,jmax
3425 2 tmp(ju(mu+j)) = tmp(ju(mu+j)) + tmpk * u(j)
3426 3 continue
3427 ! ****** solve lt x = y by back substitution **********************
3428 k = n
3429 do 6 i=1,n
3430 sum = -tmp(k)
3431 jmin = il(k)
3432 jmax = il(k+1) - 1
3433 if (jmin .gt. jmax) go to 5
3434 ml = ijl(k) - jmin
3435 do 4 j=jmin,jmax
3436 4 sum = sum + l(j) * tmp(jl(ml+j))
3437 5 tmp(k) = -sum * d(k)
3438 z(r(k)) = tmp(k)
3439 k = k - 1
3440 6 continue
3441 return
3442 end subroutine nntc
3443 subroutine nroc (n, ic, ia, ja, a, jar, ar, p, flag)
3444 !lll. optimize
3445 !
3446 ! ----------------------------------------------------------------
3447 !
3448 ! yale sparse matrix package - nonsymmetric codes
3449 ! solving the system of equations mx = b
3450 !
3451 ! i. calling sequences
3452 ! the coefficient matrix can be processed by an ordering routine
3453 ! (e.g., to reduce fillin or ensure numerical stability) before using
3454 ! the remaining subroutines. if no reordering is done, then set
3455 ! r(i) = c(i) = ic(i) = i for i=1,...,n. if an ordering subroutine
3456 ! is used, then nroc should be used to reorder the coefficient matrix
3457 ! the calling sequence is --
3458 ! ( (matrix ordering))
3459 ! (nroc (matrix reordering))
3460 ! nsfc (symbolic factorization to determine where fillin will
3461 ! occur during numeric factorization)
3462 ! nnfc (numeric factorization into product ldu of unit lower
3463 ! triangular matrix l, diagonal matrix d, and unit
3464 ! upper triangular matrix u, and solution of linear
3465 ! system)
3466 ! nnsc (solution of linear system for additional right-hand
3467 ! side using ldu factorization from nnfc)
3468 ! (if only one system of equations is to be solved, then the
3469 ! subroutine trk should be used.)
3470 !
3471 ! ii. storage of sparse matrices
3472 ! the nonzero entries of the coefficient matrix m are stored
3473 ! row-by-row in the array a. to identify the individual nonzero
3474 ! entries in each row, we need to know in which column each entry
3475 ! lies. the column indices which correspond to the nonzero entries
3476 ! of m are stored in the array ja. i.e., if a(k) = m(i,j), then
3477 ! ja(k) = j. in addition, we need to know where each row starts and
3478 ! how long it is. the index positions in ja and a where the rows of
3479 ! m begin are stored in the array ia. i.e., if m(i,j) is the first
3480 ! (leftmost) entry in the i-th row and a(k) = m(i,j), then
3481 ! ia(i) = k. moreover, the index in ja and a of the first location
3482 ! following the last element in the last row is stored in ia(n+1).
3483 ! thus, the number of entries in the i-th row is given by
3484 ! ia(i+1) - ia(i), the nonzero entries of the i-th row are stored
3485 ! consecutively in
3486 ! a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1),
3487 ! and the corresponding column indices are stored consecutively in
3488 ! ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1).
3489 ! for example, the 5 by 5 matrix
3490 ! ( 1. 0. 2. 0. 0.)
3491 ! ( 0. 3. 0. 0. 0.)
3492 ! m = ( 0. 4. 5. 6. 0.)
3493 ! ( 0. 0. 0. 7. 0.)
3494 ! ( 0. 0. 0. 8. 9.)
3495 ! would be stored as
3496 ! - 1 2 3 4 5 6 7 8 9
3497 ! ---+--------------------------
3498 ! ia - 1 3 4 7 8 10
3499 ! ja - 1 3 2 2 3 4 4 4 5
3500 ! a - 1. 2. 3. 4. 5. 6. 7. 8. 9. .
3501 !
3502 ! the strict upper (lower) triangular portion of the matrix
3503 ! u (l) is stored in a similar fashion using the arrays iu, ju, u
3504 ! (il, jl, l) except that an additional array iju (ijl) is used to
3505 ! compress storage of ju (jl) by allowing some sequences of column
3506 ! (row) indices to used for more than one row (column) (n.b., l is
3507 ! stored by columns). iju(k) (ijl(k)) points to the starting
3508 ! location in ju (jl) of entries for the kth row (column).
3509 ! compression in ju (jl) occurs in two ways. first, if a row
3510 ! (column) i was merged into the current row (column) k, and the
3511 ! number of elements merged in from (the tail portion of) row
3512 ! (column) i is the same as the final length of row (column) k, then
3513 ! the kth row (column) and the tail of row (column) i are identical
3514 ! and iju(k) (ijl(k)) points to the start of the tail. second, if
3515 ! some tail portion of the (k-1)st row (column) is identical to the
3516 ! head of the kth row (column), then iju(k) (ijl(k)) points to the
3517 ! start of that tail portion. for example, the nonzero structure of
3518 ! the strict upper triangular part of the matrix
3519 ! d 0 x x x
3520 ! 0 d 0 x x
3521 ! 0 0 d x 0
3522 ! 0 0 0 d x
3523 ! 0 0 0 0 d
3524 ! would be represented as
3525 ! - 1 2 3 4 5 6
3526 ! ----+------------
3527 ! iu - 1 4 6 7 8 8
3528 ! ju - 3 4 5 4
3529 ! iju - 1 2 4 3 .
3530 ! the diagonal entries of l and u are assumed to be equal to one and
3531 ! are not stored. the array d contains the reciprocals of the
3532 ! diagonal entries of the matrix d.
3533 !
3534 ! iii. additional storage savings
3535 ! in nsfc, r and ic can be the same array in the calling
3536 ! sequence if no reordering of the coefficient matrix has been done.
3537 ! in nnfc, r, c, and ic can all be the same array if no
3538 ! reordering has been done. if only the rows have been reordered,
3539 ! then c and ic can be the same array. if the row and column
3540 ! orderings are the same, then r and c can be the same array. z and
3541 ! row can be the same array.
3542 ! in nnsc or nntc, r and c can be the same array if no
3543 ! reordering has been done or if the row and column orderings are the
3544 ! same. z and b can be the same array. however, then b will be
3545 ! destroyed.
3546 !
3547 ! iv. parameters
3548 ! following is a list of parameters to the programs. names are
3549 ! uniform among the various subroutines. class abbreviations are --
3550 ! n - integer variable
3551 ! f - real variable
3552 ! v - supplies a value to a subroutine
3553 ! r - returns a result from a subroutine
3554 ! i - used internally by a subroutine
3555 ! a - array
3556 !
3557 ! class - parameter
3558 ! ------+----------
3559 ! fva - a - nonzero entries of the coefficient matrix m, stored
3560 ! - by rows.
3561 ! - size = number of nonzero entries in m.
3562 ! fva - b - right-hand side b.
3563 ! - size = n.
3564 ! nva - c - ordering of the columns of m.
3565 ! - size = n.
3566 ! fvra - d - reciprocals of the diagonal entries of the matrix d.
3567 ! - size = n.
3568 ! nr - flag - error flag. values and their meanings are --
3569 ! - 0 no errors detected
3570 ! - n+k null row in a -- row = k
3571 ! - 2n+k duplicate entry in a -- row = k
3572 ! - 3n+k insufficient storage for jl -- row = k
3573 ! - 4n+1 insufficient storage for l
3574 ! - 5n+k null pivot -- row = k
3575 ! - 6n+k insufficient storage for ju -- row = k
3576 ! - 7n+1 insufficient storage for u
3577 ! - 8n+k zero pivot -- row = k
3578 ! nva - ia - pointers to delimit the rows of a.
3579 ! - size = n+1.
3580 ! nvra - ijl - pointers to the first element in each column in jl,
3581 ! - used to compress storage in jl.
3582 ! - size = n.
3583 ! nvra - iju - pointers to the first element in each row in ju, used
3584 ! - to compress storage in ju.
3585 ! - size = n.
3586 ! nvra - il - pointers to delimit the columns of l.
3587 ! - size = n+1.
3588 ! nvra - iu - pointers to delimit the rows of u.
3589 ! - size = n+1.
3590 ! nva - ja - column numbers corresponding to the elements of a.
3591 ! - size = size of a.
3592 ! nvra - jl - row numbers corresponding to the elements of l.
3593 ! - size = jlmax.
3594 ! nv - jlmax - declared dimension of jl. jlmax must be larger than
3595 ! - the number of nonzeros in the strict lower triangle
3596 ! - of m plus fillin minus compression.
3597 ! nvra - ju - column numbers corresponding to the elements of u.
3598 ! - size = jumax.
3599 ! nv - jumax - declared dimension of ju. jumax must be larger than
3600 ! - the number of nonzeros in the strict upper triangle
3601 ! - of m plus fillin minus compression.
3602 ! fvra - l - nonzero entries in the strict lower triangular portion
3603 ! - of the matrix l, stored by columns.
3604 ! - size = lmax.
3605 ! nv - lmax - declared dimension of l. lmax must be larger than
3606 ! - the number of nonzeros in the strict lower triangle
3607 ! - of m plus fillin (il(n+1)-1 after nsfc).
3608 ! nv - n - number of variables/equations.
3609 ! nva - r - ordering of the rows of m.
3610 ! - size = n.
3611 ! fvra - u - nonzero entries in the strict upper triangular portion
3612 ! - of the matrix u, stored by rows.
3613 ! - size = umax.
3614 ! nv - umax - declared dimension of u. umax must be larger than
3615 ! - the number of nonzeros in the strict upper triangle
3616 ! - of m plus fillin (iu(n+1)-1 after nsfc).
3617 ! fra - z - solution x.
3618 ! - size = n.
3619 !
3620 ! ----------------------------------------------------------------
3621 !
3622 !*** subroutine nroc
3623 !*** reorders rows of a, leaving row order unchanged
3624 !
3625 !
3626 ! input parameters.. n, ic, ia, ja, a
3627 ! output parameters.. ja, a, flag
3628 !
3629 ! parameters used internally..
3630 ! nia - p - at the kth step, p is a linked list of the reordered
3631 ! - column indices of the kth row of a. p(n+1) points
3632 ! - to the first entry in the list.
3633 ! - size = n+1.
3634 ! nia - jar - at the kth step,jar contains the elements of the
3635 ! - reordered column indices of a.
3636 ! - size = n.
3637 ! fia - ar - at the kth step, ar contains the elements of the
3638 ! - reordered row of a.
3639 ! - size = n.
3640 !
3641 !jdf integer ic(1), ia(1), ja(1), jar(1), p(1), flag
3642 !jdf real a(1), ar(1)
3643 integer ic(*), ia(*), ja(*), jar(*), p(*), flag
3644 real a(*), ar(*)
3645 ! double precision a(1), ar(1)
3646 !
3647 ! ****** for each nonempty row *******************************
3648 do 5 k=1,n
3649 jmin = ia(k)
3650 jmax = ia(k+1) - 1
3651 if(jmin .gt. jmax) go to 5
3652 p(n+1) = n + 1
3653 ! ****** insert each element in the list *********************
3654 do 3 j=jmin,jmax
3655 newj = ic(ja(j))
3656 i = n + 1
3657 1 if(p(i) .ge. newj) go to 2
3658 i = p(i)
3659 go to 1
3660 2 if(p(i) .eq. newj) go to 102
3661 p(newj) = p(i)
3662 p(i) = newj
3663 jar(newj) = ja(j)
3664 ar(newj) = a(j)
3665 3 continue
3666 ! ****** replace old row in ja and a *************************
3667 i = n + 1
3668 do 4 j=jmin,jmax
3669 i = p(i)
3670 ja(j) = jar(i)
3671 4 a(j) = ar(i)
3672 5 continue
3673 flag = 0
3674 return
3675 !
3676 ! ** error.. duplicate entry in a
3677 102 flag = n + k
3678 return
3679 end subroutine nroc
3680 subroutine nsfc &
3681 (n, r, ic, ia,ja, jlmax,il,jl,ijl, jumax,iu,ju,iju, &
3682 q, ira,jra, irac, irl,jrl, iru,jru, flag)
3683 !lll. optimize
3684 !*** subroutine nsfc
3685 !*** symbolic ldu-factorization of nonsymmetric sparse matrix
3686 ! (compressed pointer storage)
3687 !
3688 !
3689 ! input variables.. n, r, ic, ia, ja, jlmax, jumax.
3690 ! output variables.. il, jl, ijl, iu, ju, iju, flag.
3691 !
3692 ! parameters used internally..
3693 ! nia - q - suppose m* is the result of reordering m. if
3694 ! - processing of the ith row of m* (hence the ith
3695 ! - row of u) is being done, q(j) is initially
3696 ! - nonzero if m*(i,j) is nonzero (j.ge.i). since
3697 ! - values need not be stored, each entry points to the
3698 ! - next nonzero and q(n+1) points to the first. n+1
3699 ! - indicates the end of the list. for example, if n=9
3700 ! - and the 5th row of m* is
3701 ! - 0 x x 0 x 0 0 x 0
3702 ! - then q will initially be
3703 ! - a a a a 8 a a 10 5 (a - arbitrary).
3704 ! - as the algorithm proceeds, other elements of q
3705 ! - are inserted in the list because of fillin.
3706 ! - q is used in an analogous manner to compute the
3707 ! - ith column of l.
3708 ! - size = n+1.
3709 ! nia - ira, - vectors used to find the columns of m. at the kth
3710 ! nia - jra, step of the factorization, irac(k) points to the
3711 ! nia - irac head of a linked list in jra of row indices i
3712 ! - such that i .ge. k and m(i,k) is nonzero. zero
3713 ! - indicates the end of the list. ira(i) (i.ge.k)
3714 ! - points to the smallest j such that j .ge. k and
3715 ! - m(i,j) is nonzero.
3716 ! - size of each = n.
3717 ! nia - irl, - vectors used to find the rows of l. at the kth step
3718 ! nia - jrl of the factorization, jrl(k) points to the head
3719 ! - of a linked list in jrl of column indices j
3720 ! - such j .lt. k and l(k,j) is nonzero. zero
3721 ! - indicates the end of the list. irl(j) (j.lt.k)
3722 ! - points to the smallest i such that i .ge. k and
3723 ! - l(i,j) is nonzero.
3724 ! - size of each = n.
3725 ! nia - iru, - vectors used in a manner analogous to irl and jrl
3726 ! nia - jru to find the columns of u.
3727 ! - size of each = n.
3728 !
3729 ! internal variables..
3730 ! jlptr - points to the last position used in jl.
3731 ! juptr - points to the last position used in ju.
3732 ! jmin,jmax - are the indices in a or u of the first and last
3733 ! elements to be examined in a given row.
3734 ! for example, jmin=ia(k), jmax=ia(k+1)-1.
3735 !
3736 integer cend, qm, rend, rk, vj
3737 !jdf integer ia(1), ja(1), ira(1), jra(1), il(1), jl(1), ijl(1)
3738 !jdf integer iu(1), ju(1), iju(1), irl(1), jrl(1), iru(1), jru(1)
3739 !jdf integer r(1), ic(1), q(1), irac(1), flag
3740 integer ia(*), ja(*), ira(*), jra(*), il(*), jl(*), ijl(*)
3741 integer iu(*), ju(*), iju(*), irl(*), jrl(*), iru(*), jru(*)
3742 integer r(*), ic(*), q(*), irac(*), flag
3743 !
3744 ! ****** initialize pointers ****************************************
3745 np1 = n + 1
3746 jlmin = 1
3747 jlptr = 0
3748 il(1) = 1
3749 jumin = 1
3750 juptr = 0
3751 iu(1) = 1
3752 do 1 k=1,n
3753 irac(k) = 0
3754 jra(k) = 0
3755 jrl(k) = 0
3756 1 jru(k) = 0
3757 ! ****** initialize column pointers for a ***************************
3758 do 2 k=1,n
3759 rk = r(k)
3760 iak = ia(rk)
3761 if (iak .ge. ia(rk+1)) go to 101
3762 jaiak = ic(ja(iak))
3763 if (jaiak .gt. k) go to 105
3764 jra(k) = irac(jaiak)
3765 irac(jaiak) = k
3766 2 ira(k) = iak
3767 !
3768 ! ****** for each column of l and row of u **************************
3769 do 41 k=1,n
3770 !
3771 ! ****** initialize q for computing kth column of l *****************
3772 q(np1) = np1
3773 luk = -1
3774 ! ****** by filling in kth column of a ******************************
3775 vj = irac(k)
3776 if (vj .eq. 0) go to 5
3777 3 qm = np1
3778 4 m = qm
3779 qm = q(m)
3780 if (qm .lt. vj) go to 4
3781 if (qm .eq. vj) go to 102
3782 luk = luk + 1
3783 q(m) = vj
3784 q(vj) = qm
3785 vj = jra(vj)
3786 if (vj .ne. 0) go to 3
3787 ! ****** link through jru *******************************************
3788 5 lastid = 0
3789 lasti = 0
3790 ijl(k) = jlptr
3791 i = k
3792 6 i = jru(i)
3793 if (i .eq. 0) go to 10
3794 qm = np1
3795 jmin = irl(i)
3796 jmax = ijl(i) + il(i+1) - il(i) - 1
3797 long = jmax - jmin
3798 if (long .lt. 0) go to 6
3799 jtmp = jl(jmin)
3800 if (jtmp .ne. k) long = long + 1
3801 if (jtmp .eq. k) r(i) = -r(i)
3802 if (lastid .ge. long) go to 7
3803 lasti = i
3804 lastid = long
3805 ! ****** and merge the corresponding columns into the kth column ****
3806 7 do 9 j=jmin,jmax
3807 vj = jl(j)
3808 8 m = qm
3809 qm = q(m)
3810 if (qm .lt. vj) go to 8
3811 if (qm .eq. vj) go to 9
3812 luk = luk + 1
3813 q(m) = vj
3814 q(vj) = qm
3815 qm = vj
3816 9 continue
3817 go to 6
3818 ! ****** lasti is the longest column merged into the kth ************
3819 ! ****** see if it equals the entire kth column *********************
3820 10 qm = q(np1)
3821 if (qm .ne. k) go to 105
3822 if (luk .eq. 0) go to 17
3823 if (lastid .ne. luk) go to 11
3824 ! ****** if so, jl can be compressed ********************************
3825 irll = irl(lasti)
3826 ijl(k) = irll + 1
3827 if (jl(irll) .ne. k) ijl(k) = ijl(k) - 1
3828 go to 17
3829 ! ****** if not, see if kth column can overlap the previous one *****
3830 11 if (jlmin .gt. jlptr) go to 15
3831 qm = q(qm)
3832 do 12 j=jlmin,jlptr
3833 if (jl(j) - qm) 12, 13, 15
3834 12 continue
3835 go to 15
3836 13 ijl(k) = j
3837 do 14 i=j,jlptr
3838 if (jl(i) .ne. qm) go to 15
3839 qm = q(qm)
3840 if (qm .gt. n) go to 17
3841 14 continue
3842 jlptr = j - 1
3843 ! ****** move column indices from q to jl, update vectors ***********
3844 15 jlmin = jlptr + 1
3845 ijl(k) = jlmin
3846 if (luk .eq. 0) go to 17
3847 jlptr = jlptr + luk
3848 if (jlptr .gt. jlmax) go to 103
3849 qm = q(np1)
3850 do 16 j=jlmin,jlptr
3851 qm = q(qm)
3852 16 jl(j) = qm
3853 17 irl(k) = ijl(k)
3854 il(k+1) = il(k) + luk
3855 !
3856 ! ****** initialize q for computing kth row of u ********************
3857 q(np1) = np1
3858 luk = -1
3859 ! ****** by filling in kth row of reordered a ***********************
3860 rk = r(k)
3861 jmin = ira(k)
3862 jmax = ia(rk+1) - 1
3863 if (jmin .gt. jmax) go to 20
3864 do 19 j=jmin,jmax
3865 vj = ic(ja(j))
3866 qm = np1
3867 18 m = qm
3868 qm = q(m)
3869 if (qm .lt. vj) go to 18
3870 if (qm .eq. vj) go to 102
3871 luk = luk + 1
3872 q(m) = vj
3873 q(vj) = qm
3874 19 continue
3875 ! ****** link through jrl, ******************************************
3876 20 lastid = 0
3877 lasti = 0
3878 iju(k) = juptr
3879 i = k
3880 i1 = jrl(k)
3881 21 i = i1
3882 if (i .eq. 0) go to 26
3883 i1 = jrl(i)
3884 qm = np1
3885 jmin = iru(i)
3886 jmax = iju(i) + iu(i+1) - iu(i) - 1
3887 long = jmax - jmin
3888 if (long .lt. 0) go to 21
3889 jtmp = ju(jmin)
3890 if (jtmp .eq. k) go to 22
3891 ! ****** update irl and jrl, *****************************************
3892 long = long + 1
3893 cend = ijl(i) + il(i+1) - il(i)
3894 irl(i) = irl(i) + 1
3895 if (irl(i) .ge. cend) go to 22
3896 j = jl(irl(i))
3897 jrl(i) = jrl(j)
3898 jrl(j) = i
3899 22 if (lastid .ge. long) go to 23
3900 lasti = i
3901 lastid = long
3902 ! ****** and merge the corresponding rows into the kth row **********
3903 23 do 25 j=jmin,jmax
3904 vj = ju(j)
3905 24 m = qm
3906 qm = q(m)
3907 if (qm .lt. vj) go to 24
3908 if (qm .eq. vj) go to 25
3909 luk = luk + 1
3910 q(m) = vj
3911 q(vj) = qm
3912 qm = vj
3913 25 continue
3914 go to 21
3915 ! ****** update jrl(k) and irl(k) ***********************************
3916 26 if (il(k+1) .le. il(k)) go to 27
3917 j = jl(irl(k))
3918 jrl(k) = jrl(j)
3919 jrl(j) = k
3920 ! ****** lasti is the longest row merged into the kth ***************
3921 ! ****** see if it equals the entire kth row ************************
3922 27 qm = q(np1)
3923 if (qm .ne. k) go to 105
3924 if (luk .eq. 0) go to 34
3925 if (lastid .ne. luk) go to 28
3926 ! ****** if so, ju can be compressed ********************************
3927 irul = iru(lasti)
3928 iju(k) = irul + 1
3929 if (ju(irul) .ne. k) iju(k) = iju(k) - 1
3930 go to 34
3931 ! ****** if not, see if kth row can overlap the previous one ********
3932 28 if (jumin .gt. juptr) go to 32
3933 qm = q(qm)
3934 do 29 j=jumin,juptr
3935 if (ju(j) - qm) 29, 30, 32
3936 29 continue
3937 go to 32
3938 30 iju(k) = j
3939 do 31 i=j,juptr
3940 if (ju(i) .ne. qm) go to 32
3941 qm = q(qm)
3942 if (qm .gt. n) go to 34
3943 31 continue
3944 juptr = j - 1
3945 ! ****** move row indices from q to ju, update vectors **************
3946 32 jumin = juptr + 1
3947 iju(k) = jumin
3948 if (luk .eq. 0) go to 34
3949 juptr = juptr + luk
3950 if (juptr .gt. jumax) go to 106
3951 qm = q(np1)
3952 do 33 j=jumin,juptr
3953 qm = q(qm)
3954 33 ju(j) = qm
3955 34 iru(k) = iju(k)
3956 iu(k+1) = iu(k) + luk
3957 !
3958 ! ****** update iru, jru ********************************************
3959 i = k
3960 35 i1 = jru(i)
3961 if (r(i) .lt. 0) go to 36
3962 rend = iju(i) + iu(i+1) - iu(i)
3963 if (iru(i) .ge. rend) go to 37
3964 j = ju(iru(i))
3965 jru(i) = jru(j)
3966 jru(j) = i
3967 go to 37
3968 36 r(i) = -r(i)
3969 37 i = i1
3970 if (i .eq. 0) go to 38
3971 iru(i) = iru(i) + 1
3972 go to 35
3973 !
3974 ! ****** update ira, jra, irac **************************************
3975 38 i = irac(k)
3976 if (i .eq. 0) go to 41
3977 39 i1 = jra(i)
3978 ira(i) = ira(i) + 1
3979 if (ira(i) .ge. ia(r(i)+1)) go to 40
3980 irai = ira(i)
3981 jairai = ic(ja(irai))
3982 if (jairai .gt. i) go to 40
3983 jra(i) = irac(jairai)
3984 irac(jairai) = i
3985 40 i = i1
3986 if (i .ne. 0) go to 39
3987 41 continue
3988 !
3989 ijl(n) = jlptr
3990 iju(n) = juptr
3991 flag = 0
3992 return
3993 !
3994 ! ** error.. null row in a
3995 101 flag = n + rk
3996 return
3997 ! ** error.. duplicate entry in a
3998 102 flag = 2*n + rk
3999 return
4000 ! ** error.. insufficient storage for jl
4001 103 flag = 3*n + k
4002 return
4003 ! ** error.. null pivot
4004 105 flag = 5*n + k
4005 return
4006 ! ** error.. insufficient storage for ju
4007 106 flag = 6*n + k
4008 return
4009 end subroutine nsfc
4010 subroutine odrv &
4011 (n, ia,ja,a, p,ip, nsp,isp, path, flag)
4012 !lll. optimize
4013 ! 5/2/83
4014 !***********************************************************************
4015 ! odrv -- driver for sparse matrix reordering routines
4016 !***********************************************************************
4017 !
4018 ! description
4019 !
4020 ! odrv finds a minimum degree ordering of the rows and columns
4021 ! of a matrix m stored in (ia,ja,a) format (see below). for the
4022 ! reordered matrix, the work and storage required to perform
4023 ! gaussian elimination is (usually) significantly less.
4024 !
4025 ! note.. odrv and its subordinate routines have been modified to
4026 ! compute orderings for general matrices, not necessarily having any
4027 ! symmetry. the miminum degree ordering is computed for the
4028 ! structure of the symmetric matrix m + m-transpose.
4029 ! modifications to the original odrv module have been made in
4030 ! the coding in subroutine mdi, and in the initial comments in
4031 ! subroutines odrv and md.
4032 !
4033 ! if only the nonzero entries in the upper triangle of m are being
4034 ! stored, then odrv symmetrically reorders (ia,ja,a), (optionally)
4035 ! with the diagonal entries placed first in each row. this is to
4036 ! ensure that if m(i,j) will be in the upper triangle of m with
4037 ! respect to the new ordering, then m(i,j) is stored in row i (and
4038 ! thus m(j,i) is not stored), whereas if m(i,j) will be in the
4039 ! strict lower triangle of m, then m(j,i) is stored in row j (and
4040 ! thus m(i,j) is not stored).
4041 !
4042 !
4043 ! storage of sparse matrices
4044 !
4045 ! the nonzero entries of the matrix m are stored row-by-row in the
4046 ! array a. to identify the individual nonzero entries in each row,
4047 ! we need to know in which column each entry lies. these column
4048 ! indices are stored in the array ja. i.e., if a(k) = m(i,j), then
4049 ! ja(k) = j. to identify the individual rows, we need to know where
4050 ! each row starts. these row pointers are stored in the array ia.
4051 ! i.e., if m(i,j) is the first nonzero entry (stored) in the i-th row
4052 ! and a(k) = m(i,j), then ia(i) = k. moreover, ia(n+1) points to
4053 ! the first location following the last element in the last row.
4054 ! thus, the number of entries in the i-th row is ia(i+1) - ia(i),
4055 ! the nonzero entries in the i-th row are stored consecutively in
4056 !
4057 ! a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1),
4058 !
4059 ! and the corresponding column indices are stored consecutively in
4060 !
4061 ! ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1).
4062 !
4063 ! when the coefficient matrix is symmetric, only the nonzero entries
4064 ! in the upper triangle need be stored. for example, the matrix
4065 !
4066 ! ( 1 0 2 3 0 )
4067 ! ( 0 4 0 0 0 )
4068 ! m = ( 2 0 5 6 0 )
4069 ! ( 3 0 6 7 8 )
4070 ! ( 0 0 0 8 9 )
4071 !
4072 ! could be stored as
4073 !
4074 ! - 1 2 3 4 5 6 7 8 9 10 11 12 13
4075 ! ---+--------------------------------------
4076 ! ia - 1 4 5 8 12 14
4077 ! ja - 1 3 4 2 1 3 4 1 3 4 5 4 5
4078 ! a - 1 2 3 4 2 5 6 3 6 7 8 8 9
4079 !
4080 ! or (symmetrically) as
4081 !
4082 ! - 1 2 3 4 5 6 7 8 9
4083 ! ---+--------------------------
4084 ! ia - 1 4 5 7 9 10
4085 ! ja - 1 3 4 2 3 4 4 5 5
4086 ! a - 1 2 3 4 5 6 7 8 9 .
4087 !
4088 !
4089 ! parameters
4090 !
4091 ! n - order of the matrix
4092 !
4093 ! ia - integer one-dimensional array containing pointers to delimit
4094 ! rows in ja and a. dimension = n+1
4095 !
4096 ! ja - integer one-dimensional array containing the column indices
4097 ! corresponding to the elements of a. dimension = number of
4098 ! nonzero entries in (the upper triangle of) m
4099 !
4100 ! a - real one-dimensional array containing the nonzero entries in
4101 ! (the upper triangle of) m, stored by rows. dimension =
4102 ! number of nonzero entries in (the upper triangle of) m
4103 !
4104 ! p - integer one-dimensional array used to return the permutation
4105 ! of the rows and columns of m corresponding to the minimum
4106 ! degree ordering. dimension = n
4107 !
4108 ! ip - integer one-dimensional array used to return the inverse of
4109 ! the permutation returned in p. dimension = n
4110 !
4111 ! nsp - declared dimension of the one-dimensional array isp. nsp
4112 ! must be at least 3n+4k, where k is the number of nonzeroes
4113 ! in the strict upper triangle of m
4114 !
4115 ! isp - integer one-dimensional array used for working storage.
4116 ! dimension = nsp
4117 !
4118 ! path - integer path specification. values and their meanings are -
4119 ! 1 find minimum degree ordering only
4120 ! 2 find minimum degree ordering and reorder symmetrically
4121 ! stored matrix (used when only the nonzero entries in
4122 ! the upper triangle of m are being stored)
4123 ! 3 reorder symmetrically stored matrix as specified by
4124 ! input permutation (used when an ordering has already
4125 ! been determined and only the nonzero entries in the
4126 ! upper triangle of m are being stored)
4127 ! 4 same as 2 but put diagonal entries at start of each row
4128 ! 5 same as 3 but put diagonal entries at start of each row
4129 !
4130 ! flag - integer error flag. values and their meanings are -
4131 ! 0 no errors detected
4132 ! 9n+k insufficient storage in md
4133 ! 10n+1 insufficient storage in odrv
4134 ! 11n+1 illegal path specification
4135 !
4136 !
4137 ! conversion from real to double precision
4138 !
4139 ! change the real declarations in odrv and sro to double precision
4140 ! declarations.
4141 !
4142 !-----------------------------------------------------------------------
4143 !
4144 !jdf integer ia(1), ja(1), p(1), ip(1), isp(1), path, flag,
4145 !jdf * v, l, head, tmp, q
4146 !jdf real a(1)
4147 integer ia(*), ja(*), p(*), ip(*), isp(*), path, flag, &
4148 v, l, head, tmp, q
4149 real a(*)
4150 !... double precision a(1)
4151 logical dflag
4152 !
4153 !----initialize error flag and validate path specification
4154 flag = 0
4155 if (path.lt.1 .or. 5.lt.path) go to 111
4156 !
4157 !----allocate storage and find minimum degree ordering
4158 if ((path-1) * (path-2) * (path-4) .ne. 0) go to 1
4159 max = (nsp-n)/2
4160 v = 1
4161 l = v + max
4162 head = l + max
4163 next = head + n
4164 if (max.lt.n) go to 110
4165 !
4166 call md &
4167 (n, ia,ja, max,isp(v),isp(l), isp(head),p,ip, isp(v), flag)
4168 if (flag.ne.0) go to 100
4169 !
4170 !----allocate storage and symmetrically reorder matrix
4171 1 if ((path-2) * (path-3) * (path-4) * (path-5) .ne. 0) go to 2
4172 tmp = (nsp+1) - n
4173 q = tmp - (ia(n+1)-1)
4174 if (q.lt.1) go to 110
4175 !
4176 dflag = path.eq.4 .or. path.eq.5
4177 call sro &
4178 (n, ip, ia, ja, a, isp(tmp), isp(q), dflag)
4179 !
4180 2 return
4181 !
4182 ! ** error -- error detected in md
4183 100 return
4184 ! ** error -- insufficient storage
4185 110 flag = 10*n + 1
4186 return
4187 ! ** error -- illegal path specified
4188 111 flag = 11*n + 1
4189 return
4190 end subroutine odrv
4191
4192
4193
4194 subroutine prjs (neq,y,yh,nyh,ewt,ftem,savf,wk,iwk,f,jac, &
4195 ruserpar, nruserpar, iuserpar, niuserpar )
4196 !lll. optimize
4197 external f,jac
4198 integer neq, nyh, iwk
4199 integer iownd, iowns, &
4200 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
4201 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
4202 integer iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp, &
4203 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa, &
4204 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj, &
4205 nslj, ngp, nlu, nnz, nsp, nzl, nzu
4206 integer i, imul, j, jj, jok, jmax, jmin, k, kmax, kmin, ng
4207 integer nruserpar, iuserpar, niuserpar
4208 real y, yh, ewt, ftem, savf, wk
4209 real rowns, &
4210 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
4211 real con0, conmin, ccmxj, psmall, rbig, seth
4212 !rce real con, di, fac, hl0, pij, r, r0, rcon, rcont, &
4213 !rce srur, vnorm
4214 real con, di, fac, hl0, pij, r, r0, rcon, rcont, &
4215 srur
4216 real ruserpar
4217 !jdf dimension neq(1), y(1), yh(nyh,1), ewt(1), ftem(1), savf(1),
4218 !jdf 1 wk(1), iwk(1)
4219 dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*), &
4220 wk(*), iwk(*)
4221 dimension ruserpar(nruserpar), iuserpar(niuserpar)
4222 common /ls0001/ rowns(209), &
4223 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, &
4224 iownd(14), iowns(6), &
4225 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
4226 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
4227 common /lss001/ con0, conmin, ccmxj, psmall, rbig, seth, &
4228 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp, &
4229 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa, &
4230 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj, &
4231 nslj, ngp, nlu, nnz, nsp, nzl, nzu
4232 !-----------------------------------------------------------------------
4233 ! prjs is called to compute and process the matrix
4234 ! p = i - h*el(1)*j , where j is an approximation to the jacobian.
4235 ! j is computed by columns, either by the user-supplied routine jac
4236 ! if miter = 1, or by finite differencing if miter = 2.
4237 ! if miter = 3, a diagonal approximation to j is used.
4238 ! if miter = 1 or 2, and if the existing value of the jacobian
4239 ! (as contained in p) is considered acceptable, then a new value of
4240 ! p is reconstructed from the old value. in any case, when miter
4241 ! is 1 or 2, the p matrix is subjected to lu decomposition in cdrv.
4242 ! p and its lu decomposition are stored (separately) in wk.
4243 !
4244 ! in addition to variables described previously, communication
4245 ! with prjs uses the following..
4246 ! y = array containing predicted values on entry.
4247 ! ftem = work array of length n (acor in stode).
4248 ! savf = array containing f evaluated at predicted y.
4249 ! wk = real work space for matrices. on output it contains the
4250 ! inverse diagonal matrix if miter = 3, and p and its sparse
4251 ! lu decomposition if miter is 1 or 2.
4252 ! storage of matrix elements starts at wk(3).
4253 ! wk also contains the following matrix-related data..
4254 ! wk(1) = sqrt(uround), used in numerical jacobian increments.
4255 ! wk(2) = h*el0, saved for later use if miter = 3.
4256 ! iwk = integer work space for matrix-related data, assumed to
4257 ! be equivalenced to wk. in addition, wk(iprsp) and iwk(ipisp)
4258 ! are assumed to have identical locations.
4259 ! el0 = el(1) (input).
4260 ! ierpj = output error flag (in common).
4261 ! = 0 if no error.
4262 ! = 1 if zero pivot found in cdrv.
4263 ! = 2 if a singular matrix arose with miter = 3.
4264 ! = -1 if insufficient storage for cdrv (should not occur here).
4265 ! = -2 if other error found in cdrv (should not occur here).
4266 ! jcur = output flag = 1 to indicate that the jacobian matrix
4267 ! (or approximation) is now current.
4268 ! this routine also uses other variables in common.
4269 !-----------------------------------------------------------------------
4270 hl0 = h*el0
4271 con = -hl0
4272 if (miter .eq. 3) go to 300
4273 ! see whether j should be reevaluated (jok = 0) or not (jok = 1). ------
4274 jok = 1
4275 if (nst .eq. 0 .or. nst .ge. nslj+msbj) jok = 0
4276 if (icf .eq. 1 .and. abs(rc - 1.0e0) .lt. ccmxj) jok = 0
4277 if (icf .eq. 2) jok = 0
4278 if (jok .eq. 1) go to 250
4279 !
4280 ! miter = 1 or 2, and the jacobian is to be reevaluated. ---------------
4281 20 jcur = 1
4282 nje = nje + 1
4283 nslj = nst
4284 iplost = 0
4285 conmin = abs(con)
4286 go to (100, 200), miter
4287 !
4288 ! if miter = 1, call jac, multiply by scalar, and add identity. --------
4289 100 continue
4290 kmin = iwk(ipian)
4291 do 130 j = 1, n
4292 kmax = iwk(ipian+j) - 1
4293 do 110 i = 1,n
4294 110 ftem(i) = 0.0e0
4295 call jac (neq, tn, y, j, iwk(ipian), iwk(ipjan), ftem, &
4296 ruserpar, nruserpar, iuserpar, niuserpar)
4297 do 120 k = kmin, kmax
4298 i = iwk(ibjan+k)
4299 wk(iba+k) = ftem(i)*con
4300 if (i .eq. j) wk(iba+k) = wk(iba+k) + 1.0e0
4301 120 continue
4302 kmin = kmax + 1
4303 130 continue
4304 go to 290
4305 !
4306 ! if miter = 2, make ngp calls to f to approximate j and p. ------------
4307 200 continue
4308 fac = vnorm(n, savf, ewt)
4309 r0 = 1000.0e0 * abs(h) * uround * float(n) * fac
4310 if (r0 .eq. 0.0e0) r0 = 1.0e0
4311 srur = wk(1)
4312 jmin = iwk(ipigp)
4313 do 240 ng = 1,ngp
4314 jmax = iwk(ipigp+ng) - 1
4315 do 210 j = jmin,jmax
4316 jj = iwk(ibjgp+j)
4317 r = amax1(srur*abs(y(jj)),r0/ewt(jj))
4318 210 y(jj) = y(jj) + r
4319 call f (neq, tn, y, ftem, &
4320 ruserpar, nruserpar, iuserpar, niuserpar)
4321 do 230 j = jmin,jmax
4322 jj = iwk(ibjgp+j)
4323 y(jj) = yh(jj,1)
4324 r = amax1(srur*abs(y(jj)),r0/ewt(jj))
4325 fac = -hl0/r
4326 kmin =iwk(ibian+jj)
4327 kmax =iwk(ibian+jj+1) - 1
4328 do 220 k = kmin,kmax
4329 i = iwk(ibjan+k)
4330 wk(iba+k) = (ftem(i) - savf(i))*fac
4331 if (i .eq. jj) wk(iba+k) = wk(iba+k) + 1.0e0
4332 220 continue
4333 230 continue
4334 jmin = jmax + 1
4335 240 continue
4336 nfe = nfe + ngp
4337 go to 290
4338 !
4339 ! if jok = 1, reconstruct new p from old p. ----------------------------
4340 250 jcur = 0
4341 rcon = con/con0
4342 rcont = abs(con)/conmin
4343 if (rcont .gt. rbig .and. iplost .eq. 1) go to 20
4344 kmin = iwk(ipian)
4345 do 275 j = 1,n
4346 kmax = iwk(ipian+j) - 1
4347 do 270 k = kmin,kmax
4348 i = iwk(ibjan+k)
4349 pij = wk(iba+k)
4350 if (i .ne. j) go to 260
4351 pij = pij - 1.0e0
4352 if (abs(pij) .ge. psmall) go to 260
4353 iplost = 1
4354 conmin = amin1(abs(con0),conmin)
4355 260 pij = pij*rcon
4356 if (i .eq. j) pij = pij + 1.0e0
4357 wk(iba+k) = pij
4358 270 continue
4359 kmin = kmax + 1
4360 275 continue
4361 !
4362 ! do numerical factorization of p matrix. ------------------------------
4363 290 nlu = nlu + 1
4364 con0 = con
4365 ierpj = 0
4366 do 295 i = 1,n
4367 295 ftem(i) = 0.0e0
4368 call cdrv (n,iwk(ipr),iwk(ipc),iwk(ipic),iwk(ipian),iwk(ipjan), &
4369 wk(ipa),ftem,ftem,nsp,iwk(ipisp),wk(iprsp),iesp,2,iys)
4370 if (iys .eq. 0) return
4371 imul = (iys - 1)/n
4372 ierpj = -2
4373 if (imul .eq. 8) ierpj = 1
4374 if (imul .eq. 10) ierpj = -1
4375 return
4376 !
4377 ! if miter = 3, construct a diagonal approximation to j and p. ---------
4378 300 continue
4379 jcur = 1
4380 nje = nje + 1
4381 wk(2) = hl0
4382 ierpj = 0
4383 r = el0*0.1e0
4384 do 310 i = 1,n
4385 310 y(i) = y(i) + r*(h*savf(i) - yh(i,2))
4386 call f (neq, tn, y, wk(3), &
4387 ruserpar, nruserpar, iuserpar, niuserpar)
4388 nfe = nfe + 1
4389 do 320 i = 1,n
4390 r0 = h*savf(i) - yh(i,2)
4391 di = 0.1e0*r0 - h*(wk(i+2) - savf(i))
4392 wk(i+2) = 1.0e0
4393 if (abs(r0) .lt. uround/ewt(i)) go to 320
4394 if (abs(di) .eq. 0.0e0) go to 330
4395 wk(i+2) = 0.1e0*r0/di
4396 320 continue
4397 return
4398 330 ierpj = 2
4399 return
4400 !----------------------- end of subroutine prjs ------------------------
4401 end subroutine prjs
4402 subroutine slss (wk, iwk, x, tem)
4403 !lll. optimize
4404 integer iwk
4405 integer iownd, iowns, &
4406 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
4407 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
4408 integer iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp, &
4409 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa, &
4410 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj, &
4411 nslj, ngp, nlu, nnz, nsp, nzl, nzu
4412 integer i
4413 real wk, x, tem
4414 real rowns, &
4415 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
4416 real rlss
4417 real di, hl0, phl0, r
4418 !jdf dimension wk(1), iwk(1), x(1), tem(1)
4419 dimension wk(*), iwk(*), x(*), tem(*)
4420
4421 common /ls0001/ rowns(209), &
4422 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, &
4423 iownd(14), iowns(6), &
4424 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
4425 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
4426 common /lss001/ rlss(6), &
4427 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp, &
4428 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa, &
4429 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj, &
4430 nslj, ngp, nlu, nnz, nsp, nzl, nzu
4431 !-----------------------------------------------------------------------
4432 ! this routine manages the solution of the linear system arising from
4433 ! a chord iteration. it is called if miter .ne. 0.
4434 ! if miter is 1 or 2, it calls cdrv to accomplish this.
4435 ! if miter = 3 it updates the coefficient h*el0 in the diagonal
4436 ! matrix, and then computes the solution.
4437 ! communication with slss uses the following variables..
4438 ! wk = real work space containing the inverse diagonal matrix if
4439 ! miter = 3 and the lu decomposition of the matrix otherwise.
4440 ! storage of matrix elements starts at wk(3).
4441 ! wk also contains the following matrix-related data..
4442 ! wk(1) = sqrt(uround) (not used here),
4443 ! wk(2) = hl0, the previous value of h*el0, used if miter = 3.
4444 ! iwk = integer work space for matrix-related data, assumed to
4445 ! be equivalenced to wk. in addition, wk(iprsp) and iwk(ipisp)
4446 ! are assumed to have identical locations.
4447 ! x = the right-hand side vector on input, and the solution vector
4448 ! on output, of length n.
4449 ! tem = vector of work space of length n, not used in this version.
4450 ! iersl = output flag (in common).
4451 ! iersl = 0 if no trouble occurred.
4452 ! iersl = -1 if cdrv returned an error flag (miter = 1 or 2).
4453 ! this should never occur and is considered fatal.
4454 ! iersl = 1 if a singular matrix arose with miter = 3.
4455 ! this routine also uses other variables in common.
4456 !-----------------------------------------------------------------------
4457 iersl = 0
4458 go to (100, 100, 300), miter
4459 100 call cdrv (n,iwk(ipr),iwk(ipc),iwk(ipic),iwk(ipian),iwk(ipjan), &
4460 wk(ipa),x,x,nsp,iwk(ipisp),wk(iprsp),iesp,4,iersl)
4461 if (iersl .ne. 0) iersl = -1
4462 return
4463 !
4464 300 phl0 = wk(2)
4465 hl0 = h*el0
4466 wk(2) = hl0
4467 if (hl0 .eq. phl0) go to 330
4468 r = hl0/phl0
4469 do 320 i = 1,n
4470 di = 1.0e0 - r*(1.0e0 - 1.0e0/wk(i+2))
4471 if (abs(di) .eq. 0.0e0) go to 390
4472 320 wk(i+2) = 1.0e0/di
4473 330 do 340 i = 1,n
4474 340 x(i) = wk(i+2)*x(i)
4475 return
4476 390 iersl = 1
4477 return
4478 !
4479 !----------------------- end of subroutine slss ------------------------
4480 end subroutine slss
4481 subroutine sro &
4482 (n, ip, ia,ja,a, q, r, dflag)
4483 !lll. optimize
4484 !***********************************************************************
4485 ! sro -- symmetric reordering of sparse symmetric matrix
4486 !***********************************************************************
4487 !
4488 ! description
4489 !
4490 ! the nonzero entries of the matrix m are assumed to be stored
4491 ! symmetrically in (ia,ja,a) format (i.e., not both m(i,j) and m(j,i)
4492 ! are stored if i ne j).
4493 !
4494 ! sro does not rearrange the order of the rows, but does move
4495 ! nonzeroes from one row to another to ensure that if m(i,j) will be
4496 ! in the upper triangle of m with respect to the new ordering, then
4497 ! m(i,j) is stored in row i (and thus m(j,i) is not stored), whereas
4498 ! if m(i,j) will be in the strict lower triangle of m, then m(j,i) is
4499 ! stored in row j (and thus m(i,j) is not stored).
4500 !
4501 !
4502 ! additional parameters
4503 !
4504 ! q - integer one-dimensional work array. dimension = n
4505 !
4506 ! r - integer one-dimensional work array. dimension = number of
4507 ! nonzero entries in the upper triangle of m
4508 !
4509 ! dflag - logical variable. if dflag = .true., then store nonzero
4510 ! diagonal elements at the beginning of the row
4511 !
4512 !-----------------------------------------------------------------------
4513 !
4514 !jdf integer ip(1), ia(1), ja(1), q(1), r(1)
4515 !jdf real a(1), ak
4516 integer ip(*), ia(*), ja(*), q(*), r(*)
4517 real a(*), ak
4518 !... double precision a(1), ak
4519 logical dflag
4520 !
4521 !
4522 !--phase 1 -- find row in which to store each nonzero
4523 !----initialize count of nonzeroes to be stored in each row
4524 do 1 i=1,n
4525 1 q(i) = 0
4526 !
4527 !----for each nonzero element a(j)
4528 do 3 i=1,n
4529 jmin = ia(i)
4530 jmax = ia(i+1) - 1
4531 if (jmin.gt.jmax) go to 3
4532 do 2 j=jmin,jmax
4533 !
4534 !--------find row (=r(j)) and column (=ja(j)) in which to store a(j) ...
4535 k = ja(j)
4536 if (ip(k).lt.ip(i)) ja(j) = i
4537 if (ip(k).ge.ip(i)) k = i
4538 r(j) = k
4539 !
4540 !--------... and increment count of nonzeroes (=q(r(j)) in that row
4541 2 q(k) = q(k) + 1
4542 3 continue
4543 !
4544 !
4545 !--phase 2 -- find new ia and permutation to apply to (ja,a)
4546 !----determine pointers to delimit rows in permuted (ja,a)
4547 do 4 i=1,n
4548 ia(i+1) = ia(i) + q(i)
4549 4 q(i) = ia(i+1)
4550 !
4551 !----determine where each (ja(j),a(j)) is stored in permuted (ja,a)
4552 !----for each nonzero element (in reverse order)
4553 ilast = 0
4554 jmin = ia(1)
4555 jmax = ia(n+1) - 1
4556 j = jmax
4557 do 6 jdummy=jmin,jmax
4558 i = r(j)
4559 if (.not.dflag .or. ja(j).ne.i .or. i.eq.ilast) go to 5
4560 !
4561 !------if dflag, then put diagonal nonzero at beginning of row
4562 r(j) = ia(i)
4563 ilast = i
4564 go to 6
4565 !
4566 !------put (off-diagonal) nonzero in last unused location in row
4567 5 q(i) = q(i) - 1
4568 r(j) = q(i)
4569 !
4570 6 j = j-1
4571 !
4572 !
4573 !--phase 3 -- permute (ja,a) to upper triangular form (wrt new ordering)
4574 do 8 j=jmin,jmax
4575 7 if (r(j).eq.j) go to 8
4576 k = r(j)
4577 r(j) = r(k)
4578 r(k) = k
4579 jak = ja(k)
4580 ja(k) = ja(j)
4581 ja(j) = jak
4582 ak = a(k)
4583 a(k) = a(j)
4584 a(j) = ak
4585 go to 7
4586 8 continue
4587 !
4588 return
4589 end subroutine sro
4590
4591
4592
4593 real function vnorm (n, v, w)
4594 !lll. optimize
4595 !-----------------------------------------------------------------------
4596 ! this function routine computes the weighted root-mean-square norm
4597 ! of the vector of length n contained in the array v, with weights
4598 ! contained in the array w of length n..
4599 ! vnorm = sqrt( (1/n) * sum( v(i)*w(i) )**2 )
4600 !-----------------------------------------------------------------------
4601 integer n, i
4602 real v, w, sum
4603 dimension v(n), w(n)
4604 sum = 0.0e0
4605 do 10 i = 1,n
4606 10 sum = sum + (v(i)*w(i))**2
4607 vnorm = sqrt(sum/float(n))
4608 return
4609 !----------------------- end of function vnorm -------------------------
4610 end function vnorm
4611 subroutine xerrwv (msg, nmes, nerr, level, ni, i1, i2, nr, r1, r2)
4612 use module_peg_util, only: peg_message, peg_error_fatal
4613 ! integer msg, nmes, nerr, level, ni, i1, i2, nr, &
4614 integer nmes, nerr, level, ni, i1, i2, nr, &
4615 i, lun, lunit, mesflg, ncpw, nch, nwds
4616 real r1, r2
4617 character(*) msg
4618 character*80 errmsg
4619 !-----------------------------------------------------------------------
4620 ! subroutines xerrwv, xsetf, and xsetun, as given here, constitute
4621 ! a simplified version of the slatec error handling package.
4622 ! written by a. c. hindmarsh at llnl. version of march 30, 1987.
4623 !
4624 ! all arguments are input arguments.
4625 !
4626 ! msg = the message (hollerith literal or integer array).
4627 ! nmes = the length of msg (number of characters).
4628 ! nerr = the error number (not used).
4629 ! level = the error level..
4630 ! 0 or 1 means recoverable (control returns to caller).
4631 ! 2 means fatal (run is aborted--see note below).
4632 ! ni = number of integers (0, 1, or 2) to be printed with message.
4633 ! i1,i2 = integers to be printed, depending on ni.
4634 ! nr = number of reals (0, 1, or 2) to be printed with message.
4635 ! r1,r2 = reals to be printed, depending on nr.
4636 !
4637 ! note.. this routine is machine-dependent and specialized for use
4638 ! in limited context, in the following ways..
4639 ! 1. the number of hollerith characters stored per word, denoted
4640 ! by ncpw below, is a data-loaded constant.
4641 ! 2. the value of nmes is assumed to be at most 60.
4642 ! (multi-line messages are generated by repeated calls.)
4643 ! 3. if level = 2, control passes to the statement stop
4644 ! to abort the run. this statement may be machine-dependent.
4645 ! 4. r1 and r2 are assumed to be in single precision and are printed
4646 ! in e21.13 format.
4647 ! 5. the common block /eh0001/ below is data-loaded (a machine-
4648 ! dependent feature) with default values.
4649 ! this block is needed for proper retention of parameters used by
4650 ! this routine which the user can reset by calling xsetf or xsetun.
4651 ! the variables in this block are as follows..
4652 ! mesflg = print control flag..
4653 ! 1 means print all messages (the default).
4654 ! 0 means no printing.
4655 ! lunit = logical unit number for messages.
4656 ! the default is 6 (machine-dependent).
4657 !-----------------------------------------------------------------------
4658 ! the following are instructions for installing this routine
4659 ! in different machine environments.
4660 !
4661 ! to change the default output unit, change the data statement below.
4662 !
4663 ! for some systems, the data statement below must be replaced
4664 ! by a separate block data subprogram.
4665 !
4666 ! for a different number of characters per word, change the
4667 ! data statement setting ncpw below, and format 10. alternatives for
4668 ! various computers are shown in comment cards.
4669 !
4670 ! for a different run-abort command, change the statement following
4671 ! statement 100 at the end.
4672 !-----------------------------------------------------------------------
4673 common /eh0001/ mesflg, lunit
4674 !
4675 !raz data mesflg/1/, lunit/6/
4676 mesflg = 1
4677 lunit = 6
4678 !-----------------------------------------------------------------------
4679 ! the following data-loaded value of ncpw is valid for the cdc-6600
4680 ! and cdc-7600 computers.
4681 ! data ncpw/10/
4682 ! the following is valid for the cray-1 computer.
4683 ! data ncpw/8/
4684 ! the following is valid for the burroughs 6700 and 7800 computers.
4685 ! data ncpw/6/
4686 ! the following is valid for the pdp-10 computer.
4687 ! data ncpw/5/
4688 ! the following is valid for the vax computer with 4 bytes per integer,
4689 ! and for the ibm-360, ibm-370, ibm-303x, and ibm-43xx computers.
4690 data ncpw/4/
4691 ! the following is valid for the pdp-11, or vax with 2-byte integers.
4692 ! data ncpw/2/
4693 !-----------------------------------------------------------------------
4694 !
4695 if (mesflg .eq. 0) go to 100
4696 ! get logical unit number. ---------------------------------------------
4697 lun = lunit
4698 ! get number of words in message. --------------------------------------
4699 nch = min0(nmes,60)
4700 nwds = nch/ncpw
4701 if (nch .ne. nwds*ncpw) nwds = nwds + 1
4702 ! write the message. ---------------------------------------------------
4703 ! write (lun, 10) (msg(i),i=1,nwds)
4704 ! write (lun, 10) msg
4705 call peg_message( lun, msg )
4706 !-----------------------------------------------------------------------
4707 ! the following format statement is to have the form
4708 ! 10 format(1x,mmann)
4709 ! where nn = ncpw and mm is the smallest integer .ge. 60/ncpw.
4710 ! the following is valid for ncpw = 10.
4711 ! 10 format(1x,6a10)
4712 ! the following is valid for ncpw = 8.
4713 ! 10 format(1x,8a8)
4714 ! the following is valid for ncpw = 6.
4715 ! 10 format(1x,10a6)
4716 ! the following is valid for ncpw = 5.
4717 ! 10 format(1x,12a5)
4718 ! the following is valid for ncpw = 4.
4719 ! 10 format(1x,15a4)
4720 10 format(1x,a)
4721 ! the following is valid for ncpw = 2.
4722 ! 10 format(1x,30a2)
4723 !-----------------------------------------------------------------------
4724 errmsg = ' '
4725 ! if (ni .eq. 1) write (lun, 20) i1
4726 if (ni .eq. 1) write (errmsg, 20) i1
4727 20 format(6x,23hin above message, i1 =,i10)
4728
4729 ! if (ni .eq. 2) write (lun, 30) i1,i2
4730 if (ni .eq. 2) write (errmsg, 30) i1,i2
4731 30 format(6x,23hin above message, i1 =,i10,3x,4hi2 =,i10)
4732
4733 ! if (nr .eq. 1) write (lun, 40) r1
4734 if (nr .eq. 1) write (errmsg, 40) r1
4735 40 format(6x,23hin above message, r1 =,e21.13)
4736
4737 ! if (nr .eq. 2) write (lun, 50) r1,r2
4738 if (nr .eq. 2) write (errmsg, 50) r1,r2
4739 50 format(6x,15hin above, r1 =,e21.13,3x,4hr2 =,e21.13)
4740
4741 if (errmsg .ne. ' ') call peg_message( lun, errmsg )
4742
4743 ! abort the run if level = 2. ------------------------------------------
4744 100 if (level .ne. 2) return
4745 call peg_error_fatal( lun, '*** subr xerrwv fatal error' )
4746 stop
4747 !----------------------- end of subroutine xerrwv ----------------------
4748 end subroutine xerrwv
4749 !-----------------------------------------------------------------------
4750 real function r1mach(i)
4751 use module_peg_util, only: peg_error_fatal
4752 !
4753 ! single-precision machine constants
4754 !
4755 ! r1mach(1) = b**(emin-1), the smallest positive magnitude.
4756 !
4757 ! r1mach(2) = b**emax*(1 - b**(-t)), the largest magnitude.
4758 !
4759 ! r1mach(3) = b**(-t), the smallest relative spacing.
4760 !
4761 ! r1mach(4) = b**(1-t), the largest relative spacing.
4762 !
4763 ! r1mach(5) = log10(b)
4764 !
4765 ! to alter this function for a particular environment,
4766 ! the desired set of data statements should be activated by
4767 ! removing the c from column 1.
4768 ! on rare machines a static statement may need to be added.
4769 ! (but probably more systems prohibit it than require it.)
4770 !
4771 ! for ieee-arithmetic machines (binary standard), the first
4772 ! set of constants below should be appropriate.
4773 !
4774 ! where possible, decimal, octal or hexadecimal constants are used
4775 ! to specify the constants exactly. sometimes this requires using
4776 ! equivalent integer arrays. if your compiler uses half-word
4777 ! integers by default (sometimes called integer*2), you may need to
4778 ! change integer to integer*4 or otherwise instruct your compiler
4779 ! to use full-word integers in the next 5 declarations.
4780 !
4781 integer small(2)
4782 integer large(2)
4783 integer right(2)
4784 integer diver(2)
4785 integer log10(2)
4786 integer sc
4787 !
4788 character*80 errmsg
4789 !
4790 real rmach(5)
4791 !
4792 equivalence (rmach(1),small(1))
4793 equivalence (rmach(2),large(1))
4794 equivalence (rmach(3),right(1))
4795 equivalence (rmach(4),diver(1))
4796 equivalence (rmach(5),log10(1))
4797 !
4798 ! machine constants for ieee arithmetic machines, such as the at&t
4799 ! 3b series, motorola 68000 based machines (e.g. sun 3 and at&t
4800 ! pc 7300), and 8087 based micros (e.g. ibm pc and at&t 6300).
4801 !
4802 data small(1) / 8388608 /
4803 data large(1) / 2139095039 /
4804 data right(1) / 864026624 /
4805 data diver(1) / 872415232 /
4806 data log10(1) / 1050288283 /, sc/987/
4807 !
4808 ! machine constants for amdahl machines.
4809 !
4810 ! data small(1) / 1048576 /
4811 ! data large(1) / 2147483647 /
4812 ! data right(1) / 990904320 /
4813 ! data diver(1) / 1007681536 /
4814 ! data log10(1) / 1091781651 /, sc/987/
4815 !
4816 ! machine constants for the burroughs 1700 system.
4817 !
4818 ! data rmach(1) / z400800000 /
4819 ! data rmach(2) / z5ffffffff /
4820 ! data rmach(3) / z4e9800000 /
4821 ! data rmach(4) / z4ea800000 /
4822 ! data rmach(5) / z500e730e8 /, sc/987/
4823 !
4824 ! machine constants for the burroughs 5700/6700/7700 systems.
4825 !
4826 ! data rmach(1) / o1771000000000000 /
4827 ! data rmach(2) / o0777777777777777 /
4828 ! data rmach(3) / o1311000000000000 /
4829 ! data rmach(4) / o1301000000000000 /
4830 ! data rmach(5) / o1157163034761675 /, sc/987/
4831 !
4832 ! machine constants for ftn4 on the cdc 6000/7000 series.
4833 !
4834 ! data rmach(1) / 00564000000000000000b /
4835 ! data rmach(2) / 37767777777777777776b /
4836 ! data rmach(3) / 16414000000000000000b /
4837 ! data rmach(4) / 16424000000000000000b /
4838 ! data rmach(5) / 17164642023241175720b /, sc/987/
4839 !
4840 ! machine constants for ftn5 on the cdc 6000/7000 series.
4841 !
4842 ! data rmach(1) / o"00564000000000000000" /
4843 ! data rmach(2) / o"37767777777777777776" /
4844 ! data rmach(3) / o"16414000000000000000" /
4845 ! data rmach(4) / o"16424000000000000000" /
4846 ! data rmach(5) / o"17164642023241175720" /, sc/987/
4847 !
4848 ! machine constants for convex c-1.
4849 !
4850 ! data rmach(1) / '00800000'x /
4851 ! data rmach(2) / '7fffffff'x /
4852 ! data rmach(3) / '34800000'x /
4853 ! data rmach(4) / '35000000'x /
4854 ! data rmach(5) / '3f9a209b'x /, sc/987/
4855 !
4856 ! machine constants for the cray 1, xmp, 2, and 3.
4857 !
4858 ! data rmach(1) / 200034000000000000000b /
4859 ! data rmach(2) / 577767777777777777776b /
4860 ! data rmach(3) / 377224000000000000000b /
4861 ! data rmach(4) / 377234000000000000000b /
4862 ! data rmach(5) / 377774642023241175720b /, sc/987/
4863 !
4864 ! machine constants for the data general eclipse s/200.
4865 !
4866 ! note - it may be appropriate to include the following line -
4867 ! static rmach(5)
4868 !
4869 ! data small/20k,0/,large/77777k,177777k/
4870 ! data right/35420k,0/,diver/36020k,0/
4871 ! data log10/40423k,42023k/, sc/987/
4872 !
4873 ! machine constants for the harris slash 6 and slash 7.
4874 !
4875 ! data small(1),small(2) / '20000000, '00000201 /
4876 ! data large(1),large(2) / '37777777, '00000177 /
4877 ! data right(1),right(2) / '20000000, '00000352 /
4878 ! data diver(1),diver(2) / '20000000, '00000353 /
4879 ! data log10(1),log10(2) / '23210115, '00000377 /, sc/987/
4880 !
4881 ! machine constants for the honeywell dps 8/70 series.
4882 !
4883 ! data rmach(1) / o402400000000 /
4884 ! data rmach(2) / o376777777777 /
4885 ! data rmach(3) / o714400000000 /
4886 ! data rmach(4) / o716400000000 /
4887 ! data rmach(5) / o776464202324 /, sc/987/
4888 !
4889 ! machine constants for the ibm 360/370 series,
4890 ! the xerox sigma 5/7/9 and the sel systems 85/86.
4891 !
4892 ! data rmach(1) / z00100000 /
4893 ! data rmach(2) / z7fffffff /
4894 ! data rmach(3) / z3b100000 /
4895 ! data rmach(4) / z3c100000 /
4896 ! data rmach(5) / z41134413 /, sc/987/
4897 !
4898 ! machine constants for the interdata 8/32
4899 ! with the unix system fortran 77 compiler.
4900 !
4901 ! for the interdata fortran vii compiler replace
4902 ! the z's specifying hex constants with y's.
4903 !
4904 ! data rmach(1) / z'00100000' /
4905 ! data rmach(2) / z'7effffff' /
4906 ! data rmach(3) / z'3b100000' /
4907 ! data rmach(4) / z'3c100000' /
4908 ! data rmach(5) / z'41134413' /, sc/987/
4909 !
4910 ! machine constants for the pdp-10 (ka or ki processor).
4911 !----------------------------------------------------------------------
4912 ! rce 2004-01-07
4913 ! The following 5 lines for rmach(1-5) each contained one
4914 ! quotation-mark character.
4915 ! The WRF preprocessor did not like this, so I changed the
4916 ! quotation-mark characters to QUOTE.
4917 !
4918 ! data rmach(1) / QUOTE000400000000 /
4919 ! data rmach(2) / QUOTE377777777777 /
4920 ! data rmach(3) / QUOTE146400000000 /
4921 ! data rmach(4) / QUOTE147400000000 /
4922 ! data rmach(5) / QUOTE177464202324 /, sc/987/
4923 !----------------------------------------------------------------------
4924 !
4925 ! machine constants for pdp-11 fortrans supporting
4926 ! 32-bit integers (expressed in integer and octal).
4927 !
4928 ! data small(1) / 8388608 /
4929 ! data large(1) / 2147483647 /
4930 ! data right(1) / 880803840 /
4931 ! data diver(1) / 889192448 /
4932 ! data log10(1) / 1067065499 /, sc/987/
4933 !
4934 ! data rmach(1) / o00040000000 /
4935 ! data rmach(2) / o17777777777 /
4936 ! data rmach(3) / o06440000000 /
4937 ! data rmach(4) / o06500000000 /
4938 ! data rmach(5) / o07746420233 /, sc/987/
4939 !
4940 ! machine constants for pdp-11 fortrans supporting
4941 ! 16-bit integers (expressed in integer and octal).
4942 !
4943 ! data small(1),small(2) / 128, 0 /
4944 ! data large(1),large(2) / 32767, -1 /
4945 ! data right(1),right(2) / 13440, 0 /
4946 ! data diver(1),diver(2) / 13568, 0 /
4947 ! data log10(1),log10(2) / 16282, 8347 /, sc/987/
4948 !
4949 ! data small(1),small(2) / o000200, o000000 /
4950 ! data large(1),large(2) / o077777, o177777 /
4951 ! data right(1),right(2) / o032200, o000000 /
4952 ! data diver(1),diver(2) / o032400, o000000 /
4953 ! data log10(1),log10(2) / o037632, o020233 /, sc/987/
4954 !
4955 ! machine constants for the sequent balance 8000.
4956 !
4957 ! data small(1) / $00800000 /
4958 ! data large(1) / $7f7fffff /
4959 ! data right(1) / $33800000 /
4960 ! data diver(1) / $34000000 /
4961 ! data log10(1) / $3e9a209b /, sc/987/
4962 !
4963 ! machine constants for the univac 1100 series.
4964 !
4965 ! data rmach(1) / o000400000000 /
4966 ! data rmach(2) / o377777777777 /
4967 ! data rmach(3) / o146400000000 /
4968 ! data rmach(4) / o147400000000 /
4969 ! data rmach(5) / o177464202324 /, sc/987/
4970 !
4971 ! machine constants for the vax unix f77 compiler.
4972 !
4973 ! data small(1) / 128 /
4974 ! data large(1) / -32769 /
4975 ! data right(1) / 13440 /
4976 ! data diver(1) / 13568 /
4977 ! data log10(1) / 547045274 /, sc/987/
4978 !
4979 ! machine constants for the vax-11 with
4980 ! fortran iv-plus compiler.
4981 !
4982 ! data rmach(1) / z00000080 /
4983 ! data rmach(2) / zffff7fff /
4984 ! data rmach(3) / z00003480 /
4985 ! data rmach(4) / z00003500 /
4986 ! data rmach(5) / z209b3f9a /, sc/987/
4987 !
4988 ! machine constants for vax/vms version 2.2.
4989 !
4990 ! data rmach(1) / '80'x /
4991 ! data rmach(2) / 'ffff7fff'x /
4992 ! data rmach(3) / '3480'x /
4993 ! data rmach(4) / '3500'x /
4994 ! data rmach(5) / '209b3f9a'x /, sc/987/
4995 !
4996 ! *** issue stop 778 if all data statements are commented...
4997 ! if (sc .ne. 987) stop 778
4998 if (sc .ne. 987) then
4999 call peg_error_fatal( -1, &
5000 '*** func r1mach fatal error -- all data statements inactive' )
5001 stop
5002 end if
5003
5004 if (i .lt. 1 .or. i .gt. 5) goto 999
5005 r1mach = rmach(i)
5006 return
5007
5008 ! 999 write(*,1999) i
5009 !1999 format(' r1mach - i out of bounds',i10)
5010 999 write(errmsg,1999) i
5011 1999 format('*** func r1mach fatal error -- i out of bounds',i10)
5012 call peg_error_fatal( -1, errmsg )
5013 stop
5014 end function r1mach
5015 !
5016 ! subroutine xsetf
5017
5018 subroutine xsetf (mflag)
5019 !
5020 ! this routine resets the print control flag mflag.
5021 !
5022 integer mflag, mesflg, lunit
5023 common /eh0001/ mesflg, lunit
5024 !
5025 if (mflag .eq. 0 .or. mflag .eq. 1) mesflg = mflag
5026 return
5027 !----------------------- end of subroutine xsetf -----------------------
5028 end subroutine xsetf
5029
5030
5031 !-----------------------------------------------------------------------
5032 subroutine set_lsodes_common_vars()
5033 !
5034 ! place various constant or initial values into lsodes common blocks
5035 !
5036 common /eh0001/ mesflg, lunit
5037 common /ls0001/ rowns(209), &
5038 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, &
5039 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm, &
5040 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6), &
5041 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
5042 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5043
5044 ! lsodes parameters
5045 illin = 0
5046 ntrep = 0
5047 mesflg = 1
5048 lunit = 6
5049
5050 return
5051 !--------------- end of subroutine set_lsodes_common_vars ---------------
5052 end subroutine set_lsodes_common_vars
5053
5054
5055 end module module_cbmz_lsodes_solver
5056
5057
5058 !----------------------------------------------------------------------
5059 ! Subr stode and prep must be outside of the module definition.
5060 ! When lsodes calls stode, the rwork array (in lsodes) is passed to
5061 ! both the wm and iwm arrays (in stode). This is treated as a
5062 ! severe error if stode is within the module.
5063 ! The same problem arises when iprep calls prep.
5064 ! They were renamed to stode_lsodes and prep_lsodes to reduce the
5065 ! chance of name conflicts.
5066 !
5067 subroutine stode_lsodes (neq, y, yh, nyh, yh1, ewt, savf, acor, &
5068 wm, iwm, f, jac, pjac, slvs, &
5069 ruserpar, nruserpar, iuserpar, niuserpar )
5070 use module_cbmz_lsodes_solver, only: cfode, prjs, slss, r1mach, vnorm
5071 !lll. optimize
5072 external f, jac, pjac, slvs
5073 integer neq, nyh, iwm
5074 integer iownd, ialth, ipup, lmax, meo, nqnyh, nslp, &
5075 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
5076 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5077 integer i, i1, iredo, iret, j, jb, m, ncf, newq
5078 integer nruserpar, iuserpar, niuserpar
5079 real y, yh, yh1, ewt, savf, acor, wm
5080 real conit, crate, el, elco, hold, rmax, tesco, &
5081 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
5082 !rce real dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup, &
5083 !rce r, rh, rhdn, rhsm, rhup, told, vnorm
5084 real dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup, &
5085 r, rh, rhdn, rhsm, rhup, told
5086 real ruserpar
5087 !jdf dimension neq(1), y(1), yh(nyh,1), yh1(1), ewt(1), savf(1),
5088 !jdf 1 acor(1), wm(1), iwm(1)
5089 dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*), &
5090 acor(*), wm(*), iwm(*)
5091 dimension ruserpar(nruserpar), iuserpar(niuserpar)
5092 common /ls0001/ conit, crate, el(13), elco(13,12), &
5093 hold, rmax, tesco(3,12), &
5094 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, iownd(14), &
5095 ialth, ipup, lmax, meo, nqnyh, nslp, &
5096 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
5097 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5098 !-----------------------------------------------------------------------
5099 ! stode performs one step of the integration of an initial value
5100 ! problem for a system of ordinary differential equations.
5101 ! note.. stode is independent of the value of the iteration method
5102 ! indicator miter, when this is .ne. 0, and hence is independent
5103 ! of the type of chord method used, or the jacobian structure.
5104 ! communication with stode is done with the following variables..
5105 !
5106 ! neq = integer array containing problem size in neq(1), and
5107 ! passed as the neq argument in all calls to f and jac.
5108 ! y = an array of length .ge. n used as the y argument in
5109 ! all calls to f and jac.
5110 ! yh = an nyh by lmax array containing the dependent variables
5111 ! and their approximate scaled derivatives, where
5112 ! lmax = maxord + 1. yh(i,j+1) contains the approximate
5113 ! j-th derivative of y(i), scaled by h**j/factorial(j)
5114 ! (j = 0,1,...,nq). on entry for the first step, the first
5115 ! two columns of yh must be set from the initial values.
5116 ! nyh = a constant integer .ge. n, the first dimension of yh.
5117 ! yh1 = a one-dimensional array occupying the same space as yh.
5118 ! ewt = an array of length n containing multiplicative weights
5119 ! for local error measurements. local errors in y(i) are
5120 ! compared to 1.0/ewt(i) in various error tests.
5121 ! savf = an array of working storage, of length n.
5122 ! also used for input of yh(*,maxord+2) when jstart = -1
5123 ! and maxord .lt. the current order nq.
5124 ! acor = a work array of length n, used for the accumulated
5125 ! corrections. on a successful return, acor(i) contains
5126 ! the estimated one-step local error in y(i).
5127 ! wm,iwm = real and integer work arrays associated with matrix
5128 ! operations in chord iteration (miter .ne. 0).
5129 ! pjac = name of routine to evaluate and preprocess jacobian matrix
5130 ! and p = i - h*el0*jac, if a chord method is being used.
5131 ! slvs = name of routine to solve linear system in chord iteration.
5132 ! ccmax = maximum relative change in h*el0 before pjac is called.
5133 ! h = the step size to be attempted on the next step.
5134 ! h is altered by the error control algorithm during the
5135 ! problem. h can be either positive or negative, but its
5136 ! sign must remain constant throughout the problem.
5137 ! hmin = the minimum absolute value of the step size h to be used.
5138 ! hmxi = inverse of the maximum absolute value of h to be used.
5139 ! hmxi = 0.0 is allowed and corresponds to an infinite hmax.
5140 ! hmin and hmxi may be changed at any time, but will not
5141 ! take effect until the next change of h is considered.
5142 ! tn = the independent variable. tn is updated on each step taken.
5143 ! jstart = an integer used for input only, with the following
5144 ! values and meanings..
5145 ! 0 perform the first step.
5146 ! .gt.0 take a new step continuing from the last.
5147 ! -1 take the next step with a new value of h, maxord,
5148 ! n, meth, miter, and/or matrix parameters.
5149 ! -2 take the next step with a new value of h,
5150 ! but with other inputs unchanged.
5151 ! on return, jstart is set to 1 to facilitate continuation.
5152 ! kflag = a completion code with the following meanings..
5153 ! 0 the step was succesful.
5154 ! -1 the requested error could not be achieved.
5155 ! -2 corrector convergence could not be achieved.
5156 ! -3 fatal error in pjac or slvs.
5157 ! a return with kflag = -1 or -2 means either
5158 ! abs(h) = hmin or 10 consecutive failures occurred.
5159 ! on a return with kflag negative, the values of tn and
5160 ! the yh array are as of the beginning of the last
5161 ! step, and h is the last step size attempted.
5162 ! maxord = the maximum order of integration method to be allowed.
5163 ! maxcor = the maximum number of corrector iterations allowed.
5164 ! msbp = maximum number of steps between pjac calls (miter .gt. 0).
5165 ! mxncf = maximum number of convergence failures allowed.
5166 ! meth/miter = the method flags. see description in driver.
5167 ! n = the number of first-order differential equations.
5168 !-----------------------------------------------------------------------
5169 kflag = 0
5170 told = tn
5171 ncf = 0
5172 ierpj = 0
5173 iersl = 0
5174 jcur = 0
5175 icf = 0
5176 delp = 0.0e0
5177 if (jstart .gt. 0) go to 200
5178 if (jstart .eq. -1) go to 100
5179 if (jstart .eq. -2) go to 160
5180 !-----------------------------------------------------------------------
5181 ! on the first call, the order is set to 1, and other variables are
5182 ! initialized. rmax is the maximum ratio by which h can be increased
5183 ! in a single step. it is initially 1.e4 to compensate for the small
5184 ! initial h, but then is normally equal to 10. if a failure
5185 ! occurs (in corrector convergence or error test), rmax is set at 2
5186 ! for the next increase.
5187 !-----------------------------------------------------------------------
5188 lmax = maxord + 1
5189 nq = 1
5190 l = 2
5191 ialth = 2
5192 rmax = 10000.0e0
5193 rc = 0.0e0
5194 el0 = 1.0e0
5195 crate = 0.7e0
5196 hold = h
5197 meo = meth
5198 nslp = 0
5199 ipup = miter
5200 iret = 3
5201 go to 140
5202 !-----------------------------------------------------------------------
5203 ! the following block handles preliminaries needed when jstart = -1.
5204 ! ipup is set to miter to force a matrix update.
5205 ! if an order increase is about to be considered (ialth = 1),
5206 ! ialth is reset to 2 to postpone consideration one more step.
5207 ! if the caller has changed meth, cfode is called to reset
5208 ! the coefficients of the method.
5209 ! if the caller has changed maxord to a value less than the current
5210 ! order nq, nq is reduced to maxord, and a new h chosen accordingly.
5211 ! if h is to be changed, yh must be rescaled.
5212 ! if h or meth is being changed, ialth is reset to l = nq + 1
5213 ! to prevent further changes in h for that many steps.
5214 !-----------------------------------------------------------------------
5215 100 ipup = miter
5216 lmax = maxord + 1
5217 if (ialth .eq. 1) ialth = 2
5218 if (meth .eq. meo) go to 110
5219 call cfode (meth, elco, tesco)
5220 meo = meth
5221 if (nq .gt. maxord) go to 120
5222 ialth = l
5223 iret = 1
5224 go to 150
5225 110 if (nq .le. maxord) go to 160
5226 120 nq = maxord
5227 l = lmax
5228 do 125 i = 1,l
5229 125 el(i) = elco(i,nq)
5230 nqnyh = nq*nyh
5231 rc = rc*el(1)/el0
5232 el0 = el(1)
5233 conit = 0.5e0/float(nq+2)
5234 ddn = vnorm (n, savf, ewt)/tesco(1,l)
5235 exdn = 1.0e0/float(l)
5236 rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
5237 rh = amin1(rhdn,1.0e0)
5238 iredo = 3
5239 if (h .eq. hold) go to 170
5240 rh = amin1(rh,abs(h/hold))
5241 h = hold
5242 go to 175
5243 !-----------------------------------------------------------------------
5244 ! cfode is called to get all the integration coefficients for the
5245 ! current meth. then the el vector and related constants are reset
5246 ! whenever the order nq is changed, or at the start of the problem.
5247 !-----------------------------------------------------------------------
5248 140 call cfode (meth, elco, tesco)
5249 150 do 155 i = 1,l
5250 155 el(i) = elco(i,nq)
5251 nqnyh = nq*nyh
5252 rc = rc*el(1)/el0
5253 el0 = el(1)
5254 conit = 0.5e0/float(nq+2)
5255 go to (160, 170, 200), iret
5256 !-----------------------------------------------------------------------
5257 ! if h is being changed, the h ratio rh is checked against
5258 ! rmax, hmin, and hmxi, and the yh array rescaled. ialth is set to
5259 ! l = nq + 1 to prevent a change of h for that many steps, unless
5260 ! forced by a convergence or error test failure.
5261 !-----------------------------------------------------------------------
5262 160 if (h .eq. hold) go to 200
5263 rh = h/hold
5264 h = hold
5265 iredo = 3
5266 go to 175
5267 170 rh = amax1(rh,hmin/abs(h))
5268 175 rh = amin1(rh,rmax)
5269 rh = rh/amax1(1.0e0,abs(h)*hmxi*rh)
5270 r = 1.0e0
5271 do 180 j = 2,l
5272 r = r*rh
5273 do 180 i = 1,n
5274 180 yh(i,j) = yh(i,j)*r
5275 h = h*rh
5276 rc = rc*rh
5277 ialth = l
5278 if (iredo .eq. 0) go to 690
5279 !-----------------------------------------------------------------------
5280 ! this section computes the predicted values by effectively
5281 ! multiplying the yh array by the pascal triangle matrix.
5282 ! rc is the ratio of new to old values of the coefficient h*el(1).
5283 ! when rc differs from 1 by more than ccmax, ipup is set to miter
5284 ! to force pjac to be called, if a jacobian is involved.
5285 ! in any case, pjac is called at least every msbp steps.
5286 !-----------------------------------------------------------------------
5287 200 if (abs(rc-1.0e0) .gt. ccmax) ipup = miter
5288 if (nst .ge. nslp+msbp) ipup = miter
5289 tn = tn + h
5290 i1 = nqnyh + 1
5291 do 215 jb = 1,nq
5292 i1 = i1 - nyh
5293 !dir$ ivdep
5294 do 210 i = i1,nqnyh
5295 210 yh1(i) = yh1(i) + yh1(i+nyh)
5296 215 continue
5297 !-----------------------------------------------------------------------
5298 ! up to maxcor corrector iterations are taken. a convergence test is
5299 ! made on the r.m.s. norm of each correction, weighted by the error
5300 ! weight vector ewt. the sum of the corrections is accumulated in the
5301 ! vector acor(i). the yh array is not altered in the corrector loop.
5302 !-----------------------------------------------------------------------
5303 220 m = 0
5304 do 230 i = 1,n
5305 230 y(i) = yh(i,1)
5306 call f (neq, tn, y, savf, &
5307 ruserpar, nruserpar, iuserpar, niuserpar)
5308 nfe = nfe + 1
5309 if (ipup .le. 0) go to 250
5310 !-----------------------------------------------------------------------
5311 ! if indicated, the matrix p = i - h*el(1)*j is reevaluated and
5312 ! preprocessed before starting the corrector iteration. ipup is set
5313 ! to 0 as an indicator that this has been done.
5314 !-----------------------------------------------------------------------
5315 call pjac (neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac, &
5316 ruserpar, nruserpar, iuserpar, niuserpar )
5317 ipup = 0
5318 rc = 1.0e0
5319 nslp = nst
5320 crate = 0.7e0
5321 if (ierpj .ne. 0) go to 430
5322 250 do 260 i = 1,n
5323 260 acor(i) = 0.0e0
5324 270 if (miter .ne. 0) go to 350
5325 !-----------------------------------------------------------------------
5326 ! in the case of functional iteration, update y directly from
5327 ! the result of the last function evaluation.
5328 !-----------------------------------------------------------------------
5329 do 290 i = 1,n
5330 savf(i) = h*savf(i) - yh(i,2)
5331 290 y(i) = savf(i) - acor(i)
5332 del = vnorm (n, y, ewt)
5333 do 300 i = 1,n
5334 y(i) = yh(i,1) + el(1)*savf(i)
5335 300 acor(i) = savf(i)
5336 go to 400
5337 !-----------------------------------------------------------------------
5338 ! in the case of the chord method, compute the corrector error,
5339 ! and solve the linear system with that as right-hand side and
5340 ! p as coefficient matrix.
5341 !-----------------------------------------------------------------------
5342 350 do 360 i = 1,n
5343 360 y(i) = h*savf(i) - (yh(i,2) + acor(i))
5344 call slvs (wm, iwm, y, savf)
5345 if (iersl .lt. 0) go to 430
5346 if (iersl .gt. 0) go to 410
5347 del = vnorm (n, y, ewt)
5348 do 380 i = 1,n
5349 acor(i) = acor(i) + y(i)
5350 380 y(i) = yh(i,1) + el(1)*acor(i)
5351 !-----------------------------------------------------------------------
5352 ! test for convergence. if m.gt.0, an estimate of the convergence
5353 ! rate constant is stored in crate, and this is used in the test.
5354 !-----------------------------------------------------------------------
5355 400 if (m .ne. 0) crate = amax1(0.2e0*crate,del/delp)
5356 dcon = del*amin1(1.0e0,1.5e0*crate)/(tesco(2,nq)*conit)
5357 if (dcon .le. 1.0e0) go to 450
5358 m = m + 1
5359 if (m .eq. maxcor) go to 410
5360 if (m .ge. 2 .and. del .gt. 2.0e0*delp) go to 410
5361 delp = del
5362 call f (neq, tn, y, savf, &
5363 ruserpar, nruserpar, iuserpar, niuserpar)
5364 nfe = nfe + 1
5365 go to 270
5366 !-----------------------------------------------------------------------
5367 ! the corrector iteration failed to converge.
5368 ! if miter .ne. 0 and the jacobian is out of date, pjac is called for
5369 ! the next try. otherwise the yh array is retracted to its values
5370 ! before prediction, and h is reduced, if possible. if h cannot be
5371 ! reduced or mxncf failures have occurred, exit with kflag = -2.
5372 !-----------------------------------------------------------------------
5373 410 if (miter .eq. 0 .or. jcur .eq. 1) go to 430
5374 icf = 1
5375 ipup = miter
5376 go to 220
5377 430 icf = 2
5378 ncf = ncf + 1
5379 rmax = 2.0e0
5380 tn = told
5381 i1 = nqnyh + 1
5382 do 445 jb = 1,nq
5383 i1 = i1 - nyh
5384 !dir$ ivdep
5385 do 440 i = i1,nqnyh
5386 440 yh1(i) = yh1(i) - yh1(i+nyh)
5387 445 continue
5388 if (ierpj .lt. 0 .or. iersl .lt. 0) go to 680
5389 if (abs(h) .le. hmin*1.00001e0) go to 670
5390 if (ncf .eq. mxncf) go to 670
5391 rh = 0.25e0
5392 ipup = miter
5393 iredo = 1
5394 go to 170
5395 !-----------------------------------------------------------------------
5396 ! the corrector has converged. jcur is set to 0
5397 ! to signal that the jacobian involved may need updating later.
5398 ! the local error test is made and control passes to statement 500
5399 ! if it fails.
5400 !-----------------------------------------------------------------------
5401 450 jcur = 0
5402 if (m .eq. 0) dsm = del/tesco(2,nq)
5403 if (m .gt. 0) dsm = vnorm (n, acor, ewt)/tesco(2,nq)
5404 if (dsm .gt. 1.0e0) go to 500
5405 !-----------------------------------------------------------------------
5406 ! after a successful step, update the yh array.
5407 ! consider changing h if ialth = 1. otherwise decrease ialth by 1.
5408 ! if ialth is then 1 and nq .lt. maxord, then acor is saved for
5409 ! use in a possible order increase on the next step.
5410 ! if a change in h is considered, an increase or decrease in order
5411 ! by one is considered also. a change in h is made only if it is by a
5412 ! factor of at least 1.1. if not, ialth is set to 3 to prevent
5413 ! testing for that many steps.
5414 !-----------------------------------------------------------------------
5415 kflag = 0
5416 iredo = 0
5417 nst = nst + 1
5418 hu = h
5419 nqu = nq
5420 do 470 j = 1,l
5421 do 470 i = 1,n
5422 470 yh(i,j) = yh(i,j) + el(j)*acor(i)
5423 ialth = ialth - 1
5424 if (ialth .eq. 0) go to 520
5425 if (ialth .gt. 1) go to 700
5426 if (l .eq. lmax) go to 700
5427 do 490 i = 1,n
5428 490 yh(i,lmax) = acor(i)
5429 go to 700
5430 !-----------------------------------------------------------------------
5431 ! the error test failed. kflag keeps track of multiple failures.
5432 ! restore tn and the yh array to their previous values, and prepare
5433 ! to try the step again. compute the optimum step size for this or
5434 ! one lower order. after 2 or more failures, h is forced to decrease
5435 ! by a factor of 0.2 or less.
5436 !-----------------------------------------------------------------------
5437 500 kflag = kflag - 1
5438 tn = told
5439 i1 = nqnyh + 1
5440 do 515 jb = 1,nq
5441 i1 = i1 - nyh
5442 !dir$ ivdep
5443 do 510 i = i1,nqnyh
5444 510 yh1(i) = yh1(i) - yh1(i+nyh)
5445 515 continue
5446 rmax = 2.0e0
5447 if (abs(h) .le. hmin*1.00001e0) go to 660
5448 if (kflag .le. -3) go to 640
5449 iredo = 2
5450 rhup = 0.0e0
5451 go to 540
5452 !-----------------------------------------------------------------------
5453 ! regardless of the success or failure of the step, factors
5454 ! rhdn, rhsm, and rhup are computed, by which h could be multiplied
5455 ! at order nq - 1, order nq, or order nq + 1, respectively.
5456 ! in the case of failure, rhup = 0.0 to avoid an order increase.
5457 ! the largest of these is determined and the new order chosen
5458 ! accordingly. if the order is to be increased, we compute one
5459 ! additional scaled derivative.
5460 !-----------------------------------------------------------------------
5461 520 rhup = 0.0e0
5462 if (l .eq. lmax) go to 540
5463 do 530 i = 1,n
5464 530 savf(i) = acor(i) - yh(i,lmax)
5465 dup = vnorm (n, savf, ewt)/tesco(3,nq)
5466 exup = 1.0e0/float(l+1)
5467 rhup = 1.0e0/(1.4e0*dup**exup + 0.0000014e0)
5468 540 exsm = 1.0e0/float(l)
5469 rhsm = 1.0e0/(1.2e0*dsm**exsm + 0.0000012e0)
5470 rhdn = 0.0e0
5471 if (nq .eq. 1) go to 560
5472 ddn = vnorm (n, yh(1,l), ewt)/tesco(1,nq)
5473 exdn = 1.0e0/float(nq)
5474 rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
5475 560 if (rhsm .ge. rhup) go to 570
5476 if (rhup .gt. rhdn) go to 590
5477 go to 580
5478 570 if (rhsm .lt. rhdn) go to 580
5479 newq = nq
5480 rh = rhsm
5481 go to 620
5482 580 newq = nq - 1
5483 rh = rhdn
5484 if (kflag .lt. 0 .and. rh .gt. 1.0e0) rh = 1.0e0
5485 go to 620
5486 590 newq = l
5487 rh = rhup
5488 if (rh .lt. 1.1e0) go to 610
5489 r = el(l)/float(l)
5490 do 600 i = 1,n
5491 600 yh(i,newq+1) = acor(i)*r
5492 go to 630
5493 610 ialth = 3
5494 go to 700
5495 620 if ((kflag .eq. 0) .and. (rh .lt. 1.1e0)) go to 610
5496 if (kflag .le. -2) rh = amin1(rh,0.2e0)
5497 !-----------------------------------------------------------------------
5498 ! if there is a change of order, reset nq, l, and the coefficients.
5499 ! in any case h is reset according to rh and the yh array is rescaled.
5500 ! then exit from 690 if the step was ok, or redo the step otherwise.
5501 !-----------------------------------------------------------------------
5502 if (newq .eq. nq) go to 170
5503 630 nq = newq
5504 l = nq + 1
5505 iret = 2
5506 go to 150
5507 !-----------------------------------------------------------------------
5508 ! control reaches this section if 3 or more failures have occured.
5509 ! if 10 failures have occurred, exit with kflag = -1.
5510 ! it is assumed that the derivatives that have accumulated in the
5511 ! yh array have errors of the wrong order. hence the first
5512 ! derivative is recomputed, and the order is set to 1. then
5513 ! h is reduced by a factor of 10, and the step is retried,
5514 ! until it succeeds or h reaches hmin.
5515 !-----------------------------------------------------------------------
5516 640 if (kflag .eq. -10) go to 660
5517 rh = 0.1e0
5518 rh = amax1(hmin/abs(h),rh)
5519 h = h*rh
5520 do 645 i = 1,n
5521 645 y(i) = yh(i,1)
5522 call f (neq, tn, y, savf, &
5523 ruserpar, nruserpar, iuserpar, niuserpar)
5524 nfe = nfe + 1
5525 do 650 i = 1,n
5526 650 yh(i,2) = h*savf(i)
5527 ipup = miter
5528 ialth = 5
5529 if (nq .eq. 1) go to 200
5530 nq = 1
5531 l = 2
5532 iret = 3
5533 go to 150
5534 !-----------------------------------------------------------------------
5535 ! all returns are made through this section. h is saved in hold
5536 ! to allow the caller to change h on the next step.
5537 !-----------------------------------------------------------------------
5538 660 kflag = -1
5539 go to 720
5540 670 kflag = -2
5541 go to 720
5542 680 kflag = -3
5543 go to 720
5544 690 rmax = 10.0e0
5545 700 r = 1.0e0/tesco(2,nqu)
5546 do 710 i = 1,n
5547 710 acor(i) = acor(i)*r
5548 720 hold = h
5549 jstart = 1
5550 return
5551 !----------------------- end of subroutine stode_lsodes -----------------------
5552 end subroutine stode_lsodes
5553
5554
5555
5556 subroutine prep_lsodes (neq, y, yh, savf, ewt, ftem, ia, ja, &
5557 wk, iwk, ipper, f, jac, &
5558 ruserpar, nruserpar, iuserpar, niuserpar )
5559 use module_cbmz_lsodes_solver, only: adjlr, cdrv, cntnzu, jgroup, &
5560 odrv
5561 !lll. optimize
5562 external f,jac
5563 integer neq, ia, ja, iwk, ipper
5564 integer iownd, iowns, &
5565 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
5566 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5567 integer iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp, &
5568 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa, &
5569 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj, &
5570 nslj, ngp, nlu, nnz, nsp, nzl, nzu
5571 integer i, ibr, ier, ipil, ipiu, iptt1, iptt2, j, jfound, k, &
5572 knew, kmax, kmin, ldif, lenigp, liwk, maxg, np1, nzsut
5573 integer nruserpar, iuserpar, niuserpar
5574 real y, yh, savf, ewt, ftem, wk
5575 real rowns, &
5576 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
5577 real con0, conmin, ccmxj, psmall, rbig, seth
5578 real dq, dyj, erwt, fac, yj
5579 real ruserpar
5580 !jdf dimension neq(1), y(1), yh(1), savf(1), ewt(1), ftem(1),
5581 !jdf 1 ia(1), ja(1), wk(1), iwk(1)
5582 dimension neq(*), y(*), yh(*), savf(*), ewt(*), ftem(*), &
5583 ia(*), ja(*), wk(*), iwk(*)
5584 dimension ruserpar(nruserpar), iuserpar(niuserpar)
5585 common /ls0001/ rowns(209), &
5586 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, &
5587 iownd(14), iowns(6), &
5588 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, &
5589 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
5590 common /lss001/ con0, conmin, ccmxj, psmall, rbig, seth, &
5591 iplost, iesp, istatc, iys, iba, ibian, ibjan, ibjgp, &
5592 ipian, ipjan, ipjgp, ipigp, ipr, ipc, ipic, ipisp, iprsp, ipa, &
5593 lenyh, lenyhm, lenwk, lreq, lrat, lrest, lwmin, moss, msbj, &
5594 nslj, ngp, nlu, nnz, nsp, nzl, nzu
5595 !-----------------------------------------------------------------------
5596 ! this routine performs preprocessing related to the sparse linear
5597 ! systems that must be solved if miter = 1 or 2.
5598 ! the operations that are performed here are..
5599 ! * compute sparseness structure of jacobian according to moss,
5600 ! * compute grouping of column indices (miter = 2),
5601 ! * compute a new ordering of rows and columns of the matrix,
5602 ! * reorder ja corresponding to the new ordering,
5603 ! * perform a symbolic lu factorization of the matrix, and
5604 ! * set pointers for segments of the iwk/wk array.
5605 ! in addition to variables described previously, prep uses the
5606 ! following for communication..
5607 ! yh = the history array. only the first column, containing the
5608 ! current y vector, is used. used only if moss .ne. 0.
5609 ! savf = a work array of length neq, used only if moss .ne. 0.
5610 ! ewt = array of length neq containing (inverted) error weights.
5611 ! used only if moss = 2 or if istate = moss = 1.
5612 ! ftem = a work array of length neq, identical to acor in the driver,
5613 ! used only if moss = 2.
5614 ! wk = a real work array of length lenwk, identical to wm in
5615 ! the driver.
5616 ! iwk = integer work array, assumed to occupy the same space as wk.
5617 ! lenwk = the length of the work arrays wk and iwk.
5618 ! istatc = a copy of the driver input argument istate (= 1 on the
5619 ! first call, = 3 on a continuation call).
5620 ! iys = flag value from odrv or cdrv.
5621 ! ipper = output error flag with the following values and meanings..
5622 ! 0 no error.
5623 ! -1 insufficient storage for internal structure pointers.
5624 ! -2 insufficient storage for jgroup.
5625 ! -3 insufficient storage for odrv.
5626 ! -4 other error flag from odrv (should never occur).
5627 ! -5 insufficient storage for cdrv.
5628 ! -6 other error flag from cdrv.
5629 !-----------------------------------------------------------------------
5630 ibian = lrat*2
5631 ipian = ibian + 1
5632 np1 = n + 1
5633 ipjan = ipian + np1
5634 ibjan = ipjan - 1
5635 liwk = lenwk*lrat
5636 if (ipjan+n-1 .gt. liwk) go to 210
5637 if (moss .eq. 0) go to 30
5638 !
5639 if (istatc .eq. 3) go to 20
5640 ! istate = 1 and moss .ne. 0. perturb y for structure determination. --
5641 do 10 i = 1,n
5642 erwt = 1.0e0/ewt(i)
5643 fac = 1.0e0 + 1.0e0/(float(i)+1.0e0)
5644 y(i) = y(i) + fac*sign(erwt,y(i))
5645 10 continue
5646 go to (70, 100), moss
5647 !
5648 20 continue
5649 ! istate = 3 and moss .ne. 0. load y from yh(*,1). --------------------
5650 do 25 i = 1,n
5651 25 y(i) = yh(i)
5652 go to (70, 100), moss
5653 !
5654 ! moss = 0. process user-s ia,ja. add diagonal entries if necessary. -
5655 30 knew = ipjan
5656 kmin = ia(1)
5657 iwk(ipian) = 1
5658 do 60 j = 1,n
5659 jfound = 0
5660 kmax = ia(j+1) - 1
5661 if (kmin .gt. kmax) go to 45
5662 do 40 k = kmin,kmax
5663 i = ja(k)
5664 if (i .eq. j) jfound = 1
5665 if (knew .gt. liwk) go to 210
5666 iwk(knew) = i
5667 knew = knew + 1
5668 40 continue
5669 if (jfound .eq. 1) go to 50
5670 45 if (knew .gt. liwk) go to 210
5671 iwk(knew) = j
5672 knew = knew + 1
5673 50 iwk(ipian+j) = knew + 1 - ipjan
5674 kmin = kmax + 1
5675 60 continue
5676 go to 140
5677 !
5678 ! moss = 1. compute structure from user-supplied jacobian routine jac.
5679 70 continue
5680 ! a dummy call to f allows user to create temporaries for use in jac. --
5681 call f (neq, tn, y, savf, &
5682 ruserpar, nruserpar, iuserpar, niuserpar)
5683 k = ipjan
5684 iwk(ipian) = 1
5685 do 90 j = 1,n
5686 if (k .gt. liwk) go to 210
5687 iwk(k) = j
5688 k = k + 1
5689 do 75 i = 1,n
5690 75 savf(i) = 0.0e0
5691 call jac (neq, tn, y, j, iwk(ipian), iwk(ipjan), savf, &
5692 ruserpar, nruserpar, iuserpar, niuserpar)
5693 do 80 i = 1,n
5694 if (abs(savf(i)) .le. seth) go to 80
5695 if (i .eq. j) go to 80
5696 if (k .gt. liwk) go to 210
5697 iwk(k) = i
5698 k = k + 1
5699 80 continue
5700 iwk(ipian+j) = k + 1 - ipjan
5701 90 continue
5702 go to 140
5703 !
5704 ! moss = 2. compute structure from results of n + 1 calls to f. -------
5705 100 k = ipjan
5706 iwk(ipian) = 1
5707 call f (neq, tn, y, savf, &
5708 ruserpar, nruserpar, iuserpar, niuserpar)
5709 do 120 j = 1,n
5710 if (k .gt. liwk) go to 210
5711 iwk(k) = j
5712 k = k + 1
5713 yj = y(j)
5714 erwt = 1.0e0/ewt(j)
5715 dyj = sign(erwt,yj)
5716 y(j) = yj + dyj
5717 call f (neq, tn, y, ftem, &
5718 ruserpar, nruserpar, iuserpar, niuserpar)
5719 y(j) = yj
5720 do 110 i = 1,n
5721 dq = (ftem(i) - savf(i))/dyj
5722 if (abs(dq) .le. seth) go to 110
5723 if (i .eq. j) go to 110
5724 if (k .gt. liwk) go to 210
5725 iwk(k) = i
5726 k = k + 1
5727 110 continue
5728 iwk(ipian+j) = k + 1 - ipjan
5729 120 continue
5730 !
5731 140 continue
5732 if (moss .eq. 0 .or. istatc .ne. 1) go to 150
5733 ! if istate = 1 and moss .ne. 0, restore y from yh. --------------------
5734 do 145 i = 1,n
5735 145 y(i) = yh(i)
5736 150 nnz = iwk(ipian+n) - 1
5737 lenigp = 0
5738 ipigp = ipjan + nnz
5739 if (miter .ne. 2) go to 160
5740 !
5741 ! compute grouping of column indices (miter = 2). ----------------------
5742 maxg = np1
5743 ipjgp = ipjan + nnz
5744 ibjgp = ipjgp - 1
5745 ipigp = ipjgp + n
5746 iptt1 = ipigp + np1
5747 iptt2 = iptt1 + n
5748 lreq = iptt2 + n - 1
5749 if (lreq .gt. liwk) go to 220
5750 call jgroup (n, iwk(ipian), iwk(ipjan), maxg, ngp, iwk(ipigp), &
5751 iwk(ipjgp), iwk(iptt1), iwk(iptt2), ier)
5752 if (ier .ne. 0) go to 220
5753 lenigp = ngp + 1
5754 !
5755 ! compute new ordering of rows/columns of jacobian. --------------------
5756 160 ipr = ipigp + lenigp
5757 ipc = ipr
5758 ipic = ipc + n
5759 ipisp = ipic + n
5760 iprsp = (ipisp - 2)/lrat + 2
5761 iesp = lenwk + 1 - iprsp
5762 if (iesp .lt. 0) go to 230
5763 ibr = ipr - 1
5764 do 170 i = 1,n
5765 170 iwk(ibr+i) = i
5766 nsp = liwk + 1 - ipisp
5767 call odrv (n, iwk(ipian), iwk(ipjan), wk, iwk(ipr), iwk(ipic), &
5768 nsp, iwk(ipisp), 1, iys)
5769 if (iys .eq. 11*n+1) go to 240
5770 if (iys .ne. 0) go to 230
5771 !
5772 ! reorder jan and do symbolic lu factorization of matrix. --------------
5773 ipa = lenwk + 1 - nnz
5774 nsp = ipa - iprsp
5775 lreq = max0(12*n/lrat, 6*n/lrat+2*n+nnz) + 3
5776 lreq = lreq + iprsp - 1 + nnz
5777 if (lreq .gt. lenwk) go to 250
5778 iba = ipa - 1
5779 do 180 i = 1,nnz
5780 180 wk(iba+i) = 0.0e0
5781 ipisp = lrat*(iprsp - 1) + 1
5782 call cdrv (n,iwk(ipr),iwk(ipc),iwk(ipic),iwk(ipian),iwk(ipjan), &
5783 wk(ipa),wk(ipa),wk(ipa),nsp,iwk(ipisp),wk(iprsp),iesp,5,iys)
5784 lreq = lenwk - iesp
5785 if (iys .eq. 10*n+1) go to 250
5786 if (iys .ne. 0) go to 260
5787 ipil = ipisp
5788 ipiu = ipil + 2*n + 1
5789 nzu = iwk(ipil+n) - iwk(ipil)
5790 nzl = iwk(ipiu+n) - iwk(ipiu)
5791 if (lrat .gt. 1) go to 190
5792 call adjlr (n, iwk(ipisp), ldif)
5793 lreq = lreq + ldif
5794 190 continue
5795 if (lrat .eq. 2 .and. nnz .eq. n) lreq = lreq + 1
5796 nsp = nsp + lreq - lenwk
5797 ipa = lreq + 1 - nnz
5798 iba = ipa - 1
5799 ipper = 0
5800 return
5801 !
5802 210 ipper = -1
5803 lreq = 2 + (2*n + 1)/lrat
5804 lreq = max0(lenwk+1,lreq)
5805 return
5806 !
5807 220 ipper = -2
5808 lreq = (lreq - 1)/lrat + 1
5809 return
5810 !
5811 230 ipper = -3
5812 call cntnzu (n, iwk(ipian), iwk(ipjan), nzsut)
5813 lreq = lenwk - iesp + (3*n + 4*nzsut - 1)/lrat + 1
5814 return
5815 !
5816 240 ipper = -4
5817 return
5818 !
5819 250 ipper = -5
5820 return
5821 !
5822 260 ipper = -6
5823 lreq = lenwk
5824 return
5825 !----------------------- end of subroutine prep_lsodes ------------------------
5826 end subroutine prep_lsodes