?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
In documentation for lreduce and rreduce, supply second argument as an explicit list
[maxima.git] / share / odepack / src / dstoka.lisp
blobef753564431cd9444e9cd7a8956c769578c2b79c
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 $")
9
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))
16
17 (in-package "ODEPACK")
18
19
20 (defun dstoka (neq y yh nyh yh1 ewt savf savx acor wm iwm f jac psol)
21 (declare (type (f2cl-lib:integer4) nyh)
22 (type (array double-float (*)) wm acor savx savf ewt yh1 yh y)
23 (type (array f2cl-lib:integer4 (*)) iwm neq))
24 (let ((dls001-el
25 (make-array 13
26 :element-type 'double-float
27 :displaced-to (dls001-part-0 *dls001-common-block*)
28 :displaced-index-offset 2))
29 (dls001-elco
30 (make-array 156
31 :element-type 'double-float
32 :displaced-to (dls001-part-0 *dls001-common-block*)
33 :displaced-index-offset 15))
34 (dls001-tesco
35 (make-array 36
36 :element-type 'double-float
37 :displaced-to (dls001-part-0 *dls001-common-block*)
38 :displaced-index-offset 173)))
39 (symbol-macrolet ((conit (aref (dls001-part-0 *dls001-common-block*) 0))
40 (crate (aref (dls001-part-0 *dls001-common-block*) 1))
41 (el dls001-el)
42 (elco dls001-elco)
43 (hold (aref (dls001-part-0 *dls001-common-block*) 171))
44 (rmax (aref (dls001-part-0 *dls001-common-block*) 172))
45 (tesco dls001-tesco)
46 (ccmax (aref (dls001-part-0 *dls001-common-block*) 209))
47 (el0 (aref (dls001-part-0 *dls001-common-block*) 210))
48 (h (aref (dls001-part-0 *dls001-common-block*) 211))
49 (hmin (aref (dls001-part-0 *dls001-common-block*) 212))
50 (hmxi (aref (dls001-part-0 *dls001-common-block*) 213))
51 (hu (aref (dls001-part-0 *dls001-common-block*) 214))
52 (rc (aref (dls001-part-0 *dls001-common-block*) 215))
53 (tn (aref (dls001-part-0 *dls001-common-block*) 216))
54 (ialth (aref (dls001-part-1 *dls001-common-block*) 6))
55 (ipup (aref (dls001-part-1 *dls001-common-block*) 7))
56 (lmax (aref (dls001-part-1 *dls001-common-block*) 8))
57 (meo (aref (dls001-part-1 *dls001-common-block*) 9))
58 (nqnyh (aref (dls001-part-1 *dls001-common-block*) 10))
59 (nslp (aref (dls001-part-1 *dls001-common-block*) 11))
60 (icf (aref (dls001-part-1 *dls001-common-block*) 12))
61 (ierpj (aref (dls001-part-1 *dls001-common-block*) 13))
62 (iersl (aref (dls001-part-1 *dls001-common-block*) 14))
63 (jcur (aref (dls001-part-1 *dls001-common-block*) 15))
64 (jstart (aref (dls001-part-1 *dls001-common-block*) 16))
65 (kflag (aref (dls001-part-1 *dls001-common-block*) 17))
66 (l (aref (dls001-part-1 *dls001-common-block*) 18))
67 (meth (aref (dls001-part-1 *dls001-common-block*) 25))
68 (miter (aref (dls001-part-1 *dls001-common-block*) 26))
69 (maxord (aref (dls001-part-1 *dls001-common-block*) 27))
70 (maxcor (aref (dls001-part-1 *dls001-common-block*) 28))
71 (msbp (aref (dls001-part-1 *dls001-common-block*) 29))
72 (mxncf (aref (dls001-part-1 *dls001-common-block*) 30))
73 (n (aref (dls001-part-1 *dls001-common-block*) 31))
74 (nq (aref (dls001-part-1 *dls001-common-block*) 32))
75 (nst (aref (dls001-part-1 *dls001-common-block*) 33))
76 (nfe (aref (dls001-part-1 *dls001-common-block*) 34))
77 (nqu (aref (dls001-part-1 *dls001-common-block*) 36))
78 (stifr (aref (dls002-part-0 *dls002-common-block*) 0))
79 (newt (aref (dls002-part-1 *dls002-common-block*) 0))
80 (nsfi (aref (dls002-part-1 *dls002-common-block*) 1))
81 (nslj (aref (dls002-part-1 *dls002-common-block*) 2))
82 (njev (aref (dls002-part-1 *dls002-common-block*) 3))
83 (epcon (aref (dlpk01-part-0 *dlpk01-common-block*) 1))
84 (jacflg (aref (dlpk01-part-1 *dlpk01-common-block*) 1))
85 (mnewt (aref (dlpk01-part-1 *dlpk01-common-block*) 7))
86 (ncfn (aref (dlpk01-part-1 *dlpk01-common-block*) 11)))
87 (f2cl-lib:with-multi-array-data
88 ((neq f2cl-lib:integer4 neq-%data% neq-%offset%)
89 (iwm f2cl-lib:integer4 iwm-%data% iwm-%offset%)
90 (y double-float y-%data% y-%offset%)
91 (yh double-float yh-%data% yh-%offset%)
92 (yh1 double-float yh1-%data% yh1-%offset%)
93 (ewt double-float ewt-%data% ewt-%offset%)
94 (savf double-float savf-%data% savf-%offset%)
95 (savx double-float savx-%data% savx-%offset%)
96 (acor double-float acor-%data% acor-%offset%)
97 (wm double-float wm-%data% wm-%offset%))
98 (prog ((nslow 0) (newq 0) (ncf 0) (m 0) (jok 0) (jb 0) (j 0) (iret 0)
99 (iredo 0) (i1 0) (i 0) (told 0.0d0) (stiff 0.0d0) (roc 0.0d0)
100 (rhup 0.0d0) (rhsm 0.0d0) (rhdn 0.0d0) (rh 0.0d0) (r 0.0d0)
101 (dfnorm 0.0d0) (exup 0.0d0) (exsm 0.0d0) (exdn 0.0d0)
102 (dup 0.0d0) (dsm 0.0d0) (drc 0.0d0) (delp 0.0d0) (del 0.0d0)
103 (ddn 0.0d0) (dcon 0.0d0))
104 (declare (type (double-float) dcon ddn del delp drc dsm dup exdn exsm
105 exup dfnorm r rh rhdn rhsm rhup roc
106 stiff told)
107 (type (f2cl-lib:integer4) i i1 iredo iret j jb jok m ncf
108 newq nslow))
109 (setf kflag 0)
110 (setf told tn)
111 (setf ncf 0)
112 (setf ierpj 0)
113 (setf iersl 0)
114 (setf jcur 0)
115 (setf icf 0)
116 (setf delp 0.0d0)
117 (if (> jstart 0) (go label200))
118 (if (= jstart -1) (go label100))
119 (if (= jstart -2) (go label160))
120 (setf lmax (f2cl-lib:int-add maxord 1))
121 (setf nq 1)
122 (setf l 2)
123 (setf ialth 2)
124 (setf rmax 10000.0d0)
125 (setf rc 0.0d0)
126 (setf el0 1.0d0)
127 (setf crate 0.7d0)
128 (setf hold h)
129 (setf meo meth)
130 (setf nslp 0)
131 (setf nslj 0)
132 (setf ipup 0)
133 (setf iret 3)
134 (setf newt 0)
135 (setf stifr 0.0d0)
136 (go label140)
137 label100
138 (setf ipup miter)
139 (setf lmax (f2cl-lib:int-add maxord 1))
140 (if (= ialth 1) (setf ialth 2))
141 (if (= meth meo) (go label110))
142 (dcfode meth elco tesco)
143 (setf meo meth)
144 (if (> nq maxord) (go label120))
145 (setf ialth l)
146 (setf iret 1)
147 (go label150)
148 label110
149 (if (<= nq maxord) (go label160))
150 label120
151 (setf nq maxord)
152 (setf l lmax)
153 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
154 ((> i l) nil)
155 (tagbody
156 label125
157 (setf (f2cl-lib:fref el (i) ((1 13)))
158 (f2cl-lib:fref elco (i nq) ((1 13) (1 12))))))
159 (setf nqnyh (f2cl-lib:int-mul nq nyh))
160 (setf rc (/ (* rc (f2cl-lib:fref el (1) ((1 13)))) el0))
161 (setf el0 (f2cl-lib:fref el (1) ((1 13))))
162 (setf conit (/ 0.5d0 (f2cl-lib:int-add nq 2)))
163 (setf epcon (* conit (f2cl-lib:fref tesco (2 nq) ((1 3) (1 12)))))
164 (setf ddn
165 (/ (dvnorm n savf ewt)
166 (f2cl-lib:fref tesco (1 l) ((1 3) (1 12)))))
167 (setf exdn (/ 1.0d0 l))
168 (setf rhdn (/ 1.0d0 (+ (* 1.3d0 (expt ddn exdn)) 1.3d-6)))
169 (setf rh (min rhdn 1.0d0))
170 (setf iredo 3)
171 (if (= h hold) (go label170))
172 (setf rh (min rh (abs (/ h hold))))
173 (setf h hold)
174 (go label175)
175 label140
176 (dcfode meth elco tesco)
177 label150
178 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
179 ((> i l) nil)
180 (tagbody
181 label155
182 (setf (f2cl-lib:fref el (i) ((1 13)))
183 (f2cl-lib:fref elco (i nq) ((1 13) (1 12))))))
184 (setf nqnyh (f2cl-lib:int-mul nq nyh))
185 (setf rc (/ (* rc (f2cl-lib:fref el (1) ((1 13)))) el0))
186 (setf el0 (f2cl-lib:fref el (1) ((1 13))))
187 (setf conit (/ 0.5d0 (f2cl-lib:int-add nq 2)))
188 (setf epcon (* conit (f2cl-lib:fref tesco (2 nq) ((1 3) (1 12)))))
189 (f2cl-lib:computed-goto (label160 label170 label200) iret)
190 label160
191 (if (= h hold) (go label200))
192 (setf rh (/ h hold))
193 (setf h hold)
194 (setf iredo 3)
195 (go label175)
196 label170
197 (setf rh (max rh (/ hmin (abs h))))
198 label175
199 (setf rh (min rh rmax))
200 (setf rh (/ rh (max 1.0d0 (* (abs h) hmxi rh))))
201 (setf r 1.0d0)
202 (f2cl-lib:fdo (j 2 (f2cl-lib:int-add j 1))
203 ((> j l) nil)
204 (tagbody
205 (setf r (* r rh))
206 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
207 ((> i n) nil)
208 (tagbody
209 (setf (f2cl-lib:fref yh-%data%
210 (i j)
211 ((1 nyh) (1 *))
212 yh-%offset%)
213 (*
214 (f2cl-lib:fref yh-%data%
215 (i j)
216 ((1 nyh) (1 *))
217 yh-%offset%)
218 r))))))
219 label180
220 (setf h (* h rh))
221 (setf rc (* rc rh))
222 (setf ialth l)
223 (if (= iredo 0) (go label690))
224 label200
225 (cond
226 ((or (= newt 0) (= jacflg 0))
227 (setf drc 0.0d0)
228 (setf ipup 0)
229 (setf crate 0.7d0))
230 (t
231 (setf drc (abs (- rc 1.0d0)))
232 (if (> drc ccmax) (setf ipup miter))
233 (if (>= nst (f2cl-lib:int-add nslp msbp)) (setf ipup miter))))
234 (setf tn (+ tn h))
235 (setf i1 (f2cl-lib:int-add nqnyh 1))
236 (f2cl-lib:fdo (jb 1 (f2cl-lib:int-add jb 1))
237 ((> jb nq) nil)
238 (tagbody
239 (setf i1 (f2cl-lib:int-sub i1 nyh))
240 (f2cl-lib:fdo (i i1 (f2cl-lib:int-add i 1))
241 ((> i nqnyh) nil)
242 (tagbody
243 label210
244 (setf (f2cl-lib:fref yh1-%data% (i) ((1 *)) yh1-%offset%)
245 (+
246 (f2cl-lib:fref yh1-%data% (i) ((1 *)) yh1-%offset%)
247 (f2cl-lib:fref yh1-%data%
248 ((f2cl-lib:int-add i nyh))
249 ((1 *))
250 yh1-%offset%)))))
251 label215))
252 label220
253 (setf m 0)
254 (setf mnewt 0)
255 (setf stiff 0.0d0)
256 (setf roc 0.05d0)
257 (setf nslow 0)
258 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
259 ((> i n) nil)
260 (tagbody
261 label230
262 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
263 (f2cl-lib:fref yh-%data%
264 (i 1)
265 ((1 nyh) (1 *))
266 yh-%offset%))))
267 (multiple-value-bind (var-0 var-1 var-2 var-3)
268 (funcall f neq tn y savf)
269 (declare (ignore var-0 var-2 var-3))
270 (when var-1
271 (setf tn var-1)))
272 (setf nfe (f2cl-lib:int-add nfe 1))
273 (if (or (= newt 0) (<= ipup 0)) (go label250))
274 (setf jok 1)
275 (if (or (= nst 0) (> nst (f2cl-lib:int-add nslj 50))) (setf jok -1))
276 (if (and (= icf 1) (< drc 0.2d0)) (setf jok -1))
277 (if (= icf 2) (setf jok -1))
278 (cond
279 ((= jok (f2cl-lib:int-sub 1))
280 (setf nslj nst)
281 (setf njev (f2cl-lib:int-add njev 1))))
282 (multiple-value-bind
283 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
284 var-10)
285 (dsetpk neq y yh1 ewt acor savf jok wm iwm f jac)
286 (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-7 var-8
287 var-9 var-10))
288 (setf jok var-6))
289 (setf ipup 0)
290 (setf rc 1.0d0)
291 (setf drc 0.0d0)
292 (setf nslp nst)
293 (setf crate 0.7d0)
294 (if (/= ierpj 0) (go label430))
295 label250
296 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
297 ((> i n) nil)
298 (tagbody
299 label260
300 (setf (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
301 0.0d0)))
302 label270
303 (if (/= newt 0) (go label350))
304 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
305 ((> i n) nil)
306 (tagbody
307 (setf (f2cl-lib:fref savf-%data% (i) ((1 *)) savf-%offset%)
308 (-
309 (* h
310 (f2cl-lib:fref savf-%data%
311 (i)
312 ((1 *))
313 savf-%offset%))
314 (f2cl-lib:fref yh-%data%
315 (i 2)
316 ((1 nyh) (1 *))
317 yh-%offset%)))
318 label290
319 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
320 (- (f2cl-lib:fref savf-%data% (i) ((1 *)) savf-%offset%)
321 (f2cl-lib:fref acor-%data%
322 (i)
323 ((1 *))
324 acor-%offset%)))))
325 (setf del (dvnorm n y ewt))
326 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
327 ((> i n) nil)
328 (tagbody
329 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
330 (+
331 (f2cl-lib:fref yh-%data%
332 (i 1)
333 ((1 nyh) (1 *))
334 yh-%offset%)
335 (* (f2cl-lib:fref el (1) ((1 13)))
336 (f2cl-lib:fref savf-%data%
337 (i)
338 ((1 *))
339 savf-%offset%))))
340 label300
341 (setf (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
342 (f2cl-lib:fref savf-%data% (i) ((1 *)) savf-%offset%))))
343 (setf stiff 1.0d0)
344 (go label400)
345 label350
346 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
347 ((> i n) nil)
348 (tagbody
349 label360
350 (setf (f2cl-lib:fref savx-%data% (i) ((1 *)) savx-%offset%)
351 (-
352 (* h
353 (f2cl-lib:fref savf-%data%
354 (i)
355 ((1 *))
356 savf-%offset%))
357 (+
358 (f2cl-lib:fref yh-%data%
359 (i 2)
360 ((1 nyh) (1 *))
361 yh-%offset%)
362 (f2cl-lib:fref acor-%data%
363 (i)
364 ((1 *))
365 acor-%offset%))))))
366 (setf dfnorm (dvnorm n savx ewt))
367 (dsolpk neq y savf savx ewt wm iwm f psol)
368 (if (< iersl 0) (go label430))
369 (if (> iersl 0) (go label410))
370 (setf del (dvnorm n savx ewt))
371 (if (> del 1.0d-8) (setf stiff (max stiff (/ dfnorm del))))
372 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
373 ((> i n) nil)
374 (tagbody
375 (setf (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
376 (+ (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
377 (f2cl-lib:fref savx-%data%
378 (i)
379 ((1 *))
380 savx-%offset%)))
381 label380
382 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
383 (+
384 (f2cl-lib:fref yh-%data%
385 (i 1)
386 ((1 nyh) (1 *))
387 yh-%offset%)
388 (* (f2cl-lib:fref el (1) ((1 13)))
389 (f2cl-lib:fref acor-%data%
390 (i)
391 ((1 *))
392 acor-%offset%))))))
393 label400
394 (cond
395 ((/= m 0)
396 (setf roc (max 0.05d0 (/ del delp)))
397 (setf crate (max (* 0.2d0 crate) roc))))
398 (setf dcon (/ (* del (min 1.0d0 (* 1.5d0 crate))) epcon))
399 (if (<= dcon 1.0d0) (go label450))
400 (setf m (f2cl-lib:int-add m 1))
401 (if (= m maxcor) (go label410))
402 (if (and (>= m 2) (> del (* 2.0d0 delp))) (go label410))
403 (if (> roc 10.0d0) (go label410))
404 (if (> roc 0.8d0) (setf nslow (f2cl-lib:int-add nslow 1)))
405 (if (>= nslow 2) (go label410))
406 (setf mnewt m)
407 (setf delp del)
408 (multiple-value-bind (var-0 var-1 var-2 var-3)
409 (funcall f neq tn y savf)
410 (declare (ignore var-0 var-2 var-3))
411 (when var-1
412 (setf tn var-1)))
413 (setf nfe (f2cl-lib:int-add nfe 1))
414 (go label270)
415 label410
416 (setf icf 1)
417 (cond
418 ((= newt 0)
419 (if (= nst 0) (go label430))
420 (if (= miter 0) (go label430))
421 (setf newt miter)
422 (setf stifr 1023.0d0)
423 (setf ipup miter)
424 (go label220)))
425 (if (or (= jcur 1) (= jacflg 0)) (go label430))
426 (setf ipup miter)
427 (go label220)
428 label430
429 (setf icf 2)
430 (setf ncf (f2cl-lib:int-add ncf 1))
431 (setf ncfn (f2cl-lib:int-add ncfn 1))
432 (setf rmax 2.0d0)
433 (setf tn told)
434 (setf i1 (f2cl-lib:int-add nqnyh 1))
435 (f2cl-lib:fdo (jb 1 (f2cl-lib:int-add jb 1))
436 ((> jb nq) nil)
437 (tagbody
438 (setf i1 (f2cl-lib:int-sub i1 nyh))
439 (f2cl-lib:fdo (i i1 (f2cl-lib:int-add i 1))
440 ((> i nqnyh) nil)
441 (tagbody
442 label440
443 (setf (f2cl-lib:fref yh1-%data% (i) ((1 *)) yh1-%offset%)
444 (-
445 (f2cl-lib:fref yh1-%data% (i) ((1 *)) yh1-%offset%)
446 (f2cl-lib:fref yh1-%data%
447 ((f2cl-lib:int-add i nyh))
448 ((1 *))
449 yh1-%offset%)))))
450 label445))
451 (if (or (< ierpj 0) (< iersl 0)) (go label680))
452 (if (<= (abs h) (* hmin 1.00001d0)) (go label670))
453 (if (= ncf mxncf) (go label670))
454 (setf rh 0.5d0)
455 (setf ipup miter)
456 (setf iredo 1)
457 (go label170)
458 label450
459 (setf jcur 0)
460 (if (> newt 0) (setf stifr (* 0.5d0 (+ stifr stiff))))
461 (if (= m 0)
462 (setf dsm (/ del (f2cl-lib:fref tesco (2 nq) ((1 3) (1 12))))))
463 (if (> m 0)
464 (setf dsm
465 (/ (dvnorm n acor ewt)
466 (f2cl-lib:fref tesco (2 nq) ((1 3) (1 12))))))
467 (if (> dsm 1.0d0) (go label500))
468 (setf kflag 0)
469 (setf iredo 0)
470 (setf nst (f2cl-lib:int-add nst 1))
471 (if (= newt 0) (setf nsfi (f2cl-lib:int-add nsfi 1)))
472 (if (and (> newt 0) (< stifr 1.5d0)) (setf newt 0))
473 (setf hu h)
474 (setf nqu nq)
475 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
476 ((> j l) nil)
477 (tagbody
478 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
479 ((> i n) nil)
480 (tagbody
481 (setf (f2cl-lib:fref yh-%data%
482 (i j)
483 ((1 nyh) (1 *))
484 yh-%offset%)
485 (+
486 (f2cl-lib:fref yh-%data%
487 (i j)
488 ((1 nyh) (1 *))
489 yh-%offset%)
490 (* (f2cl-lib:fref el (j) ((1 13)))
491 (f2cl-lib:fref acor-%data%
492 (i)
493 ((1 *))
494 acor-%offset%))))))))
495 label470
496 (setf ialth (f2cl-lib:int-sub ialth 1))
497 (if (= ialth 0) (go label520))
498 (if (> ialth 1) (go label700))
499 (if (= l lmax) (go label700))
500 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
501 ((> i n) nil)
502 (tagbody
503 label490
504 (setf (f2cl-lib:fref yh-%data%
505 (i lmax)
506 ((1 nyh) (1 *))
507 yh-%offset%)
508 (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%))))
509 (go label700)
510 label500
511 (setf kflag (f2cl-lib:int-sub kflag 1))
512 (setf tn told)
513 (setf i1 (f2cl-lib:int-add nqnyh 1))
514 (f2cl-lib:fdo (jb 1 (f2cl-lib:int-add jb 1))
515 ((> jb nq) nil)
516 (tagbody
517 (setf i1 (f2cl-lib:int-sub i1 nyh))
518 (f2cl-lib:fdo (i i1 (f2cl-lib:int-add i 1))
519 ((> i nqnyh) nil)
520 (tagbody
521 label510
522 (setf (f2cl-lib:fref yh1-%data% (i) ((1 *)) yh1-%offset%)
523 (-
524 (f2cl-lib:fref yh1-%data% (i) ((1 *)) yh1-%offset%)
525 (f2cl-lib:fref yh1-%data%
526 ((f2cl-lib:int-add i nyh))
527 ((1 *))
528 yh1-%offset%)))))
529 label515))
530 (setf rmax 2.0d0)
531 (if (<= (abs h) (* hmin 1.00001d0)) (go label660))
532 (if (<= kflag -3) (go label640))
533 (setf iredo 2)
534 (setf rhup 0.0d0)
535 (go label540)
536 label520
537 (setf rhup 0.0d0)
538 (if (= l lmax) (go label540))
539 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
540 ((> i n) nil)
541 (tagbody
542 label530
543 (setf (f2cl-lib:fref savf-%data% (i) ((1 *)) savf-%offset%)
544 (- (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
545 (f2cl-lib:fref yh-%data%
546 (i lmax)
547 ((1 nyh) (1 *))
548 yh-%offset%)))))
549 (setf dup
550 (/ (dvnorm n savf ewt)
551 (f2cl-lib:fref tesco (3 nq) ((1 3) (1 12)))))
552 (setf exup (/ 1.0d0 (f2cl-lib:int-add l 1)))
553 (setf rhup (/ 1.0d0 (+ (* 1.4d0 (expt dup exup)) 1.4d-6)))
554 label540
555 (setf exsm (/ 1.0d0 l))
556 (setf rhsm (/ 1.0d0 (+ (* 1.2d0 (expt dsm exsm)) 1.2d-6)))
557 (setf rhdn 0.0d0)
558 (if (= nq 1) (go label560))
559 (setf ddn
560 (/
561 (dvnorm n
562 (f2cl-lib:array-slice yh-%data%
563 double-float
564 (1 l)
565 ((1 nyh) (1 *))
566 yh-%offset%)
567 ewt)
568 (f2cl-lib:fref tesco (1 nq) ((1 3) (1 12)))))
569 (setf exdn (/ 1.0d0 nq))
570 (setf rhdn (/ 1.0d0 (+ (* 1.3d0 (expt ddn exdn)) 1.3d-6)))
571 label560
572 (if (>= rhsm rhup) (go label570))
573 (if (> rhup rhdn) (go label590))
574 (go label580)
575 label570
576 (if (< rhsm rhdn) (go label580))
577 (setf newq nq)
578 (setf rh rhsm)
579 (go label620)
580 label580
581 (setf newq (f2cl-lib:int-sub nq 1))
582 (setf rh rhdn)
583 (if (and (< kflag 0) (> rh 1.0d0)) (setf rh 1.0d0))
584 (go label620)
585 label590
586 (setf newq l)
587 (setf rh rhup)
588 (if (< rh 1.1d0) (go label610))
589 (setf r (/ (f2cl-lib:fref el (l) ((1 13))) l))
590 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
591 ((> i n) nil)
592 (tagbody
593 label600
594 (setf (f2cl-lib:fref yh-%data%
595 (i (f2cl-lib:int-add newq 1))
596 ((1 nyh) (1 *))
597 yh-%offset%)
598 (* (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
599 r))))
600 (go label630)
601 label610
602 (setf ialth 3)
603 (go label700)
604 label620
605 (if (and (= kflag 0) (< rh 1.1d0)) (go label610))
606 (if (<= kflag -2) (setf rh (min rh 0.2d0)))
607 (if (= newq nq) (go label170))
608 label630
609 (setf nq newq)
610 (setf l (f2cl-lib:int-add nq 1))
611 (setf iret 2)
612 (go label150)
613 label640
614 (if (= kflag -10) (go label660))
615 (setf rh 0.1d0)
616 (setf rh (max (/ hmin (abs h)) rh))
617 (setf h (* h rh))
618 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
619 ((> i n) nil)
620 (tagbody
621 label645
622 (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
623 (f2cl-lib:fref yh-%data%
624 (i 1)
625 ((1 nyh) (1 *))
626 yh-%offset%))))
627 (multiple-value-bind (var-0 var-1 var-2 var-3)
628 (funcall f neq tn y savf)
629 (declare (ignore var-0 var-2 var-3))
630 (when var-1
631 (setf tn var-1)))
632 (setf nfe (f2cl-lib:int-add nfe 1))
633 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
634 ((> i n) nil)
635 (tagbody
636 label650
637 (setf (f2cl-lib:fref yh-%data% (i 2) ((1 nyh) (1 *)) yh-%offset%)
638 (* h
639 (f2cl-lib:fref savf-%data%
640 (i)
641 ((1 *))
642 savf-%offset%)))))
643 (setf ipup miter)
644 (setf ialth 5)
645 (if (= nq 1) (go label200))
646 (setf nq 1)
647 (setf l 2)
648 (setf iret 3)
649 (go label150)
650 label660
651 (setf kflag -1)
652 (go label720)
653 label670
654 (setf kflag -2)
655 (go label720)
656 label680
657 (setf kflag -3)
658 (go label720)
659 label690
660 (setf rmax 10.0d0)
661 label700
662 (setf r (/ 1.0d0 (f2cl-lib:fref tesco (2 nqu) ((1 3) (1 12)))))
663 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
664 ((> i n) nil)
665 (tagbody
666 label710
667 (setf (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
668 (* (f2cl-lib:fref acor-%data% (i) ((1 *)) acor-%offset%)
669 r))))
670 label720
671 (setf hold h)
672 (setf jstart 1)
673 (go end_label)
674 end_label
675 (return
676 (values nil
677 nil
678 nil
679 nil
680 nil
681 nil
682 nil
683 nil
684 nil
685 nil
686 nil
687 nil
688 nil
689 nil)))))))
690
691 (in-package #:cl-user)
692 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
693 (eval-when (:load-toplevel :compile-toplevel :execute)
694 (setf (gethash 'fortran-to-lisp::dstoka
695 fortran-to-lisp::*f2cl-function-info*)
696 (fortran-to-lisp::make-f2cl-finfo
697 :arg-types '((array fortran-to-lisp::integer4 (*))
698 (array double-float (*)) (array double-float (*))
699 (fortran-to-lisp::integer4) (array double-float (*))
700 (array double-float (*)) (array double-float (*))
701 (array double-float (*)) (array double-float (*))
702 (array double-float (*))
703 (array fortran-to-lisp::integer4 (*)) t t t)
704 :return-values '(nil nil nil nil nil nil nil nil nil nil nil nil nil
705 nil)
706 :calls '(fortran-to-lisp::dsolpk fortran-to-lisp::dsetpk
707 fortran-to-lisp::dvnorm fortran-to-lisp::dcfode))))
708
Maxima, a Computer Algebra System