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