?  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 / pcgds.lisp
blobd135bcc502314e5c9a4cee9c6fdb47eb058a0556
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 (defun pcgds (nn aa lenaa maxa pp start work iflag)
21 (declare (type (array f2cl-lib:integer4 (*)) maxa)
22 (type (array double-float (*)) work start pp aa)
23 (type (f2cl-lib:integer4) iflag lenaa nn))
24 (f2cl-lib:with-multi-array-data
25 ((aa double-float aa-%data% aa-%offset%)
26 (pp double-float pp-%data% pp-%offset%)
27 (start double-float start-%data% start-%offset%)
28 (work double-float work-%data% work-%offset%)
29 (maxa f2cl-lib:integer4 maxa-%data% maxa-%offset%))
30 (prog ((stillu nil) (stillb nil) (ab 0.0) (au 0.0) (bb 0.0) (bu 0.0)
31 (dznrm 0.0) (pbnprd 0.0) (punprd 0.0) (rbnprd 0.0) (rbtol 0.0)
32 (rnprd 0.0) (runprd 0.0) (rutol 0.0) (startk 0.0) (temp 0.0)
33 (unrm 0.0) (zlen 0.0) (ztol 0.0) (imax 0) (ind 0) (j 0) (k 0)
34 (np1 0) (np2 0) (n2p3 0) (n3p4 0) (n4p5 0) (n5p6 0))
35 (declare (type (f2cl-lib:integer4) n5p6 n4p5 n3p4 n2p3 np2 np1 k j ind
36 imax)
37 (type (double-float) ztol zlen unrm temp startk rutol runprd
38 rnprd rbtol rbnprd punprd pbnprd dznrm bu
39 bb au ab)
40 (type f2cl-lib:logical stillb stillu))
41 (setf np1 (f2cl-lib:int-add nn 1))
42 (setf np2 (f2cl-lib:int-add nn 2))
43 (setf n2p3 (f2cl-lib:int-add (f2cl-lib:int-mul 2 nn) 3))
44 (setf n3p4 (f2cl-lib:int-add (f2cl-lib:int-mul 3 nn) 4))
45 (setf n4p5 (f2cl-lib:int-add (f2cl-lib:int-mul 4 nn) 5))
46 (setf n5p6 (f2cl-lib:int-add (f2cl-lib:int-mul 5 nn) 6))
47 (setf k (idamax np1 start 1))
48 (setf startk
49 (f2cl-lib:fref start-%data%
50 (k)
51 ((1 (f2cl-lib:int-add nn 1)))
52 start-%offset%))
53 (dcopy lenaa aa 1
54 (f2cl-lib:array-slice work-%data%
55 double-float
56 (n5p6)
57 ((1
58 (f2cl-lib:int-add
59 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
60 lenaa
61 nn
62 1)))
63 work-%offset%)
64 1)
65 (setf (f2cl-lib:fref maxa-%data%
66 ((f2cl-lib:int-add nn 1))
67 ((1 (f2cl-lib:int-add nn 2)))
68 maxa-%offset%)
69 (f2cl-lib:int-add lenaa 1))
70 (setf (f2cl-lib:fref maxa-%data%
71 ((f2cl-lib:int-add nn 2))
72 ((1 (f2cl-lib:int-add nn 2)))
73 maxa-%offset%)
74 (f2cl-lib:int-sub (f2cl-lib:int-add lenaa nn 3) k))
75 (mfacds nn
76 (f2cl-lib:array-slice work-%data%
77 double-float
78 (n5p6)
79 ((1
80 (f2cl-lib:int-add
81 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
82 lenaa
83 nn
84 1)))
85 work-%offset%)
86 lenaa maxa)
87 (dcopy nn pp 1 work 1)
88 (if (< k np1)
89 (setf (f2cl-lib:fref work-%data%
90 (k)
91 ((1
92 (f2cl-lib:int-add
93 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
94 lenaa
95 nn
96 1)))
97 work-%offset%)
98 (+
99 (f2cl-lib:fref work-%data%
100 (k)
101 ((1
102 (f2cl-lib:int-add
103 (f2cl-lib:int-mul 5
104 (f2cl-lib:int-add nn 1))
105 lenaa
106 nn
107 1)))
108 work-%offset%)
109 1.0)))
110 (setf unrm (dnrm2 nn work 1))
111 (setf imax (f2cl-lib:int-mul 10 np1))
112 (setf stillu f2cl-lib:%true%)
113 (setf stillb f2cl-lib:%true%)
114 (setf ztol (* 100.0f0 (f2cl-lib:d1mach 4)))
115 (setf rbtol (* ztol (abs startk)))
116 (setf rutol (* ztol unrm))
117 (multds
118 (f2cl-lib:array-slice work-%data%
119 double-float
120 (n3p4)
121 ((1
122 (f2cl-lib:int-add
123 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
124 lenaa
125 nn
126 1)))
127 work-%offset%)
128 aa
129 (f2cl-lib:array-slice work-%data%
130 double-float
131 (np2)
132 ((1
133 (f2cl-lib:int-add
134 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
135 lenaa
136 nn
137 1)))
138 work-%offset%)
139 maxa nn lenaa)
140 (setf (f2cl-lib:fref work-%data%
141 ((f2cl-lib:int-add n3p4 nn))
142 ((1
143 (f2cl-lib:int-add
144 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
145 lenaa
146 nn
147 1)))
148 work-%offset%)
149 (f2cl-lib:fref work-%data%
150 ((f2cl-lib:int-sub (f2cl-lib:int-add np2 k) 1))
151 ((1
152 (f2cl-lib:int-add
153 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
154 lenaa
155 nn
156 1)))
157 work-%offset%))
158 (setf ind (f2cl-lib:int-sub (f2cl-lib:int-add n3p4 k) 1))
159 (if (< k np1)
160 (setf (f2cl-lib:fref work-%data%
161 (ind)
162 ((1
163 (f2cl-lib:int-add
164 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
165 lenaa
166 nn
167 1)))
168 work-%offset%)
169 (+
170 (f2cl-lib:fref work-%data%
171 (ind)
172 ((1
173 (f2cl-lib:int-add
174 (f2cl-lib:int-mul 5
175 (f2cl-lib:int-add nn 1))
176 lenaa
177 nn
178 1)))
179 work-%offset%)
180 (f2cl-lib:fref work-%data%
181 ((f2cl-lib:int-add np2 nn))
182 ((1
183 (f2cl-lib:int-add
184 (f2cl-lib:int-mul 5
185 (f2cl-lib:int-add nn 1))
186 lenaa
187 nn
188 1)))
189 work-%offset%))))
190 (dscal np1 -1.0
191 (f2cl-lib:array-slice work-%data%
192 double-float
193 (n3p4)
194 ((1
195 (f2cl-lib:int-add
196 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
197 lenaa
198 nn
199 1)))
200 work-%offset%)
201 1)
202 (daxpy nn -1.0 pp 1
203 (f2cl-lib:array-slice work-%data%
204 double-float
205 (n3p4)
206 ((1
207 (f2cl-lib:int-add
208 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
209 lenaa
210 nn
211 1)))
212 work-%offset%)
213 1)
214 (if (< k np1)
215 (setf (f2cl-lib:fref work-%data%
216 (ind)
217 ((1
218 (f2cl-lib:int-add
219 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
220 lenaa
221 nn
222 1)))
223 work-%offset%)
224 (-
225 (f2cl-lib:fref work-%data%
226 (ind)
227 ((1
228 (f2cl-lib:int-add
229 (f2cl-lib:int-mul 5
230 (f2cl-lib:int-add nn 1))
231 lenaa
232 nn
233 1)))
234 work-%offset%)
235 1.0)))
236 (qimuds
237 (f2cl-lib:array-slice work-%data%
238 double-float
239 (n5p6)
240 ((1
241 (f2cl-lib:int-add
242 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
243 lenaa
244 nn
245 1)))
246 work-%offset%)
247 (f2cl-lib:array-slice work-%data%
248 double-float
249 (n3p4)
250 ((1
251 (f2cl-lib:int-add
252 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
253 lenaa
254 nn
255 1)))
256 work-%offset%)
257 maxa nn lenaa)
258 (dcopy np1
259 (f2cl-lib:array-slice work-%data%
260 double-float
261 (n3p4)
262 ((1
263 (f2cl-lib:int-add
264 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
265 lenaa
266 nn
267 1)))
268 work-%offset%)
269 1 work 1)
270 (qimuds
271 (f2cl-lib:array-slice work-%data%
272 double-float
273 (n5p6)
274 ((1
275 (f2cl-lib:int-add
276 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
277 lenaa
278 nn
279 1)))
280 work-%offset%)
281 work maxa nn lenaa)
282 (multds
283 (f2cl-lib:array-slice work-%data%
284 double-float
285 (n4p5)
286 ((1
287 (f2cl-lib:int-add
288 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
289 lenaa
290 nn
291 1)))
292 work-%offset%)
293 aa work maxa nn lenaa)
294 (setf (f2cl-lib:fref work-%data%
295 ((f2cl-lib:int-add n4p5 nn))
296 ((1
297 (f2cl-lib:int-add
298 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
299 lenaa
300 nn
301 1)))
302 work-%offset%)
303 (f2cl-lib:fref work-%data%
304 (k)
305 ((1
306 (f2cl-lib:int-add
307 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
308 lenaa
309 nn
310 1)))
311 work-%offset%))
312 (if (< k np1)
313 (setf (f2cl-lib:fref work-%data%
314 ((f2cl-lib:int-sub (f2cl-lib:int-add n4p5 k) 1))
315 ((1
316 (f2cl-lib:int-add
317 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
318 lenaa
319 nn
320 1)))
321 work-%offset%)
322 (+
323 (f2cl-lib:fref work-%data%
324 ((f2cl-lib:int-sub (f2cl-lib:int-add n4p5 k)
325 1))
326 ((1
327 (f2cl-lib:int-add
328 (f2cl-lib:int-mul 5
329 (f2cl-lib:int-add nn 1))
330 lenaa
331 nn
332 1)))
333 work-%offset%)
334 (f2cl-lib:fref work-%data%
335 (np1)
336 ((1
337 (f2cl-lib:int-add
338 (f2cl-lib:int-mul 5
339 (f2cl-lib:int-add nn 1))
340 lenaa
341 nn
342 1)))
343 work-%offset%))))
344 (setf runprd
345 (ddot np1
346 (f2cl-lib:array-slice work-%data%
347 double-float
348 (n3p4)
349 ((1
350 (f2cl-lib:int-add
351 (f2cl-lib:int-mul 5
352 (f2cl-lib:int-add nn
353 1))
354 lenaa
355 nn
356 1)))
357 work-%offset%)
358 1
359 (f2cl-lib:array-slice work-%data%
360 double-float
361 (n3p4)
362 ((1
363 (f2cl-lib:int-add
364 (f2cl-lib:int-mul 5
365 (f2cl-lib:int-add nn
366 1))
367 lenaa
368 nn
369 1)))
370 work-%offset%)
371 1))
372 (setf punprd
373 (ddot np1
374 (f2cl-lib:array-slice work-%data%
375 double-float
376 (n4p5)
377 ((1
378 (f2cl-lib:int-add
379 (f2cl-lib:int-mul 5
380 (f2cl-lib:int-add nn
381 1))
382 lenaa
383 nn
384 1)))
385 work-%offset%)
386 1
387 (f2cl-lib:array-slice work-%data%
388 double-float
389 (n4p5)
390 ((1
391 (f2cl-lib:int-add
392 (f2cl-lib:int-mul 5
393 (f2cl-lib:int-add nn
394 1))
395 lenaa
396 nn
397 1)))
398 work-%offset%)
399 1))
400 (setf j 1)
401 label100
402 (if (not (and stillu (<= j imax))) (go label200))
403 (cond
404 ((> (f2cl-lib:fsqrt runprd) rutol)
405 (cond
406 ((= punprd 0.0f0)
407 (multds
408 (f2cl-lib:array-slice work-%data%
409 double-float
410 (n3p4)
411 ((1
412 (f2cl-lib:int-add
413 (f2cl-lib:int-mul 5
414 (f2cl-lib:int-add nn
415 1))
416 lenaa
417 nn
418 1)))
419 work-%offset%)
420 aa
421 (f2cl-lib:array-slice work-%data%
422 double-float
423 (np2)
424 ((1
425 (f2cl-lib:int-add
426 (f2cl-lib:int-mul 5
427 (f2cl-lib:int-add nn
428 1))
429 lenaa
430 nn
431 1)))
432 work-%offset%)
433 maxa nn lenaa)
434 (setf (f2cl-lib:fref work-%data%
435 ((f2cl-lib:int-add n3p4 nn))
436 ((1
437 (f2cl-lib:int-add
438 (f2cl-lib:int-mul 5
439 (f2cl-lib:int-add nn 1))
440 lenaa
441 nn
442 1)))
443 work-%offset%)
444 (f2cl-lib:fref work-%data%
445 ((f2cl-lib:int-sub (f2cl-lib:int-add np2 k)
446 1))
447 ((1
448 (f2cl-lib:int-add
449 (f2cl-lib:int-mul 5
450 (f2cl-lib:int-add nn
451 1))
452 lenaa
453 nn
454 1)))
455 work-%offset%))
456 (setf ind (f2cl-lib:int-sub (f2cl-lib:int-add n3p4 k) 1))
457 (if (< k np1)
458 (setf (f2cl-lib:fref work-%data%
459 (ind)
460 ((1
461 (f2cl-lib:int-add
462 (f2cl-lib:int-mul 5
463 (f2cl-lib:int-add nn
464 1))
465 lenaa
466 nn
467 1)))
468 work-%offset%)
469 (+
470 (f2cl-lib:fref work-%data%
471 (ind)
472 ((1
473 (f2cl-lib:int-add
474 (f2cl-lib:int-mul 5
475 (f2cl-lib:int-add
476 nn
477 1))
478 lenaa
479 nn
480 1)))
481 work-%offset%)
482 (f2cl-lib:fref work-%data%
483 ((f2cl-lib:int-add np2 nn))
484 ((1
485 (f2cl-lib:int-add
486 (f2cl-lib:int-mul 5
487 (f2cl-lib:int-add
488 nn
489 1))
490 lenaa
491 nn
492 1)))
493 work-%offset%))))
494 (dscal np1 -1.0
495 (f2cl-lib:array-slice work-%data%
496 double-float
497 (n3p4)
498 ((1
499 (f2cl-lib:int-add
500 (f2cl-lib:int-mul 5
501 (f2cl-lib:int-add nn
502 1))
503 lenaa
504 nn
505 1)))
506 work-%offset%)
507 1)
508 (daxpy nn -1.0 pp 1
509 (f2cl-lib:array-slice work-%data%
510 double-float
511 (n3p4)
512 ((1
513 (f2cl-lib:int-add
514 (f2cl-lib:int-mul 5
515 (f2cl-lib:int-add nn
516 1))
517 lenaa
518 nn
519 1)))
520 work-%offset%)
521 1)
522 (if (< k np1)
523 (setf (f2cl-lib:fref work-%data%
524 ((f2cl-lib:int-sub
525 (f2cl-lib:int-add n3p4 k)
526 1))
527 ((1
528 (f2cl-lib:int-add
529 (f2cl-lib:int-mul 5
530 (f2cl-lib:int-add nn
531 1))
532 lenaa
533 nn
534 1)))
535 work-%offset%)
536 (-
537 (f2cl-lib:fref work-%data%
538 ((f2cl-lib:int-sub
539 (f2cl-lib:int-add n3p4 k)
540 1))
541 ((1
542 (f2cl-lib:int-add
543 (f2cl-lib:int-mul 5
544 (f2cl-lib:int-add
545 nn
546 1))
547 lenaa
548 nn
549 1)))
550 work-%offset%)
551 1.0)))
552 (qimuds
553 (f2cl-lib:array-slice work-%data%
554 double-float
555 (n5p6)
556 ((1
557 (f2cl-lib:int-add
558 (f2cl-lib:int-mul 5
559 (f2cl-lib:int-add nn
560 1))
561 lenaa
562 nn
563 1)))
564 work-%offset%)
565 (f2cl-lib:array-slice work-%data%
566 double-float
567 (n3p4)
568 ((1
569 (f2cl-lib:int-add
570 (f2cl-lib:int-mul 5
571 (f2cl-lib:int-add nn
572 1))
573 lenaa
574 nn
575 1)))
576 work-%offset%)
577 maxa nn lenaa)
578 (dcopy np1
579 (f2cl-lib:array-slice work-%data%
580 double-float
581 (n3p4)
582 ((1
583 (f2cl-lib:int-add
584 (f2cl-lib:int-mul 5
585 (f2cl-lib:int-add nn
586 1))
587 lenaa
588 nn
589 1)))
590 work-%offset%)
591 1 work 1)
592 (qimuds
593 (f2cl-lib:array-slice work-%data%
594 double-float
595 (n5p6)
596 ((1
597 (f2cl-lib:int-add
598 (f2cl-lib:int-mul 5
599 (f2cl-lib:int-add nn
600 1))
601 lenaa
602 nn
603 1)))
604 work-%offset%)
605 work maxa nn lenaa)
606 (multds
607 (f2cl-lib:array-slice work-%data%
608 double-float
609 (n4p5)
610 ((1
611 (f2cl-lib:int-add
612 (f2cl-lib:int-mul 5
613 (f2cl-lib:int-add nn
614 1))
615 lenaa
616 nn
617 1)))
618 work-%offset%)
619 aa work maxa nn lenaa)
620 (setf (f2cl-lib:fref work-%data%
621 ((f2cl-lib:int-add n4p5 nn))
622 ((1
623 (f2cl-lib:int-add
624 (f2cl-lib:int-mul 5
625 (f2cl-lib:int-add nn 1))
626 lenaa
627 nn
628 1)))
629 work-%offset%)
630 (f2cl-lib:fref work-%data%
631 (k)
632 ((1
633 (f2cl-lib:int-add
634 (f2cl-lib:int-mul 5
635 (f2cl-lib:int-add nn
636 1))
637 lenaa
638 nn
639 1)))
640 work-%offset%))
641 (setf ind (f2cl-lib:int-sub (f2cl-lib:int-add n4p5 k) 1))
642 (if (< k np1)
643 (setf (f2cl-lib:fref work-%data%
644 (ind)
645 ((1
646 (f2cl-lib:int-add
647 (f2cl-lib:int-mul 5
648 (f2cl-lib:int-add nn
649 1))
650 lenaa
651 nn
652 1)))
653 work-%offset%)
654 (+
655 (f2cl-lib:fref work-%data%
656 (ind)
657 ((1
658 (f2cl-lib:int-add
659 (f2cl-lib:int-mul 5
660 (f2cl-lib:int-add
661 nn
662 1))
663 lenaa
664 nn
665 1)))
666 work-%offset%)
667 (f2cl-lib:fref work-%data%
668 (np1)
669 ((1
670 (f2cl-lib:int-add
671 (f2cl-lib:int-mul 5
672 (f2cl-lib:int-add
673 nn
674 1))
675 lenaa
676 nn
677 1)))
678 work-%offset%))))
679 (setf runprd
680 (ddot np1
681 (f2cl-lib:array-slice work-%data%
682 double-float
683 (n3p4)
684 ((1
685 (f2cl-lib:int-add
686 (f2cl-lib:int-mul 5
687 (f2cl-lib:int-add
688 nn
689 1))
690 lenaa
691 nn
692 1)))
693 work-%offset%)
694 1
695 (f2cl-lib:array-slice work-%data%
696 double-float
697 (n3p4)
698 ((1
699 (f2cl-lib:int-add
700 (f2cl-lib:int-mul 5
701 (f2cl-lib:int-add
702 nn
703 1))
704 lenaa
705 nn
706 1)))
707 work-%offset%)
708 1))
709 (setf punprd
710 (ddot np1
711 (f2cl-lib:array-slice work-%data%
712 double-float
713 (n4p5)
714 ((1
715 (f2cl-lib:int-add
716 (f2cl-lib:int-mul 5
717 (f2cl-lib:int-add
718 nn
719 1))
720 lenaa
721 nn
722 1)))
723 work-%offset%)
724 1
725 (f2cl-lib:array-slice work-%data%
726 double-float
727 (n4p5)
728 ((1
729 (f2cl-lib:int-add
730 (f2cl-lib:int-mul 5
731 (f2cl-lib:int-add
732 nn
733 1))
734 lenaa
735 nn
736 1)))
737 work-%offset%)
738 1))
739 (cond
740 ((<= (f2cl-lib:fsqrt runprd) rutol)
741 (setf stillu f2cl-lib:%false%)))))
742 (cond
743 (stillu
744 (setf au (/ runprd punprd))
745 (dcopy np1
746 (f2cl-lib:array-slice work-%data%
747 double-float
748 (np2)
749 ((1
750 (f2cl-lib:int-add
751 (f2cl-lib:int-mul 5
752 (f2cl-lib:int-add nn
753 1))
754 lenaa
755 nn
756 1)))
757 work-%offset%)
758 1 work 1)
759 (daxpy np1 au
760 (f2cl-lib:array-slice work-%data%
761 double-float
762 (n4p5)
763 ((1
764 (f2cl-lib:int-add
765 (f2cl-lib:int-mul 5
766 (f2cl-lib:int-add nn
767 1))
768 lenaa
769 nn
770 1)))
771 work-%offset%)
772 1
773 (f2cl-lib:array-slice work-%data%
774 double-float
775 (np2)
776 ((1
777 (f2cl-lib:int-add
778 (f2cl-lib:int-mul 5
779 (f2cl-lib:int-add nn
780 1))
781 lenaa
782 nn
783 1)))
784 work-%offset%)
785 1)
786 (daxpy np1 -1.0
787 (f2cl-lib:array-slice work-%data%
788 double-float
789 (np2)
790 ((1
791 (f2cl-lib:int-add
792 (f2cl-lib:int-mul 5
793 (f2cl-lib:int-add nn
794 1))
795 lenaa
796 nn
797 1)))
798 work-%offset%)
799 1 work 1)
800 (setf zlen
801 (dnrm2 np1
802 (f2cl-lib:array-slice work-%data%
803 double-float
804 (np2)
805 ((1
806 (f2cl-lib:int-add
807 (f2cl-lib:int-mul 5
808 (f2cl-lib:int-add
809 nn
810 1))
811 lenaa
812 nn
813 1)))
814 work-%offset%)
815 1))
816 (setf dznrm (dnrm2 np1 work 1))
817 (if (< (/ dznrm zlen) ztol) (setf stillu f2cl-lib:%false%)))))
818 (t
819 (setf stillu f2cl-lib:%false%)))
820 (cond
821 (stillu
822 (multds work aa
823 (f2cl-lib:array-slice work-%data%
824 double-float
825 (n4p5)
826 ((1
827 (f2cl-lib:int-add
828 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
829 lenaa
830 nn
831 1)))
832 work-%offset%)
833 maxa nn lenaa)
834 (setf (f2cl-lib:fref work-%data%
835 (np1)
836 ((1
837 (f2cl-lib:int-add
838 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
839 lenaa
840 nn
841 1)))
842 work-%offset%)
843 (f2cl-lib:fref work-%data%
844 ((f2cl-lib:int-sub (f2cl-lib:int-add n4p5 k)
845 1))
846 ((1
847 (f2cl-lib:int-add
848 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
849 lenaa
850 nn
851 1)))
852 work-%offset%))
853 (if (< k np1)
854 (setf (f2cl-lib:fref work-%data%
855 (k)
856 ((1
857 (f2cl-lib:int-add
858 (f2cl-lib:int-mul 5
859 (f2cl-lib:int-add nn 1))
860 lenaa
861 nn
862 1)))
863 work-%offset%)
864 (+
865 (f2cl-lib:fref work-%data%
866 (k)
867 ((1
868 (f2cl-lib:int-add
869 (f2cl-lib:int-mul 5
870 (f2cl-lib:int-add nn
871 1))
872 lenaa
873 nn
874 1)))
875 work-%offset%)
876 (f2cl-lib:fref work-%data%
877 ((f2cl-lib:int-add n4p5 nn))
878 ((1
879 (f2cl-lib:int-add
880 (f2cl-lib:int-mul 5
881 (f2cl-lib:int-add nn
882 1))
883 lenaa
884 nn
885 1)))
886 work-%offset%))))
887 (qimuds
888 (f2cl-lib:array-slice work-%data%
889 double-float
890 (n5p6)
891 ((1
892 (f2cl-lib:int-add
893 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
894 lenaa
895 nn
896 1)))
897 work-%offset%)
898 work maxa nn lenaa)
899 (daxpy np1 (- au) work 1
900 (f2cl-lib:array-slice work-%data%
901 double-float
902 (n3p4)
903 ((1
904 (f2cl-lib:int-add
905 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
906 lenaa
907 nn
908 1)))
909 work-%offset%)
910 1)
911 (setf rnprd
912 (ddot np1
913 (f2cl-lib:array-slice work-%data%
914 double-float
915 (n3p4)
916 ((1
917 (f2cl-lib:int-add
918 (f2cl-lib:int-mul 5
919 (f2cl-lib:int-add
920 nn
921 1))
922 lenaa
923 nn
924 1)))
925 work-%offset%)
926 1
927 (f2cl-lib:array-slice work-%data%
928 double-float
929 (n3p4)
930 ((1
931 (f2cl-lib:int-add
932 (f2cl-lib:int-mul 5
933 (f2cl-lib:int-add
934 nn
935 1))
936 lenaa
937 nn
938 1)))
939 work-%offset%)
940 1))
941 (setf bu (/ rnprd runprd))
942 (setf runprd rnprd)
943 (dcopy np1
944 (f2cl-lib:array-slice work-%data%
945 double-float
946 (n3p4)
947 ((1
948 (f2cl-lib:int-add
949 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
950 lenaa
951 nn
952 1)))
953 work-%offset%)
954 1 work 1)
955 (qimuds
956 (f2cl-lib:array-slice work-%data%
957 double-float
958 (n5p6)
959 ((1
960 (f2cl-lib:int-add
961 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
962 lenaa
963 nn
964 1)))
965 work-%offset%)
966 work maxa nn lenaa)
967 (multds start aa work maxa nn lenaa)
968 (setf (f2cl-lib:fref start-%data%
969 (np1)
970 ((1 (f2cl-lib:int-add nn 1)))
971 start-%offset%)
972 (f2cl-lib:fref work-%data%
973 (k)
974 ((1
975 (f2cl-lib:int-add
976 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
977 lenaa
978 nn
979 1)))
980 work-%offset%))
981 (if (< k np1)
982 (setf (f2cl-lib:fref start-%data%
983 (k)
984 ((1 (f2cl-lib:int-add nn 1)))
985 start-%offset%)
986 (+
987 (f2cl-lib:fref start-%data%
988 (k)
989 ((1 (f2cl-lib:int-add nn 1)))
990 start-%offset%)
991 (f2cl-lib:fref work-%data%
992 (np1)
993 ((1
994 (f2cl-lib:int-add
995 (f2cl-lib:int-mul 5
996 (f2cl-lib:int-add nn
997 1))
998 lenaa
999 nn
1000 1)))
1001 work-%offset%))))
1002 (daxpy np1 bu
1003 (f2cl-lib:array-slice work-%data%
1004 double-float
1005 (n4p5)
1006 ((1
1007 (f2cl-lib:int-add
1008 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1009 lenaa
1010 nn
1011 1)))
1012 work-%offset%)
1013 1 start 1)
1014 (dcopy np1 start 1
1015 (f2cl-lib:array-slice work-%data%
1016 double-float
1017 (n4p5)
1018 ((1
1019 (f2cl-lib:int-add
1020 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1021 lenaa
1022 nn
1023 1)))
1024 work-%offset%)
1025 1)
1026 (setf punprd
1027 (ddot np1
1028 (f2cl-lib:array-slice work-%data%
1029 double-float
1030 (n4p5)
1031 ((1
1032 (f2cl-lib:int-add
1033 (f2cl-lib:int-mul 5
1034 (f2cl-lib:int-add
1035 nn
1036 1))
1037 lenaa
1038 nn
1039 1)))
1040 work-%offset%)
1041 1
1042 (f2cl-lib:array-slice work-%data%
1043 double-float
1044 (n4p5)
1045 ((1
1046 (f2cl-lib:int-add
1047 (f2cl-lib:int-mul 5
1048 (f2cl-lib:int-add
1049 nn
1050 1))
1051 lenaa
1052 nn
1053 1)))
1054 work-%offset%)
1055 1))))
1056 (setf j (f2cl-lib:int-add j 1))
1057 (go label100)
1058 label200
1059 (cond
1060 ((> j imax)
1061 (setf iflag 4)
1062 (go end_label)))
1063 (multds
1064 (f2cl-lib:array-slice work-%data%
1065 double-float
1066 (n3p4)
1067 ((1
1068 (f2cl-lib:int-add
1069 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1070 lenaa
1071 nn
1072 1)))
1073 work-%offset%)
1074 aa
1075 (f2cl-lib:array-slice work-%data%
1076 double-float
1077 (n2p3)
1078 ((1
1079 (f2cl-lib:int-add
1080 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1081 lenaa
1082 nn
1083 1)))
1084 work-%offset%)
1085 maxa nn lenaa)
1086 (setf (f2cl-lib:fref work-%data%
1087 ((f2cl-lib:int-add n3p4 nn))
1088 ((1
1089 (f2cl-lib:int-add
1090 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1091 lenaa
1092 nn
1093 1)))
1094 work-%offset%)
1095 (f2cl-lib:fref work-%data%
1096 ((f2cl-lib:int-sub (f2cl-lib:int-add n2p3 k) 1))
1097 ((1
1098 (f2cl-lib:int-add
1099 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1100 lenaa
1101 nn
1102 1)))
1103 work-%offset%))
1104 (setf ind (f2cl-lib:int-sub (f2cl-lib:int-add n3p4 k) 1))
1105 (if (< k np1)
1106 (setf (f2cl-lib:fref work-%data%
1107 (ind)
1108 ((1
1109 (f2cl-lib:int-add
1110 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1111 lenaa
1112 nn
1113 1)))
1114 work-%offset%)
1115 (+
1116 (f2cl-lib:fref work-%data%
1117 (ind)
1118 ((1
1119 (f2cl-lib:int-add
1120 (f2cl-lib:int-mul 5
1121 (f2cl-lib:int-add nn 1))
1122 lenaa
1123 nn
1124 1)))
1125 work-%offset%)
1126 (f2cl-lib:fref work-%data%
1127 ((f2cl-lib:int-add n2p3 nn))
1128 ((1
1129 (f2cl-lib:int-add
1130 (f2cl-lib:int-mul 5
1131 (f2cl-lib:int-add nn 1))
1132 lenaa
1133 nn
1134 1)))
1135 work-%offset%))))
1136 (dscal np1 -1.0
1137 (f2cl-lib:array-slice work-%data%
1138 double-float
1139 (n3p4)
1140 ((1
1141 (f2cl-lib:int-add
1142 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1143 lenaa
1144 nn
1145 1)))
1146 work-%offset%)
1147 1)
1148 (setf (f2cl-lib:fref work-%data%
1149 ((f2cl-lib:int-add n3p4 nn))
1150 ((1
1151 (f2cl-lib:int-add
1152 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1153 lenaa
1154 nn
1155 1)))
1156 work-%offset%)
1157 (+ startk
1158 (f2cl-lib:fref work-%data%
1159 ((f2cl-lib:int-add n3p4 nn))
1160 ((1
1161 (f2cl-lib:int-add
1162 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1163 lenaa
1164 nn
1165 1)))
1166 work-%offset%)))
1167 (qimuds
1168 (f2cl-lib:array-slice work-%data%
1169 double-float
1170 (n5p6)
1171 ((1
1172 (f2cl-lib:int-add
1173 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1174 lenaa
1175 nn
1176 1)))
1177 work-%offset%)
1178 (f2cl-lib:array-slice work-%data%
1179 double-float
1180 (n3p4)
1181 ((1
1182 (f2cl-lib:int-add
1183 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1184 lenaa
1185 nn
1186 1)))
1187 work-%offset%)
1188 maxa nn lenaa)
1189 (dcopy np1
1190 (f2cl-lib:array-slice work-%data%
1191 double-float
1192 (n3p4)
1193 ((1
1194 (f2cl-lib:int-add
1195 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1196 lenaa
1197 nn
1198 1)))
1199 work-%offset%)
1200 1 work 1)
1201 (qimuds
1202 (f2cl-lib:array-slice work-%data%
1203 double-float
1204 (n5p6)
1205 ((1
1206 (f2cl-lib:int-add
1207 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1208 lenaa
1209 nn
1210 1)))
1211 work-%offset%)
1212 work maxa nn lenaa)
1213 (multds
1214 (f2cl-lib:array-slice work-%data%
1215 double-float
1216 (n4p5)
1217 ((1
1218 (f2cl-lib:int-add
1219 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1220 lenaa
1221 nn
1222 1)))
1223 work-%offset%)
1224 aa work maxa nn lenaa)
1225 (setf (f2cl-lib:fref work-%data%
1226 ((f2cl-lib:int-add n4p5 nn))
1227 ((1
1228 (f2cl-lib:int-add
1229 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1230 lenaa
1231 nn
1232 1)))
1233 work-%offset%)
1234 (f2cl-lib:fref work-%data%
1235 (k)
1236 ((1
1237 (f2cl-lib:int-add
1238 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1239 lenaa
1240 nn
1241 1)))
1242 work-%offset%))
1243 (if (< k np1)
1244 (setf (f2cl-lib:fref work-%data%
1245 ((f2cl-lib:int-sub (f2cl-lib:int-add n4p5 k) 1))
1246 ((1
1247 (f2cl-lib:int-add
1248 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1249 lenaa
1250 nn
1251 1)))
1252 work-%offset%)
1253 (+
1254 (f2cl-lib:fref work-%data%
1255 ((f2cl-lib:int-sub (f2cl-lib:int-add n4p5 k)
1256 1))
1257 ((1
1258 (f2cl-lib:int-add
1259 (f2cl-lib:int-mul 5
1260 (f2cl-lib:int-add nn 1))
1261 lenaa
1262 nn
1263 1)))
1264 work-%offset%)
1265 (f2cl-lib:fref work-%data%
1266 (np1)
1267 ((1
1268 (f2cl-lib:int-add
1269 (f2cl-lib:int-mul 5
1270 (f2cl-lib:int-add nn 1))
1271 lenaa
1272 nn
1273 1)))
1274 work-%offset%))))
1275 (setf rbnprd
1276 (ddot np1
1277 (f2cl-lib:array-slice work-%data%
1278 double-float
1279 (n3p4)
1280 ((1
1281 (f2cl-lib:int-add
1282 (f2cl-lib:int-mul 5
1283 (f2cl-lib:int-add nn
1284 1))
1285 lenaa
1286 nn
1287 1)))
1288 work-%offset%)
1289 1
1290 (f2cl-lib:array-slice work-%data%
1291 double-float
1292 (n3p4)
1293 ((1
1294 (f2cl-lib:int-add
1295 (f2cl-lib:int-mul 5
1296 (f2cl-lib:int-add nn
1297 1))
1298 lenaa
1299 nn
1300 1)))
1301 work-%offset%)
1302 1))
1303 (setf pbnprd
1304 (ddot np1
1305 (f2cl-lib:array-slice work-%data%
1306 double-float
1307 (n4p5)
1308 ((1
1309 (f2cl-lib:int-add
1310 (f2cl-lib:int-mul 5
1311 (f2cl-lib:int-add nn
1312 1))
1313 lenaa
1314 nn
1315 1)))
1316 work-%offset%)
1317 1
1318 (f2cl-lib:array-slice work-%data%
1319 double-float
1320 (n4p5)
1321 ((1
1322 (f2cl-lib:int-add
1323 (f2cl-lib:int-mul 5
1324 (f2cl-lib:int-add nn
1325 1))
1326 lenaa
1327 nn
1328 1)))
1329 work-%offset%)
1330 1))
1331 (setf j 1)
1332 label300
1333 (if (not (and stillb (<= j imax))) (go label400))
1334 (cond
1335 ((> (f2cl-lib:fsqrt rbnprd) rbtol)
1336 (cond
1337 ((= pbnprd 0.0f0)
1338 (multds
1339 (f2cl-lib:array-slice work-%data%
1340 double-float
1341 (n3p4)
1342 ((1
1343 (f2cl-lib:int-add
1344 (f2cl-lib:int-mul 5
1345 (f2cl-lib:int-add nn
1346 1))
1347 lenaa
1348 nn
1349 1)))
1350 work-%offset%)
1351 aa
1352 (f2cl-lib:array-slice work-%data%
1353 double-float
1354 (n2p3)
1355 ((1
1356 (f2cl-lib:int-add
1357 (f2cl-lib:int-mul 5
1358 (f2cl-lib:int-add nn
1359 1))
1360 lenaa
1361 nn
1362 1)))
1363 work-%offset%)
1364 maxa nn lenaa)
1365 (setf (f2cl-lib:fref work-%data%
1366 ((f2cl-lib:int-add n3p4 nn))
1367 ((1
1368 (f2cl-lib:int-add
1369 (f2cl-lib:int-mul 5
1370 (f2cl-lib:int-add nn 1))
1371 lenaa
1372 nn
1373 1)))
1374 work-%offset%)
1375 (f2cl-lib:fref work-%data%
1376 ((f2cl-lib:int-sub (f2cl-lib:int-add n2p3 k)
1377 1))
1378 ((1
1379 (f2cl-lib:int-add
1380 (f2cl-lib:int-mul 5
1381 (f2cl-lib:int-add nn
1382 1))
1383 lenaa
1384 nn
1385 1)))
1386 work-%offset%))
1387 (setf ind (f2cl-lib:int-sub (f2cl-lib:int-add n3p4 k) 1))
1388 (if (< k np1)
1389 (setf (f2cl-lib:fref work-%data%
1390 (ind)
1391 ((1
1392 (f2cl-lib:int-add
1393 (f2cl-lib:int-mul 5
1394 (f2cl-lib:int-add nn
1395 1))
1396 lenaa
1397 nn
1398 1)))
1399 work-%offset%)
1400 (+
1401 (f2cl-lib:fref work-%data%
1402 (ind)
1403 ((1
1404 (f2cl-lib:int-add
1405 (f2cl-lib:int-mul 5
1406 (f2cl-lib:int-add
1407 nn
1408 1))
1409 lenaa
1410 nn
1411 1)))
1412 work-%offset%)
1413 (f2cl-lib:fref work-%data%
1414 ((f2cl-lib:int-add n2p3 nn))
1415 ((1
1416 (f2cl-lib:int-add
1417 (f2cl-lib:int-mul 5
1418 (f2cl-lib:int-add
1419 nn
1420 1))
1421 lenaa
1422 nn
1423 1)))
1424 work-%offset%))))
1425 (dscal np1 -1.0
1426 (f2cl-lib:array-slice work-%data%
1427 double-float
1428 (n3p4)
1429 ((1
1430 (f2cl-lib:int-add
1431 (f2cl-lib:int-mul 5
1432 (f2cl-lib:int-add nn
1433 1))
1434 lenaa
1435 nn
1436 1)))
1437 work-%offset%)
1438 1)
1439 (setf (f2cl-lib:fref work-%data%
1440 ((f2cl-lib:int-add n3p4 nn))
1441 ((1
1442 (f2cl-lib:int-add
1443 (f2cl-lib:int-mul 5
1444 (f2cl-lib:int-add nn 1))
1445 lenaa
1446 nn
1447 1)))
1448 work-%offset%)
1449 (+ startk
1450 (f2cl-lib:fref work-%data%
1451 ((f2cl-lib:int-add n3p4 nn))
1452 ((1
1453 (f2cl-lib:int-add
1454 (f2cl-lib:int-mul 5
1455 (f2cl-lib:int-add nn
1456 1))
1457 lenaa
1458 nn
1459 1)))
1460 work-%offset%)))
1461 (qimuds
1462 (f2cl-lib:array-slice work-%data%
1463 double-float
1464 (n5p6)
1465 ((1
1466 (f2cl-lib:int-add
1467 (f2cl-lib:int-mul 5
1468 (f2cl-lib:int-add nn
1469 1))
1470 lenaa
1471 nn
1472 1)))
1473 work-%offset%)
1474 (f2cl-lib:array-slice work-%data%
1475 double-float
1476 (n3p4)
1477 ((1
1478 (f2cl-lib:int-add
1479 (f2cl-lib:int-mul 5
1480 (f2cl-lib:int-add nn
1481 1))
1482 lenaa
1483 nn
1484 1)))
1485 work-%offset%)
1486 maxa nn lenaa)
1487 (dcopy np1
1488 (f2cl-lib:array-slice work-%data%
1489 double-float
1490 (n3p4)
1491 ((1
1492 (f2cl-lib:int-add
1493 (f2cl-lib:int-mul 5
1494 (f2cl-lib:int-add nn
1495 1))
1496 lenaa
1497 nn
1498 1)))
1499 work-%offset%)
1500 1 work 1)
1501 (qimuds
1502 (f2cl-lib:array-slice work-%data%
1503 double-float
1504 (n5p6)
1505 ((1
1506 (f2cl-lib:int-add
1507 (f2cl-lib:int-mul 5
1508 (f2cl-lib:int-add nn
1509 1))
1510 lenaa
1511 nn
1512 1)))
1513 work-%offset%)
1514 work maxa nn lenaa)
1515 (multds
1516 (f2cl-lib:array-slice work-%data%
1517 double-float
1518 (n4p5)
1519 ((1
1520 (f2cl-lib:int-add
1521 (f2cl-lib:int-mul 5
1522 (f2cl-lib:int-add nn
1523 1))
1524 lenaa
1525 nn
1526 1)))
1527 work-%offset%)
1528 aa work maxa nn lenaa)
1529 (setf (f2cl-lib:fref work-%data%
1530 ((f2cl-lib:int-add n4p5 nn))
1531 ((1
1532 (f2cl-lib:int-add
1533 (f2cl-lib:int-mul 5
1534 (f2cl-lib:int-add nn 1))
1535 lenaa
1536 nn
1537 1)))
1538 work-%offset%)
1539 (f2cl-lib:fref work-%data%
1540 (k)
1541 ((1
1542 (f2cl-lib:int-add
1543 (f2cl-lib:int-mul 5
1544 (f2cl-lib:int-add nn
1545 1))
1546 lenaa
1547 nn
1548 1)))
1549 work-%offset%))
1550 (setf ind (f2cl-lib:int-sub (f2cl-lib:int-add n4p5 k) 1))
1551 (if (< k np1)
1552 (setf (f2cl-lib:fref work-%data%
1553 (ind)
1554 ((1
1555 (f2cl-lib:int-add
1556 (f2cl-lib:int-mul 5
1557 (f2cl-lib:int-add nn
1558 1))
1559 lenaa
1560 nn
1561 1)))
1562 work-%offset%)
1563 (+
1564 (f2cl-lib:fref work-%data%
1565 (ind)
1566 ((1
1567 (f2cl-lib:int-add
1568 (f2cl-lib:int-mul 5
1569 (f2cl-lib:int-add
1570 nn
1571 1))
1572 lenaa
1573 nn
1574 1)))
1575 work-%offset%)
1576 (f2cl-lib:fref work-%data%
1577 (np1)
1578 ((1
1579 (f2cl-lib:int-add
1580 (f2cl-lib:int-mul 5
1581 (f2cl-lib:int-add
1582 nn
1583 1))
1584 lenaa
1585 nn
1586 1)))
1587 work-%offset%))))
1588 (setf rbnprd
1589 (ddot np1
1590 (f2cl-lib:array-slice work-%data%
1591 double-float
1592 (n3p4)
1593 ((1
1594 (f2cl-lib:int-add
1595 (f2cl-lib:int-mul 5
1596 (f2cl-lib:int-add
1597 nn
1598 1))
1599 lenaa
1600 nn
1601 1)))
1602 work-%offset%)
1603 1
1604 (f2cl-lib:array-slice work-%data%
1605 double-float
1606 (n3p4)
1607 ((1
1608 (f2cl-lib:int-add
1609 (f2cl-lib:int-mul 5
1610 (f2cl-lib:int-add
1611 nn
1612 1))
1613 lenaa
1614 nn
1615 1)))
1616 work-%offset%)
1617 1))
1618 (setf pbnprd
1619 (ddot np1
1620 (f2cl-lib:array-slice work-%data%
1621 double-float
1622 (n4p5)
1623 ((1
1624 (f2cl-lib:int-add
1625 (f2cl-lib:int-mul 5
1626 (f2cl-lib:int-add
1627 nn
1628 1))
1629 lenaa
1630 nn
1631 1)))
1632 work-%offset%)
1633 1
1634 (f2cl-lib:array-slice work-%data%
1635 double-float
1636 (n4p5)
1637 ((1
1638 (f2cl-lib:int-add
1639 (f2cl-lib:int-mul 5
1640 (f2cl-lib:int-add
1641 nn
1642 1))
1643 lenaa
1644 nn
1645 1)))
1646 work-%offset%)
1647 1))
1648 (cond
1649 ((<= (f2cl-lib:fsqrt rbnprd) rbtol)
1650 (setf stillb f2cl-lib:%false%)))))
1651 (cond
1652 (stillb
1653 (setf ab (/ rbnprd pbnprd))
1654 (dcopy np1
1655 (f2cl-lib:array-slice work-%data%
1656 double-float
1657 (n2p3)
1658 ((1
1659 (f2cl-lib:int-add
1660 (f2cl-lib:int-mul 5
1661 (f2cl-lib:int-add nn
1662 1))
1663 lenaa
1664 nn
1665 1)))
1666 work-%offset%)
1667 1 work 1)
1668 (daxpy np1 ab
1669 (f2cl-lib:array-slice work-%data%
1670 double-float
1671 (n4p5)
1672 ((1
1673 (f2cl-lib:int-add
1674 (f2cl-lib:int-mul 5
1675 (f2cl-lib:int-add nn
1676 1))
1677 lenaa
1678 nn
1679 1)))
1680 work-%offset%)
1681 1
1682 (f2cl-lib:array-slice work-%data%
1683 double-float
1684 (n2p3)
1685 ((1
1686 (f2cl-lib:int-add
1687 (f2cl-lib:int-mul 5
1688 (f2cl-lib:int-add nn
1689 1))
1690 lenaa
1691 nn
1692 1)))
1693 work-%offset%)
1694 1)
1695 (daxpy np1 -1.0
1696 (f2cl-lib:array-slice work-%data%
1697 double-float
1698 (n2p3)
1699 ((1
1700 (f2cl-lib:int-add
1701 (f2cl-lib:int-mul 5
1702 (f2cl-lib:int-add nn
1703 1))
1704 lenaa
1705 nn
1706 1)))
1707 work-%offset%)
1708 1 work 1)
1709 (setf zlen
1710 (dnrm2 np1
1711 (f2cl-lib:array-slice work-%data%
1712 double-float
1713 (n2p3)
1714 ((1
1715 (f2cl-lib:int-add
1716 (f2cl-lib:int-mul 5
1717 (f2cl-lib:int-add
1718 nn
1719 1))
1720 lenaa
1721 nn
1722 1)))
1723 work-%offset%)
1724 1))
1725 (setf dznrm (dnrm2 np1 work 1))
1726 (if (< (/ dznrm zlen) ztol) (setf stillb f2cl-lib:%false%)))))
1727 (t
1728 (setf stillb f2cl-lib:%false%)))
1729 (cond
1730 (stillb
1731 (multds work aa
1732 (f2cl-lib:array-slice work-%data%
1733 double-float
1734 (n4p5)
1735 ((1
1736 (f2cl-lib:int-add
1737 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1738 lenaa
1739 nn
1740 1)))
1741 work-%offset%)
1742 maxa nn lenaa)
1743 (setf (f2cl-lib:fref work-%data%
1744 (np1)
1745 ((1
1746 (f2cl-lib:int-add
1747 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1748 lenaa
1749 nn
1750 1)))
1751 work-%offset%)
1752 (f2cl-lib:fref work-%data%
1753 ((f2cl-lib:int-sub (f2cl-lib:int-add n4p5 k)
1754 1))
1755 ((1
1756 (f2cl-lib:int-add
1757 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1758 lenaa
1759 nn
1760 1)))
1761 work-%offset%))
1762 (if (< k np1)
1763 (setf (f2cl-lib:fref work-%data%
1764 (k)
1765 ((1
1766 (f2cl-lib:int-add
1767 (f2cl-lib:int-mul 5
1768 (f2cl-lib:int-add nn 1))
1769 lenaa
1770 nn
1771 1)))
1772 work-%offset%)
1773 (+
1774 (f2cl-lib:fref work-%data%
1775 (k)
1776 ((1
1777 (f2cl-lib:int-add
1778 (f2cl-lib:int-mul 5
1779 (f2cl-lib:int-add nn
1780 1))
1781 lenaa
1782 nn
1783 1)))
1784 work-%offset%)
1785 (f2cl-lib:fref work-%data%
1786 ((f2cl-lib:int-add n4p5 nn))
1787 ((1
1788 (f2cl-lib:int-add
1789 (f2cl-lib:int-mul 5
1790 (f2cl-lib:int-add nn
1791 1))
1792 lenaa
1793 nn
1794 1)))
1795 work-%offset%))))
1796 (qimuds
1797 (f2cl-lib:array-slice work-%data%
1798 double-float
1799 (n5p6)
1800 ((1
1801 (f2cl-lib:int-add
1802 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1803 lenaa
1804 nn
1805 1)))
1806 work-%offset%)
1807 work maxa nn lenaa)
1808 (daxpy np1 (- ab) work 1
1809 (f2cl-lib:array-slice work-%data%
1810 double-float
1811 (n3p4)
1812 ((1
1813 (f2cl-lib:int-add
1814 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1815 lenaa
1816 nn
1817 1)))
1818 work-%offset%)
1819 1)
1820 (setf rnprd
1821 (ddot np1
1822 (f2cl-lib:array-slice work-%data%
1823 double-float
1824 (n3p4)
1825 ((1
1826 (f2cl-lib:int-add
1827 (f2cl-lib:int-mul 5
1828 (f2cl-lib:int-add
1829 nn
1830 1))
1831 lenaa
1832 nn
1833 1)))
1834 work-%offset%)
1835 1
1836 (f2cl-lib:array-slice work-%data%
1837 double-float
1838 (n3p4)
1839 ((1
1840 (f2cl-lib:int-add
1841 (f2cl-lib:int-mul 5
1842 (f2cl-lib:int-add
1843 nn
1844 1))
1845 lenaa
1846 nn
1847 1)))
1848 work-%offset%)
1849 1))
1850 (setf bb (/ rnprd rbnprd))
1851 (setf rbnprd rnprd)
1852 (dcopy np1
1853 (f2cl-lib:array-slice work-%data%
1854 double-float
1855 (n3p4)
1856 ((1
1857 (f2cl-lib:int-add
1858 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1859 lenaa
1860 nn
1861 1)))
1862 work-%offset%)
1863 1 work 1)
1864 (qimuds
1865 (f2cl-lib:array-slice work-%data%
1866 double-float
1867 (n5p6)
1868 ((1
1869 (f2cl-lib:int-add
1870 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1871 lenaa
1872 nn
1873 1)))
1874 work-%offset%)
1875 work maxa nn lenaa)
1876 (multds start aa work maxa nn lenaa)
1877 (setf (f2cl-lib:fref start-%data%
1878 (np1)
1879 ((1 (f2cl-lib:int-add nn 1)))
1880 start-%offset%)
1881 (f2cl-lib:fref work-%data%
1882 (k)
1883 ((1
1884 (f2cl-lib:int-add
1885 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1886 lenaa
1887 nn
1888 1)))
1889 work-%offset%))
1890 (if (< k np1)
1891 (setf (f2cl-lib:fref start-%data%
1892 (k)
1893 ((1 (f2cl-lib:int-add nn 1)))
1894 start-%offset%)
1895 (+
1896 (f2cl-lib:fref start-%data%
1897 (k)
1898 ((1 (f2cl-lib:int-add nn 1)))
1899 start-%offset%)
1900 (f2cl-lib:fref work-%data%
1901 (np1)
1902 ((1
1903 (f2cl-lib:int-add
1904 (f2cl-lib:int-mul 5
1905 (f2cl-lib:int-add nn
1906 1))
1907 lenaa
1908 nn
1909 1)))
1910 work-%offset%))))
1911 (daxpy np1 bb
1912 (f2cl-lib:array-slice work-%data%
1913 double-float
1914 (n4p5)
1915 ((1
1916 (f2cl-lib:int-add
1917 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1918 lenaa
1919 nn
1920 1)))
1921 work-%offset%)
1922 1 start 1)
1923 (dcopy np1 start 1
1924 (f2cl-lib:array-slice work-%data%
1925 double-float
1926 (n4p5)
1927 ((1
1928 (f2cl-lib:int-add
1929 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1930 lenaa
1931 nn
1932 1)))
1933 work-%offset%)
1934 1)
1935 (setf pbnprd
1936 (ddot np1
1937 (f2cl-lib:array-slice work-%data%
1938 double-float
1939 (n4p5)
1940 ((1
1941 (f2cl-lib:int-add
1942 (f2cl-lib:int-mul 5
1943 (f2cl-lib:int-add
1944 nn
1945 1))
1946 lenaa
1947 nn
1948 1)))
1949 work-%offset%)
1950 1
1951 (f2cl-lib:array-slice work-%data%
1952 double-float
1953 (n4p5)
1954 ((1
1955 (f2cl-lib:int-add
1956 (f2cl-lib:int-mul 5
1957 (f2cl-lib:int-add
1958 nn
1959 1))
1960 lenaa
1961 nn
1962 1)))
1963 work-%offset%)
1964 1))))
1965 (setf j (f2cl-lib:int-add j 1))
1966 (go label300)
1967 label400
1968 (cond
1969 ((> j imax)
1970 (setf iflag 4)
1971 (go end_label)))
1972 (setf temp
1973 (/
1974 (-
1975 (f2cl-lib:fref work-%data%
1976 ((f2cl-lib:int-add n2p3 nn))
1977 ((1
1978 (f2cl-lib:int-add
1979 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
1980 lenaa
1981 nn
1982 1)))
1983 work-%offset%))
1984 (+ 1.0
1985 (f2cl-lib:fref work-%data%
1986 ((f2cl-lib:int-add np2 nn))
1987 ((1
1988 (f2cl-lib:int-add
1989 (f2cl-lib:int-mul 5
1990 (f2cl-lib:int-add nn 1))
1991 lenaa
1992 nn
1993 1)))
1994 work-%offset%))))
1995 (dcopy np1
1996 (f2cl-lib:array-slice work-%data%
1997 double-float
1998 (n2p3)
1999 ((1
2000 (f2cl-lib:int-add
2001 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
2002 lenaa
2003 nn
2004 1)))
2005 work-%offset%)
2006 1 start 1)
2007 (daxpy np1 temp
2008 (f2cl-lib:array-slice work-%data%
2009 double-float
2010 (np2)
2011 ((1
2012 (f2cl-lib:int-add
2013 (f2cl-lib:int-mul 5 (f2cl-lib:int-add nn 1))
2014 lenaa
2015 nn
2016 1)))
2017 work-%offset%)
2018 1 start 1)
2019 (go end_label)
2020 end_label
2021 (return (values nil nil nil nil nil nil nil iflag)))))
2022
2023 (in-package #:cl-user)
2024 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
2025 (eval-when (:load-toplevel :compile-toplevel :execute)
2026 (setf (gethash 'fortran-to-lisp::pcgds fortran-to-lisp::*f2cl-function-info*)
2027 (fortran-to-lisp::make-f2cl-finfo
2028 :arg-types '((fortran-to-lisp::integer4) (array double-float (*))
2029 (fortran-to-lisp::integer4)
2030 (array fortran-to-lisp::integer4 (*))
2031 (array double-float (*)) (array double-float (*))
2032 (array double-float (*)) (fortran-to-lisp::integer4))
2033 :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::iflag)
2034 :calls '(fortran-to-lisp::ddot fortran-to-lisp::daxpy
2035 fortran-to-lisp::dscal fortran-to-lisp::dnrm2
2036 fortran-to-lisp::dcopy fortran-to-lisp::idamax
2037 fortran-to-lisp::qimuds fortran-to-lisp::multds
2038 fortran-to-lisp::d1mach fortran-to-lisp::mfacds))))
2039
Maxima, a Computer Algebra System