In documentation for lreduce and rreduce, supply second argument as an explicit list
[maxima.git] / share / odepack / src / dstoda.lisp
blob78f7c0883c3c21a57d4599341402af6e1b130bd6
1 ;;; Compiled by f2cl version:
2 ;;; ("f2cl1.l,v 95098eb54f13 2013/04/01 00:45:16 toy $"
3 ;;; "f2cl2.l,v 95098eb54f13 2013/04/01 00:45:16 toy $"
4 ;;; "f2cl3.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
5 ;;; "f2cl4.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
6 ;;; "f2cl5.l,v 95098eb54f13 2013/04/01 00:45:16 toy $"
7 ;;; "f2cl6.l,v 1d5cbacbb977 2008/08/24 00:56:27 rtoy $"
8 ;;; "macros.l,v 1409c1352feb 2013/03/24 20:44:50 toy $")
10 ;;; Using Lisp CMU Common Lisp snapshot-2013-11 (20E Unicode)
11 ;;;
12 ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
13 ;;; (:coerce-assigns :as-needed) (:array-type ':array)
14 ;;; (:array-slicing t) (:declare-common nil)
15 ;;; (:float-format single-float))
17 (in-package "ODEPACK")
20 (let ((sm1
21 (make-array 12
22 :element-type 'double-float
23 :initial-contents '(0.5d0 0.575d0 0.55d0 0.45d0 0.35d0
24 0.25d0 0.2d0 0.15d0 0.1d0 0.075d0 0.05d0
25 0.025d0))))
26 (declare (type (array double-float (12)) sm1))
27 (defun dstoda (neq y yh nyh yh1 ewt savf acor wm iwm f jac pjac slvs)
28 (declare (type (f2cl-lib:integer4) nyh)
29 (type (array double-float (*)) wm acor savf ewt yh1 yh y)
30 (type (array f2cl-lib:integer4 (*)) iwm neq))
31 (let ((dls001-el
32 (make-array 13
33 :element-type 'double-float
34 :displaced-to (dls001-part-0 *dls001-common-block*)
35 :displaced-index-offset 2))
36 (dls001-elco
37 (make-array 156
38 :element-type 'double-float
39 :displaced-to (dls001-part-0 *dls001-common-block*)
40 :displaced-index-offset 15))
41 (dls001-tesco
42 (make-array 36
43 :element-type 'double-float
44 :displaced-to (dls001-part-0 *dls001-common-block*)
45 :displaced-index-offset 173))
46 (dlsa01-cm1
47 (make-array 12
48 :element-type 'double-float
49 :displaced-to (dlsa01-part-0 *dlsa01-common-block*)
50 :displaced-index-offset 1))
51 (dlsa01-cm2
52 (make-array 5
53 :element-type 'double-float
54 :displaced-to (dlsa01-part-0 *dlsa01-common-block*)
55 :displaced-index-offset 13)))
56 (symbol-macrolet ((conit (aref (dls001-part-0 *dls001-common-block*) 0))
57 (crate (aref (dls001-part-0 *dls001-common-block*) 1))
58 (el dls001-el)
59 (elco dls001-elco)
60 (hold (aref (dls001-part-0 *dls001-common-block*) 171))
61 (rmax (aref (dls001-part-0 *dls001-common-block*) 172))
62 (tesco dls001-tesco)
63 (ccmax
64 (aref (dls001-part-0 *dls001-common-block*) 209))
65 (el0 (aref (dls001-part-0 *dls001-common-block*) 210))
66 (h (aref (dls001-part-0 *dls001-common-block*) 211))
67 (hmin (aref (dls001-part-0 *dls001-common-block*) 212))
68 (hmxi (aref (dls001-part-0 *dls001-common-block*) 213))
69 (hu (aref (dls001-part-0 *dls001-common-block*) 214))
70 (rc (aref (dls001-part-0 *dls001-common-block*) 215))
71 (tn (aref (dls001-part-0 *dls001-common-block*) 216))
72 (uround
73 (aref (dls001-part-0 *dls001-common-block*) 217))
74 (ialth (aref (dls001-part-1 *dls001-common-block*) 6))
75 (ipup (aref (dls001-part-1 *dls001-common-block*) 7))
76 (lmax (aref (dls001-part-1 *dls001-common-block*) 8))
77 (nqnyh (aref (dls001-part-1 *dls001-common-block*) 10))
78 (nslp (aref (dls001-part-1 *dls001-common-block*) 11))
79 (icf (aref (dls001-part-1 *dls001-common-block*) 12))
80 (ierpj (aref (dls001-part-1 *dls001-common-block*) 13))
81 (iersl (aref (dls001-part-1 *dls001-common-block*) 14))
82 (jcur (aref (dls001-part-1 *dls001-common-block*) 15))
83 (jstart
84 (aref (dls001-part-1 *dls001-common-block*) 16))
85 (kflag (aref (dls001-part-1 *dls001-common-block*) 17))
86 (l (aref (dls001-part-1 *dls001-common-block*) 18))
87 (meth (aref (dls001-part-1 *dls001-common-block*) 25))
88 (miter (aref (dls001-part-1 *dls001-common-block*) 26))
89 (maxord
90 (aref (dls001-part-1 *dls001-common-block*) 27))
91 (maxcor
92 (aref (dls001-part-1 *dls001-common-block*) 28))
93 (msbp (aref (dls001-part-1 *dls001-common-block*) 29))
94 (mxncf (aref (dls001-part-1 *dls001-common-block*) 30))
95 (n (aref (dls001-part-1 *dls001-common-block*) 31))
96 (nq (aref (dls001-part-1 *dls001-common-block*) 32))
97 (nst (aref (dls001-part-1 *dls001-common-block*) 33))
98 (nfe (aref (dls001-part-1 *dls001-common-block*) 34))
99 (nqu (aref (dls001-part-1 *dls001-common-block*) 36))
100 (cm1 dlsa01-cm1)
101 (cm2 dlsa01-cm2)
102 (pdest (aref (dlsa01-part-0 *dlsa01-common-block*) 18))
103 (pdlast
104 (aref (dlsa01-part-0 *dlsa01-common-block*) 19))
105 (ratio (aref (dlsa01-part-0 *dlsa01-common-block*) 20))
106 (pdnorm
107 (aref (dlsa01-part-0 *dlsa01-common-block*) 21))
108 (icount (aref (dlsa01-part-1 *dlsa01-common-block*) 3))
109 (irflag (aref (dlsa01-part-1 *dlsa01-common-block*) 4))
110 (jtyp (aref (dlsa01-part-1 *dlsa01-common-block*) 5))
111 (mused (aref (dlsa01-part-1 *dlsa01-common-block*) 6))
112 (mxordn (aref (dlsa01-part-1 *dlsa01-common-block*) 7))
113 (mxords (aref (dlsa01-part-1 *dlsa01-common-block*) 8)))
114 (f2cl-lib:with-multi-array-data
115 ((neq f2cl-lib:integer4 neq-%data% neq-%offset%)
116 (iwm f2cl-lib:integer4 iwm-%data% iwm-%offset%)
117 (y double-float y-%data% y-%offset%)
118 (yh double-float yh-%data% yh-%offset%)
119 (yh1 double-float yh1-%data% yh1-%offset%)
120 (ewt double-float ewt-%data% ewt-%offset%)
121 (savf double-float savf-%data% savf-%offset%)
122 (acor double-float acor-%data% acor-%offset%)
123 (wm double-float wm-%data% wm-%offset%))
124 (prog ((newq 0) (ncf 0) (m 0) (jb 0) (j 0) (iret 0) (iredo 0) (i1 0)
125 (i 0) (nqm2 0) (nqm1 0) (lm2p1 0) (lm2 0) (lm1p1 0) (lm1 0)
126 (told 0.0d0) (rhup 0.0d0) (rhsm 0.0d0) (rhdn 0.0d0) (rh 0.0d0)
127 (r 0.0d0) (exup 0.0d0) (exsm 0.0d0) (exdn 0.0d0) (dup 0.0d0)
128 (dsm 0.0d0) (delp 0.0d0) (del 0.0d0) (ddn 0.0d0) (dcon 0.0d0)
129 (rm 0.0d0) (rh2 0.0d0) (rh1it 0.0d0) (rh1 0.0d0) (rate 0.0d0)
130 (pnorm 0.0d0) (pdh 0.0d0) (exm2 0.0d0) (exm1 0.0d0)
131 (dm2 0.0d0) (dm1 0.0d0) (alpha 0.0d0))
132 (declare (type (double-float) alpha dm1 dm2 exm1 exm2 pdh pnorm
133 rate rh1 rh1it rh2 rm dcon ddn del
134 delp dsm dup exdn exsm exup r rh rhdn
135 rhsm rhup told)
136 (type (f2cl-lib:integer4) lm1 lm1p1 lm2 lm2p1 nqm1 nqm2 i
137 i1 iredo iret j jb m ncf newq))
138 (setf kflag 0)
139 (setf told tn)
140 (setf ncf 0)
141 (setf ierpj 0)
142 (setf iersl 0)
143 (setf jcur 0)
144 (setf icf 0)
145 (setf delp 0.0d0)
146 (if (> jstart 0) (go label200))
147 (if (= jstart -1) (go label100))
148 (if (= jstart -2) (go label160))
149 (setf lmax (f2cl-lib:int-add maxord 1))
150 (setf nq 1)
151 (setf l 2)
152 (setf ialth 2)
153 (setf rmax 10000.0d0)
154 (setf rc 0.0d0)
155 (setf el0 1.0d0)
156 (setf crate 0.7d0)
157 (setf hold h)
158 (setf nslp 0)
159 (setf ipup miter)
160 (setf iret 3)
161 (setf icount 20)
162 (setf irflag 0)
163 (setf pdest 0.0d0)
164 (setf pdlast 0.0d0)
165 (setf ratio 5.0d0)
166 (dcfode 2 elco tesco)
167 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
168 ((> i 5) nil)
169 (tagbody
170 label10
171 (setf (f2cl-lib:fref cm2 (i) ((1 5)))
172 (* (f2cl-lib:fref tesco (2 i) ((1 3) (1 12)))
173 (f2cl-lib:fref elco
174 ((f2cl-lib:int-add i 1) i)
175 ((1 13) (1 12)))))))
176 (dcfode 1 elco tesco)
177 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
178 ((> i 12) nil)
179 (tagbody
180 label20
181 (setf (f2cl-lib:fref cm1 (i) ((1 12)))
182 (* (f2cl-lib:fref tesco (2 i) ((1 3) (1 12)))
183 (f2cl-lib:fref elco
184 ((f2cl-lib:int-add i 1) i)
185 ((1 13) (1 12)))))))
186 (go label150)
187 label100
188 (setf ipup miter)
189 (setf lmax (f2cl-lib:int-add maxord 1))
190 (if (= ialth 1) (setf ialth 2))
191 (if (= meth mused) (go label160))
192 (dcfode meth elco tesco)
193 (setf ialth l)
194 (setf iret 1)
195 label150
196 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
197 ((> i l) nil)
198 (tagbody
199 label155
200 (setf (f2cl-lib:fref el (i) ((1 13)))
201 (f2cl-lib:fref elco (i nq) ((1 13) (1 12))))))
202 (setf nqnyh (f2cl-lib:int-mul nq nyh))
203 (setf rc (/ (* rc (f2cl-lib:fref el (1) ((1 13)))) el0))
204 (setf el0 (f2cl-lib:fref el (1) ((1 13))))
205 (setf conit (/ 0.5d0 (f2cl-lib:int-add nq 2)))
206 (f2cl-lib:computed-goto (label160 label170 label200) iret)
207 label160
208 (if (= h hold) (go label200))
209 (setf rh (/ h hold))
210 (setf h hold)
211 (setf iredo 3)
212 (go label175)
213 label170
214 (setf rh (max rh (/ hmin (abs h))))
215 label175
216 (setf rh (min rh rmax))
217 (setf rh (/ rh (max 1.0d0 (* (abs h) hmxi rh))))
218 (if (= meth 2) (go label178))
219 (setf irflag 0)
220 (setf pdh (max (* (abs h) pdlast) 1.0d-6))
221 (if (< (* rh pdh 1.00001d0) (f2cl-lib:fref sm1 (nq) ((1 12))))
222 (go label178))
223 (setf rh (/ (f2cl-lib:fref sm1 (nq) ((1 12))) pdh))
224 (setf irflag 1)
225 label178
226 (setf r 1.0d0)
227 (f2cl-lib:fdo (j 2 (f2cl-lib:int-add j 1))
228 ((> j l) nil)
229 (tagbody
230 (setf r (* r rh))
231 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
232 ((> i n) nil)
233 (tagbody
234 (setf (f2cl-lib:fref yh-%data%
235 (i j)
236 ((1 nyh) (1 *))
237 yh-%offset%)
239 (f2cl-lib:fref yh-%data%
240 (i j)
241 ((1 nyh) (1 *))
242 yh-%offset%)
243 r))))))
244 label180
245 (setf h (* h rh))
246 (setf rc (* rc rh))
247 (setf ialth l)
248 (if (= iredo 0) (go label690))
249 label200
250 (if (> (abs (- rc 1.0d0)) ccmax) (setf ipup miter))
251 (if (>= nst (f2cl-lib:int-add nslp msbp)) (setf ipup miter))
252 (setf tn (+ tn h))
253 (setf i1 (f2cl-lib:int-add nqnyh 1))
254 (f2cl-lib:fdo (jb 1 (f2cl-lib:int-add jb 1))
255 ((> jb nq) nil)
256 (tagbody
257 (setf i1 (f2cl-lib:int-sub i1 nyh))
258 (f2cl-lib:fdo (i i1 (f2cl-lib:int-add i 1))
259 ((> i nqnyh) nil)
260 (tagbody
261 label210
262 (setf (f2cl-lib:fref yh1-%data% (i) ((1 *)) yh1-%offset%)
264 (f2cl-lib:fref yh1-%data%
266 ((1 *))
267 yh1-%offset%)
268 (f2cl-lib:fref yh1-%data%
269 ((f2cl-lib:int-add i nyh))
270 ((1 *))
271 yh1-%offset%)))))
272 label215))
273 (setf pnorm (dmnorm n yh1 ewt))
274 label220
275 (setf m 0)
276 (setf rate 0.0d0)
277 (setf del 0.0d0)
278 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
279 ((> i n) nil)
280 (tagbody
281 label230
282 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
283 (f2cl-lib:fref yh-%data%
284 (i 1)
285 ((1 nyh) (1 *))
286 yh-%offset%))))
287 (multiple-value-bind (var-0 var-1 var-2 var-3)
288 (funcall f neq tn y savf)
289 (declare (ignore var-0 var-2 var-3))
290 (when var-1
291 (setf tn var-1)))
292 (setf nfe (f2cl-lib:int-add nfe 1))
293 (if (<= ipup 0) (go label250))
294 (multiple-value-bind
295 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
296 var-10)
297 (funcall pjac neq y yh nyh ewt acor savf wm iwm f jac)
298 (declare (ignore var-0 var-1 var-2 var-4 var-5 var-6 var-7 var-8
299 var-9 var-10))
300 (when var-3
301 (setf nyh var-3)))
302 (setf ipup 0)
303 (setf rc 1.0d0)
304 (setf nslp nst)
305 (setf crate 0.7d0)
306 (if (/= ierpj 0) (go label430))
307 label250
308 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
309 ((> i n) nil)
310 (tagbody
311 label260
312 (setf (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
313 0.0d0)))
314 label270
315 (if (/= miter 0) (go label350))
316 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
317 ((> i n) nil)
318 (tagbody
319 (setf (f2cl-lib:fref savf-%data% (i) ((1 *)) savf-%offset%)
321 (* h
322 (f2cl-lib:fref savf-%data%
324 ((1 *))
325 savf-%offset%))
326 (f2cl-lib:fref yh-%data%
327 (i 2)
328 ((1 nyh) (1 *))
329 yh-%offset%)))
330 label290
331 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
333 (f2cl-lib:fref savf-%data% (i) ((1 *)) savf-%offset%)
334 (f2cl-lib:fref acor-%data%
336 ((1 *))
337 acor-%offset%)))))
338 (setf del (dmnorm n y ewt))
339 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
340 ((> i n) nil)
341 (tagbody
342 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
344 (f2cl-lib:fref yh-%data%
345 (i 1)
346 ((1 nyh) (1 *))
347 yh-%offset%)
348 (* (f2cl-lib:fref el (1) ((1 13)))
349 (f2cl-lib:fref savf-%data%
351 ((1 *))
352 savf-%offset%))))
353 label300
354 (setf (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
355 (f2cl-lib:fref savf-%data%
357 ((1 *))
358 savf-%offset%))))
359 (go label400)
360 label350
361 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
362 ((> i n) nil)
363 (tagbody
364 label360
365 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
367 (* h
368 (f2cl-lib:fref savf-%data%
370 ((1 *))
371 savf-%offset%))
373 (f2cl-lib:fref yh-%data%
374 (i 2)
375 ((1 nyh) (1 *))
376 yh-%offset%)
377 (f2cl-lib:fref acor-%data%
379 ((1 *))
380 acor-%offset%))))))
381 (funcall slvs wm iwm y savf)
382 (if (< iersl 0) (go label430))
383 (if (> iersl 0) (go label410))
384 (setf del (dmnorm n y ewt))
385 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
386 ((> i n) nil)
387 (tagbody
388 (setf (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
390 (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
391 (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)))
392 label380
393 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
395 (f2cl-lib:fref yh-%data%
396 (i 1)
397 ((1 nyh) (1 *))
398 yh-%offset%)
399 (* (f2cl-lib:fref el (1) ((1 13)))
400 (f2cl-lib:fref acor-%data%
402 ((1 *))
403 acor-%offset%))))))
404 label400
405 (if (<= del (* 100.0d0 pnorm uround)) (go label450))
406 (if (and (= m 0) (= meth 1)) (go label405))
407 (if (= m 0) (go label402))
408 (setf rm 1024.0d0)
409 (if (<= del (* 1024.0d0 delp)) (setf rm (/ del delp)))
410 (setf rate (max rate rm))
411 (setf crate (max (* 0.2d0 crate) rm))
412 label402
413 (setf dcon
414 (/ (* del (min 1.0d0 (* 1.5d0 crate)))
415 (* (f2cl-lib:fref tesco (2 nq) ((1 3) (1 12))) conit)))
416 (if (> dcon 1.0d0) (go label405))
417 (setf pdest
418 (max pdest
419 (/ rate (abs (* h (f2cl-lib:fref el (1) ((1 13))))))))
420 (if (/= pdest 0.0d0) (setf pdlast pdest))
421 (go label450)
422 label405
423 (setf m (f2cl-lib:int-add m 1))
424 (if (= m maxcor) (go label410))
425 (if (and (>= m 2) (> del (* 2.0d0 delp))) (go label410))
426 (setf delp del)
427 (multiple-value-bind (var-0 var-1 var-2 var-3)
428 (funcall f neq tn y savf)
429 (declare (ignore var-0 var-2 var-3))
430 (when var-1
431 (setf tn var-1)))
432 (setf nfe (f2cl-lib:int-add nfe 1))
433 (go label270)
434 label410
435 (if (or (= miter 0) (= jcur 1)) (go label430))
436 (setf icf 1)
437 (setf ipup miter)
438 (go label220)
439 label430
440 (setf icf 2)
441 (setf ncf (f2cl-lib:int-add ncf 1))
442 (setf rmax 2.0d0)
443 (setf tn told)
444 (setf i1 (f2cl-lib:int-add nqnyh 1))
445 (f2cl-lib:fdo (jb 1 (f2cl-lib:int-add jb 1))
446 ((> jb nq) nil)
447 (tagbody
448 (setf i1 (f2cl-lib:int-sub i1 nyh))
449 (f2cl-lib:fdo (i i1 (f2cl-lib:int-add i 1))
450 ((> i nqnyh) nil)
451 (tagbody
452 label440
453 (setf (f2cl-lib:fref yh1-%data% (i) ((1 *)) yh1-%offset%)
455 (f2cl-lib:fref yh1-%data%
457 ((1 *))
458 yh1-%offset%)
459 (f2cl-lib:fref yh1-%data%
460 ((f2cl-lib:int-add i nyh))
461 ((1 *))
462 yh1-%offset%)))))
463 label445))
464 (if (or (< ierpj 0) (< iersl 0)) (go label680))
465 (if (<= (abs h) (* hmin 1.00001d0)) (go label670))
466 (if (= ncf mxncf) (go label670))
467 (setf rh 0.25d0)
468 (setf ipup miter)
469 (setf iredo 1)
470 (go label170)
471 label450
472 (setf jcur 0)
473 (if (= m 0)
474 (setf dsm (/ del (f2cl-lib:fref tesco (2 nq) ((1 3) (1 12))))))
475 (if (> m 0)
476 (setf dsm
477 (/ (dmnorm n acor ewt)
478 (f2cl-lib:fref tesco (2 nq) ((1 3) (1 12))))))
479 (if (> dsm 1.0d0) (go label500))
480 (setf kflag 0)
481 (setf iredo 0)
482 (setf nst (f2cl-lib:int-add nst 1))
483 (setf hu h)
484 (setf nqu nq)
485 (setf mused meth)
486 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
487 ((> j l) nil)
488 (tagbody
489 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
490 ((> i n) nil)
491 (tagbody
492 (setf (f2cl-lib:fref yh-%data%
493 (i j)
494 ((1 nyh) (1 *))
495 yh-%offset%)
497 (f2cl-lib:fref yh-%data%
498 (i j)
499 ((1 nyh) (1 *))
500 yh-%offset%)
501 (* (f2cl-lib:fref el (j) ((1 13)))
502 (f2cl-lib:fref acor-%data%
504 ((1 *))
505 acor-%offset%))))))))
506 label460
507 (setf icount (f2cl-lib:int-sub icount 1))
508 (if (>= icount 0) (go label488))
509 (if (= meth 2) (go label480))
510 (if (> nq 5) (go label488))
511 (if (and (> dsm (* 100.0d0 pnorm uround)) (/= pdest 0.0d0))
512 (go label470))
513 (if (= irflag 0) (go label488))
514 (setf rh2 2.0d0)
515 (setf nqm2
516 (min (the f2cl-lib:integer4 nq)
517 (the f2cl-lib:integer4 mxords)))
518 (go label478)
519 label470
520 (setf exsm (/ 1.0d0 l))
521 (setf rh1 (/ 1.0d0 (+ (* 1.2d0 (expt dsm exsm)) 1.2d-6)))
522 (setf rh1it (* 2.0d0 rh1))
523 (setf pdh (* pdlast (abs h)))
524 (if (> (* pdh rh1) 1.0d-5)
525 (setf rh1it (/ (f2cl-lib:fref sm1 (nq) ((1 12))) pdh)))
526 (setf rh1 (min rh1 rh1it))
527 (if (<= nq mxords) (go label474))
528 (setf nqm2 mxords)
529 (setf lm2 (f2cl-lib:int-add mxords 1))
530 (setf exm2 (/ 1.0d0 lm2))
531 (setf lm2p1 (f2cl-lib:int-add lm2 1))
532 (setf dm2
534 (dmnorm n
535 (f2cl-lib:array-slice yh-%data%
536 double-float
537 (1 lm2p1)
538 ((1 nyh) (1 *))
539 yh-%offset%)
540 ewt)
541 (f2cl-lib:fref cm2 (mxords) ((1 5)))))
542 (setf rh2 (/ 1.0d0 (+ (* 1.2d0 (expt dm2 exm2)) 1.2d-6)))
543 (go label476)
544 label474
545 (setf dm2
546 (* dsm
547 (/ (f2cl-lib:fref cm1 (nq) ((1 12)))
548 (f2cl-lib:fref cm2 (nq) ((1 5))))))
549 (setf rh2 (/ 1.0d0 (+ (* 1.2d0 (expt dm2 exsm)) 1.2d-6)))
550 (setf nqm2 nq)
551 label476
552 (if (< rh2 (* ratio rh1)) (go label488))
553 label478
554 (setf rh rh2)
555 (setf icount 20)
556 (setf meth 2)
557 (setf miter jtyp)
558 (setf pdlast 0.0d0)
559 (setf nq nqm2)
560 (setf l (f2cl-lib:int-add nq 1))
561 (go label170)
562 label480
563 (setf exsm (/ 1.0d0 l))
564 (if (>= mxordn nq) (go label484))
565 (setf nqm1 mxordn)
566 (setf lm1 (f2cl-lib:int-add mxordn 1))
567 (setf exm1 (/ 1.0d0 lm1))
568 (setf lm1p1 (f2cl-lib:int-add lm1 1))
569 (setf dm1
571 (dmnorm n
572 (f2cl-lib:array-slice yh-%data%
573 double-float
574 (1 lm1p1)
575 ((1 nyh) (1 *))
576 yh-%offset%)
577 ewt)
578 (f2cl-lib:fref cm1 (mxordn) ((1 12)))))
579 (setf rh1 (/ 1.0d0 (+ (* 1.2d0 (expt dm1 exm1)) 1.2d-6)))
580 (go label486)
581 label484
582 (setf dm1
583 (* dsm
584 (/ (f2cl-lib:fref cm2 (nq) ((1 5)))
585 (f2cl-lib:fref cm1 (nq) ((1 12))))))
586 (setf rh1 (/ 1.0d0 (+ (* 1.2d0 (expt dm1 exsm)) 1.2d-6)))
587 (setf nqm1 nq)
588 (setf exm1 exsm)
589 label486
590 (setf rh1it (* 2.0d0 rh1))
591 (setf pdh (* pdnorm (abs h)))
592 (if (> (* pdh rh1) 1.0d-5)
593 (setf rh1it (/ (f2cl-lib:fref sm1 (nqm1) ((1 12))) pdh)))
594 (setf rh1 (min rh1 rh1it))
595 (setf rh2 (/ 1.0d0 (+ (* 1.2d0 (expt dsm exsm)) 1.2d-6)))
596 (if (< (* rh1 ratio) (* 5.0d0 rh2)) (go label488))
597 (setf alpha (max 0.001d0 rh1))
598 (setf dm1 (* (expt alpha exm1) dm1))
599 (if (<= dm1 (* 1000.0d0 uround pnorm)) (go label488))
600 (setf rh rh1)
601 (setf icount 20)
602 (setf meth 1)
603 (setf miter 0)
604 (setf pdlast 0.0d0)
605 (setf nq nqm1)
606 (setf l (f2cl-lib:int-add nq 1))
607 (go label170)
608 label488
609 (setf ialth (f2cl-lib:int-sub ialth 1))
610 (if (= ialth 0) (go label520))
611 (if (> ialth 1) (go label700))
612 (if (= l lmax) (go label700))
613 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
614 ((> i n) nil)
615 (tagbody
616 label490
617 (setf (f2cl-lib:fref yh-%data%
618 (i lmax)
619 ((1 nyh) (1 *))
620 yh-%offset%)
621 (f2cl-lib:fref acor-%data%
623 ((1 *))
624 acor-%offset%))))
625 (go label700)
626 label500
627 (setf kflag (f2cl-lib:int-sub kflag 1))
628 (setf tn told)
629 (setf i1 (f2cl-lib:int-add nqnyh 1))
630 (f2cl-lib:fdo (jb 1 (f2cl-lib:int-add jb 1))
631 ((> jb nq) nil)
632 (tagbody
633 (setf i1 (f2cl-lib:int-sub i1 nyh))
634 (f2cl-lib:fdo (i i1 (f2cl-lib:int-add i 1))
635 ((> i nqnyh) nil)
636 (tagbody
637 label510
638 (setf (f2cl-lib:fref yh1-%data% (i) ((1 *)) yh1-%offset%)
640 (f2cl-lib:fref yh1-%data%
642 ((1 *))
643 yh1-%offset%)
644 (f2cl-lib:fref yh1-%data%
645 ((f2cl-lib:int-add i nyh))
646 ((1 *))
647 yh1-%offset%)))))
648 label515))
649 (setf rmax 2.0d0)
650 (if (<= (abs h) (* hmin 1.00001d0)) (go label660))
651 (if (<= kflag -3) (go label640))
652 (setf iredo 2)
653 (setf rhup 0.0d0)
654 (go label540)
655 label520
656 (setf rhup 0.0d0)
657 (if (= l lmax) (go label540))
658 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
659 ((> i n) nil)
660 (tagbody
661 label530
662 (setf (f2cl-lib:fref savf-%data% (i) ((1 *)) savf-%offset%)
664 (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
665 (f2cl-lib:fref yh-%data%
666 (i lmax)
667 ((1 nyh) (1 *))
668 yh-%offset%)))))
669 (setf dup
670 (/ (dmnorm n savf ewt)
671 (f2cl-lib:fref tesco (3 nq) ((1 3) (1 12)))))
672 (setf exup (/ 1.0d0 (f2cl-lib:int-add l 1)))
673 (setf rhup (/ 1.0d0 (+ (* 1.4d0 (expt dup exup)) 1.4d-6)))
674 label540
675 (setf exsm (/ 1.0d0 l))
676 (setf rhsm (/ 1.0d0 (+ (* 1.2d0 (expt dsm exsm)) 1.2d-6)))
677 (setf rhdn 0.0d0)
678 (if (= nq 1) (go label550))
679 (setf ddn
681 (dmnorm n
682 (f2cl-lib:array-slice yh-%data%
683 double-float
684 (1 l)
685 ((1 nyh) (1 *))
686 yh-%offset%)
687 ewt)
688 (f2cl-lib:fref tesco (1 nq) ((1 3) (1 12)))))
689 (setf exdn (/ 1.0d0 nq))
690 (setf rhdn (/ 1.0d0 (+ (* 1.3d0 (expt ddn exdn)) 1.3d-6)))
691 label550
692 (if (= meth 2) (go label560))
693 (setf pdh (max (* (abs h) pdlast) 1.0d-6))
694 (if (< l lmax)
695 (setf rhup
696 (min rhup (/ (f2cl-lib:fref sm1 (l) ((1 12))) pdh))))
697 (setf rhsm (min rhsm (/ (f2cl-lib:fref sm1 (nq) ((1 12))) pdh)))
698 (if (> nq 1)
699 (setf rhdn
700 (min rhdn
702 (f2cl-lib:fref sm1
703 ((f2cl-lib:int-sub nq 1))
704 ((1 12)))
705 pdh))))
706 (setf pdest 0.0d0)
707 label560
708 (if (>= rhsm rhup) (go label570))
709 (if (> rhup rhdn) (go label590))
710 (go label580)
711 label570
712 (if (< rhsm rhdn) (go label580))
713 (setf newq nq)
714 (setf rh rhsm)
715 (go label620)
716 label580
717 (setf newq (f2cl-lib:int-sub nq 1))
718 (setf rh rhdn)
719 (if (and (< kflag 0) (> rh 1.0d0)) (setf rh 1.0d0))
720 (go label620)
721 label590
722 (setf newq l)
723 (setf rh rhup)
724 (if (< rh 1.1d0) (go label610))
725 (setf r (/ (f2cl-lib:fref el (l) ((1 13))) l))
726 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
727 ((> i n) nil)
728 (tagbody
729 label600
730 (setf (f2cl-lib:fref yh-%data%
731 (i (f2cl-lib:int-add newq 1))
732 ((1 nyh) (1 *))
733 yh-%offset%)
735 (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
736 r))))
737 (go label630)
738 label610
739 (setf ialth 3)
740 (go label700)
741 label620
742 (if (= meth 2) (go label622))
743 (if (>= (* rh pdh 1.00001d0) (f2cl-lib:fref sm1 (newq) ((1 12))))
744 (go label625))
745 label622
746 (if (and (= kflag 0) (< rh 1.1d0)) (go label610))
747 label625
748 (if (<= kflag -2) (setf rh (min rh 0.2d0)))
749 (if (= newq nq) (go label170))
750 label630
751 (setf nq newq)
752 (setf l (f2cl-lib:int-add nq 1))
753 (setf iret 2)
754 (go label150)
755 label640
756 (if (= kflag -10) (go label660))
757 (setf rh 0.1d0)
758 (setf rh (max (/ hmin (abs h)) rh))
759 (setf h (* h rh))
760 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
761 ((> i n) nil)
762 (tagbody
763 label645
764 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
765 (f2cl-lib:fref yh-%data%
766 (i 1)
767 ((1 nyh) (1 *))
768 yh-%offset%))))
769 (multiple-value-bind (var-0 var-1 var-2 var-3)
770 (funcall f neq tn y savf)
771 (declare (ignore var-0 var-2 var-3))
772 (when var-1
773 (setf tn var-1)))
774 (setf nfe (f2cl-lib:int-add nfe 1))
775 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
776 ((> i n) nil)
777 (tagbody
778 label650
779 (setf (f2cl-lib:fref yh-%data%
780 (i 2)
781 ((1 nyh) (1 *))
782 yh-%offset%)
783 (* h
784 (f2cl-lib:fref savf-%data%
786 ((1 *))
787 savf-%offset%)))))
788 (setf ipup miter)
789 (setf ialth 5)
790 (if (= nq 1) (go label200))
791 (setf nq 1)
792 (setf l 2)
793 (setf iret 3)
794 (go label150)
795 label660
796 (setf kflag -1)
797 (go label720)
798 label670
799 (setf kflag -2)
800 (go label720)
801 label680
802 (setf kflag -3)
803 (go label720)
804 label690
805 (setf rmax 10.0d0)
806 label700
807 (setf r (/ 1.0d0 (f2cl-lib:fref tesco (2 nqu) ((1 3) (1 12)))))
808 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
809 ((> i n) nil)
810 (tagbody
811 label710
812 (setf (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
814 (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
815 r))))
816 label720
817 (setf hold h)
818 (setf jstart 1)
819 (go end_label)
820 end_label
821 (return
822 (values nil
835 nil))))))))
837 (in-package #:cl-user)
838 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
839 (eval-when (:load-toplevel :compile-toplevel :execute)
840 (setf (gethash 'fortran-to-lisp::dstoda
841 fortran-to-lisp::*f2cl-function-info*)
842 (fortran-to-lisp::make-f2cl-finfo
843 :arg-types '((array fortran-to-lisp::integer4 (*))
844 (array double-float (*)) (array double-float (*))
845 (fortran-to-lisp::integer4) (array double-float (*))
846 (array double-float (*)) (array double-float (*))
847 (array double-float (*)) (array double-float (*))
848 (array fortran-to-lisp::integer4 (*)) t t t t)
849 :return-values '(nil nil nil fortran-to-lisp::nyh nil nil nil nil
850 nil nil nil nil nil nil)
851 :calls '(fortran-to-lisp::dmnorm fortran-to-lisp::dcfode))))