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