?  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 / fixpdf.lisp
blob85f4efbb92de79afc867803167df2c7ca7f73a5b
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* ((limitd 1000))
21 (declare (type (f2cl-lib:integer4 1000 1000) limitd) (ignorable limitd))
22 (let ((y1sout 0.0)
23 (xold 0.0)
24 (sqnp1 0.0)
25 (sout 0.0)
26 (sb 0.0)
27 (sa 0.0)
28 (s99 0.0)
29 (s 0.0)
30 (hold 0.0)
31 (h 0.0)
32 (epst 0.0)
33 (epsstp 0.0)
34 (curtol 0.0)
35 (cursw 0.0)
36 (aold 0.0)
37 (np1 0)
38 (nfec 0)
39 (limit 0)
40 (lcode 0)
41 (ksteps 0)
42 (kold 0)
43 (kgi 0)
44 (k 0)
45 (jw 0)
46 (judy 0)
47 (j 0)
48 (ivc 0)
49 (iter 0)
50 (iflagc 0)
51 (st99 nil)
52 (crash nil)
53 (start nil)
54 (gi (make-array 11 :element-type 'double-float))
55 (g (make-array 13 :element-type 'double-float))
56 (w (make-array 12 :element-type 'double-float))
57 (alphas (make-array 12 :element-type 'double-float))
58 (iv (make-array 10 :element-type 'f2cl-lib:integer4)))
59 (declare (type (double-float) y1sout xold sqnp1 sout sb sa s99 s hold h
60 epst epsstp curtol cursw aold)
61 (type (f2cl-lib:integer4) np1 nfec limit lcode ksteps kold kgi k
62 jw judy j ivc iter iflagc)
63 (type f2cl-lib:logical st99 crash start)
64 (type (array double-float (11)) gi)
65 (type (array double-float (13)) g)
66 (type (array double-float (12)) w alphas)
67 (type (array f2cl-lib:integer4 (10)) iv))
68 (defun fixpdf
69 (n y iflag arctol eps trace$ a ndima nfe arclen yp ypold qr alpha tz
70 pivot wt phi p par ipar)
71 (declare (type (array f2cl-lib:integer4 (*)) ipar)
72 (type (array double-float (*)) par)
73 (type (array f2cl-lib:integer4 (*)) pivot)
74 (type (double-float) arclen eps arctol)
75 (type (array double-float (*)) p phi wt tz alpha qr ypold yp a
76 y)
77 (type (f2cl-lib:integer4) nfe ndima trace$ iflag n))
78 (f2cl-lib:with-multi-array-data
79 ((y double-float y-%data% y-%offset%)
80 (a double-float a-%data% a-%offset%)
81 (yp double-float yp-%data% yp-%offset%)
82 (ypold double-float ypold-%data% ypold-%offset%)
83 (qr double-float qr-%data% qr-%offset%)
84 (alpha double-float alpha-%data% alpha-%offset%)
85 (tz double-float tz-%data% tz-%offset%)
86 (wt double-float wt-%data% wt-%offset%)
87 (phi double-float phi-%data% phi-%offset%)
88 (p double-float p-%data% p-%offset%)
89 (pivot f2cl-lib:integer4 pivot-%data% pivot-%offset%)
90 (par double-float par-%data% par-%offset%)
91 (ipar f2cl-lib:integer4 ipar-%data% ipar-%offset%))
92 (prog ()
93 (declare)
94 (if (or (<= n 0) (<= eps 0.0f0)) (setf iflag 7))
95 (if (and (>= iflag -2) (<= iflag 0)) (go label10))
96 (if (= iflag 2) (go label35))
97 (if (= iflag 3) (go label30))
98 (setf iflag 7)
99 (go end_label)
100 label10
101 (setf arclen (coerce 0.0f0 'double-float))
102 (setf s (coerce 0.0f0 'double-float))
103 (if (<= arctol 0.0f0) (setf arctol (* 0.5f0 (f2cl-lib:fsqrt eps))))
104 (setf nfec 0)
105 (setf iflagc iflag)
106 (setf np1 (f2cl-lib:int-add n 1))
107 (setf sqnp1 (f2cl-lib:fsqrt (f2cl-lib:dble np1)))
108 (setf cursw (coerce 10.0f0 'double-float))
109 (setf st99 f2cl-lib:%false%)
110 (setf start f2cl-lib:%true%)
111 (setf crash f2cl-lib:%false%)
112 (setf hold (coerce 1.0f0 'double-float))
113 (setf h (coerce 0.1f0 'double-float))
114 (setf epsstp arctol)
115 (setf ksteps 0)
116 (setf (f2cl-lib:fref ypold-%data%
117 (1)
118 ((1 (f2cl-lib:int-add n 1)))
119 ypold-%offset%)
120 (coerce 1.0f0 'double-float))
121 (setf (f2cl-lib:fref yp-%data%
122 (1)
123 ((1 (f2cl-lib:int-add n 1)))
124 yp-%offset%)
125 (coerce 1.0f0 'double-float))
126 (setf (f2cl-lib:fref y-%data%
127 (1)
128 ((1 (f2cl-lib:int-add n 1)))
129 y-%offset%)
130 (coerce 0.0f0 'double-float))
131 (f2cl-lib:fdo (j 2 (f2cl-lib:int-add j 1))
132 ((> j np1) nil)
133 (tagbody
134 (setf (f2cl-lib:fref ypold-%data%
135 (j)
136 ((1 (f2cl-lib:int-add n 1)))
137 ypold-%offset%)
138 (coerce 0.0f0 'double-float))
139 (setf (f2cl-lib:fref yp-%data%
140 (j)
141 ((1 (f2cl-lib:int-add n 1)))
142 yp-%offset%)
143 (coerce 0.0f0 'double-float))
144 label20))
145 (cond
146 ((>= iflagc (f2cl-lib:int-sub 1))
147 (f2cl-lib:fdo (j 2 (f2cl-lib:int-add j 1))
148 ((> j np1) nil)
149 (tagbody
150 (setf (f2cl-lib:fref a-%data%
151 ((f2cl-lib:int-sub j 1))
152 ((1 n))
153 a-%offset%)
154 (f2cl-lib:fref y-%data%
155 (j)
156 ((1 (f2cl-lib:int-add n 1)))
157 y-%offset%))
158 label23))))
159 label30
160 (setf limit limitd)
161 label35
162 (f2cl-lib:fdo (iter 1 (f2cl-lib:int-add iter 1))
163 ((> iter limit) nil)
164 (tagbody
165 (cond
166 ((< (f2cl-lib:fref y (1) ((1 (f2cl-lib:int-add n 1)))) 0.0f0)
167 (tagbody
168 label40
169 (setf arclen (+ arclen s))
170 (setf iflag 5)
171 (go end_label))))
172 label50
173 (if (<= s (* 7.0f0 sqnp1)) (go label80))
174 (setf arclen (+ arclen s))
175 (setf s (coerce 0.0f0 'double-float))
176 label60
177 (setf start f2cl-lib:%true%)
178 (setf crash f2cl-lib:%false%)
179 (cond
180 ((= iflagc (f2cl-lib:int-sub 2))
181 (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
182 ((> jw ndima) nil)
183 (tagbody
184 (setf (f2cl-lib:fref qr-%data%
185 (jw 1)
186 ((1 n) (1 (f2cl-lib:int-add n 1)))
187 qr-%offset%)
188 (f2cl-lib:fref a-%data% (jw) ((1 n)) a-%offset%))
189 label63))
190 (rhoa a
191 (f2cl-lib:fref y-%data%
192 (1)
193 ((1 (f2cl-lib:int-add n 1)))
194 y-%offset%)
195 (f2cl-lib:array-slice y-%data%
196 double-float
197 (2)
198 ((1 (f2cl-lib:int-add n 1)))
199 y-%offset%)
200 par ipar)
201 (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
202 ((> jw ndima) nil)
203 (tagbody
204 (setf aold
205 (f2cl-lib:fref qr-%data%
206 (jw 1)
207 ((1 n) (1 (f2cl-lib:int-add n 1)))
208 qr-%offset%))
209 (cond
210 ((> (abs (+ (f2cl-lib:fref a (jw) ((1 n))) (- aold)))
211 (+ 1.0f0 (abs aold)))
212 (setf arclen (+ arclen s))
213 (setf iflag 5)
214 (go end_label)))
215 label65)))
216 (t
217 (f
218 (f2cl-lib:array-slice y-%data%
219 double-float
220 (2)
221 ((1 (f2cl-lib:int-add n 1)))
222 y-%offset%)
223 yp)
224 (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
225 ((> jw n) nil)
226 (tagbody
227 (setf aold
228 (f2cl-lib:fref a-%data% (jw) ((1 n)) a-%offset%))
229 (cond
230 ((= iflagc (f2cl-lib:int-sub 1))
231 (setf (f2cl-lib:fref a-%data% (jw) ((1 n)) a-%offset%)
232 (+
233 (/
234 (*
235 (f2cl-lib:fref y-%data%
236 (1)
237 ((1 (f2cl-lib:int-add n 1)))
238 y-%offset%)
239 (f2cl-lib:fref yp-%data%
240 (jw)
241 ((1 (f2cl-lib:int-add n 1)))
242 yp-%offset%))
243 (- 1.0f0
244 (f2cl-lib:fref y-%data%
245 (1)
246 ((1
247 (f2cl-lib:int-add n 1)))
248 y-%offset%)))
249 (f2cl-lib:fref y-%data%
250 ((f2cl-lib:int-add jw 1))
251 ((1 (f2cl-lib:int-add n 1)))
252 y-%offset%))))
253 (t
254 (setf (f2cl-lib:fref a-%data% (jw) ((1 n)) a-%offset%)
255 (/
256 (-
257 (f2cl-lib:fref y-%data%
258 ((f2cl-lib:int-add jw 1))
259 ((1 (f2cl-lib:int-add n 1)))
260 y-%offset%)
261 (*
262 (f2cl-lib:fref y-%data%
263 (1)
264 ((1 (f2cl-lib:int-add n 1)))
265 y-%offset%)
266 (f2cl-lib:fref yp-%data%
267 (jw)
268 ((1 (f2cl-lib:int-add n 1)))
269 yp-%offset%)))
270 (- 1.0f0
271 (f2cl-lib:fref y-%data%
272 (1)
273 ((1 (f2cl-lib:int-add n 1)))
274 y-%offset%))))))
275 (cond
276 ((> (abs (+ (f2cl-lib:fref a (jw) ((1 n))) (- aold)))
277 (+ 1.0f0 (abs aold)))
278 (setf arclen (+ arclen s))
279 (setf iflag 5)
280 (go end_label)))
281 label70))))
282 (go label100)
283 label80
284 (if
285 (or
286 (<=
287 (f2cl-lib:fref y-%data%
288 (1)
289 ((1 (f2cl-lib:int-add n 1)))
290 y-%offset%)
291 0.99f0)
292 st99)
293 (go label100))
294 label90
295 (setf st99 f2cl-lib:%true%)
296 (setf epsstp eps)
297 (setf arctol eps)
298 (go label60)
299 label100
300 (setf curtol (* cursw hold))
301 (setf epst (/ eps epsstp))
302 (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
303 ((> jw np1) nil)
304 (tagbody
305 (cond
306 ((<=
307 (abs
308 (+ (f2cl-lib:fref yp (jw) ((1 (f2cl-lib:int-add n 1))))
309 (-
310 (f2cl-lib:fref ypold
311 (jw)
312 ((1 (f2cl-lib:int-add n 1)))))))
313 curtol)
314 (setf (f2cl-lib:fref wt-%data%
315 (jw)
316 ((1 (f2cl-lib:int-add n 1)))
317 wt-%offset%)
318 (+
319 (abs
320 (f2cl-lib:fref y-%data%
321 (jw)
322 ((1 (f2cl-lib:int-add n 1)))
323 y-%offset%))
324 1.0f0)))
325 (t
326 (setf (f2cl-lib:fref wt-%data%
327 (jw)
328 ((1 (f2cl-lib:int-add n 1)))
329 wt-%offset%)
330 (*
331 (+
332 (abs
333 (f2cl-lib:fref y-%data%
334 (jw)
335 ((1 (f2cl-lib:int-add n 1)))
336 y-%offset%))
337 1.0f0)
338 epst))))
339 label110))
340 (multiple-value-bind
341 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
342 var-9 var-10 var-11 var-12 var-13 var-14 var-15 var-16
343 var-17 var-18 var-19 var-20 var-21 var-22 var-23 var-24
344 var-25 var-26 var-27 var-28 var-29 var-30 var-31 var-32
345 var-33)
346 (steps #'fode np1 y s h epsstp wt start hold k kold crash phi
347 p yp alphas w g ksteps xold ivc iv kgi gi ypold a qr alpha
348 tz pivot nfec iflagc par ipar)
349 (declare (ignore var-0 var-1 var-2 var-6 var-12 var-13 var-14
350 var-15 var-16 var-17 var-21 var-23 var-24
351 var-25 var-26 var-27 var-28 var-29 var-32
352 var-33))
353 (setf s var-3)
354 (setf h var-4)
355 (setf epsstp var-5)
356 (setf start var-7)
357 (setf hold var-8)
358 (setf k var-9)
359 (setf kold var-10)
360 (setf crash var-11)
361 (setf ksteps var-18)
362 (setf xold var-19)
363 (setf ivc var-20)
364 (setf kgi var-22)
365 (setf nfec var-30)
366 (setf iflagc var-31))
367 (cond
368 ((> trace$ 0)
369 (f2cl-lib:fformat trace$
370 ("~%" " STEP" 1 (("~5D")) "~3@T" "NFE =" 1
371 (("~5D")) "~3@T" "ARC LENGTH =" 1
372 (("~9,4,0,'*,F")) "~3@T" "LAMBDA =" 1
373 (("~7,4,0,'*,F")) "~5@T" "X vector:" "~%" t
374 ("~1@T" 6 (("~12,4,2,0,'*,,'EE"))) "~%")
375 iter
376 nfec
377 s
378 (f2cl-lib:fref y-%data%
379 (1)
380 ((1 (f2cl-lib:int-add n 1)))
381 y-%offset%)
382 (do ((jw 2 (f2cl-lib:int-add jw 1))
383 (%ret nil))
384 ((> jw np1) (nreverse %ret))
385 (declare (type f2cl-lib:integer4 jw))
386 (push
387 (f2cl-lib:fref y-%data%
388 (jw)
389 ((1
390 (f2cl-lib:int-add n 1)))
391 y-%offset%)
392 %ret)))))
393 (setf nfe nfec)
394 (cond
395 ((= iflagc 4)
396 (setf arclen (+ arclen s))
397 (setf iflag 4)
398 (go end_label)))
399 label120
400 (cond
401 (crash
402 (setf iflag 2)
403 (setf eps epsstp)
404 (if (< arctol eps) (setf arctol eps))
405 (setf limit (f2cl-lib:int-sub limit iter))
406 (go end_label)))
407 label130
408 (cond
409 ((>= (f2cl-lib:fref y (1) ((1 (f2cl-lib:int-add n 1)))) 1.0f0)
410 (tagbody
411 (if st99 (go label160))
412 (setf s99 (- s (* 0.5f0 hold)))
413 label135
414 (sintrp s y s99 wt yp np1 kold phi ivc iv kgi gi alphas g w
415 xold p)
416 (if
417 (<
418 (f2cl-lib:fref wt-%data%
419 (1)
420 ((1 (f2cl-lib:int-add n 1)))
421 wt-%offset%)
422 1.0f0)
423 (go label140))
424 (setf s99 (* 0.5f0 (+ (- s hold) s99)))
425 (go label135)
426 label140
427 (f2cl-lib:fdo (judy 1 (f2cl-lib:int-add judy 1))
428 ((> judy np1) nil)
429 (tagbody
430 (setf (f2cl-lib:fref y-%data%
431 (judy)
432 ((1 (f2cl-lib:int-add n 1)))
433 y-%offset%)
434 (f2cl-lib:fref wt-%data%
435 (judy)
436 ((1 (f2cl-lib:int-add n 1)))
437 wt-%offset%))
438 (setf (f2cl-lib:fref ypold-%data%
439 (judy)
440 ((1 (f2cl-lib:int-add n 1)))
441 ypold-%offset%)
442 (f2cl-lib:fref yp-%data%
443 (judy)
444 ((1 (f2cl-lib:int-add n 1)))
445 yp-%offset%))
446 label144))
447 (setf s s99)
448 (go label90))))
449 label150))
450 (setf iflag 3)
451 (go end_label)
452 label160
453 (setf sa (- s hold))
454 (setf sb s)
455 (setf lcode 1)
456 label170
457 (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
458 (root sout y1sout sa sb eps eps lcode)
459 (declare (ignore var-1 var-4 var-5))
460 (setf sout var-0)
461 (setf sa var-2)
462 (setf sb var-3)
463 (setf lcode var-6))
464 (if (> lcode 0) (go label190))
465 (sintrp s y sout wt yp np1 kold phi ivc iv kgi gi alphas g w xold p)
466 (setf y1sout
467 (-
468 (f2cl-lib:fref wt-%data%
469 (1)
470 ((1 (f2cl-lib:int-add n 1)))
471 wt-%offset%)
472 1.0f0))
473 (go label170)
474 label190
475 (setf iflag 1)
476 (if (> lcode 2) (setf iflag 6))
477 (setf arclen (+ arclen sa))
478 (sintrp s y sa wt yp np1 kold phi ivc iv kgi gi alphas g w xold p)
479 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
480 ((> j np1) nil)
481 (tagbody
482 label210
483 (setf (f2cl-lib:fref y-%data%
484 (j)
485 ((1 (f2cl-lib:int-add n 1)))
486 y-%offset%)
487 (f2cl-lib:fref wt-%data%
488 (j)
489 ((1 (f2cl-lib:int-add n 1)))
490 wt-%offset%))))
491 (go end_label)
492 end_label
493 (return
494 (values nil
495 nil
496 iflag
497 arctol
498 eps
499 nil
500 nil
501 nil
502 nfe
503 arclen
504 nil
505 nil
506 nil
507 nil
508 nil
509 nil
510 nil
511 nil
512 nil
513 nil
514 nil)))))))
515
516 (in-package #:cl-user)
517 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
518 (eval-when (:load-toplevel :compile-toplevel :execute)
519 (setf (gethash 'fortran-to-lisp::fixpdf
520 fortran-to-lisp::*f2cl-function-info*)
521 (fortran-to-lisp::make-f2cl-finfo
522 :arg-types '((fortran-to-lisp::integer4) (array double-float (*))
523 (fortran-to-lisp::integer4) (double-float)
524 (double-float) (fortran-to-lisp::integer4)
525 (array double-float (*)) (fortran-to-lisp::integer4)
526 (fortran-to-lisp::integer4) (double-float)
527 (array double-float (*)) (array double-float (*))
528 (array double-float (*)) (array double-float (*))
529 (array double-float (*))
530 (array fortran-to-lisp::integer4 (*))
531 (array double-float (*)) (array double-float (*))
532 (array double-float (*)) (array double-float (*))
533 (array fortran-to-lisp::integer4 (*)))
534 :return-values '(nil nil fortran-to-lisp::iflag
535 fortran-to-lisp::arctol fortran-to-lisp::eps nil
536 nil nil fortran-to-lisp::nfe
537 fortran-to-lisp::arclen nil nil nil nil nil nil nil
538 nil nil nil nil)
539 :calls '(fortran-to-lisp::f fortran-to-lisp::rhoa
540 fortran-to-lisp::root fortran-to-lisp::sintrp
541 fortran-to-lisp::steps))))
542
Maxima, a Computer Algebra System