?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Remove the obsolete DEFMTRFUN-EXTERNAL macro
[maxima.git] / share / hompack / lisp / steps.lisp
blob50cf3f8bc18fc7c7ebf6c4767886a8f32a4255ee
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-2020-04 (21D 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 double-float))
16
17 (in-package "HOMPACK")
18
19
20 (let ((two
21 (make-array 13
22 :element-type 'double-float
23 :initial-contents '(2.0 4.0 8.0 16.0 32.0 64.0 128.0 256.0
24 512.0 1024.0 2048.0 4096.0 8192.0)))
25 (gstr
26 (make-array 13
27 :element-type 'double-float
28 :initial-contents '(0.5 0.0833 0.0417 0.0264 0.0188 0.0143
29 0.0114 0.00936 0.00789 0.00679 0.00592
30 0.00524 0.00468)))
31 (psi (make-array 12 :element-type 'double-float))
32 (beta (make-array 12 :element-type 'double-float))
33 (sig (make-array 13 :element-type 'double-float))
34 (v (make-array 12 :element-type 'double-float))
35 (phase1 nil)
36 (nornd nil)
37 (i 0)
38 (ifail 0)
39 (im1 0)
40 (ip1 0)
41 (iq 0)
42 (j 0)
43 (jv 0)
44 (km1 0)
45 (km2 0)
46 (knew 0)
47 (kp1 0)
48 (kp2 0)
49 (kprev 0)
50 (l 0)
51 (limit1 0)
52 (limit2 0)
53 (ns 0)
54 (nsm2 0)
55 (nsp1 0)
56 (nsp2 0)
57 (absh 0.0)
58 (erk 0.0)
59 (erkm1 0.0)
60 (erkm2 0.0)
61 (erkp1 0.0)
62 (err 0.0)
63 (fouru 0.0)
64 (hnew 0.0)
65 (p5eps 0.0)
66 (r 0.0)
67 (reali 0.0)
68 (realns 0.0)
69 (rho 0.0)
70 (round$ 0.0)
71 (sum 0.0)
72 (tau 0.0)
73 (temp1 0.0)
74 (temp2 0.0)
75 (temp3 0.0)
76 (temp4 0.0)
77 (temp5 0.0)
78 (temp6 0.0)
79 (twou 0.0))
80 (declare (type f2cl-lib:logical phase1 nornd)
81 (type (f2cl-lib:integer4) i ifail im1 ip1 iq j jv km1 km2 knew kp1
82 kp2 kprev l limit1 limit2 ns nsm2 nsp1
83 nsp2)
84 (type (double-float) absh erk erkm1 erkm2 erkp1 err fouru hnew p5eps
85 r reali realns rho round$ sum tau temp1 temp2
86 temp3 temp4 temp5 temp6 twou)
87 (type (array double-float (12)) psi beta v)
88 (type (array double-float (13)) two gstr sig))
89 (defun steps
90 (f neqn y x h eps wt start hold k kold crash phi p yp alpha w g ksteps
91 xold ivc iv kgi gi fpwa1 fpwa2 fpwa3 fpwa4 fpwa5 ifpwa1 ifpc1 ifpc2
92 par ipar)
93 (declare (type (array f2cl-lib:integer4 (*)) ipar)
94 (type (array double-float (*)) par)
95 (type (array f2cl-lib:integer4 (*)) ifpwa1)
96 (type (array double-float (*)) gi)
97 (type (array f2cl-lib:integer4 (*)) iv)
98 (type (array double-float (*)) g)
99 (type (array double-float (*)) w alpha)
100 (type f2cl-lib:logical crash start)
101 (type (double-float) xold hold eps h x)
102 (type (array double-float (*)) fpwa5 fpwa4 fpwa3 fpwa2 fpwa1 yp p
103 phi wt y)
104 (type (f2cl-lib:integer4) ifpc2 ifpc1 kgi ivc ksteps kold k neqn))
105 (f2cl-lib:with-multi-array-data
106 ((y double-float y-%data% y-%offset%)
107 (wt double-float wt-%data% wt-%offset%)
108 (phi double-float phi-%data% phi-%offset%)
109 (p double-float p-%data% p-%offset%)
110 (yp double-float yp-%data% yp-%offset%)
111 (fpwa1 double-float fpwa1-%data% fpwa1-%offset%)
112 (fpwa2 double-float fpwa2-%data% fpwa2-%offset%)
113 (fpwa3 double-float fpwa3-%data% fpwa3-%offset%)
114 (fpwa4 double-float fpwa4-%data% fpwa4-%offset%)
115 (fpwa5 double-float fpwa5-%data% fpwa5-%offset%)
116 (alpha double-float alpha-%data% alpha-%offset%)
117 (w double-float w-%data% w-%offset%)
118 (g double-float g-%data% g-%offset%)
119 (iv f2cl-lib:integer4 iv-%data% iv-%offset%)
120 (gi double-float gi-%data% gi-%offset%)
121 (ifpwa1 f2cl-lib:integer4 ifpwa1-%data% ifpwa1-%offset%)
122 (par double-float par-%data% par-%offset%)
123 (ipar f2cl-lib:integer4 ipar-%data% ipar-%offset%))
124 (prog ()
125 (declare)
126 (setf twou (* 2.0f0 (f2cl-lib:d1mach 4)))
127 (setf fouru (+ twou twou))
128 (setf crash f2cl-lib:%true%)
129 (if (>= (abs h) (* fouru (abs x))) (go label5))
130 (setf h (f2cl-lib:sign (* fouru (abs x)) h))
131 (go end_label)
132 label5
133 (setf p5eps (* 0.5f0 eps))
134 (setf round$ (coerce 0.0f0 'double-float))
135 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
136 ((> l neqn) nil)
137 (tagbody
138 label10
139 (setf round$
140 (+ round$
141 (expt
142 (/ (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
143 (f2cl-lib:fref wt-%data%
144 (l)
145 ((1 neqn))
146 wt-%offset%))
147 2)))))
148 (setf round$ (* twou (f2cl-lib:fsqrt round$)))
149 (if (>= p5eps round$) (go label15))
150 (setf eps (* 2.0f0 round$ (+ 1.0f0 fouru)))
151 (go end_label)
152 label15
153 (setf crash f2cl-lib:%false%)
154 (setf (f2cl-lib:fref g-%data% (1) ((1 13)) g-%offset%)
155 (coerce 1.0f0 'double-float))
156 (setf (f2cl-lib:fref g-%data% (2) ((1 13)) g-%offset%)
157 (coerce 0.5f0 'double-float))
158 (setf (f2cl-lib:fref sig (1) ((1 13))) (coerce 1.0f0 'double-float))
159 (if (not start) (go label99))
160 (multiple-value-bind
161 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
162 var-10 var-11 var-12 var-13)
163 (funcall f
164 x
165 y
166 yp
167 fpwa1
168 fpwa2
169 fpwa3
170 fpwa4
171 fpwa5
172 ifpwa1
173 ifpc1
174 (f2cl-lib:int-sub neqn 1)
175 ifpc2
176 par
177 ipar)
178 (declare (ignore var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
179 var-10 var-12 var-13))
180 (when var-0
181 (setf x var-0))
182 (when var-9
183 (setf ifpc1 var-9))
184 (when var-11
185 (setf ifpc2 var-11)))
186 (if (> ifpc2 0) (go end_label))
187 (setf sum (coerce 0.0f0 'double-float))
188 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
189 ((> l neqn) nil)
190 (tagbody
191 (setf (f2cl-lib:fref phi-%data%
192 (l 1)
193 ((1 neqn) (1 16))
194 phi-%offset%)
195 (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%))
196 (setf (f2cl-lib:fref phi-%data%
197 (l 2)
198 ((1 neqn) (1 16))
199 phi-%offset%)
200 (coerce 0.0f0 'double-float))
201 label20
202 (setf sum
203 (+ sum
204 (expt
205 (/ (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%)
206 (f2cl-lib:fref wt-%data%
207 (l)
208 ((1 neqn))
209 wt-%offset%))
210 2)))))
211 (setf sum (f2cl-lib:fsqrt sum))
212 (setf absh (abs h))
213 (if (< eps (* 16.0f0 sum h h))
214 (setf absh (* 0.25f0 (f2cl-lib:fsqrt (/ eps sum)))))
215 (setf h (f2cl-lib:sign (max absh (* fouru (abs x))) h))
216 (setf hold (coerce 0.0f0 'double-float))
217 (setf k 1)
218 (setf kold 0)
219 (setf kprev 0)
220 (setf start f2cl-lib:%false%)
221 (setf phase1 f2cl-lib:%true%)
222 (setf nornd f2cl-lib:%true%)
223 (if (> p5eps (* 100.0f0 round$)) (go label99))
224 (setf nornd f2cl-lib:%false%)
225 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
226 ((> l neqn) nil)
227 (tagbody
228 label25
229 (setf (f2cl-lib:fref phi-%data%
230 (l 15)
231 ((1 neqn) (1 16))
232 phi-%offset%)
233 (coerce 0.0f0 'double-float))))
234 label99
235 (setf ifail 0)
236 label100
237 (setf kp1 (f2cl-lib:int-add k 1))
238 (setf kp2 (f2cl-lib:int-add k 2))
239 (setf km1 (f2cl-lib:int-sub k 1))
240 (setf km2 (f2cl-lib:int-sub k 2))
241 (if (/= h hold) (setf ns 0))
242 (if (<= ns kold) (setf ns (f2cl-lib:int-add ns 1)))
243 (setf nsp1 (f2cl-lib:int-add ns 1))
244 (if (< k ns) (go label199))
245 (setf (f2cl-lib:fref beta (ns) ((1 12))) (coerce 1.0f0 'double-float))
246 (setf realns (coerce (the f2cl-lib:integer4 ns) 'double-float))
247 (setf (f2cl-lib:fref alpha-%data% (ns) ((1 12)) alpha-%offset%)
248 (/ 1.0f0 realns))
249 (setf temp1 (* h realns))
250 (setf (f2cl-lib:fref sig (nsp1) ((1 13))) (coerce 1.0f0 'double-float))
251 (if (< k nsp1) (go label110))
252 (f2cl-lib:fdo (i nsp1 (f2cl-lib:int-add i 1))
253 ((> i k) nil)
254 (tagbody
255 (setf im1 (f2cl-lib:int-sub i 1))
256 (setf temp2 (f2cl-lib:fref psi (im1) ((1 12))))
257 (setf (f2cl-lib:fref psi (im1) ((1 12))) temp1)
258 (setf (f2cl-lib:fref beta (i) ((1 12)))
259 (/
260 (* (f2cl-lib:fref beta (im1) ((1 12)))
261 (f2cl-lib:fref psi (im1) ((1 12))))
262 temp2))
263 (setf temp1 (+ temp2 h))
264 (setf (f2cl-lib:fref alpha-%data% (i) ((1 12)) alpha-%offset%)
265 (/ h temp1))
266 (setf reali (coerce (the f2cl-lib:integer4 i) 'double-float))
267 label105
268 (setf (f2cl-lib:fref sig ((f2cl-lib:int-add i 1)) ((1 13)))
269 (* reali
270 (f2cl-lib:fref alpha-%data% (i) ((1 12)) alpha-%offset%)
271 (f2cl-lib:fref sig (i) ((1 13)))))))
272 label110
273 (setf (f2cl-lib:fref psi (k) ((1 12))) temp1)
274 (if (> ns 1) (go label120))
275 (f2cl-lib:fdo (iq 1 (f2cl-lib:int-add iq 1))
276 ((> iq k) nil)
277 (tagbody
278 (setf temp3
279 (coerce
280 (the f2cl-lib:integer4
281 (f2cl-lib:int-mul iq (f2cl-lib:int-add iq 1)))
282 'double-float))
283 (setf (f2cl-lib:fref v (iq) ((1 12))) (/ 1.0f0 temp3))
284 label115
285 (setf (f2cl-lib:fref w-%data% (iq) ((1 12)) w-%offset%)
286 (f2cl-lib:fref v (iq) ((1 12))))))
287 (setf ivc 0)
288 (setf kgi 0)
289 (if (= k 1) (go label140))
290 (setf kgi 1)
291 (setf (f2cl-lib:fref gi-%data% (1) ((1 11)) gi-%offset%)
292 (f2cl-lib:fref w-%data% (2) ((1 12)) w-%offset%))
293 (go label140)
294 label120
295 (if (<= k kprev) (go label130))
296 (if (= ivc 0) (go label122))
297 (setf jv
298 (f2cl-lib:int-sub kp1
299 (f2cl-lib:fref iv-%data%
300 (ivc)
301 ((1 10))
302 iv-%offset%)))
303 (setf ivc (f2cl-lib:int-sub ivc 1))
304 (go label123)
305 label122
306 (setf jv 1)
307 (setf temp4
308 (coerce (the f2cl-lib:integer4 (f2cl-lib:int-mul k kp1))
309 'double-float))
310 (setf (f2cl-lib:fref v (k) ((1 12))) (/ 1.0f0 temp4))
311 (setf (f2cl-lib:fref w-%data% (k) ((1 12)) w-%offset%)
312 (f2cl-lib:fref v (k) ((1 12))))
313 (if (/= k 2) (go label123))
314 (setf kgi 1)
315 (setf (f2cl-lib:fref gi-%data% (1) ((1 11)) gi-%offset%)
316 (f2cl-lib:fref w-%data% (2) ((1 12)) w-%offset%))
317 label123
318 (setf nsm2 (f2cl-lib:int-sub ns 2))
319 (if (< nsm2 jv) (go label130))
320 (f2cl-lib:fdo (j jv (f2cl-lib:int-add j 1))
321 ((> j nsm2) nil)
322 (tagbody
323 (setf i (f2cl-lib:int-sub k j))
324 (setf (f2cl-lib:fref v (i) ((1 12)))
325 (- (f2cl-lib:fref v (i) ((1 12)))
326 (*
327 (f2cl-lib:fref alpha-%data%
328 ((f2cl-lib:int-add j 1))
329 ((1 12))
330 alpha-%offset%)
331 (f2cl-lib:fref v ((f2cl-lib:int-add i 1)) ((1 12))))))
332 label125
333 (setf (f2cl-lib:fref w-%data% (i) ((1 12)) w-%offset%)
334 (f2cl-lib:fref v (i) ((1 12))))))
335 (if (/= i 2) (go label130))
336 (setf kgi (f2cl-lib:int-sub ns 1))
337 (setf (f2cl-lib:fref gi-%data% (kgi) ((1 11)) gi-%offset%)
338 (f2cl-lib:fref w-%data% (2) ((1 12)) w-%offset%))
339 label130
340 (setf limit1 (f2cl-lib:int-sub kp1 ns))
341 (setf temp5 (f2cl-lib:fref alpha-%data% (ns) ((1 12)) alpha-%offset%))
342 (f2cl-lib:fdo (iq 1 (f2cl-lib:int-add iq 1))
343 ((> iq limit1) nil)
344 (tagbody
345 (setf (f2cl-lib:fref v (iq) ((1 12)))
346 (- (f2cl-lib:fref v (iq) ((1 12)))
347 (* temp5
348 (f2cl-lib:fref v
349 ((f2cl-lib:int-add iq 1))
350 ((1 12))))))
351 label135
352 (setf (f2cl-lib:fref w-%data% (iq) ((1 12)) w-%offset%)
353 (f2cl-lib:fref v (iq) ((1 12))))))
354 (setf (f2cl-lib:fref g-%data% (nsp1) ((1 13)) g-%offset%)
355 (f2cl-lib:fref w-%data% (1) ((1 12)) w-%offset%))
356 (if (= limit1 1) (go label137))
357 (setf kgi ns)
358 (setf (f2cl-lib:fref gi-%data% (kgi) ((1 11)) gi-%offset%)
359 (f2cl-lib:fref w-%data% (2) ((1 12)) w-%offset%))
360 label137
361 (setf (f2cl-lib:fref w-%data%
362 ((f2cl-lib:int-add limit1 1))
363 ((1 12))
364 w-%offset%)
365 (f2cl-lib:fref v ((f2cl-lib:int-add limit1 1)) ((1 12))))
366 (if (>= k kold) (go label140))
367 (setf ivc (f2cl-lib:int-add ivc 1))
368 (setf (f2cl-lib:fref iv-%data% (ivc) ((1 10)) iv-%offset%)
369 (f2cl-lib:int-add limit1 2))
370 label140
371 (setf nsp2 (f2cl-lib:int-add ns 2))
372 (setf kprev k)
373 (if (< kp1 nsp2) (go label199))
374 (f2cl-lib:fdo (i nsp2 (f2cl-lib:int-add i 1))
375 ((> i kp1) nil)
376 (tagbody
377 (setf limit2 (f2cl-lib:int-sub kp2 i))
378 (setf temp6
379 (f2cl-lib:fref alpha-%data%
380 ((f2cl-lib:int-sub i 1))
381 ((1 12))
382 alpha-%offset%))
383 (f2cl-lib:fdo (iq 1 (f2cl-lib:int-add iq 1))
384 ((> iq limit2) nil)
385 (tagbody
386 label145
387 (setf (f2cl-lib:fref w-%data% (iq) ((1 12)) w-%offset%)
388 (- (f2cl-lib:fref w-%data% (iq) ((1 12)) w-%offset%)
389 (* temp6
390 (f2cl-lib:fref w-%data%
391 ((f2cl-lib:int-add iq 1))
392 ((1 12))
393 w-%offset%))))))
394 label150
395 (setf (f2cl-lib:fref g-%data% (i) ((1 13)) g-%offset%)
396 (f2cl-lib:fref w-%data% (1) ((1 12)) w-%offset%))))
397 label199
398 (setf ksteps (f2cl-lib:int-add ksteps 1))
399 (if (< k nsp1) (go label215))
400 (f2cl-lib:fdo (i nsp1 (f2cl-lib:int-add i 1))
401 ((> i k) nil)
402 (tagbody
403 (setf temp1 (f2cl-lib:fref beta (i) ((1 12))))
404 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
405 ((> l neqn) nil)
406 (tagbody
407 label205
408 (setf (f2cl-lib:fref phi-%data%
409 (l i)
410 ((1 neqn) (1 16))
411 phi-%offset%)
412 (* temp1
413 (f2cl-lib:fref phi-%data%
414 (l i)
415 ((1 neqn) (1 16))
416 phi-%offset%)))))
417 label210))
418 label215
419 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
420 ((> l neqn) nil)
421 (tagbody
422 (setf (f2cl-lib:fref phi-%data%
423 (l kp2)
424 ((1 neqn) (1 16))
425 phi-%offset%)
426 (f2cl-lib:fref phi-%data%
427 (l kp1)
428 ((1 neqn) (1 16))
429 phi-%offset%))
430 (setf (f2cl-lib:fref phi-%data%
431 (l kp1)
432 ((1 neqn) (1 16))
433 phi-%offset%)
434 (coerce 0.0f0 'double-float))
435 label220
436 (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
437 (coerce 0.0f0 'double-float))))
438 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
439 ((> j k) nil)
440 (tagbody
441 (setf i (f2cl-lib:int-sub kp1 j))
442 (setf ip1 (f2cl-lib:int-add i 1))
443 (setf temp2 (f2cl-lib:fref g-%data% (i) ((1 13)) g-%offset%))
444 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
445 ((> l neqn) nil)
446 (tagbody
447 (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
448 (+ (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
449 (* temp2
450 (f2cl-lib:fref phi-%data%
451 (l i)
452 ((1 neqn) (1 16))
453 phi-%offset%))))
454 label225
455 (setf (f2cl-lib:fref phi-%data%
456 (l i)
457 ((1 neqn) (1 16))
458 phi-%offset%)
459 (+
460 (f2cl-lib:fref phi-%data%
461 (l i)
462 ((1 neqn) (1 16))
463 phi-%offset%)
464 (f2cl-lib:fref phi-%data%
465 (l ip1)
466 ((1 neqn) (1 16))
467 phi-%offset%)))))
468 label230))
469 (if nornd (go label240))
470 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
471 ((> l neqn) nil)
472 (tagbody
473 (setf tau
474 (- (* h (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%))
475 (f2cl-lib:fref phi-%data%
476 (l 15)
477 ((1 neqn) (1 16))
478 phi-%offset%)))
479 (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
480 (+ (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%) tau))
481 label235
482 (setf (f2cl-lib:fref phi-%data%
483 (l 16)
484 ((1 neqn) (1 16))
485 phi-%offset%)
486 (- (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
487 (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
488 tau))))
489 (go label250)
490 label240
491 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
492 ((> l neqn) nil)
493 (tagbody
494 label245
495 (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
496 (+ (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
497 (* h
498 (f2cl-lib:fref p-%data%
499 (l)
500 ((1 neqn))
501 p-%offset%))))))
502 label250
503 (setf xold x)
504 (setf x (+ x h))
505 (setf absh (abs h))
506 (multiple-value-bind
507 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
508 var-10 var-11 var-12 var-13)
509 (funcall f
510 x
511 p
512 yp
513 fpwa1
514 fpwa2
515 fpwa3
516 fpwa4
517 fpwa5
518 ifpwa1
519 ifpc1
520 (f2cl-lib:int-sub neqn 1)
521 ifpc2
522 par
523 ipar)
524 (declare (ignore var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
525 var-10 var-12 var-13))
526 (when var-0
527 (setf x var-0))
528 (when var-9
529 (setf ifpc1 var-9))
530 (when var-11
531 (setf ifpc2 var-11)))
532 (if (> ifpc2 0) (go end_label))
533 (setf erkm2 (coerce 0.0f0 'double-float))
534 (setf erkm1 (coerce 0.0f0 'double-float))
535 (setf erk (coerce 0.0f0 'double-float))
536 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
537 ((> l neqn) nil)
538 (tagbody
539 (setf temp3
540 (/ 1.0f0
541 (f2cl-lib:fref wt-%data% (l) ((1 neqn)) wt-%offset%)))
542 (setf temp4
543 (- (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%)
544 (f2cl-lib:fref phi-%data%
545 (l 1)
546 ((1 neqn) (1 16))
547 phi-%offset%)))
548 (f2cl-lib:arithmetic-if km2
549 (go label265)
550 (go label260)
551 (go label255))
552 label255
553 (setf erkm2
554 (+ erkm2
555 (expt
556 (*
557 (+
558 (f2cl-lib:fref phi-%data%
559 (l km1)
560 ((1 neqn) (1 16))
561 phi-%offset%)
562 temp4)
563 temp3)
564 2)))
565 label260
566 (setf erkm1
567 (+ erkm1
568 (expt
569 (*
570 (+
571 (f2cl-lib:fref phi-%data%
572 (l k)
573 ((1 neqn) (1 16))
574 phi-%offset%)
575 temp4)
576 temp3)
577 2)))
578 label265
579 (setf erk (+ erk (expt (* temp4 temp3) 2)))))
580 (f2cl-lib:arithmetic-if km2 (go label280) (go label275) (go label270))
581 label270
582 (setf erkm2
583 (* absh
584 (f2cl-lib:fref sig (km1) ((1 13)))
585 (f2cl-lib:fref gstr (km2) ((1 13)))
586 (f2cl-lib:fsqrt erkm2)))
587 label275
588 (setf erkm1
589 (* absh
590 (f2cl-lib:fref sig (k) ((1 13)))
591 (f2cl-lib:fref gstr (km1) ((1 13)))
592 (f2cl-lib:fsqrt erkm1)))
593 label280
594 (setf temp5 (* absh (f2cl-lib:fsqrt erk)))
595 (setf err
596 (* temp5
597 (- (f2cl-lib:fref g-%data% (k) ((1 13)) g-%offset%)
598 (f2cl-lib:fref g-%data% (kp1) ((1 13)) g-%offset%))))
599 (setf erk
600 (* temp5
601 (f2cl-lib:fref sig (kp1) ((1 13)))
602 (f2cl-lib:fref gstr (k) ((1 13)))))
603 (setf knew k)
604 (f2cl-lib:arithmetic-if km2 (go label299) (go label290) (go label285))
605 label285
606 (if (<= (max erkm1 erkm2) erk) (setf knew km1))
607 (go label299)
608 label290
609 (if (<= erkm1 (* 0.5f0 erk)) (setf knew km1))
610 label299
611 (if (<= err eps) (go label400))
612 (setf phase1 f2cl-lib:%false%)
613 (setf x xold)
614 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
615 ((> i k) nil)
616 (tagbody
617 (setf temp1 (/ 1.0f0 (f2cl-lib:fref beta (i) ((1 12)))))
618 (setf ip1 (f2cl-lib:int-add i 1))
619 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
620 ((> l neqn) nil)
621 (tagbody
622 label305
623 (setf (f2cl-lib:fref phi-%data%
624 (l i)
625 ((1 neqn) (1 16))
626 phi-%offset%)
627 (* temp1
628 (-
629 (f2cl-lib:fref phi-%data%
630 (l i)
631 ((1 neqn) (1 16))
632 phi-%offset%)
633 (f2cl-lib:fref phi-%data%
634 (l ip1)
635 ((1 neqn) (1 16))
636 phi-%offset%))))))
637 label310))
638 (if (< k 2) (go label320))
639 (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
640 ((> i k) nil)
641 (tagbody
642 label315
643 (setf (f2cl-lib:fref psi ((f2cl-lib:int-sub i 1)) ((1 12)))
644 (- (f2cl-lib:fref psi (i) ((1 12))) h))))
645 label320
646 (setf ifail (f2cl-lib:int-add ifail 1))
647 (setf temp2 (coerce 0.5f0 'double-float))
648 (f2cl-lib:arithmetic-if (f2cl-lib:int-sub ifail 3)
649 (go label335)
650 (go label330)
651 (go label325))
652 label325
653 (if (< p5eps (* 0.25f0 erk))
654 (setf temp2 (f2cl-lib:fsqrt (/ p5eps erk))))
655 label330
656 (setf knew 1)
657 label335
658 (setf h (* temp2 h))
659 (setf k knew)
660 (setf ns 0)
661 (if (>= (abs h) (* fouru (abs x))) (go label340))
662 (setf crash f2cl-lib:%true%)
663 (setf h (f2cl-lib:sign (* fouru (abs x)) h))
664 (setf eps (+ eps eps))
665 (go end_label)
666 label340
667 (go label100)
668 label400
669 (setf kold k)
670 (setf hold h)
671 (setf temp1 (* h (f2cl-lib:fref g-%data% (kp1) ((1 13)) g-%offset%)))
672 (if nornd (go label410))
673 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
674 ((> l neqn) nil)
675 (tagbody
676 (setf temp3 (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%))
677 (setf rho
678 (-
679 (* temp1
680 (- (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%)
681 (f2cl-lib:fref phi-%data%
682 (l 1)
683 ((1 neqn) (1 16))
684 phi-%offset%)))
685 (f2cl-lib:fref phi-%data%
686 (l 16)
687 ((1 neqn) (1 16))
688 phi-%offset%)))
689 (setf (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
690 (+ (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%) rho))
691 (setf (f2cl-lib:fref phi-%data%
692 (l 15)
693 ((1 neqn) (1 16))
694 phi-%offset%)
695 (- (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
696 (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
697 rho))
698 label405
699 (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%) temp3)))
700 (go label420)
701 label410
702 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
703 ((> l neqn) nil)
704 (tagbody
705 (setf temp3 (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%))
706 (setf (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
707 (+ (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
708 (* temp1
709 (-
710 (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%)
711 (f2cl-lib:fref phi-%data%
712 (l 1)
713 ((1 neqn) (1 16))
714 phi-%offset%)))))
715 label415
716 (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%) temp3)))
717 label420
718 (multiple-value-bind
719 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
720 var-10 var-11 var-12 var-13)
721 (funcall f
722 x
723 y
724 yp
725 fpwa1
726 fpwa2
727 fpwa3
728 fpwa4
729 fpwa5
730 ifpwa1
731 ifpc1
732 (f2cl-lib:int-sub neqn 1)
733 ifpc2
734 par
735 ipar)
736 (declare (ignore var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
737 var-10 var-12 var-13))
738 (when var-0
739 (setf x var-0))
740 (when var-9
741 (setf ifpc1 var-9))
742 (when var-11
743 (setf ifpc2 var-11)))
744 (if (> ifpc2 0) (go end_label))
745 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
746 ((> l neqn) nil)
747 (tagbody
748 (setf (f2cl-lib:fref phi-%data%
749 (l kp1)
750 ((1 neqn) (1 16))
751 phi-%offset%)
752 (- (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%)
753 (f2cl-lib:fref phi-%data%
754 (l 1)
755 ((1 neqn) (1 16))
756 phi-%offset%)))
757 label425
758 (setf (f2cl-lib:fref phi-%data%
759 (l kp2)
760 ((1 neqn) (1 16))
761 phi-%offset%)
762 (-
763 (f2cl-lib:fref phi-%data%
764 (l kp1)
765 ((1 neqn) (1 16))
766 phi-%offset%)
767 (f2cl-lib:fref phi-%data%
768 (l kp2)
769 ((1 neqn) (1 16))
770 phi-%offset%)))))
771 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
772 ((> i k) nil)
773 (tagbody
774 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
775 ((> l neqn) nil)
776 (tagbody
777 label430
778 (setf (f2cl-lib:fref phi-%data%
779 (l i)
780 ((1 neqn) (1 16))
781 phi-%offset%)
782 (+
783 (f2cl-lib:fref phi-%data%
784 (l i)
785 ((1 neqn) (1 16))
786 phi-%offset%)
787 (f2cl-lib:fref phi-%data%
788 (l kp1)
789 ((1 neqn) (1 16))
790 phi-%offset%)))))
791 label435))
792 (setf erkp1 (coerce 0.0f0 'double-float))
793 (if (or (= knew km1) (= k 12)) (setf phase1 f2cl-lib:%false%))
794 (if phase1 (go label450))
795 (if (= knew km1) (go label455))
796 (if (> kp1 ns) (go label460))
797 (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
798 ((> l neqn) nil)
799 (tagbody
800 label440
801 (setf erkp1
802 (+ erkp1
803 (expt
804 (/
805 (f2cl-lib:fref phi-%data%
806 (l kp2)
807 ((1 neqn) (1 16))
808 phi-%offset%)
809 (f2cl-lib:fref wt-%data% (l) ((1 neqn)) wt-%offset%))
810 2)))))
811 (setf erkp1
812 (* absh
813 (f2cl-lib:fref gstr (kp1) ((1 13)))
814 (f2cl-lib:fsqrt erkp1)))
815 (if (> k 1) (go label445))
816 (if (>= erkp1 (* 0.5f0 erk)) (go label460))
817 (go label450)
818 label445
819 (if (<= erkm1 (min erk erkp1)) (go label455))
820 (if (or (>= erkp1 erk) (= k 12)) (go label460))
821 label450
822 (setf k kp1)
823 (setf erk erkp1)
824 (go label460)
825 label455
826 (setf k km1)
827 (setf erk erkm1)
828 label460
829 (setf hnew (+ h h))
830 (if phase1 (go label465))
831 (if
832 (>= p5eps
833 (* erk (f2cl-lib:fref two ((f2cl-lib:int-add k 1)) ((1 13)))))
834 (go label465))
835 (setf hnew h)
836 (if (>= p5eps erk) (go label465))
837 (setf temp2
838 (coerce (the f2cl-lib:integer4 (f2cl-lib:int-add k 1))
839 'double-float))
840 (setf r (expt (/ p5eps erk) (/ 1.0f0 temp2)))
841 (setf hnew (* absh (max 0.5 (min 0.9 r))))
842 (setf hnew (f2cl-lib:sign (max hnew (* fouru (abs x))) h))
843 label465
844 (setf h hnew)
845 (go end_label)
846 end_label
847 (return
848 (values nil
849 nil
850 nil
851 x
852 h
853 eps
854 nil
855 start
856 hold
857 k
858 kold
859 crash
860 nil
861 nil
862 nil
863 nil
864 nil
865 nil
866 ksteps
867 xold
868 ivc
869 nil
870 kgi
871 nil
872 nil
873 nil
874 nil
875 nil
876 nil
877 nil
878 ifpc1
879 ifpc2
880 nil
881 nil))))))
882
883 (in-package #:cl-user)
884 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
885 (eval-when (:load-toplevel :compile-toplevel :execute)
886 (setf (gethash 'fortran-to-lisp::steps fortran-to-lisp::*f2cl-function-info*)
887 (fortran-to-lisp::make-f2cl-finfo
888 :arg-types '(t (fortran-to-lisp::integer4) (array double-float (*))
889 (double-float) (double-float) (double-float)
890 (array double-float (*)) fortran-to-lisp::logical
891 (double-float) (fortran-to-lisp::integer4)
892 (fortran-to-lisp::integer4) fortran-to-lisp::logical
893 (array double-float (*)) (array double-float (*))
894 (array double-float (*)) (array double-float (*))
895 (array double-float (*)) (array double-float (*))
896 (fortran-to-lisp::integer4) (double-float)
897 (fortran-to-lisp::integer4)
898 (array fortran-to-lisp::integer4 (*))
899 (fortran-to-lisp::integer4) (array double-float (*))
900 (array double-float (*)) (array double-float (*))
901 (array double-float (*)) (array double-float (*))
902 (array double-float (*))
903 (array fortran-to-lisp::integer4 (*))
904 (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
905 (array double-float (*))
906 (array fortran-to-lisp::integer4 (*)))
907 :return-values '(nil nil nil fortran-to-lisp::x fortran-to-lisp::h
908 fortran-to-lisp::eps nil fortran-to-lisp::start
909 fortran-to-lisp::hold fortran-to-lisp::k
910 fortran-to-lisp::kold fortran-to-lisp::crash nil
911 nil nil nil nil nil fortran-to-lisp::ksteps
912 fortran-to-lisp::xold fortran-to-lisp::ivc nil
913 fortran-to-lisp::kgi nil nil nil nil nil nil nil
914 fortran-to-lisp::ifpc1 fortran-to-lisp::ifpc2 nil
915 nil)
916 :calls '(fortran-to-lisp::d1mach))))
917
Maxima, a Computer Algebra System